28 char blackhole_eq_spher_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_eq_spher.C,v 1.5 2014/10/13 08:52:46 j_novak Exp $" ;
60 #include "blackhole.h" 66 #include "utilitaires.h" 67 #include "graphique.h" 71 double spin_omega,
double precis,
72 double precis_shift) {
98 rr.std_spectral_base() ;
114 ll.set(1) = st * cp ;
115 ll.
set(2) = st * sp ;
127 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
139 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
145 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
172 lapse_bh = 1. /
sqrt(1. + 2. * mass / rr) ;
177 for (
int i=1; i<=3; i++) {
178 shift_bh.
set(i) = 2. * lapse_bh * lapse_bh * mass * ll(i) / rr ;
190 r_are =
r_coord(neumann, first) ;
208 + cc*cc*
pow(mass/r_are/rr,4.))
219 + cc*cc*
pow(mass/r_are/rr,4.)) ;
230 for (
int i=1; i<=3; i++) {
231 shift.
set(i) = mass * mass * cc * ll(i) / rr / rr
237 for (
int i=1; i<=3; i++) {
278 double diff_lp = 1. ;
279 double diff_cf = 1. ;
280 double diff_sx = 1. ;
281 double diff_sy = 1. ;
282 double diff_sz = 1. ;
297 for (
int mer=0; (diff_lp > precis) && (diff_cf > precis) && (diff_sx > precis) && (diff_sy > precis) && (diff_sz > precis) && (mer < mermax); mer++) {
299 cout <<
"--------------------------------------------------" << endl ;
300 cout <<
"step: " << mer << endl ;
301 cout <<
"diff_lapconf = " << diff_lp << endl ;
302 cout <<
"diff_confo = " << diff_cf << endl ;
303 cout <<
"diff_shift : x = " << diff_sx
304 <<
" y = " << diff_sy <<
" z = " << diff_sz << endl ;
332 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
338 source_lapconf = 7. * lapconf_jm1 %
taij_quad 339 /
pow(confo_jm1, 8.) / 8. ;
362 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
386 lapconf_m1.set_etat_qcq() ;
416 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
426 source_confo = tmp1 ;
443 confo = confo_m1 + 1. ;
455 confo7 =
pow(confo_jm1, 7.) ;
459 for (
int i=1; i<=3; i++)
460 dlappsi.
set(i) = (lapconf_jm1.
deriv(i)
465 dlappsi.annule_domain(0) ;
470 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
480 source_shift.annule_domain(0) ;
490 source_shift.std_spectral_base() ;
492 for (
int i=1; i<=3; i++) {
493 if (source_shift(i).get_etat() != ETATZERO)
494 source_shift.set(i).filtre(4) ;
510 for (
int i=0; i<3; i++) {
511 source_p.set(i) =
Cmp(source_shift(i+1)) ;
517 for (
int i=0; i<3; i++) {
518 resu_p.set(i) =
Cmp(shift_jm1(i+1)) ;
532 poisson_vect_frontiere(1./3., source_p, resu_p,
533 bc_shif_x, bc_shif_y, bc_shif_z,
534 0, precis_shift, 14) ;
538 for (
int i=1; i<=3; i++)
541 for (
int i=1; i<=3; i++)
547 for (
int i=1; i<=3; i++)
553 for (
int i=1; i<=3; i++)
588 Tbl diff_lapconf(nz) ;
612 cout <<
"Relative difference in the lapconf function : " << endl ;
613 for (
int l=0; l<nz; l++) {
614 cout << diff_lapconf(l) <<
" " ;
618 diff_lp = diff_lapconf(1) ;
619 for (
int l=2; l<nz; l++) {
620 diff_lp += diff_lapconf(l) ;
628 cout <<
"Relative difference in the conformal factor : " << endl ;
629 for (
int l=0; l<nz; l++) {
630 cout << diff_confo(l) <<
" " ;
634 diff_cf = diff_confo(1) ;
635 for (
int l=2; l<nz; l++) {
636 diff_cf += diff_confo(l) ;
642 Tbl diff_shift_x(nz) ;
643 Tbl diff_shift_y(nz) ;
644 Tbl diff_shift_z(nz) ;
657 cout <<
"Relative difference in the shift vector (x) : " << endl ;
658 for (
int l=0; l<nz; l++) {
659 cout << diff_shift_x(l) <<
" " ;
663 diff_sx = diff_shift_x(1) ;
664 for (
int l=2; l<nz; l++) {
665 diff_sx += diff_shift_x(l) ;
680 cout <<
"Relative difference in the shift vector (y) : " << endl ;
681 for (
int l=0; l<nz; l++) {
682 cout << diff_shift_y(l) <<
" " ;
686 diff_sy = diff_shift_y(1) ;
687 for (
int l=2; l<nz; l++) {
688 diff_sy += diff_shift_y(l) ;
703 cout <<
"Relative difference in the shift vector (z) : " << endl ;
704 for (
int l=0; l<nz; l++) {
705 cout << diff_shift_z(l) <<
" " ;
709 diff_sz = diff_shift_z(1) ;
710 for (
int l=2; l<nz; l++) {
711 diff_sz += diff_shift_z(l) ;
718 cout <<
"Mass_bh : " << mass_bh / msol <<
" [M_sol]" << endl ;
719 double rad_apphor =
rad_ah() ;
720 cout <<
" : " << 0.5 * rad_apphor / ggrav / msol
721 <<
" [M_sol]" << endl ;
726 cout <<
"Mass_bh : " << mass_bh / msol
727 <<
" [M_sol]" << endl ;
734 double irr_gm, adm_gm, kom_gm ;
735 irr_gm =
mass_irr() / mass_bh - 1. ;
736 adm_gm =
mass_adm() / mass_bh - 1. ;
737 kom_gm =
mass_kom() / mass_bh - 1. ;
738 cout <<
"Irreducible mass : " <<
mass_irr() / msol << endl ;
739 cout <<
"Gravitaitonal mass : " << mass_bh / msol << endl ;
740 cout <<
"ADM mass : " <<
mass_adm() / msol << endl ;
741 cout <<
"Komar mass : " <<
mass_kom() / msol << endl ;
748 cout <<
"Diff. (Mirr-Mg)/Mg : " << irr_gm << endl ;
749 cout <<
"Diff. (Madm-Mg)/Mg : " << adm_gm << endl ;
750 cout <<
"Diff. (Mkom-Mg)/Mg : " << kom_gm << endl ;
790 lapconf_exact = 1./
sqrt(1.+2.*mass/rr) ;
794 for (
int i=1; i<=3; i++)
796 2.*mass*lapconf_exact%lapconf_exact%ll(i)/rr ;
802 rare =
r_coord(neumann, first) ;
805 lapconf_exact =
sqrt(1. - 2.*mass/rare/rr
806 + cc*cc*
pow(mass/rare/rr,4.)) *
sqrt(rare) ;
808 confo_exact =
sqrt(rare) ;
810 for (
int i=1; i<=3; i++) {
811 shift_exact.set(i) = mass * mass * cc * ll(i) / rr / rr
825 shift_exact.annule_domain(0) ;
826 shift_exact.std_spectral_base() ;
827 for (
int i=1; i<=3; i++)
828 shift_exact.set(i).raccord(1) ;
877 diff_lapconf_exact.
set(0) = 0. ;
878 cout <<
"Relative difference in the lapconf function : " << endl ;
879 for (
int l=0; l<nz; l++) {
880 cout << diff_lapconf_exact(l) <<
" " ;
886 diff_confo_exact.
set(0) = 0. ;
887 cout <<
"Relative difference in the conformal factor : " << endl ;
888 for (
int l=0; l<nz; l++) {
889 cout << diff_confo_exact(l) <<
" " ;
898 cout <<
"Relative difference in the shift vector (x) : " << endl ;
899 for (
int l=0; l<nz; l++) {
900 cout << diff_shift_exact_x(l) <<
" " ;
903 cout <<
"Relative difference in the shift vector (y) : " << endl ;
904 for (
int l=0; l<nz; l++) {
905 cout << diff_shift_exact_y(l) <<
" " ;
908 cout <<
"Relative difference in the shift vector (z) : " << endl ;
909 for (
int l=0; l<nz; l++) {
910 cout << diff_shift_exact_z(l) <<
" " ;
Map & mp
Mapping associated with the black hole.
virtual double mass_irr() const
Irreducible mass of the black hole.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Scalar lapconf_rs
Part of lapconf from the numerical computation.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Scalar taij_quad
Part of the scalar generated by .
const Valeur bc_lapconf(bool neumann, bool first) const
Boundary condition on the apparent horizon of the black hole for the lapse function: 2-D Valeur...
Cmp sqrt(const Cmp &)
Square root.
virtual void del_deriv() const
Deletes all the derived quantities.
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
double mass_bh
Gravitational mass of BH.
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
virtual double mass_adm() const
ADM mass.
Values and coefficients of a (real-value) function.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tensor field of valence 1.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
void set_dzpuis(int)
Modifies the dzpuis flag.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Scalar confo
Conformal factor generated by the black hole.
const Valeur bc_shift_x(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Scalar lapconf_bh
Part of lapconf from the analytic background.
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Scalar poisson_neumann(const Valeur &, int) const
Idem as Scalar::poisson_dirichlet , the boundary condition being on the radial derivative of the solu...
Vector shift_bh
Part of the shift vector from the analytic background.
int get_dzpuis() const
Returns dzpuis.
void equilibrium_spher(bool neumann, bool first, double spin_omega, double precis=1.e-14, double precis_shift=1.e-8)
Computes a spherical, static black-hole by giving boundary conditions on the apparent horizon...
const Valeur bc_shift_y(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
int get_nzone() const
Returns the number of domains.
virtual double rad_ah() const
Radius of the apparent horizon.
const Scalar & deriv(int i) const
Returns of *this , where .
Vector shift_rs
Part of the shift vector from the numerical computation.
Cmp pow(const Cmp &, int)
Power .
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual double mass_kom() const
Komar mass.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Scalar lapse
Lapse function generated by the black hole.
Vector shift
Shift vector generated by the black hole.
Sym_tensor taij
Trace of the extrinsic curvature.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Scalar & set(int)
Read/write access to a component.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
const Valeur bc_confo() const
Boundary condition on the apparent horizon of the black hole for the conformal factor: 2-D Valeur...
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Coord r
r coordinate centered on the grid