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
28char 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
117namespace Lorene {
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)
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
161 A_hata() ;
162 B_hata() ;
163
164 set_der_0x0() ;
165
166}
167
168
169// Constructor from physical metric
170// --------------------------------
171
173 const Sym_tensor& gamma_in, const Sym_tensor& kk_in,
174 const Metric_flat& ff_in, int 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 ) ;
188 p_psi4 = new Scalar( pow( p_gamma->determinant() / ff.determinant(),
189 0.3333333333333333) ) ;
191
192 Scalar tmp = pow(*p_psi4, 0.25) ;
193 tmp.std_spectral_base() ;
195
196 hh_evol.update( (*p_psi4) * p_gamma->con() - ff.con(),
197 jtime, time_init ) ;
198
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() ;
226
227 // N Psi identically one:
229
230 // h^{ij} identically zero:
231 Sym_tensor stmp(mp, CON, triad) ;
232 stmp.set_etat_zero() ;
234
235 // \hat{A}^{ij} identically zero:
237
238 tmp.set_etat_zero() ;
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) ;
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) ;
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) ;
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) ;
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),
327 A_hata_evol(tin.A_hata_evol),
328 B_hata_evol(tin.B_hata_evol){
329
330 set_der_0x0() ;
331
332}
333 //--------------//
334 // Destructor //
335 //--------------//
336
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
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
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
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 }
425
426}
427
429
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 }
449
450}
451
452
454
456
457 // Reset of quantities depending on N Psi:
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 }
474
475}
476
477
479
481
482 // Reset of quantities depending on N Psi:
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
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 }
530
531}
532
533
535
537
538 // Reset of quantities depending on A^{ij}:
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
553
554 tmp = hata_tt.compute_tilde_B_tt(true) ;
555 if (tmp.get_dzpuis() == 3)
556 tmp.dec_dzpuis() ; // dzpuis 3->2
557
559
561 // Reset of quantities depending on A^{ij}:
564}
565
567
568 assert (p_vec_X != 0x0) ;
571
572 const Map& mp = p_vec_X->get_mp() ;
573
576 hata_tt.inc_dzpuis(2) ;
577
579
580 // Reset of quantities depending on A^{ij}:
583
584}
585
586
587 //-----------------------------------------------//
588 // Update of fields from base class Time_slice //
589 //-----------------------------------------------//
590
592
593 if (!( n_evol.is_known(jtime) ) ) {
594
597
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
669 Scalar tmp = hata_evol[jtime].compute_A(true) ;
670 if (tmp.get_dzpuis() == 3)
671 tmp.dec_dzpuis() ; // dzpuis 3->2
672
674 }
675 return A_hata_evol[jtime] ;
676}
677
679
680 if ( !( B_hata_evol.is_known(jtime) ) ) {
681
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
688 }
689 return A_hata_evol[jtime] ;
690}
691
692
694
695 if (!( psi_evol.is_known(jtime) ) ) {
696
699
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
738
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
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
777
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
788
789 }
790
791 return hata_evol[jtime] ;
792
793}
794
795
797
798 if( !(trk_evol.is_known(jtime)) ) {
799
800 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
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
826
827 if (p_vec_X == 0x0) {
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
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
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
886
889}
890
891void Time_slice_conf::set_AB_hata(const Scalar& A_in, const Scalar& B_in) {
892
895
897 // Reset of quantities depending on A^{ij}:
900
901}
902
903
904void 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 {
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
941ostream& Time_slice_conf::operator>>(ostream& flux) const {
942
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
974void 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}
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
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
TyT time_derive(int j, int n=2) const
Computes the time derivative at time step j by means of a n-th order scheme, from the values at steps...
Definition evolution.C:504
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 Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
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
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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
Transverse and traceless symmetric tensors of rank 2.
Definition sym_tensor.h:938
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
void dec_dzpuis()
dzpuis -= 1 ;
Definition tenseur.C:1104
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Spacelike time slice of a 3+1 spacetime with conformal decomposition.
Definition time_slice.h:498
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition time_slice.h:530
virtual ~Time_slice_conf()
Destructor.
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
virtual const Sym_tensor & gam_dd() const
Induced metric (covariant components ) at the current time step (jtime )
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
virtual void set_psi_del_n(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
virtual void set_hata_TT(const Sym_tensor_tt &hata_tt)
Sets the TT part of (see member hata_evol ).
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
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 const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual void set_psi_del_npsi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime )
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition time_slice.h:542
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
virtual const Scalar & A_hata() const
Returns the potential A of .
virtual const Vector & vec_X(int method_poisson=6) const
Vector representing the longitudinal part of .
void operator=(const Time_slice_conf &)
Assignment to another Time_slice_conf.
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Metric * p_tgamma
Pointer on the conformal metric at the current time step (jtime)
Definition time_slice.h:560
virtual const Scalar & B_hata() const
Returns the potential of .
virtual void del_deriv() const
Deletes all the derived quantities.
Evolution_std< Scalar > npsi_evol
Values at successive time steps of the factor .
Definition time_slice.h:522
Vector * p_vec_X
Pointer on the vector representing the longitudinal part of .
Definition time_slice.h:577
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: .
Vector * p_hdirac
Pointer on the vector (which vanishes in Dirac gauge), at the current time step (jtime).
Definition time_slice.h:571
const Scalar & psi4() const
Factor at the current time step (jtime ).
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.
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_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...
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
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.
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition time_slice.h:547
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition time_slice.h:507
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.
Scalar * p_ln_psi
Pointer on the logarithm of at the current time step (jtime)
Definition time_slice.h:566
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition time_slice.h:552
virtual void set_hh(const Sym_tensor &hh_in)
Sets the deviation of the conformal metric from the flat metric : .
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
Scalar * p_psi4
Pointer on the factor at the current time step (jtime)
Definition time_slice.h:563
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
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 void del_deriv() const
Deletes all the derived quantities.
Definition time_slice.C:367
const Metric & gam() const
Induced metric at the current time step (jtime )
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime )
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
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
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Lorene prototypes.
Definition app_hor.h:64