LORENE
binary_helical.C
1 /*
2  * Methods of Bin_star::helical
3  *
4  * (see file star.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2006 Francois Limousin
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: binary_helical.C,v 1.7 2016/12/05 16:17:47 j_novak Exp $
33  * $Log: binary_helical.C,v $
34  * Revision 1.7 2016/12/05 16:17:47 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.6 2014/10/13 08:52:45 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.5 2014/10/06 15:12:59 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.4 2008/08/19 06:41:59 j_novak
44  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
45  * cast-type operations, and constant strings that must be defined as const char*
46  *
47  * Revision 1.3 2006/08/01 14:26:50 f_limousin
48  * Small changes
49  *
50  * Revision 1.2 2006/06/05 17:05:57 f_limousin
51  * *** empty log message ***
52  *
53  * Revision 1.1 2006/04/11 14:25:15 f_limousin
54  * New version of the code : improvement of the computation of some
55  * critical sources, estimation of the dirac gauge, helical symmetry...
56  *
57 
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Binary/binary_helical.C,v 1.7 2016/12/05 16:17:47 j_novak Exp $ *
60  */
61 
62 
63 // C headers
64 #include <cmath>
65 
66 // Lorene headers
67 #include "cmp.h"
68 #include "tenseur.h"
69 #include "metrique.h"
70 #include "binary.h"
71 #include "param.h"
72 #include "graphique.h"
73 #include "utilitaires.h"
74 #include "tensor.h"
75 #include "nbr_spx.h"
76 #include "unites.h"
77 
78 namespace Lorene {
80 
81  // Fundamental constants and units
82  // -------------------------------
83  using namespace Unites ;
84 
85 
86  Sym_tensor lie_aij_1 (star1.mp, CON, star1.mp.get_bvect_cart()) ;
87  Sym_tensor lie_aij_2 (star2.mp, CON, star2.mp.get_bvect_cart()) ;
88 
89  Scalar lie_K_1 (star1.mp) ;
90  Scalar lie_K_2 (star2.mp) ;
91 
92  for (int ll=1; ll<=2; ll++) {
93 
94  Star_bin star_i (*et[ll-1]) ;
95 
96  Map& mp = star_i.mp ;
97  const Mg3d* mg = mp.get_mg() ;
98  int nz = mg->get_nzone() ; // total number of domains
99 
100  Metric_flat flat = star_i.flat ;
101  Metric gtilde = star_i.gtilde ;
102  Scalar nn = star_i.nn ;
103  Scalar psi4 = star_i.psi4 ;
104 
105  // -------------------------------
106  // AUXILIARY QUANTITIES
107  // -------------------------------
108 
109  // Derivatives of N and logN
110  //--------------------------
111 
112  const Vector dcov_logn_auto = star_i.logn_auto.derive_cov(flat) ;
113 
114  Tensor dcovdcov_logn_auto = (star_i.logn_auto.derive_cov(flat))
115  .derive_cov(flat) ;
116  dcovdcov_logn_auto.inc_dzpuis() ;
117 
118  // Derivatives of lnq, phi and Q
119  //-------------------------------
120 
121  const Scalar phi (0.5 * (star_i.lnq - star_i.logn)) ;
122  const Scalar phi_auto (0.5 * (star_i.lnq_auto - star_i.logn_auto)) ;
123 
124  const Vector dcov_phi_auto = phi_auto.derive_cov(flat) ;
125 
126  const Vector dcov_lnq = 2*star_i.dcov_phi + star_i.dcov_logn ;
127  const Vector dcon_lnq = 2*star_i.dcon_phi + star_i.dcon_logn ;
128  const Vector dcov_lnq_auto = star_i.lnq_auto.derive_cov(flat) ;
129  Tensor dcovdcov_lnq_auto = dcov_lnq_auto.derive_cov(flat) ;
130  dcovdcov_lnq_auto.inc_dzpuis() ;
131 
132  Scalar qq = exp(star_i.lnq) ;
133  qq.std_spectral_base() ;
134  const Vector& dcov_qq = qq.derive_cov(flat) ;
135 
136  Tensor dcovdcov_beta_auto = star_i.beta_auto.derive_cov(flat)
137  .derive_cov(flat) ;
138  dcovdcov_beta_auto.inc_dzpuis() ;
139 
140 
141  // Derivatives of hij, gtilde...
142  //------------------------------
143 
144  Scalar psi2 (pow(star_i.psi4, 0.5)) ;
145  psi2.std_spectral_base() ;
146 
147  const Tensor& dcov_hij = star_i.hij.derive_cov(flat) ;
148  const Tensor& dcon_hij = star_i.hij.derive_con(flat) ;
149  const Tensor& dcov_hij_auto = star_i.hij_auto.derive_cov(flat) ;
150 
151  const Sym_tensor& gtilde_cov = star_i.gtilde.cov() ;
152  const Sym_tensor& gtilde_con = star_i.gtilde.con() ;
153  const Tensor& dcov_gtilde = gtilde_cov.derive_cov(flat) ;
154 
155  // H^i and its derivatives ( = O in Dirac gauge)
156  // ---------------------------------------------
157 
158  double lambda_dirac = 0. ;
159 
160  const Vector hdirac = lambda_dirac * star_i.hij.divergence(flat) ;
161  const Vector hdirac_auto = lambda_dirac *
162  star_i.hij_auto.divergence(flat) ;
163 
164  Tensor dcov_hdirac = lambda_dirac * hdirac.derive_cov(flat) ;
165  dcov_hdirac.inc_dzpuis() ;
166  Tensor dcov_hdirac_auto = lambda_dirac * hdirac_auto.derive_cov(flat) ;
167  dcov_hdirac_auto.inc_dzpuis() ;
168  Tensor dcon_hdirac_auto = lambda_dirac * hdirac_auto.derive_con(flat) ;
169  dcon_hdirac_auto.inc_dzpuis() ;
170 
171 
172  // Function exp(-(r-r_0)^2/sigma^2)
173  // --------------------------------
174 
175  double r0 = mp.val_r(nz-2, 1, 0, 0) ;
176  double sigma = 1.*r0 ;
177  double om = omega ;
178 
179  Scalar rr (mp) ;
180  rr = mp.r ;
181 
182  Scalar ff (mp) ;
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() ;
188 
189 
190  // omdsdp
191  // ---------------------------------
192 
193  Vector omdsdp (mp, CON, mp.get_bvect_cart()) ;
194  Scalar yya (mp) ;
195  yya = mp.ya ;
196  Scalar xxa (mp) ;
197  xxa = mp.xa ;
198  Scalar zza (mp) ;
199  zza = mp.za ;
200 
201  if (fabs(mp.get_rot_phi()) < 1e-10){
202  omdsdp.set(1) = - om * yya * ff ;
203  omdsdp.set(2) = om * xxa * ff ;
204  omdsdp.set(3).annule_hard() ;
205  }
206  else{
207  omdsdp.set(1) = om * yya * ff ;
208  omdsdp.set(2) = - om * xxa * ff ;
209  omdsdp.set(3).annule_hard() ;
210  }
211 
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() ;
215 
216 
217  // Computation of helical A^{ij}
218  // ------------------------------
219 
220  const Tensor& dbeta = star_i.beta_auto.derive_con(gtilde) ;
221  Scalar div_beta = star_i.beta_auto.divergence(gtilde) ;
222 
223  Sym_tensor tkij_a (star_i.tkij_auto) ;
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) ;
228  }
229 
230  tkij_a = tkij_a - star_i.hij_auto.derive_lie(omdsdp) ;
231  tkij_a = 0.5 * tkij_a / nn ;
232 
233  Sym_tensor tkij_auto_cov = tkij_a.up_down(gtilde) ;
234 
235  Scalar a2_auto = contract(tkij_auto_cov, 0, 1, tkij_a, 0, 1,true) ;
236 
237  // COMP
238  const Tensor& dbeta_comp = star_i.beta_comp.derive_con(gtilde) ;
239  Scalar divbeta_comp = star_i.beta_comp.divergence(gtilde) ;
240 
241  Sym_tensor tkij_c (star_i.tkij_comp) ;
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) ;
246  }
247 
248  tkij_c = tkij_c - star_i.hij_comp.derive_lie(omdsdp) ;
249  tkij_c = 0.5 * tkij_c / nn ;
250 
251  Scalar a2_comp = contract(tkij_auto_cov, 0, 1, tkij_c, 0, 1,true) ;
252 
253 
254 // tkij_a = star_i.tkij_auto ;
255 // tkij_c = star_i.tkij_comp ;
256 
257 
258  // Sources for N
259  // ---------------
260 
261  Scalar source1(mp) ;
262  Scalar source2(mp) ;
263  Scalar source3(mp) ;
264  Scalar source4(mp) ;
265  Scalar source_tot(mp) ;
266 
267  source1 = qpig * star_i.psi4 % (star_i.ener_euler + star_i.s_euler) ;
268 
269  source2 = star_i.psi4 % (a2_auto + a2_comp) ;
270 
271  source3 = - contract(dcov_logn_auto, 0, star_i.dcon_logn, 0, true)
272  - 2. * contract(contract(gtilde_con, 0, star_i.dcov_phi, 0),
273  0, dcov_logn_auto, 0, true) ;
274 
275  source4 = - contract(star_i.hij, 0, 1, dcovdcov_logn_auto +
276  dcov_logn_auto*star_i.dcov_logn, 0, 1) ;
277 
278  source_tot = source1 + source2 + source3 + source4 ;
279 
280  if (ll == 1) {
281  lie_K_1 = nn / star_i.psi4 * (source_tot - star_i.logn_auto
282  .laplacian()) ;
283  lie_K_1.dec_dzpuis(4) ;
284  }
285  if (ll == 2) {
286  lie_K_2 = nn / star_i.psi4 * (source_tot - star_i.logn_auto
287  .laplacian()) ;
288  lie_K_2.dec_dzpuis(4) ;
289  }
290 
291 
292  // Sources for hij
293  // --------------
294 
295  Scalar source_tot_hij(mp) ;
296  Tensor source_Sij(mp, 2, CON, mp.get_bvect_cart()) ;
297  Tensor source_Rij(mp, 2, CON, mp.get_bvect_cart()) ;
298  Tensor tens_temp(mp, 2, CON, mp.get_bvect_cart()) ;
299 
300  Tensor source_1 (mp, 2, CON, mp.get_bvect_cart()) ;
301  Tensor source_2 (mp, 2, CON, mp.get_bvect_cart()) ;
302  Tensor source_3a (mp, 2, CON, mp.get_bvect_cart()) ;
303  Tensor source_3b (mp, 2, CON, mp.get_bvect_cart()) ;
304  Tensor source_4 (mp, 2, CON, mp.get_bvect_cart()) ;
305  Tensor source_5 (mp, 2, CON, mp.get_bvect_cart()) ;
306  Tensor source_6 (mp, 2, CON, mp.get_bvect_cart()) ;
307 
308 
309  source_1 = contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
310 
311  source_2 = - contract(dcon_hij, 2, dcov_lnq_auto, 0)
312  - 2./3. * contract(hdirac, 0, dcov_lnq_auto, 0) * flat.con() ;
313 
314  // Lie derivative of A^{ij}
315  // --------------------------
316 
317  Scalar decouple_logn = (star_i.logn_auto - 1.e-8)/
318  (star_i.logn - 2.e-8) ;
319 
320  // Construction of Omega d/dphi
321  // ----------------------------
322 
323  // Construction of D_k \Phi^i
324  Itbl type (2) ;
325  type.set(0) = CON ;
326  type.set(1) = COV ;
327 
328  Tensor dcov_omdsdphi (mp, 2, type, mp.get_bvect_cart()) ;
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() ;
339 
340  source_3a = contract(tkij_a, 0, dcov_omdsdphi, 1) ;
341  source_3a.inc_dzpuis() ;
342 
343  // Source 3b
344  // ------------
345 
346  source_3b = - contract(omdsdp, 0, tkij_a.derive_cov(flat), 2) ;
347 
348 
349  // Source 4
350  // ---------
351 
352  source_4 = - tkij_a.derive_lie(star_i.beta) ;
353  source_4.inc_dzpuis() ;
354  source_4 += - 2./3. * star_i.beta.divergence(flat) * tkij_a ;
355 
356  source_5 = dcon_hdirac_auto ;
357 
358  source_6 = - 2./3. * hdirac_auto.divergence(flat) * flat.con() ;
359  source_6.inc_dzpuis() ;
360 
361  // Source terms for Sij
362  //---------------------
363 
364  source_Sij = 8. * nn / psi4 * phi_auto.derive_con(gtilde) *
365  contract(gtilde_con, 0, star_i.dcov_phi, 0) ;
366 
367  source_Sij += 4. / psi4 * phi_auto.derive_con(gtilde) *
368  nn * contract(gtilde_con, 0, star_i.dcov_logn, 0) +
369  4. / psi4 * nn * contract(gtilde_con, 0, star_i.dcov_logn, 0) *
370  phi_auto.derive_con(gtilde) ;
371 
372  source_Sij += - nn / (3.*psi4) * gtilde_con *
373  ( 0.25 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
374  dcov_gtilde, 0, 1), 0, 1)
375  - 0.5 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
376  dcov_gtilde, 0, 2), 0, 1)) ;
377 
378  source_Sij += - 8.*nn / (3.*psi4) * gtilde_con *
379  contract(dcov_phi_auto, 0, contract(gtilde_con, 0, star_i.dcov_phi, 0), 0) ;
380 
381  tens_temp = nn / (3.*psi4) * hdirac.divergence(flat)*star_i.hij_auto ;
382  tens_temp.inc_dzpuis() ;
383 
384  source_Sij += tens_temp ;
385 
386  source_Sij += - 8./(3.*psi4) * contract(dcov_phi_auto, 0,
387  nn*contract(gtilde_con, 0, star_i.dcov_logn, 0), 0) * gtilde_con ;
388 
389  source_Sij += 2.*nn* contract(gtilde_cov, 0, 1, tkij_a *
390  (tkij_a+tkij_c), 1, 3) ;
391 
392  source_Sij += - 2. * qpig * nn * ( psi4 * star_i.stress_euler
393  - 0.33333333333333333 * star_i.s_euler * gtilde_con ) ;
394 
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) ;
398 
399  source_Sij += - 0.5/(psi4*psi2) * contract(contract(star_i.hij, 1,
400  dcov_hij_auto, 2), 1, dcov_qq, 0) -
401  0.5/(psi4*psi2) * contract(contract(dcov_hij_auto, 2,
402  star_i.hij, 1), 1, dcov_qq, 0) ;
403 
404  source_Sij += 0.5/(psi4*psi2) * contract(contract(star_i.hij, 0,
405  dcov_hij_auto, 2), 0, dcov_qq, 0) ;
406 
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)
409  *gtilde_con ;
410 
411  source_Sij += 1./(3.*psi4*psi2) * contract(hdirac, 0,
412  dcov_qq, 0)*star_i.hij_auto ;
413 
414  // Source terms for Rij
415  //---------------------
416 
417  source_Rij = contract(star_i.hij, 0, 1, dcov_hij_auto.derive_cov(flat),
418  2, 3) ;
419  source_Rij.inc_dzpuis() ;
420 
421 
422  source_Rij += - contract(star_i.hij_auto, 1, dcov_hdirac, 1) -
423  contract(dcov_hdirac, 1, star_i.hij_auto, 1) ;
424 
425  source_Rij += contract(hdirac, 0, dcov_hij_auto, 2) ;
426 
427  source_Rij += - contract(contract(dcov_hij_auto, 1, dcov_hij, 2),
428  1, 3) ;
429 
430  source_Rij += - contract(gtilde_cov, 0, 1, contract(contract(
431  gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
432 
433  source_Rij += contract(contract(contract(contract(gtilde_cov, 0,
434  dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
435  contract(contract(contract(contract(gtilde_cov, 0,
436  dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
437 
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) ;
440 
441  source_Rij = source_Rij * 0.5 ;
442 
443  for(int i=1; i<=3; i++)
444  for(int j=1; j<=i; j++) {
445 
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) ;
451  source_tot_hij.dec_dzpuis() ;
452 
453  source3 = 2.*psi4/nn * (source_3a(i,j) + source_3a(j,i)
454  + source_3b(i,j)) ;
455 
456  source_tot_hij = source_tot_hij + source3 ;
457 
458 
459  if (ll == 1){
460  if(i==1 && j==1) {
461  Scalar lapl (star_i.hij_auto(1,1).laplacian()) ;
462  lapl.dec_dzpuis() ;
463  lie_aij_1.set(1,1) = nn / (2.*psi4) *
464  (lapl-source_tot_hij) ;
465  }
466  if(i==2 && j==1) {
467  Scalar lapl (star_i.hij_auto(2,1).laplacian()) ;
468  lapl.dec_dzpuis() ;
469  lie_aij_1.set(2,1) = nn / (2.*psi4) *
470  (lapl-source_tot_hij) ;
471  }
472  if(i==3 && j==1) {
473  Scalar lapl (star_i.hij_auto(3,1).laplacian()) ;
474  lapl.dec_dzpuis() ;
475  lie_aij_1.set(3,1) = nn / (2.*psi4) *
476  (lapl-source_tot_hij) ;
477  }
478  if(i==2 && j==2) {
479  Scalar lapl (star_i.hij_auto(2,2).laplacian()) ;
480  lapl.dec_dzpuis() ;
481  lie_aij_1.set(2,2) = nn / (2.*psi4) *
482  (lapl-source_tot_hij) ;
483  }
484  if(i==3 && j==2) {
485  Scalar lapl (star_i.hij_auto(3,2).laplacian()) ;
486  lapl.dec_dzpuis() ;
487  lie_aij_1.set(3,2) = nn / (2.*psi4) *
488  (lapl-source_tot_hij) ;
489  }
490  if(i==3 && j==3) {
491  Scalar lapl (star_i.hij_auto(3,3).laplacian()) ;
492  lapl.dec_dzpuis() ;
493  lie_aij_1.set(3,3) = nn / (2.*psi4) *
494  (lapl-source_tot_hij) ;
495  }
496  }
497 
498  if (ll == 2){
499  if(i==1 && j==1) {
500  Scalar lapl (star_i.hij_auto(1,1).laplacian()) ;
501  lapl.dec_dzpuis() ;
502  lie_aij_2.set(1,1) = nn / (2.*psi4) *
503  (lapl-source_tot_hij) ;
504  }
505  if(i==2 && j==1) {
506  Scalar lapl (star_i.hij_auto(2,1).laplacian()) ;
507  lapl.dec_dzpuis() ;
508  lie_aij_2.set(2,1) = nn / (2.*psi4) *
509  (lapl-source_tot_hij) ;
510  }
511  if(i==3 && j==1) {
512  Scalar lapl (star_i.hij_auto(3,1).laplacian()) ;
513  lapl.dec_dzpuis() ;
514  lie_aij_2.set(3,1) = nn / (2.*psi4) *
515  (lapl-source_tot_hij) ;
516  }
517  if(i==2 && j==2) {
518  Scalar lapl (star_i.hij_auto(2,2).laplacian()) ;
519  lapl.dec_dzpuis() ;
520  lie_aij_2.set(2,2) = nn / (2.*psi4) *
521  (lapl-source_tot_hij) ;
522  }
523  if(i==3 && j==2) {
524  Scalar lapl (star_i.hij_auto(3,2).laplacian()) ;
525  lapl.dec_dzpuis() ;
526  lie_aij_2.set(3,2) = nn / (2.*psi4) *
527  (lapl-source_tot_hij) ;
528  }
529  if(i==3 && j==3) {
530  Scalar lapl (star_i.hij_auto(3,3).laplacian()) ;
531  lapl.dec_dzpuis() ;
532  lie_aij_2.set(3,3) = nn / (2.*psi4) *
533  (lapl-source_tot_hij) ;
534  }
535  }
536  }
537  }
538 
539  lie_aij_1.dec_dzpuis(3) ;
540  lie_aij_2.dec_dzpuis(3) ;
541 
542  int nz = star1.mp.get_mg()->get_nzone() ;
543 
544  // Construction of an auxiliar grid and mapping (Last domain is at lambda)
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 ;
551  bornes[0] = 0 ;
552 
553  Map_af mapping (*(star1.mp.get_mg()), bornes) ;
554 
555  delete [] bornes ;
556 
557 
558  Sym_tensor lie_aij2_1 (star1.mp, CON, star1.mp.get_bvect_cart()) ;
559  Sym_tensor lie_aij_tot_1 (star1.mp, CON, star1.mp.get_bvect_cart()) ;
560  Sym_tensor lie_aij_tot (mapping, CON, mapping.get_bvect_cart()) ;
561 
562  Scalar lie_K2_1 (star1.mp) ;
563  lie_K2_1.set_etat_qcq() ;
564  Scalar lie_K_tot_1 (star1.mp) ;
565 
566 
567  // Importation on the mapping 1
568  // -------------------------------
569 
570  lie_K2_1.import(lie_K_2) ;
571  lie_K_tot_1 = lie_K_1 + lie_K2_1 ;
572  lie_K_tot_1.inc_dzpuis(2) ;
573 
574  lie_aij_2.change_triad(star1.mp.get_bvect_cart()) ;
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()) ;
580  }
581 
582  lie_aij_tot_1 = lie_aij_1 + lie_aij2_1 ;
583  lie_aij_tot_1.inc_dzpuis(2) ;
584 
585 
586  Sym_tensor lie_kij_tot (star1.mp, CON, star1.mp.get_bvect_cart()) ;
587  lie_kij_tot = lie_aij_tot_1/star1.psi4 + 1./3.*star1.gamma.con()*
588  lie_K_tot_1 ;
589 
590 
591  cout << " IN THE CENTER OF STAR 1 " << endl
592  << " ----------------------- " << endl ;
593  /*
594  cout << " components xx, xy, yy, xz, yz, zz" << endl ;
595  for(int i=1; i<=3; i++)
596  for(int j=1; j<=i; j++) {
597  Scalar resu(lie_kij_tot(i,j)*lie_kij_tot(i,j)) ;
598  cout << "i = " << i << ", j = " << j << endl ;
599  cout << "norme de la diff " << endl
600  << norme(resu)/(nr*nt*np) << endl ;
601 
602  // Computation of the integral
603  // -----------------------------
604 
605  Tbl integral (nz) ;
606  integral.annule_hard() ;
607  Tbl integ (resu.integrale_domains()) ;
608  for (int mm=0; mm<nz; mm++)
609  for (int pp=0; pp<=mm; pp++)
610  integral.set(mm) += integ(pp) ;
611  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
612  // dimensionless quantity
613  }
614  */
615 
616  cout << "L2 norm of L_k K^{ab} " << endl ;
617  Scalar determinant (pow(star1.get_gamma().determinant(), 0.5)) ;
618  determinant.std_spectral_base() ;
619  Scalar resu(2.*contract(lie_kij_tot, 0, 1,
620  lie_kij_tot.up_down(star1.gamma), 0, 1)
621  *determinant) ;
622  Tbl integral (nz) ;
623  integral.annule_hard() ;
624  Tbl integ (resu.integrale_domains()) ;
625  for (int mm=0; mm<nz; mm++)
626  for (int pp=0; pp<=mm; pp++)
627  integral.set(mm) += integ(pp) ;
628  cout << sqrt(integral) * sqrt(mass_adm()*ggrav) << endl ;
629  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
630  // dimensionless quantity
631 
632 
633  cout << "omega = " << omega << endl ;
634  cout << "mass_adm = " << mass_adm() << endl ;
635 
636 
637  lie_kij_tot.dec_dzpuis(2) ;
638 
639  cout << "Position du centre de l'etoile x/lambda = "
640  << -star1.get_mp().get_ori_x() * omega / M_PI << " in Lorene units"
641  << endl ;
642 
643 
644 
645  // Importation on the mapping defined in the center of mass
646  // -------------------------------------------------------------
647 /*
648  lie_aij_tot_1.change_triad(mapping.get_bvect_cart()) ;
649  for(int i=1; i<=3; i++)
650  for(int j=1; j<=i; j++) {
651  lie_aij_tot.set(i,j).import(lie_aij_tot_1(i,j)) ;
652  lie_aij_tot.set(i,j).set_spectral_va().set_base(lie_aij_tot_1(i,j).
653  get_spectral_va().get_base()) ;
654  }
655 
656  lie_aij_tot.inc_dzpuis(2) ;
657 
658  cout << " IN THE CENTER OF MASS : " << endl
659  << " ----------------------- " << endl ;
660 
661  cout << " components xx, xy, yy, xz, yz, zz" << endl ;
662  for(int i=1; i<=3; i++)
663  for(int j=1; j<=i; j++) {
664  Scalar resu(lie_aij_tot(i,j)*lie_aij_tot(i,j)) ;
665  cout << "i = " << i << ", j = " << j << endl ;
666  cout << "norme de la diff " << endl
667  << norme(resu)/(nr*nt*np) << endl ;
668 
669  Tbl integral (nz) ;
670  integral.annule_hard() ;
671  Tbl integ (resu.integrale_domains()) ;
672  for (int mm=0; mm<nz; mm++)
673  for (int pp=0; pp<=mm; pp++)
674  integral.set(mm) += integ(pp) ;
675  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
676  // dimensionless quantity
677  }
678 
679  cout << "L2 norm of L_k K^{ab} " << endl ;
680  resu = contract(lie_aij_tot, 0, 1,
681  lie_aij_tot.up_down(star1.gtilde), 0, 1) ;
682  integral.annule_hard() ;
683  integ = resu.integrale_domains() ;
684  for (int mm=0; mm<nz; mm++)
685  for (int pp=0; pp<=mm; pp++)
686  integral.set(mm) += integ(pp) ;
687  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
688  // dimensionless quantity
689  */
690 
691  cout << "Omega M = " << omega * mass_adm()*ggrav << endl ;
692 
693 }
694 }
Sym_tensor hij_comp
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition: star.h:594
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Definition: star.h:555
Coord xa
Absolute x coordinate.
Definition: map.h:748
Metric for tensor calculation.
Definition: metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Map & mp
Mapping associated with the star.
Definition: star.h:180
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition: star.h:588
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Metric gamma
3-metric
Definition: star.h:235
Star_bin star1
First star of the system.
Definition: binary.h:80
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition: star.h:212
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Flat metric for tensor calculation.
Definition: metric.h:261
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binary.h:94
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
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.
Definition: itbl.h:122
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
const Metric & get_gamma() const
Returns the 3-metric .
Definition: star.h:409
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Vector beta
Shift vector.
Definition: star.h:228
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:137
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:527
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 ;
Definition: tenseur.C:1133
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
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 .
Definition: scalar_integ.C:82
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition: star.h:581
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) ...
Definition: star.h:562
Scalar lnq_auto
Scalar field generated principally by the star.
Definition: star.h:543
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
Sym_tensor tkij_comp
Part of the extrinsic curvature tensor generated by beta_comp.
Definition: star.h:606
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tensor handling.
Definition: tensor.h:294
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Definition: star.h:535
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 .
Definition: tensor.C:1011
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Class for stars in binary system.
Definition: star.h:483
virtual const Scalar & determinant() const
Returns the determinant.
Definition: metric.C:395
Metric gtilde
Conformal metric .
Definition: star.h:565
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:825
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Definition: star.h:538
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]...
Definition: binary.h:89
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Coord ya
Absolute y coordinate.
Definition: map.h:749
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1023
Multi-domain grid.
Definition: grilles.h:279
Affine radial mapping.
Definition: map.h:2048
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Scalar nn
Lapse function N .
Definition: star.h:225
Vector beta_comp
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:575
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition: star.h:600
Coord za
Absolute z coordinate.
Definition: map.h:750
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:570
Basic array class.
Definition: tbl.h:164
void helical()
Function testing the helical symmetry.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
double mass_adm() const
Total ADM mass.
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375
Star_bin star2
Second star of the system.
Definition: binary.h:83
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
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.
Definition: sym_tensor.C:363
Scalar psi4
Conformal factor .
Definition: star.h:552
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
Definition: star.h:557
Coord r
r coordinate centered on the grid
Definition: map.h:736