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) {
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 ;
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)
470 cout <<
"!!!!! WARNING: Not yet available !!!!!" << endl ;
492 for (
int i=1; i<=3; i++) {
493 if (source_shift(i).get_etat() != ETATZERO)
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 ;
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
827 for (
int i=1; i<=3; i++)
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) <<
" " ;
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Scalar poisson_neumann(const Valeur &, int) const
Idem as Scalar::poisson_dirichlet , the boundary condition being on the radial derivative of the solu...
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
double & set(int i)
Read/write of a particular element (index i) (1D case)
virtual double mass_irr() const
Irreducible mass of the black hole.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Values and coefficients of a (real-value) function.
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...
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.
Scalar lapconf_bh
Part of lapconf from the analytic background.
Vector shift_rs
Part of the shift vector from the numerical computation.
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 \fx\f direct...
Vector shift_bh
Part of the shift vector from the analytic background.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
virtual double mass_kom() const
Komar mass.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
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 \fy\f direct...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Scalar lapse
Lapse function generated by the black hole.
const Scalar & deriv(int i) const
Returns of *this , where .
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
double mass_bh
Gravitational mass of BH.
Coord r
r coordinate centered on the grid
Map & mp
Mapping associated with the black hole.
Sym_tensor taij
Trace of the extrinsic curvature.
Tensor field of valence 0 (or component of a tensorial field).
Cmp pow(const Cmp &, int)
Power .
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Standard units of space, time and mass.
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
Scalar & set(int)
Read/write access to a component.
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fz\f direct...
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
virtual double mass_adm() const
ADM mass.
int get_nzone() const
Returns the number of domains.
Scalar taij_quad
Part of the scalar generated by .
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Cmp sqrt(const Cmp &)
Square root.
void set_dzpuis(int)
Modifies the dzpuis flag.
Tensor field of valence 1.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Scalar lapconf_rs
Part of lapconf from the numerical computation.
virtual double rad_ah() const
Radius of the apparent horizon.
virtual void del_deriv() const
Deletes all the derived quantities.
int get_dzpuis() const
Returns dzpuis.
Scalar confo
Conformal factor 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.
Vector shift
Shift vector generated by the black hole.
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
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.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.