LORENE
sources_hor.C
1  /*
2  * Methods of class Iso_hor to compute sources for Psi, N y beta
3  *
4  * (see file hor_isol.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004-2005 Jose Luis Jaramillo
10  * Francois Limousin
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: sources_hor.C,v 1.17 2016/12/05 16:17:56 j_novak Exp $
33  * $Log: sources_hor.C,v $
34  * Revision 1.17 2016/12/05 16:17:56 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.16 2014/10/13 08:53:01 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.15 2014/10/06 15:13:11 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.14 2005/06/09 08:05:32 f_limousin
44  * Implement a new function sol_elliptic_boundary() and
45  * Vector::poisson_boundary(...) which solve the vectorial poisson
46  * equation (method 6) with an inner boundary condition.
47  *
48  * Revision 1.13 2005/04/02 15:49:21 f_limousin
49  * New choice (Lichnerowicz) for aaquad. New member data nz.
50  *
51  * Revision 1.12 2005/03/28 19:42:39 f_limousin
52  * Implement the metric and A^{ij}A_{ij} of Sergio for pertubations
53  * of Kerr black holes.
54  *
55  * Revision 1.11 2005/03/22 13:25:36 f_limousin
56  * Small changes. The angular velocity and A^{ij} are computed
57  * with a differnet sign.
58  *
59  * Revision 1.10 2005/03/03 10:11:57 f_limousin
60  * The function source_psi(), source_nn() and source_beta() are
61  * now const and return an object const.
62  *
63  * Revision 1.9 2005/02/07 10:37:21 f_limousin
64  * Minor changes.
65  *
66  * Revision 1.8 2004/12/31 15:35:18 f_limousin
67  * Modifications to avoid warnings.
68  *
69  * Revision 1.7 2004/12/22 18:16:43 f_limousin
70  * Many different changes.
71  *
72  * Revision 1.6 2004/11/05 10:59:07 f_limousin
73  * Delete ener_dens, mom_dens and trace stress in functions
74  * source_nn, source_psi and source_beta.
75  * And some modification to avoid warnings (source_nn change to source...).
76  *
77  * Revision 1.5 2004/11/03 17:16:44 f_limousin
78  * Delete argument trk_point for source_nn()
79  *
80  * Revision 1.4 2004/10/29 15:45:08 jl_jaramillo
81  * Change name of functions
82  *
83  * Revision 1.3 2004/09/28 16:08:26 f_limousin
84  * Minor modifs.
85  *
86  * Revision 1.2 2004/09/09 16:55:08 f_limousin
87  * Add the 2 lines $Id: sources_hor.C,v 1.17 2016/12/05 16:17:56 j_novak Exp $Log: for CVS
88  *
89  *
90  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/sources_hor.C,v 1.17 2016/12/05 16:17:56 j_novak Exp $
91  *
92  */
93 
94 // C++ headers
95 #include "headcpp.h"
96 
97 // C headers
98 #include <cstdlib>
99 #include <cassert>
100 
101 // Lorene headers
102 #include "time_slice.h"
103 #include "isol_hor.h"
104 #include "metric.h"
105 #include "evolution.h"
106 #include "unites.h"
107 #include "graphique.h"
108 #include "utilitaires.h"
109 
110 namespace Lorene {
112 
113  using namespace Unites ;
114 
115  // Initialisations
116  // ---------------
117 
118  const Base_vect& triad = *(ff.get_triad()) ;
119 
120  Scalar tmp(mp) ;
121  Scalar tmp_scal(mp) ;
122  Sym_tensor tmp_sym(mp, CON, triad) ;
123 
124  Scalar source(mp) ;
125 
126  //===============================================
127  // Computations of the source for Psi
128  //===============================================
129 
130  const Vector& d_psi = psi().derive_cov(ff) ; // D_i Psi
131 
132  // Source for Psi
133  // --------------
134  tmp = 0.125* psi() * met_gamt.ricci_scal()
135  - contract(hh(), 0, 1, d_psi.derive_cov(ff), 0, 1 ) ;
136  tmp.inc_dzpuis() ; // dzpuis : 3 -> 4
137 
138  tmp -= contract(hdirac(), 0, d_psi, 0) ;
139 
140  source = tmp - 0.125* aa_quad() / psi4() / psi()/ psi()/ psi()
141  + psi()*psi4() * 8.33333333333333e-2* trK*trK ; //##
142  source.annule_domain(0) ;
143 
144  return source ;
145 }
146 
147 
149 
150  using namespace Unites ;
151 
152  // Initialisations
153  // ---------------
154 
155  const Base_vect& triad = *(ff.get_triad()) ;
156 
157  Scalar tmp(mp) ;
158  Scalar tmp_scal(mp) ;
159  Sym_tensor tmp_sym(mp, CON, triad) ;
160 
161  Scalar source(mp) ;
162 
163  //===============================================
164  // Computations of the source for NN
165  //===============================================
166 
167  const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
168  const Vector& dnnn = nn().derive_cov(ff) ; // D_i N
169 
170  // Source for N
171  // ------------
172 
173  source = aa_quad() / psi4() / psi4() * nn() +
174  psi4()*( nn()* 0.3333333333333333* trK*trK - trK_point )
175  - 2.* contract(dln_psi, 0, nn().derive_con(met_gamt), 0)
176  - contract(hdirac(), 0, dnnn, 0) ;
177 
178  tmp = psi4()* contract(beta(), 0, trK.derive_cov(ff), 0)
179  - contract( hh(), 0, 1, dnnn.derive_cov(ff), 0, 1 ) ;
180 
181  tmp.inc_dzpuis() ; // dzpuis: 3 -> 4
182 
183  source += tmp ;
184 
185  source.annule_domain(0) ;
186 
187  return source ;
188 
189 }
190 
191 
192 
194 
195  using namespace Unites ;
196 
197  // Initialisations
198  // ---------------
199 
200  const Base_vect& triad = *(ff.get_triad()) ;
201 
202  Scalar tmp(mp) ;
203  Scalar tmp_scal(mp) ;
204  Sym_tensor tmp_sym(mp, CON, triad) ;
205  Vector tmp_vect(mp, CON, triad) ;
206 
207  Vector source(mp, CON, triad) ;
208 
209  //===============================================
210  // Computations of the source for beta
211  //===============================================
212 
213  const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
214  const Vector& dnnn = nn().derive_cov(ff) ; // D_i N
215 
216  // Source for beta (dzpuis = 4)
217  // ----------------------------
218 
219  source = 2.* contract(aa(), 1,
220  dnnn - 6.*nn() * dln_psi, 0) ;
221 
222  tmp_vect = 0.66666666666666666* trK.derive_con(met_gamt) ;
223  tmp_vect.inc_dzpuis() ;
224 
225  source += 2.* nn() * ( tmp_vect
226  - contract(met_gamt.connect().get_delta(), 1, 2,
227  aa(), 0, 1) ) ;
228 
229  Vector vtmp = contract(hh(), 0, 1,
230  beta().derive_cov(ff).derive_cov(ff), 1, 2)
231  + 0.3333333333333333*
232  contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0)
233  - hdirac().derive_lie(beta())
234  + gamt_point.divergence(ff) ; // zero in the Dirac gauge
235  vtmp.inc_dzpuis() ; // dzpuis: 3 -> 4
236 
237  source -= vtmp ;
238 
239  source += 0.66666666666666666* beta().divergence(ff) * hdirac() ;
240 
241  source.annule_domain(0) ;
242 
243  /*
244  // Source for beta (dzpuis = 3)
245  // ----------------------------
246 
247  source = 2.* contract(aa(), 1,
248  dnnn - 6.*nn() * dln_psi, 0) ;
249  source.dec_dzpuis() ;
250 
251  tmp_vect = 0.66666666666666666* trK.derive_con(met_gamt) ;
252  Vector tmp_vect2 (- contract(met_gamt.connect().get_delta(), 1, 2,
253  aa(), 0, 1)) ;
254  tmp_vect2.dec_dzpuis() ;
255 
256  source += 2.* nn() * ( tmp_vect + tmp_vect2 ) ;
257 
258  Vector vtmp = contract(hh(), 0, 1,
259  beta().derive_cov(ff).derive_cov(ff), 1, 2)
260  + 0.3333333333333333*
261  contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0)
262  - hdirac().derive_lie(beta())
263  + gamt_point.divergence(ff) ; // zero in the Dirac gauge
264 
265  source -= vtmp ;
266 
267  tmp_vect = 0.66666666666666666* beta().divergence(ff) * hdirac() ;
268  tmp_vect.dec_dzpuis() ;
269  source += tmp_vect ;
270 
271  source.annule_domain(0) ;
272  */
273 
274  return source ;
275 
276 }
277 
278 
279 
280 
281 
282 
283 
284 
285 }
virtual const Vector & beta() const
shift vector at the current time step (jtime )
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
const Scalar & psi4() const
Factor at the current time step (jtime ).
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition: isol_hor.h:329
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
const Scalar source_psi() const
Source for .
Definition: sources_hor.C:111
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:332
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
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 : .
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
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 & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:353
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
virtual const Connection & connect() const
Returns the connection.
Definition: metric.C:304
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
const Scalar source_nn() const
Source for N.
Definition: sources_hor.C:148
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
const Vector source_beta() const
Source for .
Definition: sources_hor.C:193
virtual const Scalar & aa_quad() const
Conformal representation .
Definition: isol_hor.C:887
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
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:825
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition: isol_hor.h:335
Vector derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: vector.C:398
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )