29 char binary_helical_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Binary/binary_helical.C,v 1.6 2014/10/13 08:52:45 j_novak Exp $" ;
69 #include "graphique.h" 70 #include "utilitaires.h" 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++)
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 ;
209 omdsdp.set(1).set_outer_boundary(nz-1, 0) ;
210 omdsdp.set(2).set_outer_boundary(nz-1, 0) ;
211 omdsdp.std_spectral_base() ;
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) ;
338 source_3a.inc_dzpuis() ;
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) ;
537 lie_aij_2.dec_dzpuis(3) ;
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 ;
560 lie_K2_1.set_etat_qcq() ;
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)) ;
575 lie_aij2_1.set(i,j).set_spectral_va().set_base(lie_aij_2(i,j).
576 get_spectral_va().get_base()) ;
579 lie_aij_tot_1 = lie_aij_1 + lie_aij2_1 ;
580 lie_aij_tot_1.inc_dzpuis(2) ;
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 ;
634 lie_kij_tot.dec_dzpuis(2) ;
636 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 ;
Tbl & set_domain(int l)
Read/write of the value in a given domain.
void annule_hard()
Sets the Scalar to zero in a hard way.
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.
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
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