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 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 $" ;
30 
31 /*
32  * $Id: binary_helical.C,v 1.6 2014/10/13 08:52:45 j_novak Exp $
33  * $Log: binary_helical.C,v $
34  * Revision 1.6 2014/10/13 08:52:45 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.5 2014/10/06 15:12:59 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.4 2008/08/19 06:41:59 j_novak
41  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
42  * cast-type operations, and constant strings that must be defined as const char*
43  *
44  * Revision 1.3 2006/08/01 14:26:50 f_limousin
45  * Small changes
46  *
47  * Revision 1.2 2006/06/05 17:05:57 f_limousin
48  * *** empty log message ***
49  *
50  * Revision 1.1 2006/04/11 14:25:15 f_limousin
51  * New version of the code : improvement of the computation of some
52  * critical sources, estimation of the dirac gauge, helical symmetry...
53  *
54 
55  *
56  * $Header: /cvsroot/Lorene/C++/Source/Binary/binary_helical.C,v 1.6 2014/10/13 08:52:45 j_novak Exp $ *
57  */
58 
59 
60 // C headers
61 #include <cmath>
62 
63 // Lorene headers
64 #include "cmp.h"
65 #include "tenseur.h"
66 #include "metrique.h"
67 #include "binary.h"
68 #include "param.h"
69 #include "graphique.h"
70 #include "utilitaires.h"
71 #include "tensor.h"
72 #include "nbr_spx.h"
73 #include "unites.h"
74 
75 namespace Lorene {
77 
78  // Fundamental constants and units
79  // -------------------------------
80  using namespace Unites ;
81 
82 
83  Sym_tensor lie_aij_1 (star1.mp, CON, star1.mp.get_bvect_cart()) ;
84  Sym_tensor lie_aij_2 (star2.mp, CON, star2.mp.get_bvect_cart()) ;
85 
86  Scalar lie_K_1 (star1.mp) ;
87  Scalar lie_K_2 (star2.mp) ;
88 
89  for (int ll=1; ll<=2; ll++) {
90 
91  Star_bin star_i (*et[ll-1]) ;
92 
93  Map& mp = star_i.mp ;
94  const Mg3d* mg = mp.get_mg() ;
95  int nz = mg->get_nzone() ; // total number of domains
96 
97  Metric_flat flat = star_i.flat ;
98  Metric gtilde = star_i.gtilde ;
99  Scalar nn = star_i.nn ;
100  Scalar psi4 = star_i.psi4 ;
101 
102  // -------------------------------
103  // AUXILIARY QUANTITIES
104  // -------------------------------
105 
106  // Derivatives of N and logN
107  //--------------------------
108 
109  const Vector dcov_logn_auto = star_i.logn_auto.derive_cov(flat) ;
110 
111  Tensor dcovdcov_logn_auto = (star_i.logn_auto.derive_cov(flat))
112  .derive_cov(flat) ;
113  dcovdcov_logn_auto.inc_dzpuis() ;
114 
115  // Derivatives of lnq, phi and Q
116  //-------------------------------
117 
118  const Scalar phi (0.5 * (star_i.lnq - star_i.logn)) ;
119  const Scalar phi_auto (0.5 * (star_i.lnq_auto - star_i.logn_auto)) ;
120 
121  const Vector dcov_phi_auto = phi_auto.derive_cov(flat) ;
122 
123  const Vector dcov_lnq = 2*star_i.dcov_phi + star_i.dcov_logn ;
124  const Vector dcon_lnq = 2*star_i.dcon_phi + star_i.dcon_logn ;
125  const Vector dcov_lnq_auto = star_i.lnq_auto.derive_cov(flat) ;
126  Tensor dcovdcov_lnq_auto = dcov_lnq_auto.derive_cov(flat) ;
127  dcovdcov_lnq_auto.inc_dzpuis() ;
128 
129  Scalar qq = exp(star_i.lnq) ;
130  qq.std_spectral_base() ;
131  const Vector& dcov_qq = qq.derive_cov(flat) ;
132 
133  Tensor dcovdcov_beta_auto = star_i.beta_auto.derive_cov(flat)
134  .derive_cov(flat) ;
135  dcovdcov_beta_auto.inc_dzpuis() ;
136 
137 
138  // Derivatives of hij, gtilde...
139  //------------------------------
140 
141  Scalar psi2 (pow(star_i.psi4, 0.5)) ;
142  psi2.std_spectral_base() ;
143 
144  const Tensor& dcov_hij = star_i.hij.derive_cov(flat) ;
145  const Tensor& dcon_hij = star_i.hij.derive_con(flat) ;
146  const Tensor& dcov_hij_auto = star_i.hij_auto.derive_cov(flat) ;
147 
148  const Sym_tensor& gtilde_cov = star_i.gtilde.cov() ;
149  const Sym_tensor& gtilde_con = star_i.gtilde.con() ;
150  const Tensor& dcov_gtilde = gtilde_cov.derive_cov(flat) ;
151 
152  // H^i and its derivatives ( = O in Dirac gauge)
153  // ---------------------------------------------
154 
155  double lambda_dirac = 0. ;
156 
157  const Vector hdirac = lambda_dirac * star_i.hij.divergence(flat) ;
158  const Vector hdirac_auto = lambda_dirac *
159  star_i.hij_auto.divergence(flat) ;
160 
161  Tensor dcov_hdirac = lambda_dirac * hdirac.derive_cov(flat) ;
162  dcov_hdirac.inc_dzpuis() ;
163  Tensor dcov_hdirac_auto = lambda_dirac * hdirac_auto.derive_cov(flat) ;
164  dcov_hdirac_auto.inc_dzpuis() ;
165  Tensor dcon_hdirac_auto = lambda_dirac * hdirac_auto.derive_con(flat) ;
166  dcon_hdirac_auto.inc_dzpuis() ;
167 
168 
169  // Function exp(-(r-r_0)^2/sigma^2)
170  // --------------------------------
171 
172  double r0 = mp.val_r(nz-2, 1, 0, 0) ;
173  double sigma = 1.*r0 ;
174  double om = omega ;
175 
176  Scalar rr (mp) ;
177  rr = mp.r ;
178 
179  Scalar ff (mp) ;
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() ;
185 
186 
187  // omdsdp
188  // ---------------------------------
189 
190  Vector omdsdp (mp, CON, mp.get_bvect_cart()) ;
191  Scalar yya (mp) ;
192  yya = mp.ya ;
193  Scalar xxa (mp) ;
194  xxa = mp.xa ;
195  Scalar zza (mp) ;
196  zza = mp.za ;
197 
198  if (fabs(mp.get_rot_phi()) < 1e-10){
199  omdsdp.set(1) = - om * yya * ff ;
200  omdsdp.set(2) = om * xxa * ff ;
201  omdsdp.set(3).annule_hard() ;
202  }
203  else{
204  omdsdp.set(1) = om * yya * ff ;
205  omdsdp.set(2) = - om * xxa * ff ;
206  omdsdp.set(3).annule_hard() ;
207  }
208 
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() ;
212 
213 
214  // Computation of helical A^{ij}
215  // ------------------------------
216 
217  const Tensor& dbeta = star_i.beta_auto.derive_con(gtilde) ;
218  Scalar div_beta = star_i.beta_auto.divergence(gtilde) ;
219 
220  Sym_tensor tkij_a (star_i.tkij_auto) ;
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) ;
225  }
226 
227  tkij_a = tkij_a - star_i.hij_auto.derive_lie(omdsdp) ;
228  tkij_a = 0.5 * tkij_a / nn ;
229 
230  Sym_tensor tkij_auto_cov = tkij_a.up_down(gtilde) ;
231 
232  Scalar a2_auto = contract(tkij_auto_cov, 0, 1, tkij_a, 0, 1,true) ;
233 
234  // COMP
235  const Tensor& dbeta_comp = star_i.beta_comp.derive_con(gtilde) ;
236  Scalar divbeta_comp = star_i.beta_comp.divergence(gtilde) ;
237 
238  Sym_tensor tkij_c (star_i.tkij_comp) ;
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) ;
243  }
244 
245  tkij_c = tkij_c - star_i.hij_comp.derive_lie(omdsdp) ;
246  tkij_c = 0.5 * tkij_c / nn ;
247 
248  Scalar a2_comp = contract(tkij_auto_cov, 0, 1, tkij_c, 0, 1,true) ;
249 
250 
251 // tkij_a = star_i.tkij_auto ;
252 // tkij_c = star_i.tkij_comp ;
253 
254 
255  // Sources for N
256  // ---------------
257 
258  Scalar source1(mp) ;
259  Scalar source2(mp) ;
260  Scalar source3(mp) ;
261  Scalar source4(mp) ;
262  Scalar source_tot(mp) ;
263 
264  source1 = qpig * star_i.psi4 % (star_i.ener_euler + star_i.s_euler) ;
265 
266  source2 = star_i.psi4 % (a2_auto + a2_comp) ;
267 
268  source3 = - contract(dcov_logn_auto, 0, star_i.dcon_logn, 0, true)
269  - 2. * contract(contract(gtilde_con, 0, star_i.dcov_phi, 0),
270  0, dcov_logn_auto, 0, true) ;
271 
272  source4 = - contract(star_i.hij, 0, 1, dcovdcov_logn_auto +
273  dcov_logn_auto*star_i.dcov_logn, 0, 1) ;
274 
275  source_tot = source1 + source2 + source3 + source4 ;
276 
277  if (ll == 1) {
278  lie_K_1 = nn / star_i.psi4 * (source_tot - star_i.logn_auto
279  .laplacian()) ;
280  lie_K_1.dec_dzpuis(4) ;
281  }
282  if (ll == 2) {
283  lie_K_2 = nn / star_i.psi4 * (source_tot - star_i.logn_auto
284  .laplacian()) ;
285  lie_K_2.dec_dzpuis(4) ;
286  }
287 
288 
289  // Sources for hij
290  // --------------
291 
292  Scalar source_tot_hij(mp) ;
293  Tensor source_Sij(mp, 2, CON, mp.get_bvect_cart()) ;
294  Tensor source_Rij(mp, 2, CON, mp.get_bvect_cart()) ;
295  Tensor tens_temp(mp, 2, CON, mp.get_bvect_cart()) ;
296 
297  Tensor source_1 (mp, 2, CON, mp.get_bvect_cart()) ;
298  Tensor source_2 (mp, 2, CON, mp.get_bvect_cart()) ;
299  Tensor source_3a (mp, 2, CON, mp.get_bvect_cart()) ;
300  Tensor source_3b (mp, 2, CON, mp.get_bvect_cart()) ;
301  Tensor source_4 (mp, 2, CON, mp.get_bvect_cart()) ;
302  Tensor source_5 (mp, 2, CON, mp.get_bvect_cart()) ;
303  Tensor source_6 (mp, 2, CON, mp.get_bvect_cart()) ;
304 
305 
306  source_1 = contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
307 
308  source_2 = - contract(dcon_hij, 2, dcov_lnq_auto, 0)
309  - 2./3. * contract(hdirac, 0, dcov_lnq_auto, 0) * flat.con() ;
310 
311  // Lie derivative of A^{ij}
312  // --------------------------
313 
314  Scalar decouple_logn = (star_i.logn_auto - 1.e-8)/
315  (star_i.logn - 2.e-8) ;
316 
317  // Construction of Omega d/dphi
318  // ----------------------------
319 
320  // Construction of D_k \Phi^i
321  Itbl type (2) ;
322  type.set(0) = CON ;
323  type.set(1) = COV ;
324 
325  Tensor dcov_omdsdphi (mp, 2, type, mp.get_bvect_cart()) ;
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. ;
335  dcov_omdsdphi.std_spectral_base() ;
336 
337  source_3a = contract(tkij_a, 0, dcov_omdsdphi, 1) ;
338  source_3a.inc_dzpuis() ;
339 
340  // Source 3b
341  // ------------
342 
343  source_3b = - contract(omdsdp, 0, tkij_a.derive_cov(flat), 2) ;
344 
345 
346  // Source 4
347  // ---------
348 
349  source_4 = - tkij_a.derive_lie(star_i.beta) ;
350  source_4.inc_dzpuis() ;
351  source_4 += - 2./3. * star_i.beta.divergence(flat) * tkij_a ;
352 
353  source_5 = dcon_hdirac_auto ;
354 
355  source_6 = - 2./3. * hdirac_auto.divergence(flat) * flat.con() ;
356  source_6.inc_dzpuis() ;
357 
358  // Source terms for Sij
359  //---------------------
360 
361  source_Sij = 8. * nn / psi4 * phi_auto.derive_con(gtilde) *
362  contract(gtilde_con, 0, star_i.dcov_phi, 0) ;
363 
364  source_Sij += 4. / psi4 * phi_auto.derive_con(gtilde) *
365  nn * contract(gtilde_con, 0, star_i.dcov_logn, 0) +
366  4. / psi4 * nn * contract(gtilde_con, 0, star_i.dcov_logn, 0) *
367  phi_auto.derive_con(gtilde) ;
368 
369  source_Sij += - nn / (3.*psi4) * gtilde_con *
370  ( 0.25 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
371  dcov_gtilde, 0, 1), 0, 1)
372  - 0.5 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
373  dcov_gtilde, 0, 2), 0, 1)) ;
374 
375  source_Sij += - 8.*nn / (3.*psi4) * gtilde_con *
376  contract(dcov_phi_auto, 0, contract(gtilde_con, 0, star_i.dcov_phi, 0), 0) ;
377 
378  tens_temp = nn / (3.*psi4) * hdirac.divergence(flat)*star_i.hij_auto ;
379  tens_temp.inc_dzpuis() ;
380 
381  source_Sij += tens_temp ;
382 
383  source_Sij += - 8./(3.*psi4) * contract(dcov_phi_auto, 0,
384  nn*contract(gtilde_con, 0, star_i.dcov_logn, 0), 0) * gtilde_con ;
385 
386  source_Sij += 2.*nn* contract(gtilde_cov, 0, 1, tkij_a *
387  (tkij_a+tkij_c), 1, 3) ;
388 
389  source_Sij += - 2. * qpig * nn * ( psi4 * star_i.stress_euler
390  - 0.33333333333333333 * star_i.s_euler * gtilde_con ) ;
391 
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) ;
395 
396  source_Sij += - 0.5/(psi4*psi2) * contract(contract(star_i.hij, 1,
397  dcov_hij_auto, 2), 1, dcov_qq, 0) -
398  0.5/(psi4*psi2) * contract(contract(dcov_hij_auto, 2,
399  star_i.hij, 1), 1, dcov_qq, 0) ;
400 
401  source_Sij += 0.5/(psi4*psi2) * contract(contract(star_i.hij, 0,
402  dcov_hij_auto, 2), 0, dcov_qq, 0) ;
403 
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)
406  *gtilde_con ;
407 
408  source_Sij += 1./(3.*psi4*psi2) * contract(hdirac, 0,
409  dcov_qq, 0)*star_i.hij_auto ;
410 
411  // Source terms for Rij
412  //---------------------
413 
414  source_Rij = contract(star_i.hij, 0, 1, dcov_hij_auto.derive_cov(flat),
415  2, 3) ;
416  source_Rij.inc_dzpuis() ;
417 
418 
419  source_Rij += - contract(star_i.hij_auto, 1, dcov_hdirac, 1) -
420  contract(dcov_hdirac, 1, star_i.hij_auto, 1) ;
421 
422  source_Rij += contract(hdirac, 0, dcov_hij_auto, 2) ;
423 
424  source_Rij += - contract(contract(dcov_hij_auto, 1, dcov_hij, 2),
425  1, 3) ;
426 
427  source_Rij += - contract(gtilde_cov, 0, 1, contract(contract(
428  gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
429 
430  source_Rij += contract(contract(contract(contract(gtilde_cov, 0,
431  dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
432  contract(contract(contract(contract(gtilde_cov, 0,
433  dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
434 
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) ;
437 
438  source_Rij = source_Rij * 0.5 ;
439 
440  for(int i=1; i<=3; i++)
441  for(int j=1; j<=i; j++) {
442 
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) ;
448  source_tot_hij.dec_dzpuis() ;
449 
450  source3 = 2.*psi4/nn * (source_3a(i,j) + source_3a(j,i)
451  + source_3b(i,j)) ;
452 
453  source_tot_hij = source_tot_hij + source3 ;
454 
455 
456  if (ll == 1){
457  if(i==1 && j==1) {
458  Scalar lapl (star_i.hij_auto(1,1).laplacian()) ;
459  lapl.dec_dzpuis() ;
460  lie_aij_1.set(1,1) = nn / (2.*psi4) *
461  (lapl-source_tot_hij) ;
462  }
463  if(i==2 && j==1) {
464  Scalar lapl (star_i.hij_auto(2,1).laplacian()) ;
465  lapl.dec_dzpuis() ;
466  lie_aij_1.set(2,1) = nn / (2.*psi4) *
467  (lapl-source_tot_hij) ;
468  }
469  if(i==3 && j==1) {
470  Scalar lapl (star_i.hij_auto(3,1).laplacian()) ;
471  lapl.dec_dzpuis() ;
472  lie_aij_1.set(3,1) = nn / (2.*psi4) *
473  (lapl-source_tot_hij) ;
474  }
475  if(i==2 && j==2) {
476  Scalar lapl (star_i.hij_auto(2,2).laplacian()) ;
477  lapl.dec_dzpuis() ;
478  lie_aij_1.set(2,2) = nn / (2.*psi4) *
479  (lapl-source_tot_hij) ;
480  }
481  if(i==3 && j==2) {
482  Scalar lapl (star_i.hij_auto(3,2).laplacian()) ;
483  lapl.dec_dzpuis() ;
484  lie_aij_1.set(3,2) = nn / (2.*psi4) *
485  (lapl-source_tot_hij) ;
486  }
487  if(i==3 && j==3) {
488  Scalar lapl (star_i.hij_auto(3,3).laplacian()) ;
489  lapl.dec_dzpuis() ;
490  lie_aij_1.set(3,3) = nn / (2.*psi4) *
491  (lapl-source_tot_hij) ;
492  }
493  }
494 
495  if (ll == 2){
496  if(i==1 && j==1) {
497  Scalar lapl (star_i.hij_auto(1,1).laplacian()) ;
498  lapl.dec_dzpuis() ;
499  lie_aij_2.set(1,1) = nn / (2.*psi4) *
500  (lapl-source_tot_hij) ;
501  }
502  if(i==2 && j==1) {
503  Scalar lapl (star_i.hij_auto(2,1).laplacian()) ;
504  lapl.dec_dzpuis() ;
505  lie_aij_2.set(2,1) = nn / (2.*psi4) *
506  (lapl-source_tot_hij) ;
507  }
508  if(i==3 && j==1) {
509  Scalar lapl (star_i.hij_auto(3,1).laplacian()) ;
510  lapl.dec_dzpuis() ;
511  lie_aij_2.set(3,1) = nn / (2.*psi4) *
512  (lapl-source_tot_hij) ;
513  }
514  if(i==2 && j==2) {
515  Scalar lapl (star_i.hij_auto(2,2).laplacian()) ;
516  lapl.dec_dzpuis() ;
517  lie_aij_2.set(2,2) = nn / (2.*psi4) *
518  (lapl-source_tot_hij) ;
519  }
520  if(i==3 && j==2) {
521  Scalar lapl (star_i.hij_auto(3,2).laplacian()) ;
522  lapl.dec_dzpuis() ;
523  lie_aij_2.set(3,2) = nn / (2.*psi4) *
524  (lapl-source_tot_hij) ;
525  }
526  if(i==3 && j==3) {
527  Scalar lapl (star_i.hij_auto(3,3).laplacian()) ;
528  lapl.dec_dzpuis() ;
529  lie_aij_2.set(3,3) = nn / (2.*psi4) *
530  (lapl-source_tot_hij) ;
531  }
532  }
533  }
534  }
535 
536  lie_aij_1.dec_dzpuis(3) ;
537  lie_aij_2.dec_dzpuis(3) ;
538 
539  int nz = star1.mp.get_mg()->get_nzone() ;
540 
541  // Construction of an auxiliar grid and mapping (Last domain is at lambda)
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 ;
548  bornes[0] = 0 ;
549 
550  Map_af mapping (*(star1.mp.get_mg()), bornes) ;
551 
552  delete [] bornes ;
553 
554 
555  Sym_tensor lie_aij2_1 (star1.mp, CON, star1.mp.get_bvect_cart()) ;
556  Sym_tensor lie_aij_tot_1 (star1.mp, CON, star1.mp.get_bvect_cart()) ;
557  Sym_tensor lie_aij_tot (mapping, CON, mapping.get_bvect_cart()) ;
558 
559  Scalar lie_K2_1 (star1.mp) ;
560  lie_K2_1.set_etat_qcq() ;
561  Scalar lie_K_tot_1 (star1.mp) ;
562 
563 
564  // Importation on the mapping 1
565  // -------------------------------
566 
567  lie_K2_1.import(lie_K_2) ;
568  lie_K_tot_1 = lie_K_1 + lie_K2_1 ;
569  lie_K_tot_1.inc_dzpuis(2) ;
570 
571  lie_aij_2.change_triad(star1.mp.get_bvect_cart()) ;
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()) ;
577  }
578 
579  lie_aij_tot_1 = lie_aij_1 + lie_aij2_1 ;
580  lie_aij_tot_1.inc_dzpuis(2) ;
581 
582 
583  Sym_tensor lie_kij_tot (star1.mp, CON, star1.mp.get_bvect_cart()) ;
584  lie_kij_tot = lie_aij_tot_1/star1.psi4 + 1./3.*star1.gamma.con()*
585  lie_K_tot_1 ;
586 
587 
588  cout << " IN THE CENTER OF STAR 1 " << endl
589  << " ----------------------- " << endl ;
590  /*
591  cout << " components xx, xy, yy, xz, yz, zz" << endl ;
592  for(int i=1; i<=3; i++)
593  for(int j=1; j<=i; j++) {
594  Scalar resu(lie_kij_tot(i,j)*lie_kij_tot(i,j)) ;
595  cout << "i = " << i << ", j = " << j << endl ;
596  cout << "norme de la diff " << endl
597  << norme(resu)/(nr*nt*np) << endl ;
598 
599  // Computation of the integral
600  // -----------------------------
601 
602  Tbl integral (nz) ;
603  integral.annule_hard() ;
604  Tbl integ (resu.integrale_domains()) ;
605  for (int mm=0; mm<nz; mm++)
606  for (int pp=0; pp<=mm; pp++)
607  integral.set(mm) += integ(pp) ;
608  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
609  // dimensionless quantity
610  }
611  */
612 
613  cout << "L2 norm of L_k K^{ab} " << endl ;
614  Scalar determinant (pow(star1.get_gamma().determinant(), 0.5)) ;
615  determinant.std_spectral_base() ;
616  Scalar resu(2.*contract(lie_kij_tot, 0, 1,
617  lie_kij_tot.up_down(star1.gamma), 0, 1)
618  *determinant) ;
619  Tbl integral (nz) ;
620  integral.annule_hard() ;
621  Tbl integ (resu.integrale_domains()) ;
622  for (int mm=0; mm<nz; mm++)
623  for (int pp=0; pp<=mm; pp++)
624  integral.set(mm) += integ(pp) ;
625  cout << sqrt(integral) * sqrt(mass_adm()*ggrav) << endl ;
626  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
627  // dimensionless quantity
628 
629 
630  cout << "omega = " << omega << endl ;
631  cout << "mass_adm = " << mass_adm() << endl ;
632 
633 
634  lie_kij_tot.dec_dzpuis(2) ;
635 
636  cout << "Position du centre de l'etoile x/lambda = "
637  << -star1.get_mp().get_ori_x() * omega / M_PI << " in Lorene units"
638  << endl ;
639 
640 
641 
642  // Importation on the mapping defined in the center of mass
643  // -------------------------------------------------------------
644 /*
645  lie_aij_tot_1.change_triad(mapping.get_bvect_cart()) ;
646  for(int i=1; i<=3; i++)
647  for(int j=1; j<=i; j++) {
648  lie_aij_tot.set(i,j).import(lie_aij_tot_1(i,j)) ;
649  lie_aij_tot.set(i,j).set_spectral_va().set_base(lie_aij_tot_1(i,j).
650  get_spectral_va().get_base()) ;
651  }
652 
653  lie_aij_tot.inc_dzpuis(2) ;
654 
655  cout << " IN THE CENTER OF MASS : " << endl
656  << " ----------------------- " << endl ;
657 
658  cout << " components xx, xy, yy, xz, yz, zz" << endl ;
659  for(int i=1; i<=3; i++)
660  for(int j=1; j<=i; j++) {
661  Scalar resu(lie_aij_tot(i,j)*lie_aij_tot(i,j)) ;
662  cout << "i = " << i << ", j = " << j << endl ;
663  cout << "norme de la diff " << endl
664  << norme(resu)/(nr*nt*np) << endl ;
665 
666  Tbl integral (nz) ;
667  integral.annule_hard() ;
668  Tbl integ (resu.integrale_domains()) ;
669  for (int mm=0; mm<nz; mm++)
670  for (int pp=0; pp<=mm; pp++)
671  integral.set(mm) += integ(pp) ;
672  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
673  // dimensionless quantity
674  }
675 
676  cout << "L2 norm of L_k K^{ab} " << endl ;
677  resu = contract(lie_aij_tot, 0, 1,
678  lie_aij_tot.up_down(star1.gtilde), 0, 1) ;
679  integral.annule_hard() ;
680  integ = resu.integrale_domains() ;
681  for (int mm=0; mm<nz; mm++)
682  for (int pp=0; pp<=mm; pp++)
683  integral.set(mm) += integ(pp) ;
684  cout << sqrt(integral) / sqrt(omega) << endl ; // To get
685  // dimensionless quantity
686  */
687 
688  cout << "Omega M = " << omega * mass_adm()*ggrav << endl ;
689 
690 }
691 }
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:730
Metric for tensor calculation.
Definition: metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
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:220
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:64
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:765
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:281
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:387
Base class for coordinate mappings.
Definition: map.h:670
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
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:784
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:775
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:134
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:1117
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
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:79
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
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
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition: scalar_manip.C:315
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
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:348
Tensor handling.
Definition: tensor.h:288
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:1002
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:392
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:816
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:280
Coord ya
Absolute y coordinate.
Definition: map.h:731
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:1014
Multi-domain grid.
Definition: grilles.h:273
Affine radial mapping.
Definition: map.h:2027
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:791
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
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:732
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:161
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:372
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:360
Scalar psi4
Conformal factor .
Definition: star.h:552
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
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:718