LORENE
et_bin_bhns_extr_velpot.C
1/*
2 * Method of class Et_bin_bhns_extr to compute the velocity scalar
3 * potential $\psi$ by solving the continuity equation
4 * in the Kerr-Schild background metric or in the conformally flat one
5 * with extreme mass ratio
6 *
7 * (see file et_bin_bhns_extr.h for documentation).
8 *
9 */
10
11/*
12 * Copyright (c) 2004-2005 Keisuke Taniguchi
13 *
14 * This file is part of LORENE.
15 *
16 * LORENE is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License version 2
18 * as published by the Free Software Foundation.
19 *
20 * LORENE is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
24 *
25 * You should have received a copy of the GNU General Public License
26 * along with LORENE; if not, write to the Free Software
27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 *
29 */
30
31
32
33/*
34 * $Id: et_bin_bhns_extr_velpot.C,v 1.4 2016/12/05 16:17:52 j_novak Exp $
35 * $Log: et_bin_bhns_extr_velpot.C,v $
36 * Revision 1.4 2016/12/05 16:17:52 j_novak
37 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38 *
39 * Revision 1.3 2014/10/13 08:52:55 j_novak
40 * Lorene classes and functions now belong to the namespace Lorene.
41 *
42 * Revision 1.2 2005/02/28 23:17:18 k_taniguchi
43 * Modification to include the case of the conformally flat background metric
44 *
45 * Revision 1.1 2004/11/30 20:51:56 k_taniguchi
46 * *** empty log message ***
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_velpot.C,v 1.4 2016/12/05 16:17:52 j_novak Exp $
50 *
51 */
52
53// C headers
54//#include <math.h>
55
56// Lorene headers
57#include "et_bin_bhns_extr.h"
58#include "scalar.h"
59#include "metrique.h"
60#include "etoile.h"
61#include "eos.h"
62#include "param.h"
63#include "coord.h"
64#include "unites.h"
65
66// Local prototype
67namespace Lorene {
68Cmp raccord_c1(const Cmp& uu, int l1) ;
69
70double Et_bin_bhns_extr::velocity_pot_extr(const double& mass,
71 const double& sepa,
72 int mermax, double precis,
73 double relax) {
74
75 using namespace Unites ;
76
77 if (kerrschild) {
78
79 if (eos.identify() == 5 || eos.identify() == 4 ||
80 eos.identify() == 3) {
81
82 cout << "Etoile_bin::velocity_pot_extr" << endl ;
83 cout << "The code has not prepared for this kind of EOS!!!"
84 << endl;
85 abort() ;
86
87 } // End of strange stars case
88 else {
89
90 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
91
92 //--------------------------------
93 // Specific relativistic enthalpy ---> hhh
94 //--------------------------------
95
96 Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
97 hhh.set_std_base() ;
98
99 //--------------------------------------------
100 // Computation of W^i = - A^2 h Gamma_n B^i/N
101 //--------------------------------------------
102
103 Tenseur www = - a_car * hhh * gam_euler * bsn ;
104
105 www.change_triad( mp.get_bvect_cart() ) ;
106 // components on the mapping Cartesian basis
107
108 //-------------------------------------------------
109 // Constant value of W^i at the center of the star
110 //-------------------------------------------------
111
112 Tenseur v_orb(mp, 1, CON, mp.get_bvect_cart()) ;
113
114 v_orb.set_etat_qcq() ;
115 for (int i=0; i<3; i++) {
116 v_orb.set(i) = www(i)(0, 0, 0, 0) ;
117 }
118
119 v_orb.set_triad( *(www.get_triad()) ) ;
120 v_orb.set_std_base() ;
121
122 // Some auxiliary terms
123 // --------------------
124
125 const Coord& xx = mp.x ;
126 const Coord& yy = mp.y ;
127 const Coord& zz = mp.z ;
128
129 Tenseur r_bh(mp) ;
130 r_bh.set_etat_qcq() ;
131 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
132 r_bh.set_std_base() ;
133
134 Tenseur msr(mp) ;
135 msr = ggrav * mass / r_bh ;
136 msr.set_std_base() ;
137
138 Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
139 lapse_bh2 = 1. / (1.+2.*msr) ;
140 lapse_bh2.set_std_base() ;
141
142 Tenseur lapse_bh3(mp) ; // lapse_bh * lapse_bh * lapse_bh
143 lapse_bh3 = 1./pow(1.+2.*msr, 1.5) ;
144 lapse_bh3.set_std_base() ;
145
146 Tenseur xx_con(mp, 1, CON, mp.get_bvect_cart()) ;
147 xx_con.set_etat_qcq() ;
148 xx_con.set(0) = xx + sepa ;
149 xx_con.set(1) = yy ;
150 xx_con.set(2) = zz ;
151 xx_con.set_std_base() ;
152
153 Tenseur xsr_con(mp, 1, CON, mp.get_bvect_cart()) ;
154 xsr_con = xx_con / r_bh ;
155 xsr_con.set_std_base() ;
156
157 Tenseur xx_cov(mp, 1, COV, mp.get_bvect_cart()) ;
158 xx_cov.set_etat_qcq() ;
159 xx_cov.set(0) = xx + sepa ;
160 xx_cov.set(1) = yy ;
161 xx_cov.set(2) = zz ;
162 xx_cov.set_std_base() ;
163
164 Tenseur xsr_cov(mp, 1, COV, mp.get_bvect_cart()) ;
165 xsr_cov = xx_cov / r_bh ;
166 xsr_cov.set_std_base() ;
167
168 // X_i/r_bh in the spherical coordinate
169 const Coord& rr = mp.r ;
170 const Coord& st = mp.sint ;
171 const Coord& ct = mp.cost ;
172 const Coord& sp = mp.sinp ;
173 const Coord& cp = mp.cosp ;
174
175 Tenseur rr_spher_cov(mp, 1, COV, mp.get_bvect_spher()) ;
176 rr_spher_cov.set_etat_qcq() ;
177 rr_spher_cov.set(0) = rr + sepa*st*cp ;
178 rr_spher_cov.set(1) = sepa*ct*cp ;
179 rr_spher_cov.set(2) = - sepa*sp ;
180 rr_spher_cov.set_std_base() ;
181
182 Tenseur xsr_spher_cov(mp, 1, COV, mp.get_bvect_spher()) ;
183 xsr_spher_cov = rr_spher_cov / r_bh ;
184 xsr_spher_cov.set_std_base() ;
185
186 // l^j D_j H_ent
187 Tenseur ldent = flat_scalar_prod(xsr_con, ent.gradient()) ;
188
189 // l^j D_j beta_auto
190 Tenseur ldbeta = flat_scalar_prod(xsr_con, beta_auto.gradient()) ;
191
192 // l^j D_j psi0
193 Tenseur ldpsi0 = flat_scalar_prod(xsr_con, psi0.gradient()) ;
194
195 // l^i D_i (l^j D_j psi0)
196 Tenseur ldldpsi0 = flat_scalar_prod(xsr_con, ldpsi0.gradient()) ;
197
198 //-------------------------------------------------
199 // Source and coefficients a,b for poisson_compact
200 // (idenpendent from psi0)
201 //-------------------------------------------------
202
203 Cmp dndh_log = eos.der_nbar_ent(ent(), nzet) ;
204
205 // In order to avoid any division by zero in the computation of
206 // zeta_h the value of dndh_log is set to 1
207 // in the external domains:
208 for (int l=nzet; l <= nzm1; l++) {
209 dndh_log.set(l) = 1 ;
210 }
211
212 double erreur ;
213
214 Tenseur zeta_h( ent() / dndh_log ) ;
215 zeta_h.set_std_base() ;
216
217 Tenseur tmp_zeta = 1 - unsurc2 * zeta_h ;
218 tmp_zeta.set_std_base() ;
219
220 Tenseur bb = tmp_zeta * (ent.gradient_spher()
221 - 2.*lapse_bh2 * msr * ldent
222 * xsr_spher_cov)
223 + unsurc2 * zeta_h * (beta_auto.gradient_spher()
224 - 2.*lapse_bh2 * msr * ldbeta
225 * xsr_spher_cov)
226 - unsurc2 * 2. * zeta_h * lapse_bh2 * lapse_bh2 * msr / r_bh
227 * (1.+4.*msr) * xsr_spher_cov ;
228
229 Tenseur entmb = ent - beta_auto ;
230
231 Tenseur source = flat_scalar_prod(www - v_orb, ent.gradient())
232 + unsurc2*zeta_h*( flat_scalar_prod(v_orb, entmb.gradient())
233 + flat_scalar_prod(www, gam_euler.gradient())
234 / gam_euler )
235 + 2.*lapse_bh2 * msr * flat_scalar_prod(xsr_cov, v_orb)
236 * flat_scalar_prod(xsr_con, ent.gradient())
237 + unsurc2 * 2. * zeta_h
238 * (lapse_bh2*msr*(ldldpsi0
239 - flat_scalar_prod(xsr_cov, v_orb)
240 * (flat_scalar_prod(xsr_con, entmb.gradient())
241 - lapse_bh2 * (1.+4.*msr) / r_bh))
242 + a_car * hhh * gam_euler * lapse_bh3 * msr * (1.+3.*msr)
243 / r_bh) ;
244
245 source.annule(nzet, nzm1) ;
246
247 //----------------------------------------------------
248 // Resolution by means of Map_radial::poisson_compact
249 //----------------------------------------------------
250
251 Param par ;
252 int niter ;
253 par.add_int(mermax) ;
254 par.add_double(precis, 0) ;
255 par.add_double(relax, 1) ;
256 par.add_int_mod(niter) ;
257
258 if (psi0.get_etat() == ETATZERO) {
259 psi0.set_etat_qcq() ;
260 psi0.set() = 0 ;
261 }
262
263 source.set().va.ylm() ;
264
265 mp.poisson_compact(source(), zeta_h(), bb, par, psi0.set() ) ;
266
267 //-----------------------
268 // Check of the solution
269 //-----------------------
270
271 Tenseur bb_dpsi0 = flat_scalar_prod( bb, psi0.gradient_spher() ) ;
272
273 Cmp oper = zeta_h() * psi0().laplacien() + bb_dpsi0() ;
274
275 source.set().va.ylm_i() ;
276
277 erreur = diffrel(oper, source())(0) ;
278
279 cout << "Check of the resolution of the continuity equation : "
280 << endl ;
281 cout << "norme(source) : " << norme(source())(0)
282 << " diff oper/source : " << erreur << endl ;
283
284 //--------------------------
285 // Computation of grad(psi)
286 //--------------------------
287
288 // The computation is done component by component because
289 // psi0.gradient() is a covariant vector, whereas v_orb is a
290 // contravariant one.
291
292 d_psi.set_etat_qcq() ;
293
294 for (int i=0; i<3; i++) {
295 d_psi.set(i) = (psi0.gradient())(i) + v_orb(i) ;
296 }
297
298 d_psi.set_triad( *(v_orb.get_triad()) ) ;
299
300 // C^1 continuation of d_psi outside the star
301 // (to ensure a smooth enthalpy field accross the stellar surface)
302 // ----------------------------------------------------------------
303
304 d_psi.annule(nzet, nzm1) ;
305 for (int i=0; i<3; i++) {
306 d_psi.set(i) = raccord_c1(d_psi(i), nzet) ;
307 }
308
309 assert( d_psi.get_triad() == &(mp.get_bvect_cart()) ) ;
310
311 d_psi.change_triad(ref_triad) ;
312
313 return erreur ;
314
315 }
316
317 }
318 else {
319
320 if (eos.identify() == 5 || eos.identify() == 4 ||
321 eos.identify() == 3) {
322
323 cout << "Etoile_bin::velocity_pot_extr" << endl ;
324 cout << "The code has not prepared for this kind of EOS!!!"
325 << endl;
326 abort() ;
327
328 } // End of strange stars case
329 else {
330
331 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
332
333 //--------------------------------
334 // Specific relativistic enthalpy ---> hhh
335 //--------------------------------
336
337 Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
338 hhh.set_std_base() ;
339
340 //--------------------------------------------
341 // Computation of W^i = - A^2 h Gamma_n B^i/N
342 //--------------------------------------------
343
344 Tenseur www = - a_car * hhh * gam_euler * bsn ;
345
346 www.change_triad( mp.get_bvect_cart() ) ;
347 // components on the mapping Cartesian basis
348
349 //-------------------------------------------------
350 // Constant value of W^i at the center of the star
351 //-------------------------------------------------
352
353 Tenseur v_orb(mp, 1, CON, mp.get_bvect_cart()) ;
354
355 v_orb.set_etat_qcq() ;
356 for (int i=0; i<3; i++) {
357 v_orb.set(i) = www(i)(0, 0, 0, 0) ;
358 }
359
360 v_orb.set_triad( *(www.get_triad()) ) ;
361 v_orb.set_std_base() ;
362
363 // Some auxiliary terms
364 // --------------------
365
366 const Coord& xx = mp.x ;
367 const Coord& yy = mp.y ;
368 const Coord& zz = mp.z ;
369
370 Tenseur r_bh(mp) ;
371 r_bh.set_etat_qcq() ;
372 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
373 r_bh.set_std_base() ;
374
375 Tenseur msr(mp) ;
376 msr = ggrav * mass / r_bh ;
377 msr.set_std_base() ;
378
379 Tenseur xx_cov(mp, 1, COV, mp.get_bvect_cart()) ;
380 xx_cov.set_etat_qcq() ;
381 xx_cov.set(0) = xx + sepa ;
382 xx_cov.set(1) = yy ;
383 xx_cov.set(2) = zz ;
384 xx_cov.set_std_base() ;
385
386 // X_i in the spherical coordinate
387 const Coord& rr = mp.r ;
388 const Coord& st = mp.sint ;
389 const Coord& ct = mp.cost ;
390 const Coord& sp = mp.sinp ;
391 const Coord& cp = mp.cosp ;
392
393 Tenseur rr_spher_cov(mp, 1, COV, mp.get_bvect_spher()) ;
394 rr_spher_cov.set_etat_qcq() ;
395 rr_spher_cov.set(0) = rr + sepa*st*cp ;
396 rr_spher_cov.set(1) = sepa*ct*cp ;
397 rr_spher_cov.set(2) = - sepa*sp ;
398 rr_spher_cov.set_std_base() ;
399
400 Tenseur tmp_bh(mp) ;
401 tmp_bh = 0.5 * msr * msr / (1.-0.25*msr*msr) / r_bh / r_bh ;
402 tmp_bh.set_std_base() ;
403
404 //-------------------------------------------------
405 // Source and coefficients a,b for poisson_compact
406 // (idenpendent from psi0)
407 //-------------------------------------------------
408
409 Cmp dndh_log = eos.der_nbar_ent(ent(), nzet) ;
410
411 // In order to avoid any division by zero in the computation of
412 // zeta_h the value of dndh_log is set to 1
413 // in the external domains:
414 for (int l=nzet; l <= nzm1; l++) {
415 dndh_log.set(l) = 1 ;
416 }
417
418 double erreur ;
419
420 Tenseur zeta_h( ent() / dndh_log ) ;
421 zeta_h.set_std_base() ;
422
423 Tenseur tmp_zeta = 1 - unsurc2 * zeta_h ;
424 tmp_zeta.set_std_base() ;
425
426 Tenseur bb = tmp_zeta * ent.gradient_spher()
427 + unsurc2 * zeta_h * (beta_auto.gradient_spher()
428 + tmp_bh * rr_spher_cov) ;
429
430 Tenseur entmb = ent - beta_auto ;
431
432 Tenseur source = flat_scalar_prod(www - v_orb, ent.gradient())
433 + unsurc2*zeta_h*( flat_scalar_prod(v_orb, entmb.gradient())
434 - tmp_bh * flat_scalar_prod(v_orb, xx_cov)
435 + flat_scalar_prod(www, gam_euler.gradient())
436 / gam_euler ) ;
437
438 source.annule(nzet, nzm1) ;
439
440 //----------------------------------------------------
441 // Resolution by means of Map_radial::poisson_compact
442 //----------------------------------------------------
443
444 Param par ;
445 int niter ;
446 par.add_int(mermax) ;
447 par.add_double(precis, 0) ;
448 par.add_double(relax, 1) ;
449 par.add_int_mod(niter) ;
450
451 if (psi0.get_etat() == ETATZERO) {
452 psi0.set_etat_qcq() ;
453 psi0.set() = 0 ;
454 }
455
456 source.set().va.ylm() ;
457
458 mp.poisson_compact(source(), zeta_h(), bb, par, psi0.set() ) ;
459
460 //-----------------------
461 // Check of the solution
462 //-----------------------
463
464 Tenseur bb_dpsi0 = flat_scalar_prod( bb, psi0.gradient_spher() ) ;
465
466 Cmp oper = zeta_h() * psi0().laplacien() + bb_dpsi0() ;
467
468 source.set().va.ylm_i() ;
469
470 erreur = diffrel(oper, source())(0) ;
471
472 cout << "Check of the resolution of the continuity equation : "
473 << endl ;
474 cout << "norme(source) : " << norme(source())(0)
475 << " diff oper/source : " << erreur << endl ;
476
477 //--------------------------
478 // Computation of grad(psi)
479 //--------------------------
480
481 // The computation is done component by component because
482 // psi0.gradient() is a covariant vector, whereas v_orb is a
483 // contravariant one.
484
485 d_psi.set_etat_qcq() ;
486
487 for (int i=0; i<3; i++) {
488 d_psi.set(i) = (psi0.gradient())(i) + v_orb(i) ;
489 }
490
491 d_psi.set_triad( *(v_orb.get_triad()) ) ;
492
493 // C^1 continuation of d_psi outside the star
494 // (to ensure a smooth enthalpy field accross the stellar surface)
495 // ----------------------------------------------------------------
496
497 d_psi.annule(nzet, nzm1) ;
498 for (int i=0; i<3; i++) {
499 d_psi.set(i) = raccord_c1(d_psi(i), nzet) ;
500 }
501
502 assert( d_psi.get_triad() == &(mp.get_bvect_cart()) ) ;
503
504 d_psi.change_triad(ref_triad) ;
505
506 return erreur ;
507
508 }
509
510 }
511
512}
513}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::set
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Lorene::Coord
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Lorene::Et_bin_bhns_extr::kerrschild
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
Definition et_bin_bhns_extr.h:78
Lorene::Et_bin_bhns_extr::velocity_pot_extr
double velocity_pot_extr(const double &mass, const double &sepa, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Definition et_bin_bhns_extr_velpot.C:70
Lorene::Etoile_bin::psi0
Tenseur psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case).
Definition etoile.h:836
Lorene::Etoile_bin::bsn
Tenseur bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition etoile.h:953
Lorene::Etoile_bin::d_psi
Tenseur d_psi
Gradient of (in the irrotational case) (Cartesian components with respect to ref_triad ).
Definition etoile.h:841
Lorene::Etoile_bin::ref_triad
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition etoile.h:831
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:435
Lorene::Etoile::eos
const Eos & eos
Equation of state of the stellar matter.
Definition etoile.h:454
Lorene::Etoile::gam_euler
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition etoile.h:474
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:432
Lorene::Etoile::ent
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case).
Definition etoile.h:460
Lorene::Etoile::beta_auto
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition etoile.h:509
Lorene::Etoile::a_car
Tenseur a_car
Total conformal factor .
Definition etoile.h:518
Lorene::Etoile::unsurc2
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition etoile.h:445
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:318
Lorene::Param::add_int_mod
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:388
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:249
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:304
Lorene::Tenseur::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:707
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:830
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:642
Lorene::Tenseur::annule
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition tenseur.C:906
Lorene::Tenseur::gradient
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates).
Definition tenseur.C:1548
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1176
Lorene::Tenseur::change_triad
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tenseur.C:674
Lorene::Tenseur::set_triad
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tenseur.C:680
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:141
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm().
Definition valeur_ylm_i.C:134
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:273
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:507
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:484
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::flat_scalar_prod
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Definition tenseur_operateur.C:656
Lorene
Lorene prototypes.
Definition app_hor.h:67
Unites
Standard units of space, time and mass.