72 #include "graphique.h" 73 #include "utilitaires.h" 92 for (
int ll=1; ll<=2; ll++) {
121 const Scalar phi (0.5 * (star_i.lnq - star_i.
logn)) ;
158 double lambda_dirac = 0. ;
161 const Vector hdirac_auto = lambda_dirac *
175 double r0 = mp.
val_r(nz-2, 1, 0, 0) ;
176 double sigma = 1.*r0 ;
183 ff =
exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
184 for (
int ii=0; ii<nz-1; ii++)
185 ff.set_domain(ii) = 1. ;
186 ff.set_outer_boundary(nz-1, 0) ;
187 ff.std_spectral_base() ;
202 omdsdp.
set(1) = - om * yya * ff ;
203 omdsdp.set(2) = om * xxa * ff ;
204 omdsdp.set(3).annule_hard() ;
207 omdsdp.set(1) = om * yya * ff ;
208 omdsdp.set(2) = - om * xxa * ff ;
209 omdsdp.set(3).annule_hard() ;
212 omdsdp.set(1).set_outer_boundary(nz-1, 0) ;
213 omdsdp.set(2).set_outer_boundary(nz-1, 0) ;
214 omdsdp.std_spectral_base() ;
224 for (
int i=1; i<=3; i++)
225 for (
int j=1; j<=i; j++) {
226 tkij_a.
set(i, j) = dbeta(i, j) + dbeta(j, i) -
227 double(2) /double(3) * div_beta * (gtilde.
con())(i,j) ;
231 tkij_a = 0.5 * tkij_a / nn ;
235 Scalar a2_auto =
contract(tkij_auto_cov, 0, 1, tkij_a, 0, 1,
true) ;
242 for (
int i=1; i<=3; i++)
243 for (
int j=i; j<=3; j++) {
244 tkij_c.
set(i, j) = dbeta_comp(i, j) + dbeta_comp(j, i) -
245 double(2) /double(3) * divbeta_comp * (gtilde.
con())(i,j) ;
249 tkij_c = 0.5 * tkij_c / nn ;
251 Scalar a2_comp =
contract(tkij_auto_cov, 0, 1, tkij_c, 0, 1,
true) ;
269 source2 = star_i.
psi4 % (a2_auto + a2_comp) ;
273 0, dcov_logn_auto, 0,
true) ;
275 source4 = -
contract(star_i.
hij, 0, 1, dcovdcov_logn_auto +
278 source_tot = source1 + source2 + source3 + source4 ;
295 Scalar source_tot_hij(mp) ;
309 source_1 =
contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
311 source_2 = -
contract(dcon_hij, 2, dcov_lnq_auto, 0)
312 - 2./3. *
contract(hdirac, 0, dcov_lnq_auto, 0) * flat.
con() ;
318 (star_i.
logn - 2.e-8) ;
329 dcov_omdsdphi.
set(1,1) = 0. ;
330 dcov_omdsdphi.
set(2,1) = om * ff ;
331 dcov_omdsdphi.set(3,1) = 0. ;
332 dcov_omdsdphi.set(2,2) = 0. ;
333 dcov_omdsdphi.set(3,2) = 0. ;
334 dcov_omdsdphi.set(3,3) = 0. ;
335 dcov_omdsdphi.set(1,2) = -om *ff ;
336 dcov_omdsdphi.set(1,3) = 0. ;
337 dcov_omdsdphi.set(2,3) = 0. ;
338 dcov_omdsdphi.std_spectral_base() ;
340 source_3a =
contract(tkij_a, 0, dcov_omdsdphi, 1) ;
341 source_3a.inc_dzpuis() ;
356 source_5 = dcon_hdirac_auto ;
358 source_6 = - 2./3. * hdirac_auto.
divergence(flat) * flat.
con() ;
364 source_Sij = 8. * nn / psi4 * phi_auto.
derive_con(gtilde) *
367 source_Sij += 4. / psi4 * phi_auto.
derive_con(gtilde) *
372 source_Sij += - nn / (3.*psi4) * gtilde_con *
374 dcov_gtilde, 0, 1), 0, 1)
376 dcov_gtilde, 0, 2), 0, 1)) ;
378 source_Sij += - 8.*nn / (3.*psi4) * gtilde_con *
384 source_Sij += tens_temp ;
386 source_Sij += - 8./(3.*psi4) *
contract(dcov_phi_auto, 0,
389 source_Sij += 2.*nn*
contract(gtilde_cov, 0, 1, tkij_a *
390 (tkij_a+tkij_c), 1, 3) ;
392 source_Sij += - 2. * qpig * nn * ( psi4 * star_i.
stress_euler 393 - 0.33333333333333333 * star_i.
s_euler * gtilde_con ) ;
395 source_Sij += - 1./(psi4*psi2) *
contract(gtilde_con, 1,
396 contract(gtilde_con, 1, qq*dcovdcov_lnq_auto +
397 qq*dcov_lnq_auto*dcov_lnq, 0), 1) ;
400 dcov_hij_auto, 2), 1, dcov_qq, 0) -
402 star_i.
hij, 1), 1, dcov_qq, 0) ;
405 dcov_hij_auto, 2), 0, dcov_qq, 0) ;
407 source_Sij += 1./(3.*psi4*psi2)*
contract(gtilde_con, 0, 1,
408 qq*dcovdcov_lnq_auto + qq*dcov_lnq_auto*dcov_lnq, 0, 1)
411 source_Sij += 1./(3.*psi4*psi2) *
contract(hdirac, 0,
425 source_Rij +=
contract(hdirac, 0, dcov_hij_auto, 2) ;
431 gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
434 dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
436 dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
438 source_Rij += 0.5 *
contract(gtilde_con*gtilde_con, 1, 3,
439 contract(dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
441 source_Rij = source_Rij * 0.5 ;
443 for(
int i=1; i<=3; i++)
444 for(
int j=1; j<=i; j++) {
446 source_tot_hij = source_1(i,j) + source_1(j,i)
447 + source_2(i,j) + 2.*psi4/nn * (
448 source_4(i,j) - source_Sij(i,j))
449 - 2.* source_Rij(i,j) +
450 source_5(i,j) + source_5(j,i) + source_6(i,j) ;
453 source3 = 2.*psi4/nn * (source_3a(i,j) + source_3a(j,i)
456 source_tot_hij = source_tot_hij + source3 ;
463 lie_aij_1.set(1,1) = nn / (2.*psi4) *
464 (lapl-source_tot_hij) ;
469 lie_aij_1.set(2,1) = nn / (2.*psi4) *
470 (lapl-source_tot_hij) ;
475 lie_aij_1.set(3,1) = nn / (2.*psi4) *
476 (lapl-source_tot_hij) ;
481 lie_aij_1.set(2,2) = nn / (2.*psi4) *
482 (lapl-source_tot_hij) ;
487 lie_aij_1.set(3,2) = nn / (2.*psi4) *
488 (lapl-source_tot_hij) ;
493 lie_aij_1.set(3,3) = nn / (2.*psi4) *
494 (lapl-source_tot_hij) ;
502 lie_aij_2.set(1,1) = nn / (2.*psi4) *
503 (lapl-source_tot_hij) ;
508 lie_aij_2.set(2,1) = nn / (2.*psi4) *
509 (lapl-source_tot_hij) ;
514 lie_aij_2.set(3,1) = nn / (2.*psi4) *
515 (lapl-source_tot_hij) ;
520 lie_aij_2.set(2,2) = nn / (2.*psi4) *
521 (lapl-source_tot_hij) ;
526 lie_aij_2.set(3,2) = nn / (2.*psi4) *
527 (lapl-source_tot_hij) ;
532 lie_aij_2.set(3,3) = nn / (2.*psi4) *
533 (lapl-source_tot_hij) ;
540 lie_aij_2.dec_dzpuis(3) ;
545 double* bornes =
new double [6] ;
546 bornes[nz] = __infinity ;
547 bornes[4] = M_PI /
omega ;
548 bornes[3] = M_PI /
omega * 0.5 ;
549 bornes[2] = M_PI /
omega * 0.2 ;
550 bornes[1] = M_PI /
omega * 0.1 ;
563 lie_K2_1.set_etat_qcq() ;
570 lie_K2_1.import(lie_K_2) ;
571 lie_K_tot_1 = lie_K_1 + lie_K2_1 ;
575 for(
int i=1; i<=3; i++)
576 for(
int j=1; j<=i; j++) {
577 lie_aij2_1.set(i,j).import(lie_aij_2(i,j)) ;
578 lie_aij2_1.set(i,j).set_spectral_va().set_base(lie_aij_2(i,j).
579 get_spectral_va().get_base()) ;
582 lie_aij_tot_1 = lie_aij_1 + lie_aij2_1 ;
583 lie_aij_tot_1.inc_dzpuis(2) ;
591 cout <<
" IN THE CENTER OF STAR 1 " << endl
592 <<
" ----------------------- " << endl ;
616 cout <<
"L2 norm of L_k K^{ab} " << endl ;
625 for (
int mm=0; mm<nz; mm++)
626 for (
int pp=0; pp<=mm; pp++)
627 integral.
set(mm) += integ(pp) ;
633 cout <<
"omega = " <<
omega << endl ;
634 cout <<
"mass_adm = " <<
mass_adm() << endl ;
637 lie_kij_tot.dec_dzpuis(2) ;
639 cout <<
"Position du centre de l'etoile x/lambda = " Sym_tensor hij_comp
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Coord xa
Absolute x coordinate.
Metric for tensor calculation.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Cmp exp(const Cmp &)
Exponential.
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Map & mp
Mapping associated with the star.
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Cmp sqrt(const Cmp &)
Square root.
Star_bin star1
First star of the system.
Standard units of space, time and mass.
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Flat metric for tensor calculation.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Tensor field of valence 0 (or component of a tensorial field).
Base class for coordinate mappings.
double get_ori_x() const
Returns the x coordinate of the origin.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
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...
Basic integer array class.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
const Metric & get_gamma() const
Returns the 3-metric .
const Map & get_mp() const
Returns the mapping.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Tensor field of valence 1.
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
void inc_dzpuis()
dzpuis += 1 ;
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Tbl & integrale_domains() const
Computes the integral in each domain of *this .
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
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...
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) ...
Scalar lnq_auto
Scalar field generated principally by the star.
int get_nzone() const
Returns the number of domains.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Sym_tensor tkij_comp
Part of the extrinsic curvature tensor generated by beta_comp.
Cmp pow(const Cmp &, int)
Power .
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Scalar logn
Logarithm of the lapse N .
Class for stars in binary system.
virtual const Scalar & determinant() const
Returns the determinant.
Metric gtilde
Conformal 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...
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Star_bin * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Coord ya
Absolute y coordinate.
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Scalar nn
Lapse function N .
Vector beta_comp
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Coord za
Absolute z coordinate.
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
void helical()
Function testing the helical symmetry.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
double mass_adm() const
Total ADM mass.
void annule_hard()
Sets the Tbl to zero in a hard way.
Star_bin star2
Second star of the system.
Scalar ener_euler
Total energy density in the Eulerian frame.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Scalar psi4
Conformal factor .
Class intended to describe valence-2 symmetric tensors.
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
Coord r
r coordinate centered on the grid