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 char source_hor_C[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/sources_hor.C,v 1.16 2014/10/13 08:53:01 j_novak Exp $" ;
30 
31 /*
32  * $Id: sources_hor.C,v 1.16 2014/10/13 08:53:01 j_novak Exp $
33  * $Log: sources_hor.C,v $
34  * Revision 1.16 2014/10/13 08:53:01 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.15 2014/10/06 15:13:11 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.14 2005/06/09 08:05:32 f_limousin
41  * Implement a new function sol_elliptic_boundary() and
42  * Vector::poisson_boundary(...) which solve the vectorial poisson
43  * equation (method 6) with an inner boundary condition.
44  *
45  * Revision 1.13 2005/04/02 15:49:21 f_limousin
46  * New choice (Lichnerowicz) for aaquad. New member data nz.
47  *
48  * Revision 1.12 2005/03/28 19:42:39 f_limousin
49  * Implement the metric and A^{ij}A_{ij} of Sergio for pertubations
50  * of Kerr black holes.
51  *
52  * Revision 1.11 2005/03/22 13:25:36 f_limousin
53  * Small changes. The angular velocity and A^{ij} are computed
54  * with a differnet sign.
55  *
56  * Revision 1.10 2005/03/03 10:11:57 f_limousin
57  * The function source_psi(), source_nn() and source_beta() are
58  * now const and return an object const.
59  *
60  * Revision 1.9 2005/02/07 10:37:21 f_limousin
61  * Minor changes.
62  *
63  * Revision 1.8 2004/12/31 15:35:18 f_limousin
64  * Modifications to avoid warnings.
65  *
66  * Revision 1.7 2004/12/22 18:16:43 f_limousin
67  * Many different changes.
68  *
69  * Revision 1.6 2004/11/05 10:59:07 f_limousin
70  * Delete ener_dens, mom_dens and trace stress in functions
71  * source_nn, source_psi and source_beta.
72  * And some modification to avoid warnings (source_nn change to source...).
73  *
74  * Revision 1.5 2004/11/03 17:16:44 f_limousin
75  * Delete argument trk_point for source_nn()
76  *
77  * Revision 1.4 2004/10/29 15:45:08 jl_jaramillo
78  * Change name of functions
79  *
80  * Revision 1.3 2004/09/28 16:08:26 f_limousin
81  * Minor modifs.
82  *
83  * Revision 1.2 2004/09/09 16:55:08 f_limousin
84  * Add the 2 lines $Id: sources_hor.C,v 1.16 2014/10/13 08:53:01 j_novak Exp $Log: for CVS
85  *
86  *
87  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/sources_hor.C,v 1.16 2014/10/13 08:53:01 j_novak Exp $
88  *
89  */
90 
91 // C++ headers
92 #include "headcpp.h"
93 
94 // C headers
95 #include <cstdlib>
96 #include <cassert>
97 
98 // Lorene headers
99 #include "time_slice.h"
100 #include "isol_hor.h"
101 #include "metric.h"
102 #include "evolution.h"
103 #include "unites.h"
104 #include "graphique.h"
105 #include "utilitaires.h"
106 
107 namespace Lorene {
109 
110  using namespace Unites ;
111 
112  // Initialisations
113  // ---------------
114 
115  const Base_vect& triad = *(ff.get_triad()) ;
116 
117  Scalar tmp(mp) ;
118  Scalar tmp_scal(mp) ;
119  Sym_tensor tmp_sym(mp, CON, triad) ;
120 
121  Scalar source(mp) ;
122 
123  //===============================================
124  // Computations of the source for Psi
125  //===============================================
126 
127  const Vector& d_psi = psi().derive_cov(ff) ; // D_i Psi
128 
129  // Source for Psi
130  // --------------
131  tmp = 0.125* psi() * met_gamt.ricci_scal()
132  - contract(hh(), 0, 1, d_psi.derive_cov(ff), 0, 1 ) ;
133  tmp.inc_dzpuis() ; // dzpuis : 3 -> 4
134 
135  tmp -= contract(hdirac(), 0, d_psi, 0) ;
136 
137  source = tmp - 0.125* aa_quad() / psi4() / psi()/ psi()/ psi()
138  + psi()*psi4() * 8.33333333333333e-2* trK*trK ; //##
139  source.annule_domain(0) ;
140 
141  return source ;
142 }
143 
144 
146 
147  using namespace Unites ;
148 
149  // Initialisations
150  // ---------------
151 
152  const Base_vect& triad = *(ff.get_triad()) ;
153 
154  Scalar tmp(mp) ;
155  Scalar tmp_scal(mp) ;
156  Sym_tensor tmp_sym(mp, CON, triad) ;
157 
158  Scalar source(mp) ;
159 
160  //===============================================
161  // Computations of the source for NN
162  //===============================================
163 
164  const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
165  const Vector& dnnn = nn().derive_cov(ff) ; // D_i N
166 
167  // Source for N
168  // ------------
169 
170  source = aa_quad() / psi4() / psi4() * nn() +
171  psi4()*( nn()* 0.3333333333333333* trK*trK - trK_point )
172  - 2.* contract(dln_psi, 0, nn().derive_con(met_gamt), 0)
173  - contract(hdirac(), 0, dnnn, 0) ;
174 
175  tmp = psi4()* contract(beta(), 0, trK.derive_cov(ff), 0)
176  - contract( hh(), 0, 1, dnnn.derive_cov(ff), 0, 1 ) ;
177 
178  tmp.inc_dzpuis() ; // dzpuis: 3 -> 4
179 
180  source += tmp ;
181 
182  source.annule_domain(0) ;
183 
184  return source ;
185 
186 }
187 
188 
189 
191 
192  using namespace Unites ;
193 
194  // Initialisations
195  // ---------------
196 
197  const Base_vect& triad = *(ff.get_triad()) ;
198 
199  Scalar tmp(mp) ;
200  Scalar tmp_scal(mp) ;
201  Sym_tensor tmp_sym(mp, CON, triad) ;
202  Vector tmp_vect(mp, CON, triad) ;
203 
204  Vector source(mp, CON, triad) ;
205 
206  //===============================================
207  // Computations of the source for beta
208  //===============================================
209 
210  const Vector& dln_psi = ln_psi().derive_cov(ff) ; // D_i ln(Psi)
211  const Vector& dnnn = nn().derive_cov(ff) ; // D_i N
212 
213  // Source for beta (dzpuis = 4)
214  // ----------------------------
215 
216  source = 2.* contract(aa(), 1,
217  dnnn - 6.*nn() * dln_psi, 0) ;
218 
219  tmp_vect = 0.66666666666666666* trK.derive_con(met_gamt) ;
220  tmp_vect.inc_dzpuis() ;
221 
222  source += 2.* nn() * ( tmp_vect
223  - contract(met_gamt.connect().get_delta(), 1, 2,
224  aa(), 0, 1) ) ;
225 
226  Vector vtmp = contract(hh(), 0, 1,
227  beta().derive_cov(ff).derive_cov(ff), 1, 2)
228  + 0.3333333333333333*
229  contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0)
230  - hdirac().derive_lie(beta())
231  + gamt_point.divergence(ff) ; // zero in the Dirac gauge
232  vtmp.inc_dzpuis() ; // dzpuis: 3 -> 4
233 
234  source -= vtmp ;
235 
236  source += 0.66666666666666666* beta().divergence(ff) * hdirac() ;
237 
238  source.annule_domain(0) ;
239 
240  /*
241  // Source for beta (dzpuis = 3)
242  // ----------------------------
243 
244  source = 2.* contract(aa(), 1,
245  dnnn - 6.*nn() * dln_psi, 0) ;
246  source.dec_dzpuis() ;
247 
248  tmp_vect = 0.66666666666666666* trK.derive_con(met_gamt) ;
249  Vector tmp_vect2 (- contract(met_gamt.connect().get_delta(), 1, 2,
250  aa(), 0, 1)) ;
251  tmp_vect2.dec_dzpuis() ;
252 
253  source += 2.* nn() * ( tmp_vect + tmp_vect2 ) ;
254 
255  Vector vtmp = contract(hh(), 0, 1,
256  beta().derive_cov(ff).derive_cov(ff), 1, 2)
257  + 0.3333333333333333*
258  contract(hh(), 1, beta().divergence(ff).derive_cov(ff), 0)
259  - hdirac().derive_lie(beta())
260  + gamt_point.divergence(ff) ; // zero in the Dirac gauge
261 
262  source -= vtmp ;
263 
264  tmp_vect = 0.66666666666666666* beta().divergence(ff) * hdirac() ;
265  tmp_vect.dec_dzpuis() ;
266  source += tmp_vect ;
267 
268  source.annule_domain(0) ;
269  */
270 
271  return source ;
272 
273 }
274 
275 
276 
277 
278 
279 
280 
281 
282 }
Lorene::Metric::ricci_scal
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:350
Lorene::Time_slice_conf::hdirac
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Definition: time_slice_conf.C:815
Lorene::Isol_hor::source_beta
const Vector source_beta() const
Source for .
Definition: sources_hor.C:190
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::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Definition: time_slice_conf.C:591
Lorene::Time_slice_conf::ln_psi
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Definition: time_slice_conf.C:719
Lorene::Isol_hor::source_psi
const Scalar source_psi() const
Source for .
Definition: sources_hor.C:108
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Time_slice_conf::hh
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
Definition: time_slice_conf.C:758
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::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Lorene::Isol_hor::mp
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
Lorene::Isol_hor::trK
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:332
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::Isol_hor::gamt_point
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition: isol_hor.h:329
Unites
Standard units of space, time and mass.
Lorene::Isol_hor::met_gamt
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
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::Isol_hor::trK_point
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition: isol_hor.h:335
Lorene::Connection::get_delta
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
Lorene::Isol_hor::source_nn
const Scalar source_nn() const
Source for N.
Definition: sources_hor.C:145
Lorene::Time_slice_conf::ff
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:507
Lorene::Isol_hor::aa_quad
virtual const Scalar & aa_quad() const
Conformal representation .
Definition: isol_hor.C:884
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::Vector::derive_lie
Vector derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: vector.C:392
Lorene::Scalar::derive_con
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
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::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition: time_slice_access.C:87
Lorene::Metric_flat::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
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::Base_vect
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
Lorene::Metric::connect
virtual const Connection & connect() const
Returns the connection.
Definition: metric.C:301
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