29char bound_hor_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Isol_hor/bound_hor.C,v 1.36 2014/10/13 08:53:00 j_novak Exp $" ;
159#include "time_slice.h"
162#include "evolution.h"
164#include "graphique.h"
165#include "utilitaires.h"
189 for (
int k=0 ; k<nnp ; k++)
190 for (
int j=0 ; j<nnt ; j++)
209 tmp = tmp /
beta()(1) ;
220 for (
int k=0 ; k<nnp ; k++)
221 for (
int j=0 ; j<nnt ; j++)
257 for (
int k=0 ; k<nnp ; k++)
258 for (
int j=0 ; j<nnt ; j++)
287 for (
int k=0 ; k<nnp ; k++)
288 for (
int j=0 ; j<nnt ; j++)
312 for (
int k=0 ; k<nnp ; k++)
313 for (
int j=0 ; j<nnt ; j++)
345 tmp = (tmp + rho *
psi())/(1 + rho) ;
358 for (
int k=0 ; k<nnp ; k++)
359 for (
int j=0 ; j<nnt ; j++)
391 tmp = tmp / (kk_rr + rho) - 1;
394 cout <<
"k_kerr = " << k_kerr.
val_grid_point(1, 0, 0, 0) << endl ;
409 for (
int k=0 ; k<nnp ; k++)
410 for (
int j=0 ; j<nnt ; j++)
454 Scalar tmp =
psi() *
psi() * ( k_kerr + naass + 1.* traceK)
466 cout <<
"kappa_pres = " << k_kerr.
val_grid_point(1, 0, 0, 0) << endl ;
467 cout <<
"kappa_hor = " << k_hor.
val_grid_point(1, 0, 0, 0) << endl ;
479 for (
int k=0 ; k<nnp ; k++)
480 for (
int j=0 ; j<nnt ; j++)
515 Vector hom = hom1 + hom2 ;
517 Scalar div_omega = 1.*
contract( qq_uu, 0, 1, (1.*hom1 + 1.*hom2).derive_cov(
gam()), 0, 1) ;
530 Ricci_conf = 2 / rr / rr ;
536 Ricci =
pow(
psi(), -4) * (Ricci_conf - 4*log_psi.
lapang())/rr /rr ;
558 Scalar source = div_omega +
contract( qq_uu, 0, 1, hom * hom , 0, 1) - 0.5 * Ricci - L_theta;
559 source = source / theta_k ;
562 nn().derive_cov(
gam()), 0))/(1+rho)
574 for (
int k=0 ; k<nnp ; k++)
575 for (
int j=0 ; j<nnt ; j++)
595 tmp += cc * (rho -1)*
nn() ;
596 tmp = tmp / (rho*cc) ;
610 for (
int k=0 ; k<nnp ; k++)
611 for (
int j=0 ; j<nnt ; j++)
636 for (
int k=0 ; k<nnp ; k++)
637 for (
int j=0 ; j<nnt ; j++)
672 tmp = (tmp + rho *
nn())/(1 + rho) ;
683 for (
int k=0 ; k<nnp ; k++)
684 for (
int j=0 ; j<nnt ; j++)
708 Scalar det_q = q_dd(2,2) * q_dd(3,3) - q_dd(2,3) * q_dd(3,2) ;
731 Scalar source1 = div_kqs ;
732 source1 *= square_q ;
774 Scalar Ricci_conf = 2 / rr /rr ;
784 div_omega = L_theta -
contract(qqq, 0, 1, hom * hom, 0, 1) + 0.5 * Ricci
798 Scalar source3 = div_omega ;
799 source3 *= square_q ;
805 cout <<
"Constant part of div_omega = " << corr_const <<endl ;
808 corr_div_omega = corr_const ;
812 source3 -= corr_div_omega ;
819 Scalar source = source1 + source2 + source3 ;
826 Scalar source_tmp = source ;
834 lapse = (
exp(logn)*cc) ;
842 err.open (
"err_laplBC.d", ofstream::app ) ;
851 Scalar err_lapl = div_omega_test - div_omega ;
856 for (
int k=0 ; k<nnp ; k++)
857 for (
int j=0 ; j<nnt ; j++){
864 err << mer <<
" " << max_err <<
" " << min_err << endl ;
882 for (
int k=0 ; k<nnp ; k++)
883 for (
int j=0 ; j<nnt ; j++)
912 for (
int k=0 ; k<nnp ; k++)
913 for (
int j=0 ; j<nnt ; j++)
917 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
925 cout <<
"norme de lim_vr" << endl <<
norme(bnd_beta_r) << endl ;
926 cout <<
"bases" << endl << bnd_beta_r.
base << endl ;
951 bnd_beta_theta = 1. ;
953 for (
int k=0 ; k<nnp ; k++)
954 for (
int j=0 ; j<nnt ; j++)
960 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
966 cout <<
"bases" << endl << bnd_beta_theta.
base << endl ;
969 return bnd_beta_theta ;
997 for (
int k=0 ; k<nnp ; k++)
998 for (
int j=0 ; j<nnt ; j++)
1001 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
1006 for (
int l=0 ; l<(*
mp.
get_mg()).get_nzone()-1 ; l++)
1014 return bnd_beta_phi ;
1025 assert ((orientation == 0) || (orientation == M_PI)) ;
1026 int aligne = (orientation == 0) ? 1 : -1 ;
1054 dep_phi = 0.0*sint*cosp ;
1056 for (
int k=0 ; k<nnp ; k++)
1057 for (
int j=0 ; j<nnt ; j++)
1058 lim_x.
set(0, k, j, 0) = aligne * om * ya_mtbl(1, k, j, 0) * (1 +
1060 + tmp_vect(1).val_grid_point(1, k, j, 0) ;
1075 assert ((orientation == 0) || (orientation == M_PI)) ;
1076 int aligne = (orientation == 0) ? 1 : -1 ;
1105 dep_phi = 0.0*sint*cosp ;
1107 for (
int k=0 ; k<nnp ; k++)
1108 for (
int j=0 ; j<nnt ; j++)
1109 lim_y.
set(0, k, j, 0) = - aligne * om * xa_mtbl(1, k, j, 0) *(1 +
1111 + tmp_vect(2).val_grid_point(1, k, j, 0) ;
1135 for (
int k=0 ; k<nnp ; k++)
1136 for (
int j=0 ; j<nnt ; j++)
1137 lim_z.
set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
1190 tmp = tmp / (2 * s_tilde(1) ) ;
1201 for (
int k=0 ; k<nnp ; k++)
1202 for (
int j=0 ; j<nnt ; j++)
1207 return b_tilde_bound ;
1227 tmp = tmp / hh_tilde ;
1241 for (
int k=0 ; k<nnp ; k++)
1242 for (
int j=0 ; j<nnt ; j++)
1247 return b_tilde_bound ;
1287 dep_phi = 0.0*sint*cosp ;
1291 assert ((orientation == 0) || (orientation == M_PI)) ;
1292 int aligne = (orientation == 0) ? 1 : -1 ;
1300 angular.
set(1) = - aligne * om * yya * (1 + dep_phi) ;
1301 angular.
set(2) = aligne * om * xxa * (1 + dep_phi) ;
1329 kss = - 0.15 /
nn() ;
1334 cout <<
"kappa_hor = " <<
kappa_hor() << endl ;
1363 - 3 *
nn() * kss + nn_KK + accel + div_VV ;
1365 b_tilde_new = b_tilde_new / (hh_tilde * rho) ;
1369 tmp_vect.
set(1) = b_tilde_new * s_tilde(1) + angular(1) ;
1370 tmp_vect.
set(2) = b_tilde_new * s_tilde(2) + angular(2) ;
1371 tmp_vect.
set(3) = b_tilde_new * s_tilde(3) + angular(3) ;
1506 for (
int k=0 ; k<nnp ; k++)
1507 for (
int j=0 ; j<nnt ; j++)
1534 for (
int k=0 ; k<nnp ; k++)
1535 for (
int j=0 ; j<nnt ; j++)
1560 for (
int k=0 ; k<nnp ; k++)
1561 for (
int j=0 ; j<nnt ; j++)
1585 for (
int k=0 ; k<nnp ; k++)
1586 for (
int j=0 ; j<nnt ; j++)
1612 for (
int k=0 ; k<nnp ; k++)
1613 for (
int j=0 ; j<nnt ; j++)
1638 for (
int k=0 ; k<nnp ; k++)
1639 for (
int j=0 ; j<nnt ; j++)
Bases of the spectral expansions.
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
const Valeur boundary_vv_z_bin(double om, int hole=0) const
Component z of boundary value of .
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution)
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial)
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif)
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial)
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
Metric met_gamt
3 metric tilde
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook's boundary condition.
const Valeur boundary_b_tilde_Neu() const
Neumann boundary condition for b_tilde.
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
const Valeur boundary_beta_r() const
Component r of boundary value of .
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
const Valeur boundary_vv_y_bin(double om, int hole=0) const
Component y 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 .
double boost_z
Boost velocity in z-direction.
const Valeur boundary_beta_theta() const
Component theta of boundary value of .
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
const Valeur boundary_vv_x_bin(double om, int hole=0) const
Component x of boundary value of .
const Vector vv_bound_cart(double om) const
Vector for boundary conditions in cartesian
double kappa_hor() const
Surface gravity
const Valeur boundary_beta_z() const
Component z of boundary value of .
double radius
Radius of the horizon in LORENE's units.
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_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
const Vector tradial_vect_hor() const
Vector radial normal tilde.
Map_af & mp
Affine mapping.
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial)
const Valeur boundary_beta_phi(double om) const
Component phi of boundary value of .
const Vector radial_vect_hor() const
Vector radial normal.
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif) .
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
double boost_x
Boost velocity in x-direction.
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial)
const Valeur boundary_b_tilde_Dir() const
Dirichlet boundary condition for b_tilde.
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution)
const Vector vv_bound_cart_bin(double om, int hole=0) const
Vector for boundary conditions in cartesian for binary systems.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord ya
Absolute y coordinate.
Coord r
r coordinate centered on the grid
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Coord xa
Absolute x coordinate.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Tensor field of valence 0 (or component of a tensorial field).
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
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...
Valeur & set_spectral_va()
Returns va (read/write version)
const Valeur & get_spectral_va() const
Returns va (read only version)
void annule_hard()
Sets the Scalar to zero in a hard way.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
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...
Class intended to describe valence-2 symmetric tensors.
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Metric & gam() const
Induced metric at the current time step (jtime )
Values and coefficients of a (real-value) function.
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Base_val base
Bases on which the spectral expansion is performed.
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
void set_base_t(int base_t)
Sets the expansion basis for functions in all domains.
void set_base_p(int base_p)
Sets the expansion basis for functions in all domains.
Tensor field of valence 1.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Cmp exp(const Cmp &)
Exponential.
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .