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 
28 char 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 
117 namespace 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 
172 Scalar Tslice_dirac_max::solve_psi(const Scalar* p_ener_dens) const {
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 }
Lorene::Tensor_sym::derive_cov
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor_sym_calculus.C:192
Lorene::Scalar::inc_dzpuis
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...
Definition: scalar_r_manip.C:470
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Lorene::Time_slice_conf::npsi
virtual const Scalar & npsi() const
Factor at the current time step (jtime ).
Definition: time_slice_conf.C:732
Lorene::Time_slice_conf::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Definition: time_slice_conf.C:591
Lorene::Tensor_sym
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1037
Lorene::Time_slice_conf::hata
virtual const Sym_tensor & hata() const
Conformal representation of the traceless part of the extrinsic curvature: .
Definition: time_slice_conf.C:772
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Lorene::Tensor::up_down
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Definition: tensor_calculus.C:305
Lorene::Vector::poisson
Vector poisson(double lambda, int method=6) const
Solves the vector Poisson equation with *this as a source.
Definition: vector_poisson.C:518
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:861
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Tslice_dirac_max::solve_npsi
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)
Definition: tslice_dirac_max_solve.C:118
Unites
Standard units of space, time and mass.
Lorene::Tslice_dirac_max::solve_psi
virtual Scalar solve_psi(const Scalar *ener_dens=0x0) const
Solves the elliptic equation for the conformal factor $\Psi$ (Hamiltonian constraint).
Definition: tslice_dirac_max_solve.C:172
Lorene::Time_slice_conf::psi
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Definition: time_slice_conf.C:693
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
Lorene::maxabs
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.
Definition: tensor_calculus_ext.C:606
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Time_slice_conf::psi4
const Scalar & psi4() const
Factor at the current time step (jtime ).
Definition: time_slice_conf.C:707
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Tensor::inc_dzpuis
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
Lorene::Scalar::std_spectral_base
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
Lorene::Tslice_dirac_max::solve_beta
virtual Vector solve_beta(int method=6) const
Solves the elliptic equation for the shift vector from (Eq.
Definition: tslice_dirac_max_solve.C:224
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition: time_slice_access.C:87
Lorene::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
Lorene::Sym_tensor::divergence
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Lorene::Tslice_dirac_max::hh
virtual const Sym_tensor & hh(Param *par_bc=0x0, Param *par_mat=0x0) const
Deviation of the conformal metric from the flat metric : .
Definition: tslice_dirac_max.C:474
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::Time_slice_conf::aa
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
Definition: time_slice_conf.C:765
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Time_slice_conf::tgam
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Definition: time_slice_conf.C:747