LORENE
tslice_check_einstein.C
1 /*
2  * Methods of class Time_slice to check Einstein equation solutions.
3  *
4  * (see file time_slice.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Eric Gourgoulhon, Jose Luis Jaramillo & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char tslice_check_einstein_C[] = "$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_check_einstein.C,v 1.10 2014/10/13 08:53:47 j_novak Exp $" ;
29 
30 /*
31  * $Id: tslice_check_einstein.C,v 1.10 2014/10/13 08:53:47 j_novak Exp $
32  * $Log: tslice_check_einstein.C,v $
33  * Revision 1.10 2014/10/13 08:53:47 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.9 2014/10/06 15:13:22 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.8 2010/10/20 07:58:09 j_novak
40  * Better implementation of the explicit time-integration. Not fully-tested yet.
41  *
42  * Revision 1.7 2007/11/06 12:10:56 j_novak
43  * Corrected some mistakes.
44  *
45  * Revision 1.6 2007/11/06 11:53:34 j_novak
46  * Use of contravariant version of the 3+1 equations.
47  *
48  * Revision 1.5 2007/06/28 14:40:36 j_novak
49  * Dynamical check: the fields in the last domain are set to zero to avoid dzpuis problems
50  *
51  * Revision 1.4 2004/05/12 15:24:20 e_gourgoulhon
52  * Reorganized the #include 's, taking into account that
53  * time_slice.h contains now an #include "metric.h".
54  *
55  * Revision 1.3 2004/04/07 07:58:21 e_gourgoulhon
56  * Constructor as Minkowski slice: added call to std_spectral_base()
57  * after setting the lapse to 1.
58  *
59  * Revision 1.2 2004/04/05 12:38:45 j_novak
60  * Minor modifs to prevent some warnings.
61  *
62  * Revision 1.1 2004/04/05 11:54:20 j_novak
63  * First operational (but not tested!) version of checks of Eintein equations.
64  *
65  *
66  * $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_check_einstein.C,v 1.10 2014/10/13 08:53:47 j_novak Exp $
67  *
68  */
69 
70 // C headers
71 #include <cstdlib>
72 #include <cassert>
73 
74 // Lorene headers
75 #include "time_slice.h"
76 #include "unites.h"
77 
78 namespace Lorene {
80  ostream& ost, bool verb) const {
81  using namespace Unites ;
82 
83  bool vacuum = ( energy_density == 0x0 ) ;
84 
85  Scalar field = trk() * trk() - contract( k_uu(), 0, 1, k_dd(), 0, 1 ) ;
86 
87  field.dec_dzpuis() ; // dzpuis: 4 -> 3
88 
89  field += gam().ricci_scal() ;
90 
91  const Scalar* matter ;
92  if (vacuum)
93  matter = new Scalar (0*field) ;
94  else
95  matter = energy_density ;
96 
97  if (verb) ost << endl ;
98  const char* comment = 0x0 ;
99  if (verb) comment = "Check of the Hamiltonian constraint" ;
100  Tbl resu = maxabs(field - (4*qpig) * (*matter), comment, ost,verb ) ;
101  if (verb) ost << endl ;
102 
103  if (vacuum) delete matter ;
104 
105  return resu ;
106 
107 }
108 
110  ostream& ost, bool verb) const {
111  using namespace Unites ;
112 
113  bool vacuum = ( momentum_density == 0x0 ) ;
114 
115  Vector field = k_uu().divergence(gam()) - trk().derive_con(gam()) ;
116 
117 
118  const Vector* matter ;
119  if (vacuum)
120  matter = new Vector (0*field) ;
121  else {
122  assert (momentum_density->get_index_type(0) == CON) ;
123  matter = momentum_density ;
124  }
125 
126  if (verb) ost << endl ;
127  const char* comment = 0x0 ;
128  if (verb) comment = "Check of the momentum constraint" ;
129  Tbl resu = maxabs(field - (2*qpig) * (*matter), comment, ost, verb ) ;
130 
131  if (verb) ost << endl ;
132 
133  if (vacuum) delete matter ;
134 
135  return resu ;
136 
137 }
138 
140  const Scalar* energy_density,
141  ostream& ost, bool verb) const {
142 
143  using namespace Unites ;
144 
145  bool vacuum = ( ( strain_tensor == 0x0 ) && ( energy_density == 0x0 ) ) ;
146 
147  Sym_tensor dyn_lhs = (k_dd_evol.time_derive(jtime, scheme_order)).up_down(gam()) ;
148  int nz = dyn_lhs.get_mp().get_mg()->get_nzone() ;
149  dyn_lhs.annule_domain(nz-1) ;
150  dyn_lhs = dyn_lhs - k_dd().derive_lie(beta()).up_down(gam()) ;
151  dyn_lhs.annule_domain(nz-1) ;
152 
153  const Sym_tensor* matter ;
154  if (vacuum)
155  matter = new Sym_tensor(0 * dyn_lhs) ;
156  else {
157  const Scalar* ener = 0x0 ;
158  const Sym_tensor* sij = 0x0 ;
159  bool new_e = false ;
160  bool new_s = false ;
161  if (energy_density == 0x0) {
162  ener = new Scalar ( 0 * nn() ) ;
163  new_e = true ;
164  }
165  else {
166  ener = energy_density ;
167  if (strain_tensor == 0x0) {
168  sij = new Sym_tensor(0 * dyn_lhs) ;
169  new_s = true ;
170  }
171  else
172  sij = strain_tensor ;
173  }
174  matter = new Sym_tensor( (sij->trace(gam()) - *ener)*gam().con()
175  - 2*(*sij) ) ;
176  if (new_e) delete ener ;
177  if (new_s) delete sij ;
178  }
179 
180  Sym_tensor dyn_rhs = nn()*( (gam().ricci().up_down(gam()) + trk()*k_uu()
181  + qpig * (*matter))) ;
182  dyn_rhs.annule_domain(nz-1) ;
183  dyn_rhs = dyn_rhs - 2*nn()*contract(k_uu(), 1, k_dd().up(0, gam()), 1) ;
184  dyn_rhs.annule_domain(nz-1) ;
185  dyn_rhs = dyn_rhs - nn().derive_con(gam()).derive_con(gam()) ;
186  dyn_rhs.annule_domain(nz-1) ;
187 
188  if (verb) ost << endl ;
189  const char* comment = 0x0 ;
190  if (verb) comment = "Check of the dynamical equations" ;
191  Tbl resu_tmp = maxabs(dyn_lhs - dyn_rhs, comment, ost, verb) ;
192 
193  if (verb) ost << endl ;
194 
195  if (vacuum) delete matter ;
196 
197  Tbl resu(4, 4, resu_tmp.get_dim(0)) ;
198  resu.annule_hard() ; //## To initialize resu(0,0,...) which
199  // does not correspond to a valid component
200 
201  for (int i=1; i<=3; i++) {
202  for (int j=1; j<=i; j++) {
203  Itbl idx(2) ;
204  idx.set(0) = i ;
205  idx.set(1) = j ;
206  int pos = dyn_lhs.position(idx) ;
207  assert (dyn_rhs.position(idx) == pos) ;
208 
209  for (int lz=0; lz<resu_tmp.get_dim(0); lz++)
210  resu.set(i,j,lz) = resu_tmp(pos, lz) ;
211  }
212  }
213 
214  return resu ;
215 
216 }
217 }
Lorene::Metric::ricci_scal
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:350
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Lorene::Time_slice::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Definition: time_slice_access.C:79
Lorene::Time_slice::k_uu
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime )
Definition: time_slice_access.C:157
Lorene::Time_slice::scheme_order
int scheme_order
Order of the finite-differences scheme for the computation of time derivatives.
Definition: time_slice.h:187
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Time_slice::check_hamiltonian_constraint
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.
Definition: tslice_check_einstein.C:79
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Tensor::up_down
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Definition: tensor_calculus.C:305
Lorene::Itbl
Basic integer array class.
Definition: itbl.h:122
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::Scalar::dec_dzpuis
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...
Definition: scalar_r_manip.C:418
Lorene::Time_slice::check_momentum_constraint
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.
Definition: tslice_check_einstein.C:109
Unites
Standard units of space, time and mass.
Lorene::Tensor_sym::position
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition: tensor_sym.C:245
Lorene::Time_slice::k_dd_evol
Evolution_std< Sym_tensor > k_dd_evol
Values at successive time steps of the covariant components of the extrinsic curvature tensor .
Definition: time_slice.h:208
Lorene::Tbl::get_dim
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Definition: tbl.h:403
Lorene::Time_slice::trk
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
Definition: time_slice_access.C:177
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Sym_tensor::derive_lie
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: sym_tensor.C:360
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Tensor::get_index_type
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:886
Lorene::maxabs
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.
Definition: tensor_calculus_ext.C:606
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
Lorene::Time_slice::check_dynamical_equations
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.
Definition: tslice_check_einstein.C:139
Lorene::Tensor::trace
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Definition: tensor_calculus.C:94
Lorene::Scalar::derive_con
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Lorene::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition: time_slice.h:190
Lorene::Metric::ricci
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition: metric.C:338
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition: time_slice_access.C:87
Lorene::Time_slice::k_dd
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
Definition: time_slice_access.C:135
Lorene::Sym_tensor::divergence
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Time_slice::gam
const Metric & gam() const
Induced metric at the current time step (jtime )
Definition: time_slice_access.C:95