162 #include "time_slice.h" 163 #include "isol_hor.h" 165 #include "evolution.h" 167 #include "graphique.h" 168 #include "utilitaires.h" 192 for (
int k=0 ; k<nnp ; k++)
193 for (
int j=0 ; j<nnt ; j++)
196 psi_bound.std_base_scal() ;
212 tmp = tmp /
beta()(1) ;
223 for (
int k=0 ; k<nnp ; k++)
224 for (
int j=0 ; j<nnt ; j++)
227 psi_bound.std_base_scal() ;
260 for (
int k=0 ; k<nnp ; k++)
261 for (
int j=0 ; j<nnt ; j++)
264 psi_bound.std_base_scal() ;
290 for (
int k=0 ; k<nnp ; k++)
291 for (
int j=0 ; j<nnt ; j++)
294 psi_bound.std_base_scal() ;
315 for (
int k=0 ; k<nnp ; k++)
316 for (
int j=0 ; j<nnt ; j++)
321 psi_bound.std_base_scal() ;
348 tmp = (tmp + rho *
psi())/(1 + rho) ;
361 for (
int k=0 ; k<nnp ; k++)
362 for (
int j=0 ; j<nnt ; j++)
365 psi_bound.std_base_scal() ;
394 tmp = tmp / (kk_rr + rho) - 1;
397 cout <<
"k_kerr = " << k_kerr.val_grid_point(1, 0, 0, 0) << endl ;
412 for (
int k=0 ; k<nnp ; k++)
413 for (
int j=0 ; j<nnt ; j++)
416 nn_bound.std_base_scal() ;
457 Scalar tmp =
psi() *
psi() * ( k_kerr + naass + 1.* traceK)
469 cout <<
"kappa_pres = " << k_kerr.
val_grid_point(1, 0, 0, 0) << endl ;
470 cout <<
"kappa_hor = " << k_hor.
val_grid_point(1, 0, 0, 0) << endl ;
471 cout <<
"sDN = " << sdN.val_grid_point(1, 0, 0, 0) << endl ;
482 for (
int k=0 ; k<nnp ; k++)
483 for (
int j=0 ; j<nnt ; j++)
486 nn_bound.std_base_scal() ;
518 Vector hom = hom1 + hom2 ;
520 Scalar div_omega = 1.*
contract( qq_uu, 0, 1, (1.*hom1 + 1.*hom2).derive_cov(
gam()), 0, 1) ;
533 Ricci_conf = 2 / rr / rr ;
539 Ricci =
pow(
psi(), -4) * (Ricci_conf - 4*log_psi.
lapang())/rr /rr ;
561 Scalar source = div_omega +
contract( qq_uu, 0, 1, hom * hom , 0, 1) - 0.5 * Ricci - L_theta;
562 source = source / theta_k ;
565 nn().derive_cov(
gam()), 0))/(1+rho)
577 for (
int k=0 ; k<nnp ; k++)
578 for (
int j=0 ; j<nnt ; j++)
581 nn_bound.std_base_scal() ;
598 tmp += cc * (rho -1)*
nn() ;
599 tmp = tmp / (rho*cc) ;
613 for (
int k=0 ; k<nnp ; k++)
614 for (
int j=0 ; j<nnt ; j++)
617 nn_bound.std_base_scal() ;
639 for (
int k=0 ; k<nnp ; k++)
640 for (
int j=0 ; j<nnt ; j++)
643 nn_bound.std_base_scal() ;
675 tmp = (tmp + rho *
nn())/(1 + rho) ;
686 for (
int k=0 ; k<nnp ; k++)
687 for (
int j=0 ; j<nnt ; j++)
690 nn_bound.std_base_scal() ;
711 Scalar det_q = q_dd(2,2) * q_dd(3,3) - q_dd(2,3) * q_dd(3,2) ;
734 Scalar source1 = div_kqs ;
735 source1 *= square_q ;
777 Scalar Ricci_conf = 2 / rr /rr ;
787 div_omega = L_theta -
contract(qqq, 0, 1, hom * hom, 0, 1) + 0.5 * Ricci
796 div_omega.std_spectral_base() ;
797 div_omega.inc_dzpuis(3) ;
801 Scalar source3 = div_omega ;
802 source3 *= square_q ;
808 cout <<
"Constant part of div_omega = " << corr_const <<endl ;
811 corr_div_omega = corr_const ;
815 source3 -= corr_div_omega ;
822 Scalar source = source1 + source2 + source3 ;
829 Scalar source_tmp = source ;
837 lapse = (
exp(logn)*cc) ;
845 err.open (
"err_laplBC.d", ofstream::app ) ;
854 Scalar err_lapl = div_omega_test - div_omega ;
857 double min_err = err_lapl.val_grid_point(1, 0, 0, 0) ;
859 for (
int k=0 ; k<nnp ; k++)
860 for (
int j=0 ; j<nnt ; j++){
861 if (err_lapl.val_grid_point(1, k, j, 0) > max_err)
862 max_err = err_lapl.val_grid_point(1, k, j, 0) ;
863 if (err_lapl.val_grid_point(1, k, j, 0) < min_err)
864 min_err = err_lapl.val_grid_point(1, k, j, 0) ;
867 err << mer <<
" " << max_err <<
" " << min_err << endl ;
885 for (
int k=0 ; k<nnp ; k++)
886 for (
int j=0 ; j<nnt ; j++)
889 nn_bound.std_base_scal() ;
915 for (
int k=0 ; k<nnp ; k++)
916 for (
int j=0 ; j<nnt ; j++)
919 bnd_beta_r.set_base_r(0, base_tmp.
get_base_r(0)) ;
920 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
921 bnd_beta_r.set_base_r(l, base_tmp.
get_base_r(l)) ;
928 cout <<
"norme de lim_vr" << endl <<
norme(bnd_beta_r) << endl ;
929 cout <<
"bases" << endl << bnd_beta_r.base << endl ;
954 bnd_beta_theta = 1. ;
956 for (
int k=0 ; k<nnp ; k++)
957 for (
int j=0 ; j<nnt ; j++)
962 bnd_beta_theta.set_base_r(0, base_tmp.
get_base_r(0)) ;
963 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
964 bnd_beta_theta.set_base_r(l, base_tmp.
get_base_r(l)) ;
966 bnd_beta_theta.set_base_t(base_tmp.
get_base_t(1)) ;
967 bnd_beta_theta.set_base_p(base_tmp.
get_base_p(1)) ;
969 cout <<
"bases" << endl << bnd_beta_theta.base << endl ;
972 return bnd_beta_theta ;
1000 for (
int k=0 ; k<nnp ; k++)
1001 for (
int j=0 ; j<nnt ; j++)
1004 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
1008 bnd_beta_phi.set_base_r(0, base_tmp.get_base_r(0)) ;
1009 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
1010 bnd_beta_phi.set_base_r(l, base_tmp.get_base_r(l)) ;
1011 bnd_beta_phi.set_base_r((*
mp.
get_mg()).get_nzone()-1, base_tmp.get_base_r((*
mp.
get_mg()).get_nzone()-1)) ;
1012 bnd_beta_phi.set_base_t(base_tmp.get_base_t(1)) ;
1013 bnd_beta_phi.set_base_p(base_tmp.get_base_p(1)) ;
1017 return bnd_beta_phi ;
1028 assert ((orientation == 0) || (orientation == M_PI)) ;
1029 int aligne = (orientation == 0) ? 1 : -1 ;
1057 dep_phi = 0.0*sint*cosp ;
1059 for (
int k=0 ; k<nnp ; k++)
1060 for (
int j=0 ; j<nnt ; j++)
1061 lim_x.
set(0, k, j, 0) = aligne * om * ya_mtbl(1, k, j, 0) * (1 +
1063 + tmp_vect(1).val_grid_point(1, k, j, 0) ;
1078 assert ((orientation == 0) || (orientation == M_PI)) ;
1079 int aligne = (orientation == 0) ? 1 : -1 ;
1108 dep_phi = 0.0*sint*cosp ;
1110 for (
int k=0 ; k<nnp ; k++)
1111 for (
int j=0 ; j<nnt ; j++)
1112 lim_y.
set(0, k, j, 0) = - aligne * om * xa_mtbl(1, k, j, 0) *(1 +
1114 + tmp_vect(2).val_grid_point(1, k, j, 0) ;
1138 for (
int k=0 ; k<nnp ; k++)
1139 for (
int j=0 ; j<nnt ; j++)
1140 lim_z.
set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
1193 tmp = tmp / (2 * s_tilde(1) ) ;
1204 for (
int k=0 ; k<nnp ; k++)
1205 for (
int j=0 ; j<nnt ; j++)
1206 b_tilde_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1208 b_tilde_bound.std_base_scal() ;
1210 return b_tilde_bound ;
1230 tmp = tmp / hh_tilde ;
1244 for (
int k=0 ; k<nnp ; k++)
1245 for (
int j=0 ; j<nnt ; j++)
1246 b_tilde_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1248 b_tilde_bound.std_base_scal() ;
1250 return b_tilde_bound ;
1290 dep_phi = 0.0*sint*cosp ;
1294 assert ((orientation == 0) || (orientation == M_PI)) ;
1295 int aligne = (orientation == 0) ? 1 : -1 ;
1303 angular.
set(1) = - aligne * om * yya * (1 + dep_phi) ;
1304 angular.set(2) = aligne * om * xxa * (1 + dep_phi) ;
1308 angular.set(1).set_spectral_va()
1310 angular.set(2).set_spectral_va()
1312 angular.set(3).set_spectral_va()
1332 kss = - 0.15 /
nn() ;
1337 cout <<
"kappa_hor = " <<
kappa_hor() << endl ;
1366 - 3 *
nn() * kss + nn_KK + accel + div_VV ;
1368 b_tilde_new = b_tilde_new / (hh_tilde * rho) ;
1372 tmp_vect.
set(1) = b_tilde_new * s_tilde(1) + angular(1) ;
1373 tmp_vect.
set(2) = b_tilde_new * s_tilde(2) + angular(2) ;
1374 tmp_vect.
set(3) = b_tilde_new * s_tilde(3) + angular(3) ;
1509 for (
int k=0 ; k<nnp ; k++)
1510 for (
int j=0 ; j<nnt ; j++)
1537 for (
int k=0 ; k<nnp ; k++)
1538 for (
int j=0 ; j<nnt ; j++)
1563 for (
int k=0 ; k<nnp ; k++)
1564 for (
int j=0 ; j<nnt ; j++)
1588 for (
int k=0 ; k<nnp ; k++)
1589 for (
int j=0 ; j<nnt ; j++)
1615 for (
int k=0 ; k<nnp ; k++)
1616 for (
int j=0 ; j<nnt ; j++)
1641 for (
int k=0 ; k<nnp ; k++)
1642 for (
int j=0 ; j<nnt ; j++)
Coord xa
Absolute x coordinate.
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Vector radial_vect_hor() const
Vector radial normal.
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution)
Cmp log(const Cmp &)
Neperian logarithm.
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial)
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
Cmp exp(const Cmp &)
Exponential.
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
const Valeur boundary_vv_z_bin(double om, int hole=0) const
Component z of boundary value of .
const Valeur boundary_nn_Dir_lapl(int mer=1) const
Dirichlet boundary condition for N fixing the divergence of the connection form . ...
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif) .
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
double radius
Radius of the horizon in LORENE's units.
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tensor field of valence 0 (or component of a tensorial field).
const Valeur boundary_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
const Vector vv_bound_cart(double om) const
Vector for boundary conditions in cartesian.
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
double kappa_hor() const
Surface gravity.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial)
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Values and coefficients of a (real-value) function.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial)
Tensor field of valence 1.
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
void set_dzpuis(int)
Modifies the dzpuis flag.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
void annule_hard()
Sets the Scalar to zero in a hard way.
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Map_af & mp
Affine mapping.
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif) .
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
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...
double boost_x
Boost velocity in x-direction.
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Base_val base
Bases on which the spectral expansion is performed.
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Metric met_gamt
3 metric tilde
const Valeur boundary_beta_phi(double om) const
Component phi of boundary value of .
const Metric & gam() const
Induced metric at the current time step (jtime )
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
Cmp pow(const Cmp &, int)
Power .
const Valeur boundary_beta_z() const
Component z of boundary value of .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
const Valeur boundary_beta_r() const
Component r of boundary value of .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ) ...
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Coord ya
Absolute y coordinate.
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Bases of the spectral expansions.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
const Valeur boundary_vv_x_bin(double om, int hole=0) const
Component x of boundary value of .
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
const Valeur boundary_b_tilde_Dir() const
Dirichlet boundary condition for b_tilde.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial)
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar & set(int)
Read/write access to a component.
const Valeur boundary_beta_theta() const
Component theta of boundary value of .
const Vector tradial_vect_hor() const
Vector radial normal tilde.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
const Valeur boundary_vv_y_bin(double om, int hole=0) const
Component y of boundary value of .
const Vector vv_bound_cart_bin(double om, int hole=0) const
Vector for boundary conditions in cartesian for binary systems.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Valeur boundary_b_tilde_Neu() const
Neumann boundary condition for b_tilde.
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
Class intended to describe valence-2 symmetric tensors.
double boost_z
Boost velocity in z-direction.
const Valeur & get_spectral_va() const
Returns va (read only version)
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook's boundary condition.
Coord r
r coordinate centered on the grid
virtual const Sym_tensor & aa_auto() const
Conformal representation of the traceless part of the extrinsic curvature: Returns the value at the ...
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution)
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.