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
28char 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
78namespace 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
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) ;
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
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
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" ;
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) ;
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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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...
Definition evolution.C:504
int position(int j) const
Gives the position in the arrays step, the_time and val corresponding to the time step j.
Definition evolution.C:273
Basic integer array class.
Definition itbl.h:122
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Basic array class.
Definition tbl.h:161
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
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
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
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.
int scheme_order
Order of the finite-differences scheme for the computation of time derivatives.
Definition time_slice.h:187
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 )
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime )
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.
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.
Tensor field of valence 1.
Definition vector.h:188
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.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.