LORENE
tslice_dirac_max_setAB.C
1/*
2 * Methods of class Tslice_dirac_max dealing with the members potA and tildeB
3 *
4 * (see file time_slice.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2007 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char tslice_dirax_max_setAB_C[] = "$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_setAB.C,v 1.11 2014/10/13 08:53:48 j_novak Exp $" ;
29
30/*
31 * $Id: tslice_dirac_max_setAB.C,v 1.11 2014/10/13 08:53:48 j_novak Exp $
32 * $Log: tslice_dirac_max_setAB.C,v $
33 * Revision 1.11 2014/10/13 08:53:48 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.10 2014/10/06 15:13:22 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.9 2012/02/06 12:59:07 j_novak
40 * Correction of some errors.
41 *
42 * Revision 1.8 2011/07/22 13:21:02 j_novak
43 * Corrected an error on BC treatment.
44 *
45 * Revision 1.7 2010/10/20 07:58:10 j_novak
46 * Better implementation of the explicit time-integration. Not fully-tested yet.
47 *
48 * Revision 1.6 2008/12/04 18:22:49 j_novak
49 * Enhancement of the dzpuis treatment + various bug fixes.
50 *
51 * Revision 1.5 2008/12/02 15:02:22 j_novak
52 * Implementation of the new constrained formalism, following Cordero et al. 2009
53 * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
54 *
55 * Revision 1.4 2007/06/05 07:38:37 j_novak
56 * Better treatment of dzpuis for A and tilde(B) potentials. Some errors in the bases manipulation have been also corrected.
57 *
58 * Revision 1.3 2007/05/24 12:10:41 j_novak
59 * Update of khi_evol and mu_evol.
60 *
61 * Revision 1.2 2007/04/25 15:21:01 j_novak
62 * Corrected an error in the initialization of tildeB in
63 * Tslice_dirac_max::initial_dat_cts. + New method for solve_hij_AB.
64 *
65 * Revision 1.1 2007/03/21 14:51:50 j_novak
66 * Introduction of potentials A and tilde(B) of h^{ij} into Tslice_dirac_max.
67 *
68 *
69 * $Header $
70 *
71 */
72
73// C headers
74#include <cassert>
75
76// Lorene headers
77#include "time_slice.h"
78#include "param.h"
79#include "unites.h"
80#include "proto.h"
81#include "graphique.h"
82
83namespace Lorene {
85
88
89 // Reset of quantities depending on h^{ij}:
92 if (p_tgamma != 0x0) {
93 delete p_tgamma ;
94 p_tgamma = 0x0 ;
95 }
96 if (p_hdirac != 0x0) {
97 delete p_hdirac ;
98 p_hdirac = 0x0 ;
99 }
100 if (p_gamma != 0x0) {
101 delete p_gamma ;
102 p_gamma = 0x0 ;
103 }
107}
108
110
111 assert (A_hh_evol.is_known(j0)) ; // The starting point
112 assert (B_hh_evol.is_known(j0)) ; // of the computation
113
114 const Map& mp = A_hh_evol[j0].get_mp() ;
115
116 // The representation of h^{ij} as an object of class Sym_tensor_trans :
117 Sym_tensor_trans hij(mp, *(ff.get_triad()), ff) ;
118 const Scalar* ptrace = 0x0 ;
119 if (trh_evol.is_known(j0-1)) ptrace = &trh_evol[j0-1] ;
121
122 // Result set to trh_evol and hh_evol
123 // ----------------------------------
125
126 // The longitudinal part of h^{ij}, which is zero by virtue of Dirac gauge :
127 Vector wzero(mp, CON, *(ff.get_triad())) ;
128 wzero.set_etat_zero() ;
129
130 // Temporary Sym_tensor with longitudinal part set to zero :
131 Sym_tensor hh_new(mp, CON, *(ff.get_triad())) ;
132 hh_new.set_longit_trans(wzero, hij) ;
134
135 if (j0 == jtime) {
136 // Reset of quantities depending on h^{ij}:
137 if (p_tgamma != 0x0) {
138 delete p_tgamma ;
139 p_tgamma = 0x0 ;
140 }
141 if (p_hdirac != 0x0) {
142 delete p_hdirac ;
143 p_hdirac = 0x0 ;
144 }
145 if (p_gamma != 0x0) {
146 delete p_gamma ;
147 p_gamma = 0x0 ;
148 }
149 }
153
154#ifndef NDEBUG
155 // Test
156 if (j0 == jtime) {
157 maxabs(tgam().determinant() - 1,
158 "Max. of absolute value of deviation from det tgam = 1") ;
159 }
160 else {
161 Metric tgam_j0( ff.con() + hh_evol[j0] ) ;
162 maxabs(tgam_j0.determinant() - 1,
163 "Max. of absolute value of deviation from det tgam = 1") ;
164 }
165#endif
166}
167
169 const {
170
171 const Map& mp = hijtt.get_mp() ;
172
173 // The representation of h^{ij} as an object of class Sym_tensor_trans :
174 Sym_tensor_trans hij(mp, *(ff.get_triad()), ff) ;
175 const Scalar* ptrace = 0x0 ;
176 if ( trh_evol.is_known( jtime - 1 ) ) ptrace = &trh_evol[jtime-1] ;
178
179 // Result set to trh_evol and hh_evol
180 // ----------------------------------
182
183 // The longitudinal part of h^{ij}, which is zero by virtue of Dirac gauge :
184 Vector wzero(mp, CON, *(ff.get_triad())) ;
185 wzero.set_etat_zero() ;
186
187 // Temporary Sym_tensor with longitudinal part set to zero :
188 Sym_tensor hh_new(mp, CON, *(ff.get_triad())) ;
189 hh_new.set_longit_trans(wzero, hij) ;
191
192 // Reset of quantities depending on h^{ij}:
193 if (p_tgamma != 0x0) {
194 delete p_tgamma ;
195 p_tgamma = 0x0 ;
196 }
197 if (p_hdirac != 0x0) {
198 delete p_hdirac ;
199 p_hdirac = 0x0 ;
200 }
201 if (p_gamma != 0x0) {
202 delete p_gamma ;
203 p_gamma = 0x0 ;
204 }
208
209#ifndef NDEBUG
210 // Test
211 maxabs(tgam().determinant() - 1,
212 "Max. of absolute value of deviation from det tgam = 1") ;
213#endif
214}
215 //----------------------------------------------//
216 // Equations for h^{ij} and \hat{A}^{ij} //
217 //----------------------------------------------//
218
220
221 using namespace Unites ;
222
223 const Map& map = hh().get_mp() ;
224 const Base_vect& otriad = *hh().get_triad() ;
225 int nz = map.get_mg()->get_nzone() ;
226
227 Sym_tensor strain_tens(map, CON, otriad) ;
228 if (p_strain_tens != 0x0)
230 else
231 strain_tens.set_etat_zero() ;
232
233 Sym_tensor aij = aa() ;
234 aij.annule_domain(nz-1) ;
235
236 const Sym_tensor& tgam_dd = tgam().cov() ; // {\tilde \gamma}_{ij}
237 const Sym_tensor& tgam_uu = tgam().con() ; // {\tilde \gamma}^{ij}
238 const Tensor_sym& dtgam = tgam_dd.derive_cov(ff) ;// D_k {\tilde \gamma}_{ij}
239 const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
240 const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
241 const Vector& tdln_psi_u = ln_psi().derive_con(tgam()) ; // tD^i ln(Psi)
242 Scalar log_N = log(nn()) ;
243 log_N.std_spectral_base() ;
244 const Vector& tdlnn_u = log_N.derive_con(tgam()) ; // tD^i ln(N)
245 const Scalar& div_beta = beta().divergence(ff) ; // D_k beta^k
246 Scalar qq = nn()*psi()*psi() ;
247 qq.annule_domain(nz-1) ;
248 const Vector& dqq = qq.derive_cov(ff) ; // D_i Q
249 Sym_tensor a_hat = hata() ;
250 a_hat.annule_domain(nz-1) ;
251 Scalar psi6 = psi4()*psi()*psi() ;
252 Sym_tensor sym_tmp(map, CON, otriad) ;
253 Scalar tmp(map) ;
254
255 //==================================
256 // Source for hij
257 //==================================
258
259 Sym_tensor source_hij = hh().derive_lie(beta()) + 2*(nn()/psi6 - 1.)*a_hat
260 - beta().ope_killing_conf(ff) + 0.6666666666666667*div_beta*hh() ;
261 source_hij.annule_domain(nz-1) ;
262 for (int i=1; i<=3; i++)
263 for (int j=i; j<=3; j++)
264 source_hij.set( i, j ).set_dzpuis(0) ;
265
266 tmp = 2.*A_hata_evol[jtime] + source_hij.compute_A(true) ;
267 tmp.set_spectral_va().ylm() ;
268 tmp.annule_domain(nz-1) ;
269 tmp.set_dzpuis(0) ;
271
272 tmp = 2.*B_hata_evol[jtime] + source_hij.compute_tilde_B_tt(true) ;
273 tmp.set_spectral_va().ylm() ;
274 tmp.annule_domain(nz-1) ;
275 tmp.set_dzpuis(0) ;
277
278 //==================================
279 // Source for \hat{A}^{ij}
280 //==================================
281
282 Sym_tensor source_aij = a_hat.derive_lie(beta())
283 + 1.666666666666667*a_hat*div_beta ;
284
285 // Quadratic part of the Ricci tensor of gam_tilde
286 // ------------------------------------------------
287
288 Sym_tensor ricci_star(map, CON, otriad) ;
289
290 ricci_star = contract(hh(), 0, 1, dhh.derive_cov(ff), 2, 3) ;
291
292 ricci_star.inc_dzpuis() ; // dzpuis : 3 --> 4
293
294 for (int i=1; i<=3; i++) {
295 for (int j=1; j<=i; j++) {
296 tmp = 0 ;
297 for (int k=1; k<=3; k++) {
298 for (int l=1; l<=3; l++) {
299 tmp += dhh(i,k,l) * dhh(j,l,k) ;
300 }
301 }
302 ricci_star.set(i,j) -= tmp ;
303 }
304 }
305
306 for (int i=1; i<=3; i++) {
307 for (int j=1; j<=i; j++) {
308 tmp = 0 ;
309 for (int k=1; k<=3; k++) {
310 for (int l=1; l<=3; l++) {
311 for (int m=1; m<=3; m++) {
312 for (int n=1; n<=3; n++) {
313
314 tmp += 0.5 * tgam_uu(i,k)* tgam_uu(j,l)
315 * dhh(m,n,k) * dtgam(m,n,l)
316 + tgam_dd(n,l) * dhh(m,n,k)
317 * (tgam_uu(i,k) * dhh(j,l,m) + tgam_uu(j,k) * dhh(i,l,m) )
318 - tgam_dd(k,l) *tgam_uu(m,n) * dhh(i,k,m) * dhh(j,l,n) ;
319 }
320 }
321 }
322 }
323 sym_tmp.set(i,j) = tmp ;
324 }
325 }
326 ricci_star += sym_tmp ; // a factor 1/2 is still missing, shall be put later
327
328 // Curvature scalar of conformal metric :
329 // -------------------------------------
330
332 0.25 * contract(tgam_uu, 0, 1,
333 contract(dhh, 0, 1, dtgam, 0, 1), 0, 1 )
334 - 0.5 * contract(tgam_uu, 0, 1,
335 contract(dhh, 0, 1, dtgam, 0, 2), 0, 1 ) ;
336
337 Scalar lap_A = A_hh_evol[jtime].laplacian(2) ;
338 Scalar tilde_lap_B(map) ;
339 tilde_laplacian( B_hh_evol[jtime], tilde_lap_B) ;
340 Sym_tensor_tt laplace_h(map, otriad, ff) ;
341 laplace_h.set_A_tildeB(lap_A, tilde_lap_B) ;
342 laplace_h.annule_domain(nz-1) ;
343
344 // sym_tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
345
346 source_aij += (0.5*(qq - 1.))*laplace_h + qq*(0.5*ricci_star + 8.*tdln_psi_u * tdln_psi_u
347 + 4.*( tdln_psi_u * tdlnn_u + tdlnn_u * tdln_psi_u )
348 - 0.3333333333333333 * (tricci_scal + 8.*(contract(dln_psi, 0, tdln_psi_u, 0)
349 + contract(dln_psi, 0, tdlnn_u, 0) )
350 )*tgam_uu
351 ) ;
352
353 sym_tmp = contract(tgam_uu, 1, contract(tgam_uu, 1, dqq.derive_cov(ff), 0), 1) ;
354
355 for (int i=1; i<=3; i++) {
356 for (int j=1; j<=i; j++) {
357 tmp = 0 ;
358 for (int k=1; k<=3; k++) {
359 for (int l=1; l<=3; l++) {
360 tmp += ( tgam_uu(i,k)*dhh(l,j,k) + tgam_uu(k,j)*dhh(i,l,k)
361 - tgam_uu(k,l)*dhh(i,j,k) ) * dqq(l) ;
362 }
363 }
364 sym_tmp.set(i,j) += 0.5 * tmp ;
365 }
366 }
367
369 - ( 0.3333333333333333*qq.derive_con(tgam()).divergence(tgam()) ) *tgam_uu ;
370
371 for (int i=1; i<=3; i++) {
372 for (int j=1; j<=i; j++) {
373 tmp = 0 ;
374 for (int k=1; k<=3; k++) {
375 for (int l=1; l<=3; l++) {
376 tmp += tgam_dd(k,l) * a_hat(i,k) * aij(j,l) ;
377 }
378 }
379 sym_tmp.set(i,j) = tmp ;
380 }
381 }
382
383 tmp = psi4() * strain_tens.trace(tgam()) ; // S = S_i^i
384
385 source_aij += (2.*nn())
386 * (
387 sym_tmp - qpig*psi6*( psi4()* strain_tens - (0.3333333333333333 * tmp) * tgam_uu )
388 ) ;
389
390 source_aij.annule_domain(nz-1) ;
391 for (int i=1; i<=3; i++)
392 for (int j=i; j<=3; j++)
393 source_aij.set(i,j).set_dzpuis(0) ;
394#ifndef NDEBUG
395 maxabs(source_aij, "source_aij tot") ;
396#endif
397
398 tmp = 0.5*lap_A + source_aij.compute_A(true) ;
399 tmp.annule_domain(nz-1) ;
400 tmp.set_dzpuis(0) ;
402 tmp = 0.5*tilde_lap_B + source_aij.compute_tilde_B_tt(true) ;
403 tmp.annule_domain(nz-1) ;
404 tmp.set_dzpuis(0) ;
405 // Scalar dess = tmp - tilde_lap_B ;
406 // dess.set_spectral_va().ylm_i() ;
407 // des_profile(dess, 0, 8., 1, 1) ;
409}
410
411void Tslice_dirac_max::initialize_sources_copy() const {
412
417
422
427
428 int jtime1 = jtime - depth + 1;
429 for (int j=jtime1; j <= jtime; j++) {
434 }
435}
436
437}
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
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 Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Metric for tensor calculation.
Definition metric.h:90
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Transverse symmetric tensors of rank 2.
Definition sym_tensor.h:608
void set_AtBtt_det_one(const Scalar &a_in, const Scalar &tbtt_in, const Scalar *h_prev=0x0, Param *par_bc=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assigns the derived member A and computes from its TT-part (see Sym_tensor::compute_tilde_B_tt() ).
const Scalar & the_trace() const
Returns the trace of the tensor with respect to metric *met_div.
void set_tt_part_det_one(const Sym_tensor_tt &hijtt, const Scalar *h_prev=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assignes the TT-part of the tensor.
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
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition time_slice.h:530
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
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 Sym_tensor aa() const
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 const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Evolution_std< Scalar > A_hata_evol
Potential A associated with the symmetric tensor .
Definition time_slice.h:547
Evolution_std< Scalar > B_hata_evol
Potential associated with the symmetric tensor .
Definition time_slice.h:552
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Metric * p_gamma
Pointer on the induced metric at the current time step (jtime)
Definition time_slice.h:239
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
virtual const Vector & beta() const
shift vector at the current time step (jtime )
int depth
Number of stored time slices.
Definition time_slice.h:179
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
Evolution_std< Scalar > source_B_hh_evol
The potential of the source of equation for .
Definition time_slice.h:995
Evolution_std< Scalar > source_A_hh_evol
The A potential of the source of equation for .
Definition time_slice.h:989
Evolution_std< Scalar > A_hh_evol
The A potential of .
Definition time_slice.h:977
void hh_det_one(int j, Param *par_bc=0x0, Param *par_mat=0x0) const
Computes from the values of A and and using the condition , which fixes the trace of .
Evolution_std< Scalar > source_B_hata_evol
The potential of the source of equation for .
void compute_sources(const Sym_tensor *strain_tensor=0x0) const
Computes the sources source_A_XXX_evol and source_B_XXX_evol , for the solution of the evolution equa...
virtual void set_AB_hh(const Scalar &A_in, const Scalar &B_in)
Sets the potentials A and of the TT part of (see the documentation of Sym_tensor for details).
Evolution_std< Scalar > B_hh_evol
The potential of .
Definition time_slice.h:983
Evolution_std< Scalar > source_A_hata_evol
The potential A of the source of equation for .
Tensor field of valence 1.
Definition vector.h:188
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.