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