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