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