LORENE
time_slice.C
1/*
2 * Methods of class Time_slice
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 time_slice_C[] = "$Header: /cvsroot/Lorene/C++/Source/Time_slice/time_slice.C,v 1.16 2014/10/13 08:53:47 j_novak Exp $" ;
29
30/*
31 * $Id: time_slice.C,v 1.16 2014/10/13 08:53:47 j_novak Exp $
32 * $Log: time_slice.C,v $
33 * Revision 1.16 2014/10/13 08:53:47 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.15 2014/10/06 15:13:21 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.14 2008/12/02 15:02:21 j_novak
40 * Implementation of the new constrained formalism, following Cordero et al. 2009
41 * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
42 *
43 * Revision 1.13 2004/05/31 09:08:18 e_gourgoulhon
44 * Method sauve and constructor from binary file are now operational.
45 *
46 * Revision 1.12 2004/05/27 15:25:04 e_gourgoulhon
47 * Added constructors from binary file, as well as corresponding
48 * functions sauve and save.
49 *
50 * Revision 1.11 2004/05/12 15:24:20 e_gourgoulhon
51 * Reorganized the #include 's, taking into account that
52 * time_slice.h contains now an #include "metric.h".
53 *
54 * Revision 1.10 2004/05/10 09:08:34 e_gourgoulhon
55 * Added "adm_mass_evol.downdate(jtime)" in method del_deriv.
56 * Added printing of ADM mass in operator>>(ostream&).
57 *
58 * Revision 1.9 2004/05/09 20:57:34 e_gourgoulhon
59 * Added data member adm_mass_evol.
60 *
61 * Revision 1.8 2004/05/05 14:26:25 e_gourgoulhon
62 * Minor modif. in operator>>(ostream& ).
63 *
64 * Revision 1.7 2004/04/07 07:58:21 e_gourgoulhon
65 * Constructor as Minkowski slice: added call to std_spectral_base()
66 * after setting the lapse to 1.
67 *
68 * Revision 1.6 2004/04/01 16:09:02 j_novak
69 * Trace of K_ij is now member of Time_slice (it was member of Time_slice_conf).
70 * Added new methods for checking 3+1 Einstein equations (preliminary).
71 *
72 * Revision 1.5 2004/03/29 11:59:23 e_gourgoulhon
73 * Added operator>>.
74 *
75 * Revision 1.4 2004/03/28 21:29:45 e_gourgoulhon
76 * Evolution_std's renamed with suffix "_evol"
77 * Method gam() modified
78 * Added special constructor for derived classes.
79 *
80 * Revision 1.3 2004/03/26 13:33:02 j_novak
81 * New methods for accessing/updating members (nn(), beta(), gam_uu(), k_uu(), ...)
82 *
83 * Revision 1.2 2004/03/26 08:22:56 e_gourgoulhon
84 * Modifications to take into account the new setting of class
85 * Evolution.
86 *
87 * Revision 1.1 2004/03/24 14:57:17 e_gourgoulhon
88 * First version
89 *
90 *
91 * $Header: /cvsroot/Lorene/C++/Source/Time_slice/time_slice.C,v 1.16 2014/10/13 08:53:47 j_novak Exp $
92 *
93 */
94
95// C headers
96#include <cassert>
97
98// Lorene headers
99#include "time_slice.h"
100#include "utilitaires.h"
101
102
103
104 //--------------//
105 // Constructors //
106 //--------------//
107
108
109namespace Lorene {
110// Standard constructor (Hamiltonian-like)
111// ---------------------------------------
113 const Sym_tensor& gamma_in, const Sym_tensor& kk_in,
114 int depth_in)
115 : depth(depth_in),
116 scheme_order(depth_in-1),
117 jtime(0),
118 the_time(0., depth_in),
119 gam_dd_evol(depth_in),
120 gam_uu_evol(depth_in),
121 k_dd_evol(depth_in),
122 k_uu_evol(depth_in),
123 n_evol(lapse_in, depth_in),
124 beta_evol(shift_in, depth_in),
125 trk_evol(depth_in),
126 adm_mass_evol() {
127
128 set_der_0x0() ;
129
130 double time_init = the_time[jtime] ;
131
132 p_gamma = new Metric(gamma_in) ;
133
134 if (gamma_in.get_index_type(0) == COV) {
136 }
137 else {
139 }
140
141 if (kk_in.get_index_type(0) == COV) {
143 }
144 else {
146 }
147
149
150}
151
152
153// Standard constructor (Lagrangian-like)
154// ---------------------------------------
157 : depth(gamma_in.get_size()),
158 scheme_order(gamma_in.get_size()-1),
159 jtime(0),
160 the_time(0., gamma_in.get_size()),
161 gam_dd_evol( gamma_in.get_size() ),
162 gam_uu_evol( gamma_in.get_size() ),
163 k_dd_evol( gamma_in.get_size() ),
164 k_uu_evol( gamma_in.get_size() ),
165 n_evol(lapse_in, gamma_in.get_size() ),
166 beta_evol(shift_in, gamma_in.get_size() ),
167 trk_evol(gamma_in.get_size() ),
168 adm_mass_evol() {
169
170 cerr <<
171 "Time_slice constuctor from evolution of gamma not implemented yet !\n" ;
172 abort() ;
173
174 set_der_0x0() ;
175
176}
177
178// Constructor as standard time slice of flat spacetime (Minkowski)
179//-----------------------------------------------------------------
180Time_slice::Time_slice(const Map& mp, const Base_vect& triad, int depth_in)
181 : depth(depth_in),
182 scheme_order(depth_in-1),
183 jtime(0),
184 the_time(0., depth_in),
185 gam_dd_evol(depth_in),
186 gam_uu_evol(depth_in),
187 k_dd_evol(depth_in),
188 k_uu_evol(depth_in),
189 n_evol(depth_in),
190 beta_evol(depth_in),
191 trk_evol(depth_in),
192 adm_mass_evol() {
193
194 double time_init = the_time[jtime] ;
195
196 const Base_vect_spher* ptriad_s =
197 dynamic_cast<const Base_vect_spher*>(&triad) ;
198 bool spher = (ptriad_s != 0x0) ;
199
200 if (spher) {
202 }
203 else {
204 assert( dynamic_cast<const Base_vect_cart*>(&triad) != 0x0) ;
206 }
207
208
209 // K_ij identically zero:
210 Sym_tensor ktmp(mp, COV, triad) ;
211 ktmp.set_etat_zero() ;
213
214 // Lapse identically one:
215 Scalar tmp(mp) ;
216 tmp.set_etat_one() ;
217 tmp.std_spectral_base() ;
219
220 // shift identically zero:
221 Vector btmp(mp, CON, triad) ;
222 btmp.set_etat_zero() ;
224
225 // trace(K) identically zero:
226 tmp.set_etat_zero() ;
228
229 set_der_0x0() ;
230}
231
232// Constructor from binary file
233// ----------------------------
234
235Time_slice::Time_slice(const Map& mp, const Base_vect& triad, FILE* fich,
236 bool partial_read, int depth_in)
237 : depth(depth_in),
238 the_time(depth_in),
239 gam_dd_evol(depth_in),
240 gam_uu_evol(depth_in),
241 k_dd_evol(depth_in),
242 k_uu_evol(depth_in),
243 n_evol(depth_in),
244 beta_evol(depth_in),
245 trk_evol(depth_in),
246 adm_mass_evol() {
247
248 // Reading various integer parameters
249 // ----------------------------------
250
251 int depth_file ;
252 fread_be(&depth_file, sizeof(int), 1, fich) ;
253 if (depth_file != depth_in) {
254 cout <<
255 "Time_slice constructor from file: the depth read in file \n"
256 << " is different from that given in the argument list : \n"
257 << " depth_file = " << depth_file
258 << " <-> depth_in " << depth_in << " !" << endl ;
259 abort() ;
260 }
261 fread_be(&scheme_order, sizeof(int), 1, fich) ;
262 fread_be(&jtime, sizeof(int), 1, fich) ;
263
264 // Reading the_time
265 // ----------------
266 int jmin = jtime - depth + 1 ;
267 int indicator ;
268 for (int j=jmin; j<=jtime; j++) {
269 fread_be(&indicator, sizeof(int), 1, fich) ;
270 if (indicator == 1) {
271 double xx ;
272 fread_be(&xx, sizeof(double), 1, fich) ;
273 the_time.update(xx, j, xx) ;
274 }
275 }
276
277 // Reading of various fields
278 // -------------------------
279
280 // N
281 for (int j=jmin; j<=jtime; j++) {
282 fread_be(&indicator, sizeof(int), 1, fich) ;
283 if (indicator == 1) {
284 Scalar nn_file(mp, *(mp.get_mg()), fich) ;
286 }
287 }
288
289 // beta
290 for (int j=jmin; j<=jtime; j++) {
291 fread_be(&indicator, sizeof(int), 1, fich) ;
292 if (indicator == 1) {
293 Vector beta_file(mp, triad, fich) ;
295 }
296 }
297
298 // Case of a full reading
299 // -----------------------
300 if (!partial_read) {
301 cout <<
302 "Time_slice constructor from file: the case of full reading\n"
303 << " is not ready yet !" << endl ;
304 abort() ;
305 }
306
307 set_der_0x0() ;
308
309}
310
311
312// Copy constructor
313// ----------------
315 : depth(tin.depth),
316 scheme_order(tin.scheme_order),
317 jtime(tin.jtime),
318 the_time(tin.the_time),
319 gam_dd_evol(tin.gam_dd_evol),
320 gam_uu_evol(tin.gam_uu_evol),
321 k_dd_evol(tin.k_dd_evol),
322 k_uu_evol(tin.k_uu_evol),
323 n_evol(tin.n_evol),
324 beta_evol(tin.beta_evol),
325 trk_evol(tin.trk_evol),
326 adm_mass_evol(tin.adm_mass_evol) {
327
328 set_der_0x0() ;
329}
330
331// Special constructor for derived classes
332//----------------------------------------
334 : depth(depth_in),
335 scheme_order(depth_in-1),
336 jtime(0),
337 the_time(0., depth_in),
338 gam_dd_evol(depth_in),
339 gam_uu_evol(depth_in),
340 k_dd_evol(depth_in),
341 k_uu_evol(depth_in),
342 n_evol(depth_in),
343 beta_evol(depth_in),
344 trk_evol(depth_in),
345 adm_mass_evol() {
346
347 set_der_0x0() ;
348}
349
350
351
352
353 //--------------//
354 // Destructor //
355 //--------------//
356
362
363 //---------------------//
364 // Memory management //
365 //---------------------//
366
368
369 if (p_gamma != 0x0) delete p_gamma ;
370
371 set_der_0x0() ;
372
374}
375
376
378
379 p_gamma = 0x0 ;
380
381}
382
383
384 //-----------------------//
385 // Mutators / assignment //
386 //-----------------------//
387
389
390 depth = tin.depth;
391 scheme_order = tin.scheme_order ;
392 jtime = tin.jtime;
394 gam_dd_evol = tin.gam_dd_evol;
395 gam_uu_evol = tin.gam_uu_evol;
396 k_dd_evol = tin.k_dd_evol;
397 k_uu_evol = tin.k_uu_evol;
398 n_evol = tin.n_evol;
399 beta_evol = tin.beta_evol;
400 trk_evol = tin.trk_evol ;
401
402 del_deriv() ;
403
404}
405
406
407 //------------------//
408 // output //
409 //------------------//
410
412
413 flux << "\n------------------------------------------------------------\n"
414 << "Lorene class : " << typeid(*this).name() << '\n' ;
415 flux << "Number of stored slices : " << depth
416 << " order of time scheme : " << scheme_order << '\n'
417 << "Time label t = " << the_time[jtime]
418 << " index of time step j = " << jtime << '\n' << '\n' ;
420 flux << "ADM mass : " << adm_mass() << endl ;
421 }
422
423 flux << "Max. of absolute values of the various fields in each domain: \n" ;
425 maxabs( gam_dd_evol[jtime], "gam_{ij}", flux) ;
426 }
428 maxabs( gam_uu_evol[jtime], "gam^{ij}", flux) ;
429 }
430 if (k_dd_evol.is_known(jtime)) {
431 maxabs( k_dd_evol[jtime], "K_{ij}", flux) ;
432 }
433 if (k_uu_evol.is_known(jtime)) {
434 maxabs( k_uu_evol[jtime], "K^{ij}", flux) ;
435 }
436 if (n_evol.is_known(jtime)) {
437 maxabs( n_evol[jtime], "N", flux) ;
438 }
439 if (beta_evol.is_known(jtime)) {
440 maxabs( beta_evol[jtime], "beta^i", flux) ;
441 }
442 if (trk_evol.is_known(jtime)) {
443 maxabs( trk_evol[jtime], "tr K", flux) ;
444 }
445
446 if (p_gamma != 0x0) flux << "Metric gamma is up to date" << endl ;
447
448 return flux ;
449
450}
451
452
454
455 sigma >> flux ;
456 return flux ;
457
458}
459
460
461void Time_slice::save(const char* rootname) const {
462
463 // Opening of file
464 // ---------------
465 char* filename = new char[ strlen(rootname)+10 ] ;
467 char nomj[7] ;
468 sprintf(nomj, "%06d", jtime) ;
469 strcat(filename, nomj) ;
470 strcat(filename, ".d") ;
471
472 FILE* fich = fopen(filename, "w") ;
473 if (fich == 0x0) {
474 cout << "Problem in opening file " << filename << " ! " << endl ;
475 perror(" reason") ;
476 abort() ;
477 }
478
479 // Write grid, mapping, triad and depth
480 // ------------------------------------
481 const Map& map = nn().get_mp() ;
482 const Mg3d& mgrid = *(map.get_mg()) ;
483 const Base_vect& triad = *(beta().get_triad()) ;
484
485 mgrid.sauve(fich) ;
486 map.sauve(fich) ;
487 triad.sauve(fich) ;
488
489 fwrite_be(&depth, sizeof(int), 1, fich) ;
490
491 // Write all binary data by means of virtual function sauve
492 // --------------------------------------------------------
493 bool partial_save = false ;
495
496 // Close the file
497 // --------------
498
499 fclose(fich) ;
500
501 delete [] filename ;
502
503}
504
505
506
508
509 // Writing various integer parameters
510 // ----------------------------------
511
512 fwrite_be(&depth, sizeof(int), 1, fich) ;
513 fwrite_be(&scheme_order, sizeof(int), 1, fich) ;
514 fwrite_be(&jtime, sizeof(int), 1, fich) ;
515
516 // Writing the_time
517 // ----------------
518 int jmin = jtime - depth + 1 ;
519 for (int j=jmin; j<=jtime; j++) {
520 int indicator = (the_time.is_known(j)) ? 1 : 0 ;
521 fwrite_be(&indicator, sizeof(int), 1, fich) ;
522 if (indicator == 1) {
523 double xx = the_time[j] ;
524 fwrite_be(&xx, sizeof(double), 1, fich) ;
525 }
526 }
527
528 // Writing of various fields
529 // -------------------------
530
531 // N
532 nn() ; // forces the update at the current time step
533 for (int j=jmin; j<=jtime; j++) {
534 int indicator = (n_evol.is_known(j)) ? 1 : 0 ;
535 fwrite_be(&indicator, sizeof(int), 1, fich) ;
536 if (indicator == 1) n_evol[j].sauve(fich) ;
537 }
538
539 // beta
540 beta() ; // forces the update at the current time step
541 for (int j=jmin; j<=jtime; j++) {
542 int indicator = (beta_evol.is_known(j)) ? 1 : 0 ;
543 fwrite_be(&indicator, sizeof(int), 1, fich) ;
544 if (indicator == 1) beta_evol[j].sauve(fich) ;
545 }
546
547 // Case of a complete save
548 // -----------------------
549 if (!partial_save) {
550
551 cout << "Time_slice::sauve: the full writing is not ready yet !"
552 << endl ;
553 abort() ;
554 }
555
556
557}
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572}
Cartesian vectorial bases (triads).
Definition base_vect.h:201
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
virtual void sauve(FILE *) const
Save in a file.
Definition base_vect.C:150
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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.
Definition evolution.C:244
double * the_time
Array of values of t at the various time steps (size = size).
Definition evolution.h:134
bool is_known(int j) const
Checks whether the value a given time step has been set.
Definition evolution.C:335
Base class for coordinate mappings.
Definition map.h:670
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition map.C:331
virtual void sauve(FILE *) const
Save in a file.
Definition map.C:224
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321
Metric for tensor calculation.
Definition metric.h:90
Multi-domain grid.
Definition grilles.h:273
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
Spacelike time slice of a 3+1 spacetime.
Definition time_slice.h:173
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition time_slice.h:224
Metric * p_gamma
Pointer on the induced metric at the current time step (jtime)
Definition time_slice.h:239
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 void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
Definition time_slice.C:507
void operator=(const Time_slice &)
Assignment to another Time_slice.
Definition time_slice.C:388
Evolution_full< Tbl > adm_mass_evol
ADM mass at each time step, since the creation of the slice.
Definition time_slice.h:233
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
void save(const char *rootname) const
Saves in a binary file.
Definition time_slice.C:461
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition time_slice.C:377
Evolution_std< Sym_tensor > gam_dd_evol
Values at successive time steps of the covariant components of the induced metric .
Definition time_slice.h:198
Evolution_std< Sym_tensor > gam_uu_evol
Values at successive time steps of the contravariant components of the induced metric .
Definition time_slice.h:203
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 )
virtual ~Time_slice()
Destructor.
Definition time_slice.C:357
virtual void del_deriv() const
Deletes all the derived quantities.
Definition time_slice.C:367
int depth
Number of stored time slices.
Definition time_slice.h:179
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition time_slice.h:216
virtual double adm_mass() const
Returns the ADM mass (geometrical units) at the current step.
Time_slice(const Scalar &lapse_in, const Vector &shift_in, const Sym_tensor &gamma_in, const Sym_tensor &kk_in, int depth_in=3)
General constructor (Hamiltonian-like).
Definition time_slice.C:112
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Definition time_slice.C:411
Evolution_std< Sym_tensor > k_uu_evol
Values at successive time steps of the contravariant components of the extrinsic curvature tensor .
Definition time_slice.h:213
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition time_slice.h:219
Tensor field of valence 1.
Definition vector.h:188
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
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