89 for (
int ll=1; ll<=2; ll++) {
118 const Scalar phi (0.5 * (star_i.lnq - star_i.
logn)) ;
155 double lambda_dirac = 0. ;
158 const Vector hdirac_auto = lambda_dirac *
172 double r0 = mp.
val_r(nz-2, 1, 0, 0) ;
173 double sigma = 1.*r0 ;
180 ff =
exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
181 for (
int ii=0; ii<nz-1; ii++)
182 ff.set_domain(ii) = 1. ;
183 ff.set_outer_boundary(nz-1, 0) ;
184 ff.std_spectral_base() ;
199 omdsdp.
set(1) = - om * yya * ff ;
200 omdsdp.
set(2) = om * xxa * ff ;
204 omdsdp.
set(1) = om * yya * ff ;
205 omdsdp.
set(2) = - om * xxa * ff ;
221 for (
int i=1; i<=3; i++)
222 for (
int j=1; j<=i; j++) {
223 tkij_a.
set(i, j) = dbeta(i, j) + dbeta(j, i) -
224 double(2) /double(3) * div_beta * (gtilde.
con())(i,j) ;
228 tkij_a = 0.5 * tkij_a / nn ;
232 Scalar a2_auto =
contract(tkij_auto_cov, 0, 1, tkij_a, 0, 1,
true) ;
239 for (
int i=1; i<=3; i++)
240 for (
int j=i; j<=3; j++) {
241 tkij_c.
set(i, j) = dbeta_comp(i, j) + dbeta_comp(j, i) -
242 double(2) /double(3) * divbeta_comp * (gtilde.
con())(i,j) ;
246 tkij_c = 0.5 * tkij_c / nn ;
248 Scalar a2_comp =
contract(tkij_auto_cov, 0, 1, tkij_c, 0, 1,
true) ;
266 source2 = star_i.
psi4 % (a2_auto + a2_comp) ;
270 0, dcov_logn_auto, 0,
true) ;
272 source4 = -
contract(star_i.
hij, 0, 1, dcovdcov_logn_auto +
275 source_tot = source1 + source2 + source3 + source4 ;
292 Scalar source_tot_hij(mp) ;
306 source_1 =
contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
308 source_2 = -
contract(dcon_hij, 2, dcov_lnq_auto, 0)
309 - 2./3. *
contract(hdirac, 0, dcov_lnq_auto, 0) * flat.
con() ;
315 (star_i.
logn - 2.e-8) ;
326 dcov_omdsdphi.
set(1,1) = 0. ;
327 dcov_omdsdphi.
set(2,1) = om * ff ;
328 dcov_omdsdphi.
set(3,1) = 0. ;
329 dcov_omdsdphi.
set(2,2) = 0. ;
330 dcov_omdsdphi.
set(3,2) = 0. ;
331 dcov_omdsdphi.
set(3,3) = 0. ;
332 dcov_omdsdphi.
set(1,2) = -om *ff ;
333 dcov_omdsdphi.
set(1,3) = 0. ;
334 dcov_omdsdphi.
set(2,3) = 0. ;
337 source_3a =
contract(tkij_a, 0, dcov_omdsdphi, 1) ;
353 source_5 = dcon_hdirac_auto ;
355 source_6 = - 2./3. * hdirac_auto.
divergence(flat) * flat.
con() ;
361 source_Sij = 8. * nn / psi4 * phi_auto.
derive_con(gtilde) *
364 source_Sij += 4. / psi4 * phi_auto.
derive_con(gtilde) *
369 source_Sij += - nn / (3.*psi4) * gtilde_con *
371 dcov_gtilde, 0, 1), 0, 1)
373 dcov_gtilde, 0, 2), 0, 1)) ;
375 source_Sij += - 8.*nn / (3.*psi4) * gtilde_con *
381 source_Sij += tens_temp ;
383 source_Sij += - 8./(3.*psi4) *
contract(dcov_phi_auto, 0,
386 source_Sij += 2.*nn*
contract(gtilde_cov, 0, 1, tkij_a *
387 (tkij_a+tkij_c), 1, 3) ;
389 source_Sij += - 2. * qpig * nn * ( psi4 * star_i.
stress_euler
390 - 0.33333333333333333 * star_i.
s_euler * gtilde_con ) ;
392 source_Sij += - 1./(psi4*psi2) *
contract(gtilde_con, 1,
393 contract(gtilde_con, 1, qq*dcovdcov_lnq_auto +
394 qq*dcov_lnq_auto*dcov_lnq, 0), 1) ;
397 dcov_hij_auto, 2), 1, dcov_qq, 0) -
399 star_i.
hij, 1), 1, dcov_qq, 0) ;
402 dcov_hij_auto, 2), 0, dcov_qq, 0) ;
404 source_Sij += 1./(3.*psi4*psi2)*
contract(gtilde_con, 0, 1,
405 qq*dcovdcov_lnq_auto + qq*dcov_lnq_auto*dcov_lnq, 0, 1)
408 source_Sij += 1./(3.*psi4*psi2) *
contract(hdirac, 0,
422 source_Rij +=
contract(hdirac, 0, dcov_hij_auto, 2) ;
428 gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
431 dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
433 dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
435 source_Rij += 0.5 *
contract(gtilde_con*gtilde_con, 1, 3,
436 contract(dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
438 source_Rij = source_Rij * 0.5 ;
440 for(
int i=1; i<=3; i++)
441 for(
int j=1; j<=i; j++) {
443 source_tot_hij = source_1(i,j) + source_1(j,i)
444 + source_2(i,j) + 2.*psi4/nn * (
445 source_4(i,j) - source_Sij(i,j))
446 - 2.* source_Rij(i,j) +
447 source_5(i,j) + source_5(j,i) + source_6(i,j) ;
450 source3 = 2.*psi4/nn * (source_3a(i,j) + source_3a(j,i)
453 source_tot_hij = source_tot_hij + source3 ;
460 lie_aij_1.
set(1,1) = nn / (2.*psi4) *
461 (lapl-source_tot_hij) ;
466 lie_aij_1.
set(2,1) = nn / (2.*psi4) *
467 (lapl-source_tot_hij) ;
472 lie_aij_1.
set(3,1) = nn / (2.*psi4) *
473 (lapl-source_tot_hij) ;
478 lie_aij_1.
set(2,2) = nn / (2.*psi4) *
479 (lapl-source_tot_hij) ;
484 lie_aij_1.
set(3,2) = nn / (2.*psi4) *
485 (lapl-source_tot_hij) ;
490 lie_aij_1.
set(3,3) = nn / (2.*psi4) *
491 (lapl-source_tot_hij) ;
499 lie_aij_2.
set(1,1) = nn / (2.*psi4) *
500 (lapl-source_tot_hij) ;
505 lie_aij_2.
set(2,1) = nn / (2.*psi4) *
506 (lapl-source_tot_hij) ;
511 lie_aij_2.
set(3,1) = nn / (2.*psi4) *
512 (lapl-source_tot_hij) ;
517 lie_aij_2.
set(2,2) = nn / (2.*psi4) *
518 (lapl-source_tot_hij) ;
523 lie_aij_2.
set(3,2) = nn / (2.*psi4) *
524 (lapl-source_tot_hij) ;
529 lie_aij_2.
set(3,3) = nn / (2.*psi4) *
530 (lapl-source_tot_hij) ;
542 double* bornes =
new double [6] ;
543 bornes[nz] = __infinity ;
544 bornes[4] = M_PI /
omega ;
545 bornes[3] = M_PI /
omega * 0.5 ;
546 bornes[2] = M_PI /
omega * 0.2 ;
547 bornes[1] = M_PI /
omega * 0.1 ;
567 lie_K2_1.
import(lie_K_2) ;
568 lie_K_tot_1 = lie_K_1 + lie_K2_1 ;
572 for(
int i=1; i<=3; i++)
573 for(
int j=1; j<=i; j++) {
574 lie_aij2_1.
set(i,j).
import(lie_aij_2(i,j)) ;
576 get_spectral_va().get_base()) ;
579 lie_aij_tot_1 = lie_aij_1 + lie_aij2_1 ;
588 cout <<
" IN THE CENTER OF STAR 1 " << endl
589 <<
" ----------------------- " << endl ;
613 cout <<
"L2 norm of L_k K^{ab} " << endl ;
622 for (
int mm=0; mm<nz; mm++)
623 for (
int pp=0; pp<=mm; pp++)
624 integral.
set(mm) += integ(pp) ;
630 cout <<
"omega = " <<
omega << endl ;
631 cout <<
"mass_adm = " <<
mass_adm() << endl ;
636 cout <<
"Position du centre de l'etoile x/lambda = "