LORENE
time_slice_conf.C
1 /*
2  * Methods of class Time_slice_conf
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 
29 
30 /*
31  * $Id: time_slice_conf.C,v 1.19 2016/12/05 16:18:19 j_novak Exp $
32  * $Log: time_slice_conf.C,v $
33  * Revision 1.19 2016/12/05 16:18:19 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.18 2014/10/13 08:53:47 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.17 2014/10/06 15:13:21 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.16 2008/12/04 18:22:49 j_novak
43  * Enhancement of the dzpuis treatment + various bug fixes.
44  *
45  * Revision 1.15 2008/12/02 15:02:22 j_novak
46  * Implementation of the new constrained formalism, following Cordero et al. 2009
47  * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
48  *
49  * Revision 1.14 2004/06/24 20:36:54 e_gourgoulhon
50  * Added method check_psi_dot.
51  *
52  * Revision 1.13 2004/05/31 09:08:18 e_gourgoulhon
53  * Method sauve and constructor from binary file are now operational.
54  *
55  * Revision 1.12 2004/05/27 15:25:04 e_gourgoulhon
56  * Added constructors from binary file, as well as corresponding
57  * functions sauve and save.
58  *
59  * Revision 1.11 2004/05/12 15:24:20 e_gourgoulhon
60  * Reorganized the #include 's, taking into account that
61  * time_slice.h contains now an #include "metric.h".
62  *
63  * Revision 1.10 2004/05/10 09:10:05 e_gourgoulhon
64  * Added "adm_mass_evol.downdate(jtime)" in methods set_*.
65  *
66  * Revision 1.9 2004/05/05 14:31:14 e_gourgoulhon
67  * Method aa(): added *as a comment* annulation of hh_point in the compactified
68  * domain.
69  *
70  * Revision 1.8 2004/05/03 14:47:11 e_gourgoulhon
71  * Corrected method aa().
72  *
73  * Revision 1.7 2004/04/08 16:43:26 e_gourgoulhon
74  * Added methods set_*
75  * Added test of determinant one in constructor and set_hh.
76  *
77  * Revision 1.6 2004/04/05 21:25:02 e_gourgoulhon
78  * -- Added constructor as standard time slice of Minkowski spacetime.
79  * -- Added some calls to Scalar::std_spectral_base() after
80  * non-arithmetical operations.
81  *
82  * Revision 1.5 2004/04/05 12:38:45 j_novak
83  * Minor modifs to prevent some warnings.
84  *
85  * Revision 1.4 2004/04/01 16:09:02 j_novak
86  * Trace of K_ij is now member of Time_slice (it was member of Time_slice_conf).
87  * Added new methods for checking 3+1 Einstein equations (preliminary).
88  *
89  * Revision 1.3 2004/03/29 12:00:41 e_gourgoulhon
90  * Many modifs.
91  *
92  * Revision 1.2 2004/03/28 21:32:23 e_gourgoulhon
93  * Corrected error in method trk().
94  *
95  * Revision 1.1 2004/03/28 21:30:13 e_gourgoulhon
96  * First version.
97  *
98  *
99  * $Header: /cvsroot/Lorene/C++/Source/Time_slice/time_slice_conf.C,v 1.19 2016/12/05 16:18:19 j_novak Exp $
100  *
101  */
102 
103 // C headers
104 #include <cassert>
105 
106 // Lorene headers
107 #include "time_slice.h"
108 #include "utilitaires.h"
109 
110 
111 
112  //--------------//
113  // Constructors //
114  //--------------//
115 
116 
117 // Constructor from conformal decomposition
118 // ----------------------------------------
119 
120 namespace Lorene {
121 Time_slice_conf::Time_slice_conf(const Scalar& lapse_in, const Vector& shift_in,
122  const Metric_flat& ff_in, const Scalar& psi_in,
123  const Sym_tensor& hh_in, const Sym_tensor& hata_in,
124  const Scalar& trk_in, int depth_in)
125  : Time_slice(depth_in),
126  ff(ff_in),
127  psi_evol(psi_in, depth_in),
128  npsi_evol(depth_in),
129  hh_evol(hh_in, depth_in),
130  hata_evol(hata_in, depth_in),
131  A_hata_evol(depth_in), B_hata_evol(depth_in){
132 
133  assert(hh_in.get_index_type(0) == CON) ;
134  assert(hh_in.get_index_type(1) == CON) ;
135  assert(hata_in.get_index_type(0) == CON) ;
136  assert(hata_in.get_index_type(1) == CON) ;
137 
138  double time_init = the_time[jtime] ;
139 
140  // Check whether det tgam^{ij} = det f^{ij} :
141  // ----------------------------------------
142  Sym_tensor tgam_in = ff_in.con() + hh_in ;
143 
144  Scalar det_in = tgam_in(1, 1)*tgam_in(2, 2)*tgam_in(3, 3)
145  + tgam_in(1, 2)*tgam_in(2, 3)*tgam_in(3, 1)
146  + tgam_in(1, 3)*tgam_in(2, 1)*tgam_in(3, 2)
147  - tgam_in(3, 1)*tgam_in(2, 2)*tgam_in(1, 3)
148  - tgam_in(3, 2)*tgam_in(2, 3)*tgam_in(1, 1)
149  - tgam_in(3, 3)*tgam_in(2, 1)*tgam_in(1, 2) ;
150 
151  double diffdet = max(maxabs(det_in - 1. / ff.determinant(),
152  "Deviation of det tgam^{ij} from 1/f")) ;
153  if ( diffdet > 1.e-13 ) {
154  cerr <<
155  "Time_slice_conf::Time_slice_conf : the input h^{ij} does not"
156  << " ensure \n" << " det tgam_{ij} = f ! \n"
157  << " error = " << diffdet << endl ;
158  abort() ;
159  }
160 
161  n_evol.update(lapse_in, jtime, time_init) ;
162  beta_evol.update(shift_in, jtime, time_init) ;
163  trk_evol.update(trk_in, jtime, time_init) ;
164  A_hata() ;
165  B_hata() ;
166 
167  set_der_0x0() ;
168 
169 }
170 
171 
172 // Constructor from physical metric
173 // --------------------------------
174 
175 Time_slice_conf::Time_slice_conf(const Scalar& lapse_in, const Vector& shift_in,
176  const Sym_tensor& gamma_in, const Sym_tensor& kk_in,
177  const Metric_flat& ff_in, int depth_in)
178  : Time_slice(lapse_in, shift_in, gamma_in, kk_in, depth_in),
179  ff(ff_in),
180  psi_evol(depth_in),
181  npsi_evol(depth_in),
182  hh_evol(depth_in),
183  hata_evol(depth_in),
184  A_hata_evol(depth_in), B_hata_evol(depth_in) {
185 
186  set_der_0x0() ; // put here in order not to erase p_psi4
187 
188  double time_init = the_time[jtime] ;
189 
190  assert( p_gamma != 0x0 ) ;
192  0.3333333333333333) ) ;
194 
195  Scalar tmp = pow(*p_psi4, 0.25) ;
196  tmp.std_spectral_base() ;
197  psi_evol.update(tmp , jtime, time_init ) ;
198 
199  hh_evol.update( (*p_psi4) * p_gamma->con() - ff.con(),
200  jtime, time_init ) ;
201 
202  hata_evol.update( tmp*tmp*(*p_psi4)*(*p_psi4) *( Time_slice::k_uu()
203  - 0.3333333333333333 * trk_evol[jtime] * p_gamma->con() ),
204  jtime, time_init ) ;
205  A_hata() ;
206  B_hata() ;
207 }
208 
209 // Constructor as standard time slice of flat spacetime (Minkowski)
210 // ----------------------------------------------------------------
211 
213  const Metric_flat& ff_in, int depth_in)
214  : Time_slice(mp, triad, depth_in),
215  ff(ff_in),
216  psi_evol(depth_in),
217  npsi_evol(depth_in),
218  hh_evol(depth_in),
219  hata_evol(depth_in),
220  A_hata_evol(depth_in), B_hata_evol(depth_in) {
221 
222  double time_init = the_time[jtime] ;
223 
224  // Psi identically one:
225  Scalar tmp(mp) ;
226  tmp.set_etat_one() ;
227  tmp.std_spectral_base() ;
228  psi_evol.update(tmp, jtime, time_init) ;
229 
230  // N Psi identically one:
231  npsi_evol.update(tmp, jtime, time_init) ;
232 
233  // h^{ij} identically zero:
234  Sym_tensor stmp(mp, CON, triad) ;
235  stmp.set_etat_zero() ;
236  hh_evol.update(stmp, jtime, time_init) ;
237 
238  // \hat{A}^{ij} identically zero:
239  hata_evol.update(stmp, jtime, time_init) ;
240 
241  tmp.set_etat_zero() ;
242  A_hata_evol.update(tmp, jtime, time_init) ;
243  B_hata_evol.update(tmp, jtime, time_init) ;
244 
245  set_der_0x0() ;
246 
247 }
248 
249 
250 // Constructor from binary file
251 // ----------------------------
252 
254  const Metric_flat& ff_in, FILE* fich,
255  bool partial_read, int depth_in)
256  : Time_slice(mp, triad, fich, true, depth_in),
257  ff(ff_in),
258  psi_evol(depth_in),
259  npsi_evol(depth_in),
260  hh_evol(depth_in),
261  hata_evol(depth_in),
262  A_hata_evol(depth_in), B_hata_evol(depth_in) {
263 
264  // Put here, not to destroy p_vec_X
265  set_der_0x0() ;
266 
267  // Reading of various fields
268  // -------------------------
269 
270  int jmin = jtime - depth + 1 ;
271  int indicator ;
272 
273  // Psi
274  for (int j=jmin; j<=jtime; j++) {
275  fread_be(&indicator, sizeof(int), 1, fich) ;
276  if (indicator == 1) {
277  Scalar psi_file(mp, *(mp.get_mg()), fich) ;
278  psi_evol.update(psi_file, j, the_time[j]) ;
279  }
280  }
281  // hat{A}^{ij}
282  for (int j=jmin; j<=jtime; j++) {
283  fread_be(&indicator, sizeof(int), 1, fich) ;
284  if (indicator == 1) {
285  Sym_tensor hat_A_file(mp, triad, fich) ;
286  hata_evol.update( hat_A_file, j, the_time[j] ) ;
287  }
288  }
289 
290  //A and B...
291  for (int j=jmin; j<=jtime; j++) {
292  fread_be(&indicator, sizeof(int), 1, fich) ;
293  if (indicator == 1) {
294  Scalar A_file(mp, *(mp.get_mg()), fich) ;
295  A_hata_evol.update(A_file, j, the_time[j]) ;
296  }
297  }
298  for (int j=jmin; j<=jtime; j++) {
299  fread_be(&indicator, sizeof(int), 1, fich) ;
300  if (indicator == 1) {
301  Scalar B_file(mp, *(mp.get_mg()), fich) ;
302  B_hata_evol.update(B_file, j, the_time[j]) ;
303  }
304  }
305 
306  // Case of a full reading
307  // -----------------------
308  if (!partial_read) {
309  cout <<
310  "Time_slice constructor from file: the case of full reading\n"
311  << " is not ready yet !" << endl ;
312  abort() ;
313  }
314 
315 }
316 
317 
318 
319 
320 // Copy constructor
321 // ----------------
322 
324  : Time_slice(tin),
325  ff(tin.ff),
326  psi_evol(tin.psi_evol),
327  npsi_evol(tin.npsi_evol),
328  hh_evol(tin.hh_evol),
329  hata_evol(tin.hata_evol),
330  A_hata_evol(tin.A_hata_evol),
331  B_hata_evol(tin.B_hata_evol){
332 
333  set_der_0x0() ;
334 
335 }
336  //--------------//
337  // Destructor //
338  //--------------//
339 
341 
343 
344 }
345 
346  //---------------------//
347  // Memory management //
348  //---------------------//
349 
351 
352  if (p_tgamma != 0x0) delete p_tgamma ;
353  if (p_psi4 != 0x0) delete p_psi4 ;
354  if (p_ln_psi != 0x0) delete p_ln_psi ;
355  if (p_hdirac != 0x0) delete p_hdirac ;
356  if (p_vec_X != 0x0) delete p_vec_X ;
357 
358  set_der_0x0() ;
359 
361 }
362 
363 
365 
366  p_tgamma = 0x0 ;
367  p_psi4 = 0x0 ;
368  p_ln_psi = 0x0 ;
369  p_hdirac = 0x0 ;
370  p_vec_X = 0x0 ;
371 
372 }
373 
374 
375  //-----------------------//
376  // Mutators / assignment //
377  //-----------------------//
378 
380 
381  Time_slice::operator=(tin) ;
382 
383  psi_evol = tin.psi_evol ;
384  npsi_evol = tin.npsi_evol ;
385  hh_evol = tin.hh_evol ;
386  hata_evol = tin.hata_evol ;
387  A_hata_evol = tin.A_hata_evol ;
388  B_hata_evol = tin.B_hata_evol ;
389 
390  del_deriv() ;
391 
392 }
393 
395 
396  Time_slice::operator=(tin) ;
397 
398  cerr <<
399  "Time_slice_conf::operator=(const Time_slice& ) : not implemented yet !"
400  << endl ;
401  abort() ;
402  del_deriv() ;
403 
404 }
405 
406 
408 
409  psi_evol.update(psi_in, jtime, the_time[jtime]) ;
410 
411  // Reset of quantities depending on Psi:
412  if (p_psi4 != 0x0) {
413  delete p_psi4 ;
414  p_psi4 = 0x0 ;
415  }
416  if (p_ln_psi != 0x0) {
417  delete p_ln_psi ;
418  p_ln_psi = 0x0 ;
419  }
420  if (p_gamma != 0x0) {
421  delete p_gamma ;
422  p_gamma = 0x0 ;
423  }
424  npsi_evol.downdate(jtime) ;
425  gam_dd_evol.downdate(jtime) ;
426  gam_uu_evol.downdate(jtime) ;
427  adm_mass_evol.downdate(jtime) ;
428 
429 }
430 
432 
433  psi_evol.update(psi_in, jtime, the_time[jtime]) ;
434 
435  // Reset of quantities depending on Psi:
436  if (p_psi4 != 0x0) {
437  delete p_psi4 ;
438  p_psi4 = 0x0 ;
439  }
440  if (p_ln_psi != 0x0) {
441  delete p_ln_psi ;
442  p_ln_psi = 0x0 ;
443  }
444  if (p_gamma != 0x0) {
445  delete p_gamma ;
446  p_gamma = 0x0 ;
447  }
448  n_evol.downdate(jtime) ;
449  gam_dd_evol.downdate(jtime) ;
450  gam_uu_evol.downdate(jtime) ;
451  adm_mass_evol.downdate(jtime) ;
452 
453 }
454 
455 
457 
458  npsi_evol.update(npsi_in, jtime, the_time[jtime]) ;
459 
460  // Reset of quantities depending on N Psi:
461  psi_evol.downdate(jtime) ;
462  if (p_psi4 != 0x0) {
463  delete p_psi4 ;
464  p_psi4 = 0x0 ;
465  }
466  if (p_ln_psi != 0x0) {
467  delete p_ln_psi ;
468  p_ln_psi = 0x0 ;
469  }
470  if (p_gamma != 0x0) {
471  delete p_gamma ;
472  p_gamma = 0x0 ;
473  }
474  gam_dd_evol.downdate(jtime) ;
475  gam_uu_evol.downdate(jtime) ;
476  adm_mass_evol.downdate(jtime) ;
477 
478 }
479 
480 
482 
483  npsi_evol.update(npsi_in, jtime, the_time[jtime]) ;
484 
485  // Reset of quantities depending on N Psi:
486  n_evol.downdate(jtime) ;
487  adm_mass_evol.downdate(jtime) ;
488 
489 }
490 
491 
493 
494  // Check whether det tgam^{ij} = det f^{ij} :
495  // ----------------------------------------
496  Sym_tensor tgam_in = ff.con() + hh_in ;
497 
498  Scalar det_in = tgam_in(1, 1)*tgam_in(2, 2)*tgam_in(3, 3)
499  + tgam_in(1, 2)*tgam_in(2, 3)*tgam_in(3, 1)
500  + tgam_in(1, 3)*tgam_in(2, 1)*tgam_in(3, 2)
501  - tgam_in(3, 1)*tgam_in(2, 2)*tgam_in(1, 3)
502  - tgam_in(3, 2)*tgam_in(2, 3)*tgam_in(1, 1)
503  - tgam_in(3, 3)*tgam_in(2, 1)*tgam_in(1, 2) ;
504 
505  double diffdet = max(maxabs(det_in - 1. / ff.determinant(),
506  "Deviation of det tgam^{ij} from 1/f")) ;
507  if ( diffdet > 1.e-13 ) {
508  cerr <<
509  "Time_slice_conf::set_hh : the input h^{ij} does not"
510  << " ensure \n" << " det tgam_{ij} = f ! \n"
511  << " error = " << diffdet << endl ;
512  abort() ;
513  }
514 
515  hh_evol.update(hh_in, jtime, the_time[jtime]) ;
516 
517  // Reset of quantities depending on h^{ij}:
518  if (p_tgamma != 0x0) {
519  delete p_tgamma ;
520  p_tgamma = 0x0 ;
521  }
522  if (p_hdirac != 0x0) {
523  delete p_hdirac ;
524  p_hdirac = 0x0 ;
525  }
526  if (p_gamma != 0x0) {
527  delete p_gamma ;
528  p_gamma = 0x0 ;
529  }
530  gam_dd_evol.downdate(jtime) ;
531  gam_uu_evol.downdate(jtime) ;
532  adm_mass_evol.downdate(jtime) ;
533 
534 }
535 
536 
537 void Time_slice_conf::set_hata(const Sym_tensor& hata_in) {
538 
539  hata_evol.update(hata_in, jtime, the_time[jtime]) ;
540 
541  // Reset of quantities depending on A^{ij}:
542  A_hata_evol.downdate(jtime) ;
543  B_hata_evol.downdate(jtime) ;
544  k_dd_evol.downdate(jtime) ;
545  k_uu_evol.downdate(jtime) ;
546 
547 }
548 
550 
551  Scalar tmp = hata_tt.compute_A(true) ;
552  if (tmp.get_dzpuis() == 3)
553  tmp.dec_dzpuis() ; // dzpuis 3->2
554 
555  A_hata_evol.update( tmp, jtime, the_time[jtime] ) ;
556 
557  tmp = hata_tt.compute_tilde_B_tt(true) ;
558  if (tmp.get_dzpuis() == 3)
559  tmp.dec_dzpuis() ; // dzpuis 3->2
560 
561  B_hata_evol.update( tmp, jtime, the_time[jtime] ) ;
562 
563  hata_evol.downdate(jtime) ;
564  // Reset of quantities depending on A^{ij}:
565  k_dd_evol.downdate(jtime) ;
566  k_uu_evol.downdate(jtime) ;
567 }
568 
570 
571  assert (p_vec_X != 0x0) ;
572  assert (A_hata_evol.is_known(jtime)) ;
573  assert (B_hata_evol.is_known(jtime)) ;
574 
575  const Map& mp = p_vec_X->get_mp() ;
576 
577  Sym_tensor_tt hata_tt(mp, mp.get_bvect_spher(), ff) ;
578  hata_tt.set_A_tildeB(A_hata_evol[jtime], B_hata_evol[jtime], par_bc, par_mat) ;
579  hata_tt.inc_dzpuis(2) ;
580 
581  hata_evol.update( hata_tt + p_vec_X->ope_killing_conf(ff), jtime, the_time[jtime]) ;
582 
583  // Reset of quantities depending on A^{ij}:
584  k_dd_evol.downdate(jtime) ;
585  k_uu_evol.downdate(jtime) ;
586 
587 }
588 
589 
590  //-----------------------------------------------//
591  // Update of fields from base class Time_slice //
592  //-----------------------------------------------//
593 
594 const Scalar& Time_slice_conf::nn() const {
595 
596  if (!( n_evol.is_known(jtime) ) ) {
597 
598  assert( psi_evol.is_known(jtime) ) ;
599  assert( npsi_evol.is_known(jtime) ) ;
600 
601  n_evol.update( npsi_evol[jtime] / psi_evol[jtime] ,
602  jtime, the_time[jtime] ) ;
603  }
604 
605  return n_evol[jtime] ;
606 
607 }
608 
609 
610 
612 
613  if (!( gam_dd_evol.is_known(jtime)) ) {
614  gam_dd_evol.update( psi4() * tgam().cov(), jtime, the_time[jtime] ) ;
615  }
616 
617  return gam_dd_evol[jtime] ;
618 
619 }
620 
621 
623 
624  if (!( gam_uu_evol.is_known(jtime)) ) {
625  gam_uu_evol.update( tgam().con() / psi4() , jtime, the_time[jtime] ) ;
626  }
627 
628  return gam_uu_evol[jtime] ;
629 
630 }
631 
632 
634 
635  if ( ! (k_dd_evol.is_known(jtime)) ) {
636 
637  k_dd_evol.update( k_uu().up_down(gam()), jtime, the_time[jtime] ) ;
638 
639  }
640 
641  return k_dd_evol[jtime] ;
642 
643 }
644 
645 
647 
648  if ( ! (k_uu_evol.is_known(jtime)) ) {
649 
650  k_uu_evol.update( hata()/(psi4()*psi4()*psi()*psi())
651  + 0.3333333333333333* trk()* gam().con(),
652  jtime, the_time[jtime] ) ;
653  }
654 
655  return k_uu_evol[jtime] ;
656 
657 }
658 
659 
660 
661 
662 
663  //-----------------------------------//
664  // Update of fields from this class //
665  //-----------------------------------//
666 
668 
669  if ( !( A_hata_evol.is_known(jtime) ) ) {
670 
671  assert( hata_evol.is_known(jtime) ) ;
672  Scalar tmp = hata_evol[jtime].compute_A(true) ;
673  if (tmp.get_dzpuis() == 3)
674  tmp.dec_dzpuis() ; // dzpuis 3->2
675 
676  A_hata_evol.update( tmp, jtime, the_time[jtime] ) ;
677  }
678  return A_hata_evol[jtime] ;
679 }
680 
682 
683  if ( !( B_hata_evol.is_known(jtime) ) ) {
684 
685  assert( hata_evol.is_known(jtime) ) ;
686  Scalar tmp = hata_evol[jtime].compute_tilde_B_tt(true) ;
687  if (tmp.get_dzpuis() == 3)
688  tmp.dec_dzpuis() ; // dzpuis 3->2
689 
690  B_hata_evol.update( tmp, jtime, the_time[jtime] ) ;
691  }
692  return A_hata_evol[jtime] ;
693 }
694 
695 
696 const Scalar& Time_slice_conf::psi() const {
697 
698  if (!( psi_evol.is_known(jtime) ) ) {
699 
700  assert( n_evol.is_known(jtime) ) ;
701  assert( npsi_evol.is_known(jtime) ) ;
702 
703  psi_evol.update( npsi_evol[jtime] / n_evol[jtime] , jtime, the_time[jtime] ) ;
704  }
705 
706  return psi_evol[jtime] ;
707 
708 }
709 
711 
712  if (p_psi4 == 0x0) {
713 
714  p_psi4 = new Scalar( pow( psi(), 4.) ) ;
716  }
717 
718  return *p_psi4 ;
719 
720 }
721 
723 
724  if (p_ln_psi == 0x0) {
725 
726  p_ln_psi = new Scalar( log( psi() ) ) ;
728  }
729 
730  return *p_ln_psi ;
731 
732 }
733 
734 
736 
737  if (!( npsi_evol.is_known(jtime) ) ) {
738 
739  assert( n_evol.is_known(jtime) ) ;
740  assert( psi_evol.is_known(jtime) ) ;
741 
742  npsi_evol.update( psi_evol[jtime] * n_evol[jtime], jtime, the_time[jtime] ) ;
743  }
744 
745  return npsi_evol[jtime] ;
746 
747 }
748 
749 
751 
752  if (p_tgamma == 0x0) {
753  p_tgamma = new Metric( ff.con() + hh() ) ;
754  }
755 
756  return *p_tgamma ;
757 
758 }
759 
760 
762 
763  assert( hh_evol.is_known(jtime) ) ;
764  return hh_evol[jtime] ;
765 
766 }
767 
769 
770  return hata()/(psi4()*psi()*psi()) ;
771 
772 }
773 
774 
776 
777  if( !(hata_evol.is_known(jtime)) ) {
778 
779  assert( hh_evol.is_known(jtime) ) ;
780 
781  Sym_tensor hh_point = hh_evol.time_derive(jtime, scheme_order) ;
782  hh_point.inc_dzpuis(2) ; // dzpuis : 0 -> 2
783 
784  Sym_tensor resu = hh_point - hh().derive_lie(beta())
785  - 0.6666666666666666 * beta().divergence(ff) * hh()
786  + beta().ope_killing_conf(ff) ;
787 
788  resu = psi4()*psi()*psi() * resu / (2*nn()) ;
789 
790  hata_evol.update(resu, jtime, the_time[jtime]) ;
791 
792  }
793 
794  return hata_evol[jtime] ;
795 
796 }
797 
798 
799 const Scalar& Time_slice_conf::trk() const {
800 
801  if( !(trk_evol.is_known(jtime)) ) {
802 
803  psi() ;
804  Scalar resu = -6*psi_evol.time_derive(jtime, scheme_order) / psi() ;
805  resu.inc_dzpuis(2) ;
806  resu += beta().divergence(ff)
807  + 6.*contract( beta(), 0, ln_psi().derive_cov(ff), 0 ) ;
808  resu = resu / nn() ;
809 
810  trk_evol.update(resu, jtime, the_time[jtime]) ;
811  }
812 
813  return trk_evol[jtime] ;
814 
815 }
816 
817 
819 
820  if (p_hdirac == 0x0) {
821  p_hdirac = new Vector( hh().divergence(ff) ) ;
822  }
823 
824  return *p_hdirac ;
825 
826 }
827 
828 const Vector& Time_slice_conf::vec_X(int methode_poisson) const {
829 
830  if (p_vec_X == 0x0) {
831  assert ( hata_evol.is_known(jtime) ) ;
832  Vector source = hata_evol[jtime].divergence(ff) ;
833  source.inc_dzpuis() ; // dzpuis 3-> 4
834  p_vec_X = new Vector( source.poisson( 1./3., methode_poisson ) ) ;
835  }
836  return *p_vec_X ;
837 }
838 
840 (const Vector& hat_S, const Sym_tensor_tt& hata_tt, int iter_max, double precis,
841  double relax, int meth_poisson) {
842 
843  // Some checks
844  assert( hat_S.get_index_type(0) == CON ) ;
845 #ifndef NDEBUG
846  for (int i=1; i<=3; i++)
847  for (int j=i; j<=3; j++)
848  assert( hata_tt(i,j).get_dzpuis() == 2 ) ;
849 #endif
850  assert( hh_evol.is_known(jtime) ) ;
851 
852  // Initializations
853  //----------------
854  const Tensor_sym& delta = tgam().connect().get_delta() ;
855  if (p_vec_X != 0x0) {
856  delete p_vec_X ;
857  p_vec_X = 0x0 ;
858  }
859 
860  // Constant part of the source
861  Vector source = hat_S - contract( delta, 1, 2, hata_tt, 0, 1 )
862  - 2./3.*psi4()*psi()*psi()*contract( tgam().con(), 0, trk().derive_cov(ff), 0 );
863 
864  p_vec_X = new Vector( source.poisson( 1./3., meth_poisson ) ) ;
865 
866  // Iteration on the vector X
867  //--------------------------
868  for (int it=0; it<iter_max; it++) {
869 
870  Vector fuente = source
871  - contract( delta, 1, 2, p_vec_X->ope_killing_conf(ff), 0, 1 ) ;
872 
873  Vector X_new = fuente.poisson( 1./3., meth_poisson) ;
874 
875  // Control of the convergence
876  double diff = 0. ;
877  for (int i=1; i<=3; i++)
878  diff += max( max( abs( X_new(i) - (*p_vec_X)(i) ) ) ) ;
879 
880  // Relaxation
881  (*p_vec_X) = relax*X_new + ( 1. - relax )*(*p_vec_X) ;
882 
883  // If converged, gets out of the loop
884  if (diff < precis) break ;
885  }
886 
887  // Update of \hat{A}^{ij} and reset of related quantities
888  hata_evol.update( hata_tt + p_vec_X->ope_killing_conf(ff), jtime, the_time[jtime] ) ;
889 
890  k_dd_evol.downdate(jtime) ;
891  k_uu_evol.downdate(jtime) ;
892 }
893 
894 void Time_slice_conf::set_AB_hata(const Scalar& A_in, const Scalar& B_in) {
895 
896  A_hata_evol.update(A_in, jtime, the_time[jtime]) ;
897  B_hata_evol.update(B_in, jtime, the_time[jtime]) ;
898 
899  hata_evol.downdate(jtime) ;
900  // Reset of quantities depending on A^{ij}:
901  k_dd_evol.downdate(jtime) ;
902  k_uu_evol.downdate(jtime) ;
903 
904 }
905 
906 
907 void Time_slice_conf::check_psi_dot(Tbl& tlnpsi_dot, Tbl& tdiff, Tbl& tdiff_rel) const {
908 
909  // Computation of d/dt ln Psi :
910 
911  Scalar lnpsi_dot(psi().get_mp()) ;
912  if ( psi_evol.is_known(jtime-1) ) {
913  lnpsi_dot = psi_evol.time_derive(jtime, scheme_order) / psi() ;
914  }
915  else {
916  lnpsi_dot = npsi_evol.time_derive(jtime, scheme_order) / npsi()
917  - n_evol.time_derive(jtime, scheme_order) / nn() ;
918  }
919 
920  tlnpsi_dot = max(abs(lnpsi_dot)) ;
921 
922  // Error on the d/dt ln Psi relation :
923 
924  Scalar diff = contract(beta(),0, ln_psi().derive_cov(ff),0)
925  + 0.1666666666666666 * ( beta().divergence(ff) - nn() * trk() ) ;
926 
927  diff.dec_dzpuis(2) ;
928 
929  Tbl tref = max(abs(diff)) + tlnpsi_dot ;
930 
931  diff -= lnpsi_dot ;
932 
933  tdiff = max(abs(diff)) ;
934 
935  tdiff_rel = tdiff / tref ;
936 
937 }
938 
939 
940  //------------------//
941  // output //
942  //------------------//
943 
944 ostream& Time_slice_conf::operator>>(ostream& flux) const {
945 
946  Time_slice::operator>>(flux) ;
947 
948  flux << "Triad on which the components of the flat metric are defined:\n"
949  << *(ff.get_triad()) << '\n' ;
950 
951  if (psi_evol.is_known(jtime)) {
952  maxabs( psi_evol[jtime], "Psi", flux) ;
953  }
954  if (npsi_evol.is_known(jtime)) {
955  maxabs( npsi_evol[jtime], "N Psi", flux) ;
956  }
957  if (hh_evol.is_known(jtime)) {
958  maxabs( hh_evol[jtime], "h^{ij}", flux) ;
959  }
960  if (hata_evol.is_known(jtime)) {
961  maxabs( hata_evol[jtime], "hat{A}^{ij}", flux) ;
962  }
963 
964  if (p_tgamma != 0x0) flux <<
965  "Conformal metric tilde gamma is up to date" << endl ;
966  if (p_psi4 != 0x0) maxabs( *p_psi4, "Psi^4", flux) ;
967  if (p_ln_psi != 0x0) maxabs( *p_ln_psi, "ln(Psi)", flux) ;
968  if (p_hdirac != 0x0) maxabs( *p_hdirac, "H^i", flux) ;
969  if (p_vec_X != 0x0) maxabs( *p_vec_X, "X^i", flux) ;
970 
971  return flux ;
972 
973 }
974 
975 
976 
977 void Time_slice_conf::sauve(FILE* fich, bool partial_save) const {
978 
979  // Writing of quantities common to all derived classes of Time_slice
980  // -----------------------------------------------------------------
981 
982  Time_slice::sauve(fich, true) ;
983 
984  // Writing of quantities common to all derived classes of Time_slice_conf
985  // ----------------------------------------------------------------------
986 
987  int jmin = jtime - depth + 1 ;
988 
989  // Psi
990  psi() ; // forces the update at the current time step
991  for (int j=jmin; j<=jtime; j++) {
992  int indicator = (psi_evol.is_known(j)) ? 1 : 0 ;
993  fwrite_be(&indicator, sizeof(int), 1, fich) ;
994  if (indicator == 1) psi_evol[j].sauve(fich) ;
995  }
996 
997  // hat{A}
998  hata() ;
999  for (int j=jmin; j<=jtime; j++) {
1000  int indicator = ( hata_evol.is_known(j) ? 1 : 0 ) ;
1001  fwrite_be(&indicator, sizeof(int), 1, fich) ;
1002  if (indicator == 1) hata_evol[j].sauve(fich) ;
1003  }
1004 
1005  //A and B...
1006  A_hata() ;
1007  for (int j=jmin; j<=jtime; j++) {
1008  int indicator = (A_hata_evol.is_known(j)) ? 1 : 0 ;
1009  fwrite_be(&indicator, sizeof(int), 1, fich) ;
1010  if (indicator == 1) A_hata_evol[j].sauve(fich) ;
1011  }
1012 
1013  B_hata() ;
1014  for (int j=jmin; j<=jtime; j++) {
1015  int indicator = (B_hata_evol.is_known(j)) ? 1 : 0 ;
1016  fwrite_be(&indicator, sizeof(int), 1, fich) ;
1017  if (indicator == 1) B_hata_evol[j].sauve(fich) ;
1018  }
1019 
1020  // Case of a complete save
1021  // -----------------------
1022  if (!partial_save) {
1023  cout << "Time_slice_conf::sauve: the full writing is not ready yet !"
1024  << endl ;
1025  abort() ;
1026  }
1027 
1028 }
1029 }
void operator=(const Time_slice &)
Assignment to another Time_slice.
Definition: time_slice.C:391
virtual const Vector & beta() const
shift vector at the current time step (jtime )
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
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.
const Scalar & psi4() const
Factor at the current time step (jtime ).
Metric for tensor calculation.
Definition: metric.h:90
virtual void set_npsi_del_psi(const Scalar &npsi_in)
Sets the factor at the current time step (jtime ) and deletes the value of .
virtual void set_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
Definition: time_slice.h:525
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
void operator=(const Time_slice_conf &)
Assignment to another Time_slice_conf.
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
Definition: time_slice.h:520
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
virtual const Scalar & determinant() const
Returns the determinant.
Definition: metric_flat.C:217
virtual void del_deriv() const
Deletes all the derived quantities.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Scalar * p_ln_psi
Pointer on the logarithm of at the current time step (jtime)
Definition: time_slice.h:569
Lorene prototypes.
Definition: app_hor.h:67
Evolution_full< Tbl > adm_mass_evol
ADM mass at each time step, since the creation of the slice.
Definition: time_slice.h:236
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Flat metric for tensor calculation.
Definition: metric.h:261
int scheme_order
Order of the finite-differences scheme for the computation of time derivatives.
Definition: time_slice.h:190
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition: time_slice.h:550
Base class for coordinate mappings.
Definition: map.h:688
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
int jtime
Time step index of the latest slice.
Definition: time_slice.h:193
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
virtual void set_psi_del_n(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Evolution_std< Sym_tensor > gam_dd_evol
Values at successive time steps of the covariant components of the induced metric ...
Definition: time_slice.h:201
virtual void set_hh(const Sym_tensor &hh_in)
Sets the deviation of the conformal metric from the flat metric : .
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.
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.
Tensor field of valence 1.
Definition: vector.h:188
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
Vector * p_hdirac
Pointer on the vector (which vanishes in Dirac gauge), at the current time step (jtime).
Definition: time_slice.h:574
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition: time_slice.h:555
virtual const Sym_tensor & gam_dd() const
Induced metric (covariant components ) at the current time step (jtime )
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:216
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
void set_A_tildeB(const Scalar &a_in, const Scalar &tb_in, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A and .
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime ) ...
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: ...
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
Spacelike time slice of a 3+1 spacetime.
Definition: time_slice.h:176
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
virtual const Scalar & A_hata() const
Returns the potential A of .
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:899
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
const Scalar & compute_A(bool output_ylm=true, Param *par=0x0) const
Gives the field A (see member p_aaa ).
virtual void set_hata_TT(const Sym_tensor_tt &hata_tt)
Sets the TT part of (see member hata_evol ).
Evolution_std< Sym_tensor > gam_uu_evol
Values at successive time steps of the contravariant components of the induced metric ...
Definition: time_slice.h:206
virtual const Vector & vec_X(int method_poisson=6) const
Vector representing the longitudinal part of .
Parameter storage.
Definition: param.h:125
virtual ~Time_slice_conf()
Destructor.
void set_etat_one()
Sets the logical state to ETATUN (one).
Definition: scalar.C:340
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:211
const Metric & gam() const
Induced metric at the current time step (jtime )
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
virtual void set_hata_from_XAB(Param *par_bc=0x0, Param *par_mat=0x0)
Sets the conformal representation of the traceless part of the extrinsic curvature from its potentia...
Spacelike time slice of a 3+1 spacetime with conformal decomposition.
Definition: time_slice.h:501
Scalar * p_psi4
Pointer on the factor at the current time step (jtime)
Definition: time_slice.h:566
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:73
Metric * p_gamma
Pointer on the induced metric at the current time step (jtime)
Definition: time_slice.h:242
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Definition: time_slice.C:414
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
Time_slice_conf(const Scalar &lapse_in, const Vector &shift_in, const Metric_flat &ff_in, const Scalar &psi_in, const Sym_tensor &hh_in, const Sym_tensor &hata_in, const Scalar &trk_in, int depth_in=3)
Constructor from conformal decomposition.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Metric * p_tgamma
Pointer on the conformal metric at the current time step (jtime)
Definition: time_slice.h:563
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ) ...
virtual const Scalar & determinant() const
Returns the determinant.
Definition: metric.C:395
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
Evolution_std< double > the_time
Time label of each slice.
Definition: time_slice.h:196
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime ) ...
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:72
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Vector * p_vec_X
Pointer on the vector representing the longitudinal part of .
Definition: time_slice.h:580
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
Scalar compute_tilde_B_tt(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_b ) associated with the TT-part of the Sym_tensor ...
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition: time_slice.h:227
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition: time_slice.h:219
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: time_slice.C:370
Basic array class.
Definition: tbl.h:164
virtual const Scalar & B_hata() const
Returns the potential of .
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050
void check_psi_dot(Tbl &tlnpsi_dot, Tbl &tdiff, Tbl &tdiff_rel) const
Checks the relation.
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition: time_slice.h:545
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition: time_slice.h:222
virtual void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
int depth
Number of stored time slices.
Definition: time_slice.h:182
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...
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition: time_slice.h:533
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:363
virtual void set_AB_hata(const Scalar &A_in, const Scalar &B_in)
Sets the potentials A and of the TT part (see the documentation of Sym_tensor for details)...
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
Transverse and traceless symmetric tensors of rank 2.
Definition: sym_tensor.h:933
virtual void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
Definition: time_slice.C:510
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )