28char tslice_conf_init_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_conf_init.C,v 1.13 2014/10/13 08:53:48 j_novak Exp $" ;
85#include "time_slice.h"
88#include "utilitaires.h"
91void Time_slice_conf::initial_data_cts(
const Sym_tensor& uu,
92 const Scalar& trk_in,
const Scalar& trk_point,
93 double pdt,
double precis,
int method_poisson_vect,
94 const char* graph_device,
const Scalar* p_ener_dens,
95 const Vector* p_mom_dens,
const Scalar* p_trace_stress) {
101 double tr_uu =
max(
maxabs(uu.trace(
tgam()),
"trace tgam_{ij} u^{ij}")) ;
104 "Time_slice_conf::initial_data_cts : the trace of u^{ij} with respect\n"
105 <<
" to the conformal metric tgam_{ij} is not zero !\n"
106 <<
" error = " << tr_uu << endl ;
110 assert(trk_in.check_dzpuis(2)) ;
111 assert(trk_point.check_dzpuis(4)) ;
125 const Map& map = uu.get_mp() ;
126 const Base_vect& triad = *(uu.get_triad()) ;
130 int nz = map.get_mg()->get_nzone() ;
131 double ray_des = 1.25 * map.val_r(nz-2, 1., 0., 0.) ;
134 Scalar ener_dens(map) ;
135 if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
136 else ener_dens.set_etat_zero() ;
138 Vector mom_dens(map, CON, triad) ;
139 if (p_mom_dens != 0x0) mom_dens = *(p_mom_dens) ;
140 else mom_dens.set_etat_zero() ;
142 Scalar trace_stress(map) ;
143 if (p_trace_stress != 0x0) trace_stress = *(p_trace_stress) ;
144 else trace_stress.set_etat_zero() ;
147 Scalar source_psi(map) ;
148 Scalar source_nn(map) ;
149 Vector source_beta(map, CON, triad) ;
154 for (
int i=0; i<imax; i++) {
160 const Vector& dpsi =
psi().derive_cov(
ff) ;
161 const Vector& dln_psi =
ln_psi().derive_cov(
ff) ;
162 const Vector& dnn =
nn().derive_cov(
ff) ;
164 Sym_tensor taa =
aa().up_down(
tgam()) ;
165 Scalar aa_quad =
contract(taa, 0, 1,
aa(), 0, 1) ;
169 tmp = 0.125*
psi() *
tgam().ricci_scal()
175 source_psi = tmp -
psi()*
psi4()* ( 0.5*qpig* ener_dens
177 - 8.33333333333333e-2*
trk()*
trk() ) ;
182 source_nn =
psi4()*(
nn()*( qpig* (ener_dens + trace_stress) + aa_quad
183 - 0.3333333333333333*
trk()*
trk() )
200 dnn - 6.*
nn() * dln_psi, 0) ;
202 source_beta += 2.*
nn() * ( 2.*qpig*
psi4() * mom_dens
203 + 0.66666666666666666*
trk().derive_con(
tgam())
208 beta().derive_cov(
ff).derive_cov(
ff), 1, 2)
209 + 0.3333333333333333*
212 + uu.divergence(
ff) ;
215 source_beta -= vtmp ;
217 source_beta += 0.66666666666666666*
beta().divergence(
ff) *
hdirac() ;
227 Scalar psi_jp1 = source_psi.poisson() + 1. ;
229 if (psi_jp1.get_etat() == ETATUN) psi_jp1.std_spectral_base() ;
232 maxabs(psi_jp1.laplacian() - source_psi,
233 "Absolute error in the resolution of the equation for Psi") ;
235 des_meridian(psi_jp1, 0., ray_des,
"Psi", ngraph0, graph_device) ;
240 Scalar nn_jp1 = source_nn.poisson() + 1. ;
242 if (nn_jp1.get_etat() == ETATUN) nn_jp1.std_spectral_base() ;
245 maxabs(nn_jp1.laplacian() - source_nn,
246 "Absolute error in the resolution of the equation for N") ;
248 des_meridian(nn_jp1, 0., ray_des,
"N", ngraph0+1, graph_device) ;
253 Vector beta_jp1 = source_beta.poisson(0.3333333333333333,
ff,
254 method_poisson_vect) ;
256 des_meridian(beta_jp1(1), 0., ray_des,
"\\gb\\ur\\d", ngraph0+2,
258 des_meridian(beta_jp1(2), 0., ray_des,
"\\gb\\u\\gh\\d", ngraph0+3,
260 des_meridian(beta_jp1(3), 0., ray_des,
"\\gb\\u\\gf\\d", ngraph0+4,
264 Vector test_beta = (beta_jp1.derive_con(
ff)).divergence(
ff)
265 + 0.3333333333333333 * (beta_jp1.divergence(
ff)).derive_con(
ff) ;
266 test_beta.inc_dzpuis() ;
267 maxabs(test_beta - source_beta,
268 "Absolute error in the resolution for beta") ;
278 cout <<
"step = " << i <<
" : diff_psi = " << diff_psi
279 <<
" diff_nn = " << diff_nn
280 <<
" diff_beta = " << diff_beta << endl ;
281 if ( (diff_psi < precis) && (diff_nn < precis) && (diff_beta < precis) )
295 Sym_tensor aa_jp1 = (
beta().ope_killing_conf(
tgam()) + uu )
313 double ttime1 = ttime ;
315 for (
int j=1; j <
depth; j++) {
333 Sym_tensor uu0 = uu ;
336 for (
int j=1; j <
depth; j++) {
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
void downdate(int j)
Suppresses a stored value.
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual void set_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
virtual const Scalar & A_hata() const
Returns the potential A of .
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual const Scalar & B_hata() const
Returns the potential of .
virtual void del_deriv() const
Deletes all the derived quantities.
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: .
const Scalar & psi4() const
Factor at the current time step (jtime ).
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
int jtime
Time step index of the latest slice.
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Evolution_std< Sym_tensor > k_dd_evol
Values at successive time steps of the covariant components of the extrinsic curvature tensor .
virtual const Vector & beta() const
shift vector at the current time step (jtime )
int depth
Number of stored time slices.
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Evolution_std< double > the_time
Time label of each slice.
Evolution_std< Sym_tensor > k_uu_evol
Values at successive time steps of the contravariant components of the extrinsic curvature tensor .
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Standard units of space, time and mass.