63 #include "blackhole.h" 69 #include "utilitaires.h" 70 #include "graphique.h" 74 double spin_omega,
double precis,
75 double precis_shift) {
101 rr.std_spectral_base() ;
117 ll.set(1) = st * cp ;
118 ll.
set(2) = st * sp ;
130 cc = 2. * (
sqrt(13.) - 1.) / 3. ;
142 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
148 cout <<
"!!!!! WARNING: Not yet prepared !!!!!" << endl ;
175 lapse_bh = 1. /
sqrt(1. + 2. * mass / rr) ;
180 for (
int i=1; i<=3; i++) {
181 shift_bh.
set(i) = 2. * lapse_bh * lapse_bh * mass * ll(i) / rr ;
193 r_are =
r_coord(neumann, first) ;
211 + cc*cc*
pow(mass/r_are/rr,4.))
222 + cc*cc*
pow(mass/r_are/rr,4.)) ;
233 for (
int i=1; i<=3; i++) {
234 shift.
set(i) = mass * mass * cc * ll(i) / rr / rr
240 for (
int i=1; i<=3; i++) {
281 double diff_lp = 1. ;
282 double diff_cf = 1. ;
283 double diff_sx = 1. ;
284 double diff_sy = 1. ;
285 double diff_sz = 1. ;
300 for (
int mer=0; (diff_lp > precis) && (diff_cf > precis) && (diff_sx > precis) && (diff_sy > precis) && (diff_sz > precis) && (mer < mermax); mer++) {
302 cout <<
"--------------------------------------------------" << endl ;
303 cout <<
"step: " << mer << endl ;
304 cout <<
"diff_lapconf = " << diff_lp << endl ;
305 cout <<
"diff_confo = " << diff_cf << endl ;
306 cout <<
"diff_shift : x = " << diff_sx
307 <<
" y = " << diff_sy <<
" z = " << diff_sz << endl ;
335 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
341 source_lapconf = 7. * lapconf_jm1 %
taij_quad 342 /
pow(confo_jm1, 8.) / 8. ;
365 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
389 lapconf_m1.set_etat_qcq() ;
419 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
429 source_confo = tmp1 ;
446 confo = confo_m1 + 1. ;
458 confo7 =
pow(confo_jm1, 7.) ;
462 for (
int i=1; i<=3; i++)
463 dlappsi.
set(i) = (lapconf_jm1.
deriv(i)
468 dlappsi.annule_domain(0) ;
473 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
483 source_shift.annule_domain(0) ;
493 source_shift.std_spectral_base() ;
495 for (
int i=1; i<=3; i++) {
496 if (source_shift(i).get_etat() != ETATZERO)
497 source_shift.set(i).filtre(4) ;
513 for (
int i=0; i<3; i++) {
514 source_p.set(i) =
Cmp(source_shift(i+1)) ;
520 for (
int i=0; i<3; i++) {
521 resu_p.set(i) =
Cmp(shift_jm1(i+1)) ;
535 poisson_vect_frontiere(1./3., source_p, resu_p,
536 bc_shif_x, bc_shif_y, bc_shif_z,
537 0, precis_shift, 14) ;
541 for (
int i=1; i<=3; i++)
544 for (
int i=1; i<=3; i++)
550 for (
int i=1; i<=3; i++)
556 for (
int i=1; i<=3; i++)
591 Tbl diff_lapconf(nz) ;
615 cout <<
"Relative difference in the lapconf function : " << endl ;
616 for (
int l=0; l<nz; l++) {
617 cout << diff_lapconf(l) <<
" " ;
621 diff_lp = diff_lapconf(1) ;
622 for (
int l=2; l<nz; l++) {
623 diff_lp += diff_lapconf(l) ;
631 cout <<
"Relative difference in the conformal factor : " << endl ;
632 for (
int l=0; l<nz; l++) {
633 cout << diff_confo(l) <<
" " ;
637 diff_cf = diff_confo(1) ;
638 for (
int l=2; l<nz; l++) {
639 diff_cf += diff_confo(l) ;
645 Tbl diff_shift_x(nz) ;
646 Tbl diff_shift_y(nz) ;
647 Tbl diff_shift_z(nz) ;
660 cout <<
"Relative difference in the shift vector (x) : " << endl ;
661 for (
int l=0; l<nz; l++) {
662 cout << diff_shift_x(l) <<
" " ;
666 diff_sx = diff_shift_x(1) ;
667 for (
int l=2; l<nz; l++) {
668 diff_sx += diff_shift_x(l) ;
683 cout <<
"Relative difference in the shift vector (y) : " << endl ;
684 for (
int l=0; l<nz; l++) {
685 cout << diff_shift_y(l) <<
" " ;
689 diff_sy = diff_shift_y(1) ;
690 for (
int l=2; l<nz; l++) {
691 diff_sy += diff_shift_y(l) ;
706 cout <<
"Relative difference in the shift vector (z) : " << endl ;
707 for (
int l=0; l<nz; l++) {
708 cout << diff_shift_z(l) <<
" " ;
712 diff_sz = diff_shift_z(1) ;
713 for (
int l=2; l<nz; l++) {
714 diff_sz += diff_shift_z(l) ;
721 cout <<
"Mass_bh : " <<
mass_bh / msol <<
" [M_sol]" << endl ;
722 double rad_apphor =
rad_ah() ;
723 cout <<
" : " << 0.5 * rad_apphor / ggrav / msol
724 <<
" [M_sol]" << endl ;
729 cout <<
"Mass_bh : " <<
mass_bh / msol
730 <<
" [M_sol]" << endl ;
737 double irr_gm, adm_gm, kom_gm ;
741 cout <<
"Irreducible mass : " <<
mass_irr() / msol << endl ;
742 cout <<
"Gravitaitonal mass : " <<
mass_bh / msol << endl ;
743 cout <<
"ADM mass : " <<
mass_adm() / msol << endl ;
744 cout <<
"Komar mass : " <<
mass_kom() / msol << endl ;
751 cout <<
"Diff. (Mirr-Mg)/Mg : " << irr_gm << endl ;
752 cout <<
"Diff. (Madm-Mg)/Mg : " << adm_gm << endl ;
753 cout <<
"Diff. (Mkom-Mg)/Mg : " << kom_gm << endl ;
793 lapconf_exact = 1./
sqrt(1.+2.*mass/rr) ;
797 for (
int i=1; i<=3; i++)
799 2.*mass*lapconf_exact%lapconf_exact%ll(i)/rr ;
805 rare =
r_coord(neumann, first) ;
808 lapconf_exact =
sqrt(1. - 2.*mass/rare/rr
809 + cc*cc*
pow(mass/rare/rr,4.)) *
sqrt(rare) ;
811 confo_exact =
sqrt(rare) ;
813 for (
int i=1; i<=3; i++) {
814 shift_exact.set(i) = mass * mass * cc * ll(i) / rr / rr
828 shift_exact.annule_domain(0) ;
829 shift_exact.std_spectral_base() ;
830 for (
int i=1; i<=3; i++)
831 shift_exact.set(i).raccord(1) ;
880 diff_lapconf_exact.
set(0) = 0. ;
881 cout <<
"Relative difference in the lapconf function : " << endl ;
882 for (
int l=0; l<nz; l++) {
883 cout << diff_lapconf_exact(l) <<
" " ;
889 diff_confo_exact.
set(0) = 0. ;
890 cout <<
"Relative difference in the conformal factor : " << endl ;
891 for (
int l=0; l<nz; l++) {
892 cout << diff_confo_exact(l) <<
" " ;
901 cout <<
"Relative difference in the shift vector (x) : " << endl ;
902 for (
int l=0; l<nz; l++) {
903 cout << diff_shift_exact_x(l) <<
" " ;
906 cout <<
"Relative difference in the shift vector (y) : " << endl ;
907 for (
int l=0; l<nz; l++) {
908 cout << diff_shift_exact_y(l) <<
" " ;
911 cout <<
"Relative difference in the shift vector (z) : " << endl ;
912 for (
int l=0; l<nz; l++) {
913 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