28char tslice_dirac_max_evolve_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_evolve.C,v 1.22 2015/08/10 15:32:27 j_novak Exp $" ;
111#include "time_slice.h"
113#include "graphique.h"
114#include "utilitaires.h"
118const Tbl& monitor_scalar(
const Scalar& uu, Tbl& resu) ;
131 const Map& map =
nn().get_mp() ;
138 int nz = map.
get_mg()->get_nzone() ;
140 int nt = map.
get_mg()->get_nt(0) ;
141 int np = map.
get_mg()->get_np(0) ;
149 for (
int lz=1;
lz<nz;
lz++) {
169 for (
int i=1;
i<=3;
i++)
170 for (
int j=
i;
j<=3;
j++)
179 for (
int i=1;
i<=3;
i++)
180 for (
int j=
i;
j<=3;
j++)
194 tmp.mult_r() ;
tmp.mult_r() ;
198 tmp.mult_r() ;
tmp.mult_r() ;
203 double Rmax = map.
val_r(nz-2, 1., 0., 0.) ;
208 double an = 23./12. ;
209 double anm1 = -4./3. ;
210 double anm2 = 5./12. ;
247 xij_b.set_etat_qcq() ;
254 xij_a.set_etat_qcq() ;
399 "==============================================================\n"
402 <<
", Log of central lapse: " <<
log(
nn().val_grid_point(0,0,0,0)) <<
endl
403 <<
"==============================================================\n" ;
447 for (
int l=1;
l<nz-1;
l++) {
453 "Check of Hamiltonian constraint",
458 for (
int l=1;
l<nz-1;
l++) {
464 "Check of momentum constraint (r comp.)",
ngraph0_mon+2,
468 for (
int l=1;
l<nz-1;
l++) {
474 "Check of momentum constraint (\\gh comp.)",
ngraph0_mon+3,
478 for (
int l=1;
l<nz-1;
l++) {
484 "Check of momentum constraint (\\gf comp.)",
ngraph0_mon+4,
490 for (
int i=1;
i<=3;
i++)
491 for(
int j=1;
j<=
i;
j++) {
493 for (
int l=1;
l<nz-1;
l++) {
546 bc_A.set_spectral_va().ylm() ;
552 bc_B.set_spectral_va().ylm() ;
561 if (
sbcmu.get_etat() == ETATZERO) {
562 sbcmu.annule_hard() ;
565 sbcmu.set_spectral_va().ylm() ;
567 if (
tmp.get_etat() == ETATZERO) {
571 tmp.set_spectral_va().ylm() ;
578 if (
sbckhi.get_etat() == ETATZERO) {
582 sbckhi.set_spectral_va().ylm() ;
584 if (
tmp.get_etat() == ETATZERO) {
588 tmp.set_spectral_va().ylm() ;
595 if (
sbcmu.get_etat() == ETATZERO) {
596 sbcmu.annule_hard() ;
599 sbcmu.set_spectral_va().ylm() ;
601 if (
tmp.get_etat() == ETATZERO) {
605 tmp.set_spectral_va().ylm() ;
612 if (
sbckhi.get_etat() == ETATZERO) {
616 sbckhi.set_spectral_va().ylm() ;
618 if (
tmp.get_etat() == ETATZERO) {
622 tmp.set_spectral_va().ylm() ;
639 for (
int i=1;
i<=3;
i++)
640 for (
int j=
i;
j<=3;
j++)
644 for (
int i=1;
i<=3;
i++)
645 for (
int j=
i;
j<=3;
j++) {
660 tmp.mult_r() ;
tmp.mult_r() ;
664 tmp.mult_r() ;
tmp.mult_r() ;
710 tmp.set_spectral_va().ylm_i() ;
714 tmp.set_spectral_va().ylm_i() ;
749 resu.set_etat_qcq() ;
751 resu.set(0) =
uu.val_grid_point(0,0,0,0) ;
755 const Mg3d& mg = *(
uu.get_mp().get_mg()) ;
763 resu.set(3) =
uu.val_grid_point(
nzm1, 0, 0, nr-1) ;
764 resu.set(4) =
uu.val_grid_point(
nzm1, 0, nt-1, nr-1) ;
765 resu.set(5) =
uu.val_grid_point(
nzm1, np/2, nt-1, nr-1) ;
Bases of the spectral expansions.
Vectorial bases (triads) with respect to which the tensorial components are defined.
Time evolution with partial storage (*** under development ***).
virtual void update(const TyT &new_value, int j, double time_j)
Sets a new value at a given time step.
void save(const char *filename) const
Saves *this in a formatted file.
TyT time_derive(int j, int n=2) const
Computes the time derivative at time step j by means of a n-th order scheme, from the values at steps...
bool is_known(int j) const
Checks whether the value a given time step has been set.
Base class for coordinate mappings.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
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.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tensor field of valence 0 (or component of a tensorial field).
Transverse symmetric tensors of rank 2.
Transverse and traceless symmetric tensors of rank 2.
Class intended to describe valence-2 symmetric tensors.
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
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 void set_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
virtual const Vector & vec_X(int method_poisson=6) const
Vector representing the longitudinal part of .
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual void del_deriv() const
Deletes all the derived quantities.
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
virtual void set_npsi_del_n(const Scalar &npsi_in)
Sets the factor at the current time step (jtime ) and deletes the value of N.
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
void compute_X_from_momentum_constraint(const Vector &hat_S, const Sym_tensor_tt &hata_tt, int iter_max=200, double precis=1.e-12, double relax=0.8, int methode_poisson=6)
Computes the vector from the conformally-rescaled momentum , using the momentum constraint.
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
int jtime
Time step index of the latest slice.
void save(const char *rootname) const
Saves in a binary file.
Tbl check_hamiltonian_constraint(const Scalar *energy_density=0x0, ostream &ost=cout, bool verb=true) const
Checks the level at which the hamiltonian constraint is verified.
virtual const Vector & beta() const
shift vector at the current time step (jtime )
int depth
Number of stored time slices.
Evolution_std< double > the_time
Time label of each slice.
Tbl check_dynamical_equations(const Sym_tensor *strain_tensor=0x0, const Scalar *energy_density=0x0, ostream &ost=cout, bool verb=true) const
Checks the level at which the dynamical equations are verified.
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Tbl check_momentum_constraint(const Vector *momentum_density=0x0, ostream &ost=cout, bool verb=true) const
Checks the level at which the momentum constraints are verified.
virtual const Scalar & trh() const
Computes the trace h, with respect to the flat metric ff , of .
virtual double adm_mass() const
Returns the ADM mass at (geometrical units) the current step.
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $\Psi$ (Hamiltonian constraint).
Evolution_std< Scalar > source_B_hh_evol
The potential of the source of equation for .
Evolution_std< Scalar > source_A_hh_evol
The A potential of the source of equation for .
void evolve(double pdt, int nb_time_steps, int niter_elliptic, double relax_elliptic, int check_mod, int save_mod, int method_poisson_vect=6, int nopause=1, const char *graph_device=0x0, bool verbose=true, const Scalar *ener_euler=0x0, const Vector *mom_euler=0x0, const Scalar *s_euler=0x0, const Sym_tensor *strain_euler=0x0)
Time evolution by resolution of Einstein equations.
virtual const Scalar & B_hh() const
Returns the potential of .
Evolution_std< Scalar > A_hh_evol
The A potential of .
void hh_det_one(int j, Param *par_bc=0x0, Param *par_mat=0x0) const
Computes from the values of A and and using the condition , which fixes the trace of .
Evolution_std< Scalar > source_B_hata_evol
The potential of the source of equation for .
void compute_sources(const Sym_tensor *strain_tensor=0x0) const
Computes the sources source_A_XXX_evol and source_B_XXX_evol , for the solution of the evolution equa...
virtual Scalar solve_npsi(const Scalar *ener_dens=0x0, const Scalar *trace_stress=0x0) const
Solves the elliptic equation for (maximal slicing condition + Hamiltonian constraint)
virtual const Scalar & A_hh() const
Returns the potential A of .
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
virtual void set_AB_hh(const Scalar &A_in, const Scalar &B_in)
Sets the potentials A and of the TT part of (see the documentation of Sym_tensor for details).
Evolution_std< Scalar > B_hh_evol
The potential of .
Evolution_std< Scalar > source_A_hata_evol
The potential A of the source of equation for .
Tensor field of valence 1.
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp log(const Cmp &)
Neperian logarithm.
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 .
void des_evol(const Evolution< double > &uu, const char *nomy=0x0, const char *title=0x0, int ngraph=0, const char *device=0x0, bool closeit=false, bool show_time=true, const char *nomx=0x0)
Plots the variation of some quantity against time.
void arrete(int a=0)
Setting a stop point in a code.
Tbl maxabs_all_domains(const Tensor &aa, int l_excluded=-1, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maximum of the absolute value of each component of a tensor over all the domains.
Tbl central_value(const Tensor &aa, const char *comment=0x0, ostream &ost=cout)
Central value of each component of a tensor.