LORENE
Isol_hole/secmembre_hij_stat.C
1
2// Header Lorene:
3#include "nbr_spx.h"
4#include "utilitaires.h"
5#include "graphique.h"
6#include "math.h"
7#include "metric.h"
8#include "param.h"
9#include "param_elliptic.h"
10#include "tensor.h"
11#include "unites.h"
12#include "isol_hole.h"
13
14//Computes the rhs of hyperbolic equation for conformal metric assuming statioarity; WARNING; up to now, we are only able to handle void spacetimes.
15
16
17namespace Lorene {
19
20 // Getting: hij; hatA; lapse; conf_fact; shift;
21 // hij; aa; nn; ppsi; bb;
22
23
24 using namespace Unites;
25
26 const int nz = (*(hij.get_mp().get_mg())).get_nzone();
27
28
29 const Vector& beta = shift;
30 const Sym_tensor& hh = hij;
31
32
33 const Scalar& psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
34 Scalar ln_psi = log(conf_fact); ln_psi.std_spectral_base();
35
36 const Scalar qq = lapse*conf_fact*conf_fact;
37
38
39 Sym_tensor aa = hatA/(psi4*sqrt(psi4));
40 aa.std_spectral_base(); //(check...)
41
42
43 const Metric_flat& ff = (hij.get_mp()).flat_met_spher() ;
44
45 const Sym_tensor& tgam_uu = ff.con() + hh;
46
47 const Metric tgam(tgam_uu);
48
49 const Base_vect_spher& otriad = hij.get_mp().get_bvect_spher();
50
51 // On met a zero les quantités supposees etre de "matiere"
52
53 Sym_tensor strain_tens = hij;
54 for (int ii=1; ii<=3; ii++)
55 for(int jj=1; jj<=3; jj++)
56 { strain_tens.set(ii,jj).annule_hard();
57 }
58 strain_tens.std_spectral_base();
59
60 // On met a zero les quantités derivee temporelle
61
62 Vector beta_point = shift;
63
64 for (int ii=1; ii<=3; ii++)
65
66 { beta_point.set(ii).annule_hard();
67 }
68 beta_point.annule_domain(nz-1) ; // Pour faire passer un assert, je ne comprends pas pourquoi...
69
70
71
72 beta_point.std_spectral_base();
73 Scalar lapse_point(hij.get_mp());
74 lapse_point.annule_hard();
75 lapse_point.annule_domain(nz -1);
76
77 lapse_point.std_spectral_base();
78
79 Sym_tensor hh_point = hij;
80 for (int ii=1; ii<=3; ii++)
81 for(int jj=1; jj<=3; jj++)
82 { hh_point.set(ii,jj).annule_hard();
83 }
84 hh_point.annule_domain(nz-1);
85 hh_point.std_spectral_base();
86
87
88 // Note: Il sera probablement nécessaire de ne pas mettre a zero hh point;
89
90
91 //Sym_tensor Rrij(map, CON, map.get_bvect_spher());
92
93
94 // Rrij = 0.5*[ contract (( h_iju, 0, h_iju.derive_cov(mets).derive_cov(mets), 3), 0,3) - contract((contract(h_iju.derive_cov(mets),1, h_iju.derive_cov(mets),2)), 1,3)] ;
95
96 //==================================
97 // Source for hij
98 //==================================
99
100 const Sym_tensor& tgam_dd = tgam.cov() ; // {\tilde \gamma}_{ij}
101 // const Sym_tensor& tgam_uu = tgam().con() ; // {\tilde \gamma}^{ij}
102 const Tensor_sym& dtgam = tgam_dd.derive_cov(ff) ;// D_k {\tilde \gamma}_{ij} // ff etant la métrique plate
103 const Tensor_sym& dhh = hh.derive_cov(ff) ; // D_k h^{ij}
104 const Vector& dln_psi = ln_psi.derive_cov(ff) ; // D_i ln(Psi)
105 const Vector& tdln_psi_u = ln_psi.derive_con(tgam) ; // tD^i ln(Psi)
106 const Vector& tdnn_u = lapse.derive_con(tgam) ; // tD^i N
107 const Vector& dqq = qq.derive_cov(ff) ; // D_i Q
108 const Scalar& div_beta = beta.divergence(ff) ; // D_k beta^k
109
110 Scalar tmp(hij.get_mp()) ;
111 Sym_tensor sym_tmp(hij.get_mp(), CON, otriad) ;
112
113 // Quadratic part of the Ricci tensor of gam_tilde
114 // ------------------------------------------------
115
116 Sym_tensor ricci_star(hij.get_mp(), CON, otriad) ;
117
118 ricci_star = contract(hh, 0, 1, dhh.derive_cov(ff), 2, 3) ;
119
120 ricci_star.inc_dzpuis() ; // dzpuis : 3 --> 4
121
122 // des_profile (ricci_star(1,1), 1, 8, 1,1, "riccistar"); // A enlever
123
124 for (int i=1; i<=3; i++) {
125 for (int j=1; j<=i; j++) {
126 tmp = 0 ;
127 for (int k=1; k<=3; k++) {
128 for (int l=1; l<=3; l++) {
129 tmp += dhh(i,k,l) * dhh(j,l,k) ;
130 }
131 }
132 sym_tmp.set(i,j) = tmp ;
133 }
134 }
135 ricci_star -= sym_tmp ;
136 for (int i=1; i<=3; i++) {
137 for (int j=1; j<=i; j++) {
138 tmp = 0 ;
139 for (int k=1; k<=3; k++) {
140 for (int l=1; l<=3; l++) {
141 for (int m=1; m<=3; m++) {
142 for (int n=1; n<=3; n++) {
143
144 tmp += 0.5 * tgam_uu(i,k)* tgam_uu(j,l)
145 * dhh(m,n,k) * dtgam(m,n,l)
146 + tgam_dd(n,l) * dhh(m,n,k)
147 * (tgam_uu(i,k) * dhh(j,l,m) + tgam_uu(j,k) * dhh(i,l,m) )
148 - tgam_dd(k,l) *tgam_uu(m,n) * dhh(i,k,m) * dhh(j,l,n) ;
149 }
150 }
151 }
152 }
153 sym_tmp.set(i,j) = tmp ;
154 }
155 }
156 ricci_star += sym_tmp ;
157
158 ricci_star = 0.5 * ricci_star ;
159
160 // Curvature scalar of conformal metric :
161 // -------------------------------------
162
163 Scalar tricci_scal =
164 0.25 * contract(tgam_uu, 0, 1,
165 contract(dhh, 0, 1, dtgam, 0, 1), 0, 1 )
166 - 0.5 * contract(tgam_uu, 0, 1,
167 contract(dhh, 0, 1, dtgam, 0, 2), 0, 1 ) ;
168
169 // Full quadratic part of source for h : S^{ij}
170 // --------------------------------------------
171
172 Sym_tensor ss(hij.get_mp(), CON, otriad) ; // Source secondaire
173
174 sym_tmp = lapse * (ricci_star + 8.* tdln_psi_u * tdln_psi_u)
175 + 4.*( tdln_psi_u * tdnn_u + tdnn_u * tdln_psi_u )
176 - 0.3333333333333333 *
177 ( lapse * (tricci_scal + 8.* contract(dln_psi, 0, tdln_psi_u, 0) )
178 + 8.* contract(dln_psi, 0, tdnn_u, 0) ) *tgam_uu ;
179
180 ss = sym_tmp / psi4 ;
181
182 sym_tmp = contract(tgam_uu, 1,
183 contract(tgam_uu, 1, dqq.derive_cov(ff), 0), 1) ;
184
185 sym_tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
186 // des_profile (sym_tmp(1,1), 1, 8, 1,1, "sym_tmp"); // A enlever
187
188
189 for (int i=1; i<=3; i++) {
190 for (int j=1; j<=i; j++) {
191 tmp = 0 ;
192 for (int k=1; k<=3; k++) {
193 for (int l=1; l<=3; l++) {
194 tmp += ( hh(i,k)*dhh(l,j,k) + hh(k,j)*dhh(i,l,k)
195 - hh(k,l)*dhh(i,j,k) ) * dqq(l) ;
196 }
197 }
198 sym_tmp.set(i,j) += 0.5 * tmp ;
199 }
200 }
201
202 tmp = qq.derive_con(tgam).divergence(tgam) ;
203 tmp.inc_dzpuis() ; // dzpuis : 3 --> 4 // reverifier pourquoi
204
205 // des_profile (tmp, 1, 8, 1,1, "tmp"); // A enlever
206
207 sym_tmp -= 0.3333333333333333 * tmp *tgam_uu ;
208
209 ss -= sym_tmp / (psi4*conf_fact*conf_fact) ; // Voir dans quel sens sont construits psi et psi4 (eviter les multiplications d'erreurs)
210
211 for (int i=1; i<=3; i++) {
212 for (int j=1; j<=i; j++) {
213 tmp = 0 ;
214 for (int k=1; k<=3; k++) {
215 for (int l=1; l<=3; l++) {
216 tmp += tgam_dd(k,l) * aa(i,k) * aa(j,l) ;
217 }
218 }
219 sym_tmp.set(i,j) = tmp ;
220 }
221 }
222
223 tmp = psi4 * strain_tens.trace(tgam) ; // S = S_i^i
224
225 ss += (2.*lapse) * ( sym_tmp);// - qpig*( psi4* strain_tens
226 // - 0.3333333333333333 * tmp * tgam_uu );
227 Sym_tensor ss2 =2.*lapse*( qpig*(psi4*strain_tens - 0.33333333333333 * tmp * tgam_uu));
228 ss2.inc_dzpuis(4); // A retirer!
229
230 // des_profile (ss2(1,1), 1, 8, 1,1, "ss2"); // A enlever
231
232 // cout << zone << endl;
233
234 // ss2.annule_domain(nz-1);
235
236 ss += -ss2; // ATTENTION!!!! A RETABLIR!!!!
237
238 // maxabs(ss, "ss tot") ;
239
240 // Source for h^{ij}
241 // -----------------
242
243 Sym_tensor lbh = hh.derive_lie(beta) ;
244
245 source_hh =// (lapse*lapse/psi4 - 1.)
246 // * hh.derive_con(ff).divergence(ff)
247 + 2.* hh_point.derive_lie(beta); // - lbh.derive_lie(beta) ; // La double derivée de
248 // Lie en Beta est retirée (prise en charge dans tensorelliptic.C)
249
250 source_hh.inc_dzpuis() ;
251
252 // des_profile (source_hh(1,1), 1, 8, 1,1, "sourcehh"); // A enlever
253
254 source_hh += 2.* lapse * ss ;
255
256 //## Provisory: waiting for the Lie derivative to allow
257 // derivation with respect to a vector with dzpuis != 0
258 Vector vtmp = beta_point ;
259 // vtmp.dec_dzpuis(2) ; // A remettre si jamais beta point est non nul;
260
261 sym_tmp = hh.derive_lie(vtmp) ;
262 sym_tmp.inc_dzpuis(2) ;
263
264 // des_profile (sym_tmp(1,1), 1, 8, 1,1, "sym_tmp"); // A enlever
265
266 source_hh += sym_tmp
267 + 1.3333333333333333 * div_beta* (hh_point - lbh)
268 + 2. * (lapse_point - lapse.derive_lie(beta)) * aa ;
269
270
271 for (int i=1; i<=3; i++) {
272 for (int j=1; j<=i; j++) {
273 tmp = 0. ;
274 for (int k=1; k<=3; k++) {
275 tmp += ( hh.derive_con(ff)(k,j,i)
276 + hh.derive_con(ff)(i,k,j)
277 - hh.derive_con(ff)(i,j,k) ) * dqq(k) ;
278 }
279 sym_tmp.set(i,j) = tmp ;
280 }
281 }
282
283 source_hh -= lapse / (psi4*conf_fact*conf_fact) * sym_tmp ;
284
285 tmp = beta_point.divergence(ff) - div_beta.derive_lie(beta) ;
286 tmp.inc_dzpuis() ;
287
288 // des_profile (tmp, 1, 8, 1,1, "tmp"); // A enlever
289
290 source_hh += 0.6666666666666666*
291 ( tmp - 0.6666666666666666* div_beta * div_beta ) * hh ;
292
293
294 // Term (d/dt - Lie_beta) (L beta)^{ij}--> sym_tmp
295 // ---------------------------------------
296 Sym_tensor l_beta = beta.ope_killing_conf(ff) ; // Attention aux headers a inclure
297
298 sym_tmp = beta_point.ope_killing_conf(ff) - l_beta.derive_lie(beta) ;
299
300 sym_tmp.inc_dzpuis() ;
301
302 // Final source:
303 // ------------
304 source_hh += 0.6666666666666666* div_beta * l_beta - sym_tmp ;
305
306 // Invert it (because the right source is for a Laplace operator, not a Dalembertain operator as it is initially:
307
308 source_hh = -source_hh;
309
310
311 // Annulation de la source
312// for (int ii=1; ii<=3; ii++)
313// for (int jj=1; jj<=3; jj++){
314// source_hh.set(ii,jj).annule_hard();
315// }
316// source_hh.std_spectral_base();
317
318 return;
319
320}
321
322
323
324}
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Scalar lapse
Lapse function.
Definition isol_hole.h:76
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
Definition isol_hole.h:87
void secmembre_kerr(Sym_tensor &source_hh)
Computes the rhs of hyperbolic equation for conformal metric assuming statioarity; WARNING; up to now...
Vector shift
Shift vector.
Definition isol_hole.h:82
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: see Cordero et al.
Definition isol_hole.h:92
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Flat metric for tensor calculation.
Definition metric.h:261
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
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
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
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...
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
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_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.