LORENE
tslice_dirac_max_solve.C
1/*
2 * Methods of class Tslice_dirac_max for solving Einstein equations
3 *
4 * (see file time_slice.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2004 Eric Gourgoulhon & 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_dirac_max_solve_C[] = "$Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_solve.C,v 1.18 2014/10/13 08:53:48 j_novak Exp $" ;
29
30/*
31 * $Id: tslice_dirac_max_solve.C,v 1.18 2014/10/13 08:53:48 j_novak Exp $
32 * $Log: tslice_dirac_max_solve.C,v $
33 * Revision 1.18 2014/10/13 08:53:48 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.17 2014/10/06 15:13:22 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.16 2010/10/20 07:58:10 j_novak
40 * Better implementation of the explicit time-integration. Not fully-tested yet.
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/07/08 12:29:01 j_novak
47 * use of new method Tensor::annule_extern_cn
48 *
49 * Revision 1.13 2004/06/30 08:02:40 j_novak
50 * Added filtering in l of khi_new and mu_new. ki_source is forced to go to
51 * zero at least as r^2.
52 *
53 * Revision 1.12 2004/06/17 07:07:11 e_gourgoulhon
54 * Method solve_hij:
55 * -- replaced the attenuation of khi_source with tempo by a call to
56 * the new method Tensor::annule_extern_c2
57 * -- suppressed filtre_r on khi_source and khi_new
58 * -- added graphs of W^i and LW^{ij} (transverse decomp of S^ij).
59 *
60 * Revision 1.11 2004/06/15 09:43:36 j_novak
61 * Attenuation of the source for khi in the last shell (temporary?).
62 *
63 * Revision 1.10 2004/06/14 20:47:31 e_gourgoulhon
64 * Added argument method_poisson to method solve_hij.
65 *
66 * Revision 1.9 2004/06/03 10:02:45 j_novak
67 * Some filtering is done on source_khi and khi_new.
68 *
69 * Revision 1.8 2004/05/24 21:00:44 e_gourgoulhon
70 * Method solve_hij: khi and mu.smooth_decay(2,2) --> smooth_decay(2,1) ;
71 * added exponential_decay(khi) and exponential_decay(mu) after the
72 * call to smooth_decay. Method exponential_decay is provisory defined
73 * in this file.
74 *
75 * Revision 1.7 2004/05/17 19:56:25 e_gourgoulhon
76 * -- Method solve_beta: added argument method
77 * -- Method solve_hij: added argument graph_device
78 *
79 * Revision 1.6 2004/05/12 15:24:20 e_gourgoulhon
80 * Reorganized the #include 's, taking into account that
81 * time_slice.h contains now an #include "metric.h".
82 *
83 * Revision 1.5 2004/05/05 14:47:05 e_gourgoulhon
84 * Modified text and graphical outputs.
85 *
86 * Revision 1.4 2004/05/03 15:06:27 e_gourgoulhon
87 * Added matter source in solve_hij.
88 *
89 * Revision 1.3 2004/05/03 14:50:00 e_gourgoulhon
90 * Finished the implementation of method solve_hij.
91 *
92 * Revision 1.2 2004/04/30 14:36:15 j_novak
93 * Added the method Tslice_dirac_max::solve_hij(...)
94 * NOT READY YET!!!
95 *
96 * Revision 1.1 2004/04/30 10:52:14 e_gourgoulhon
97 * First version.
98 *
99 *
100 * $Header: /cvsroot/Lorene/C++/Source/Time_slice/tslice_dirac_max_solve.C,v 1.18 2014/10/13 08:53:48 j_novak Exp $
101 *
102 */
103
104// C headers
105#include <cassert>
106
107// Lorene headers
108#include "time_slice.h"
109#include "unites.h"
110#include "graphique.h"
111#include "proto.h"
112
113 //----------------------------//
114 // Equation for N\Psi //
115 //----------------------------//
116
117namespace Lorene {
119 const Scalar* p_trace_stress) const {
120
121 using namespace Unites ;
122
123 const Map& map = npsi().get_mp() ;
124 Scalar ener_dens(map) ;
125 if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
126 else ener_dens.set_etat_zero() ;
127
128 Scalar trace_stress(map) ;
129 if (p_trace_stress != 0x0) trace_stress = *(p_trace_stress) ;
130 else trace_stress.set_etat_zero() ;
131
132 // Source for N\Psi
133 // ----------------
134
135 const Vector& dnpsi = npsi().derive_cov(ff) ; // D_i N\Psi
136 const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
137 const Tensor_sym& dtgam = tgam().cov().derive_cov(ff) ;
138 // D_k {\tilde \gamma}_{ij}
139 Sym_tensor taa = aa().up_down(tgam()) ;
140
141 Scalar aa_quad = contract(taa, 0, 1, aa(), 0, 1) ;
142 Scalar tildeR = 0.25 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgam())
143 - 0.5 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgam()) ;
144
145 Scalar source_npsi = - contract( hh(), 0, 1, dnpsi.derive_cov(ff), 0, 1 ) ;
146 source_npsi.inc_dzpuis() ;
147 source_npsi += npsi()*( 0.5*qpig*(ener_dens + 2.*trace_stress )/( psi()*psi() )
148 + 0.875*aa_quad/( psi4()*psi4() ) + 0.125*tildeR ) ;
149
150 // Resolution of the Poisson equation for N\Psi
151 // --------------------------------------------
152
153 Scalar npsi_new = source_npsi.poisson() + 1. ;
154
155 if (npsi_new.get_etat() == ETATUN) npsi_new.std_spectral_base() ;
156
157#ifndef NDEBUG
158 // Test:
159 maxabs(npsi_new.laplacian() - source_npsi,
160 "Absolute error in the resolution of the equation for N") ;
161#endif
162 return npsi_new ;
163
164}
165
166
167 //--------------------------//
168 // Equation for \Psi //
169 //--------------------------//
170
171
173
174 using namespace Unites ;
175
176 const Map& map = psi().get_mp() ;
177 Scalar ener_dens(map) ;
178 if (p_ener_dens != 0x0) ener_dens = *(p_ener_dens) ;
179 else ener_dens.set_etat_zero() ;
180
181 // Source for \Psi
182 // ---------------
183
184 const Vector& dpsi = psi().derive_cov(ff) ; // D_i Psi
185 const Tensor_sym& dhh = hh().derive_cov(ff) ; // D_k h^{ij}
186 const Tensor_sym& dtgam = tgam().cov().derive_cov(ff) ;
187 // D_k {\tilde \gamma}_{ij}
188
189 Sym_tensor taa = hata().up_down(tgam()) ;
190
191 Scalar aa_quad = contract(taa, 0, 1, hata(), 0, 1) ;
192
193 Scalar tildeR = 0.25 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgam())
194 - 0.5 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgam()) ;
195
196 Scalar source_psi = -contract( hh(), 0, 1, dpsi.derive_cov(ff), 0, 1 ) ;
197 source_psi.inc_dzpuis() ;
198 source_psi -= 0.5*qpig*ener_dens/psi()
199 + 0.125*( aa_quad*pow(psi(), -7) - tildeR*psi() ) ;
200
201 // Resolution of the Poisson equation for Psi
202 // ------------------------------------------
203
204 Scalar psi_new = source_psi.poisson() + 1. ;
205
206 if (psi_new.get_etat() == ETATUN) psi_new.std_spectral_base() ;
207
208#ifndef NDEBUG
209 // Test:
210 maxabs(psi_new.laplacian() - source_psi,
211 "Absolute error in the resolution of the equation for Psi") ;
212#endif
213
214 return psi_new ;
215}
216
217
218
219 //--------------------------//
220 // Equation for beta //
221 //--------------------------//
222
223
225 const {
226
227 // Source for beta
228 // ---------------
229 Vector source_beta =
230 - contract(hh(), 0, 1, beta().derive_cov(ff).derive_cov(ff), 1, 2)
231 - 0.3333333333333333*contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0) ;
232
233 Sym_tensor sym_tmp = 2*nn()*aa() ;
234 source_beta += sym_tmp.divergence(ff) ;
235 source_beta.inc_dzpuis() ; // dzpuis: 3 -> 4
236
237 // Resolution of the vector Poisson equation
238 //------------------------------------------
239
240 Vector beta_new = source_beta.poisson(0.3333333333333333, ff, method) ;
241
242#ifndef NDEBUG
243 // Test:
244 Vector test_beta = (beta_new.derive_con(ff)).divergence(ff)
245 + 0.3333333333333333 * (beta_new.divergence(ff)).derive_con(ff) ;
246 test_beta.inc_dzpuis() ;
247 maxabs(test_beta - source_beta,
248 "Absolute error in the resolution of the equation for beta") ;
249#endif
250 return beta_new ;
251
252}
253
254
255}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Base class for coordinate mappings.
Definition map.h:670
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition scalar_pde.C:136
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...
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
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
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: .
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
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 ).
virtual const Vector & beta() const
shift vector at the current time step (jtime )
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $\Psi$ (Hamiltonian constraint).
virtual Scalar solve_npsi(const Scalar *ener_dens=0x0, const Scalar *trace_stress=0x0) const
Solves the elliptic equation for (maximal slicing condition + Hamiltonian constraint)
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
Tensor field of valence 1.
Definition vector.h:188
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
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.