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