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
78namespace Lorene {
79void Binary::helical(){
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}
Lorene::Binary::et
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
Lorene::Binary::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binary.h:94
Lorene::Binary::star2
Star_bin star2
Second star of the system.
Definition binary.h:83
Lorene::Binary::star1
Star_bin star1
First star of the system.
Definition binary.h:80
Lorene::Binary::mass_adm
double mass_adm() const
Total ADM mass.
Definition binary_global.C:106
Lorene::Binary::helical
void helical()
Function testing the helical symmetry.
Definition binary_helical.C:79
Lorene::Itbl
Basic integer array class.
Definition itbl.h:122
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case).
Definition itbl.h:247
Lorene::Map_af
Affine radial mapping.
Definition map.h:2042
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition metric.h:261
Lorene::Metric_flat::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric_flat.C:156
Lorene::Metric
Metric for tensor calculation.
Definition metric.h:90
Lorene::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:293
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:283
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition scalar_deriv.C:390
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition scalar_deriv.C:436
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:359
Lorene::Scalar::std_spectral_base
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
Lorene::Scalar::inc_dzpuis
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...
Definition scalar_r_manip.C:473
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version).
Definition scalar.h:610
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:386
Lorene::Scalar::set_outer_boundary
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:321
Lorene::Scalar::derive_con
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
Lorene::Scalar::integrale_domains
const Tbl & integrale_domains() const
Computes the integral in each domain of *this .
Definition scalar_integ.C:82
Lorene::Scalar::dec_dzpuis
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...
Definition scalar_r_manip.C:421
Lorene::Scalar::import
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition scalar_import.C:71
Lorene::Star_bin
Class for stars in binary system.
Definition star.h:483
Lorene::Star_bin::lnq_auto
Scalar lnq_auto
Scalar field generated principally by the star.
Definition star.h:543
Lorene::Star_bin::dcon_logn
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Definition star.h:538
Lorene::Star_bin::dcon_phi
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
Definition star.h:557
Lorene::Star_bin::hij_comp
Sym_tensor hij_comp
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition star.h:594
Lorene::Star_bin::beta_comp
Vector beta_comp
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition star.h:575
Lorene::Star_bin::dcov_logn
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Definition star.h:535
Lorene::Star_bin::logn_auto
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition star.h:527
Lorene::Star_bin::tkij_comp
Sym_tensor tkij_comp
Part of the extrinsic curvature tensor generated by beta_comp.
Definition star.h:606
Lorene::Star_bin::beta_auto
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition star.h:570
Lorene::Star_bin::hij
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition star.h:581
Lorene::Star_bin::hij_auto
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition star.h:588
Lorene::Star_bin::psi4
Scalar psi4
Conformal factor .
Definition star.h:552
Lorene::Star_bin::tkij_auto
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition star.h:600
Lorene::Star_bin::dcov_phi
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Definition star.h:555
Lorene::Star_bin::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition star.h:562
Lorene::Star_bin::gtilde
Metric gtilde
Conformal metric .
Definition star.h:565
Lorene::Star::logn
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Lorene::Star::nn
Scalar nn
Lapse function N .
Definition star.h:225
Lorene::Star::ener_euler
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
Lorene::Star::s_euler
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Lorene::Star::stress_euler
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition star.h:212
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition star.h:180
Lorene::Star::beta
Vector beta
Shift vector.
Definition star.h:228
Lorene::Sym_tensor
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:226
Lorene::Sym_tensor::derive_lie
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
Lorene::Sym_tensor::divergence
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition sym_tensor.C:352
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::annule_hard
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:375
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tensor
Tensor handling.
Definition tensor.h:294
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base ).
Definition valeur.C:813
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:319
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:384
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:299
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::Tensor_sym::derive_cov
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor_sym_calculus.C:195
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:817
Lorene::Tensor::up_down
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Definition tensor_calculus.C:308
Lorene::Tensor::derive_con
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
Lorene::Tensor::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tensor_change_triad.C:80
Lorene::Tensor::inc_dzpuis
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
Lorene::Tensor_sym::derive_con
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...
Definition tensor_sym_calculus.C:207
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1011
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:663
Lorene::Tensor::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:935
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:282
Lorene
Lorene prototypes.
Definition app_hor.h:67
Lorene::Map
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Definition map.C:142
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:732
Unites
Standard units of space, time and mass.