LORENE
et_bin_bhns_extr_hydro.C
1/*
2 * Methods of the class Et_bin_bhns_extr for computing hydro quantities
3 * in the Kerr-Schild background metric or in the conformally flat one
4 * with extreme mass ratio
5 *
6 * (see file et_bin_bhns_extr.h for documentation).
7 *
8 */
9
10/*
11 * Copyright (c) 2004-2005 Keisuke Taniguchi
12 *
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License version 2
17 * as published by the Free Software Foundation.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30char et_bin_bhns_extr_hydro_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_hydro.C,v 1.3 2014/10/13 08:52:55 j_novak Exp $" ;
31
32/*
33 * $Id: et_bin_bhns_extr_hydro.C,v 1.3 2014/10/13 08:52:55 j_novak Exp $
34 * $Log: et_bin_bhns_extr_hydro.C,v $
35 * Revision 1.3 2014/10/13 08:52:55 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.2 2005/02/28 23:14:16 k_taniguchi
39 * Modification to include the case of the conformally flat background metric
40 *
41 * Revision 1.1 2004/11/30 20:49:34 k_taniguchi
42 * *** empty log message ***
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_hydro.C,v 1.3 2014/10/13 08:52:55 j_novak Exp $
46 *
47 */
48
49// Lorene headers
50#include "et_bin_bhns_extr.h"
51#include "etoile.h"
52#include "coord.h"
53#include "unites.h"
54
55namespace Lorene {
56void Et_bin_bhns_extr::hydro_euler_extr(const double& mass,
57 const double& sepa) {
58
59 using namespace Unites ;
60
61 if (kerrschild) {
62
63 int nz = mp.get_mg()->get_nzone() ;
64 int nzm1 = nz - 1 ;
65
66 //--------------------------------
67 // Specific relativistic enthalpy ---> hhh
68 //--------------------------------
69
70 Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
71 hhh.set_std_base() ;
72
73 //----------------------------------------
74 // Lorentz factor between the co-orbiting ---> gam0
75 // observer and the Eulerian one
76 //----------------------------------------
77
78 const Coord& xx = mp.x ;
79 const Coord& yy = mp.y ;
80 const Coord& zz = mp.z ;
81
82 Tenseur r_bh(mp) ;
83 r_bh.set_etat_qcq() ;
84 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
85 r_bh.set_std_base() ;
86
87 Tenseur xx_cov(mp, 1, COV, ref_triad) ;
88 xx_cov.set_etat_qcq() ;
89 xx_cov.set(0) = xx + sepa ;
90 xx_cov.set(1) = yy ;
91 xx_cov.set(2) = zz ;
92 xx_cov.set_std_base() ;
93
94 Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
95 xsr_cov = xx_cov / r_bh ;
96 xsr_cov.set_std_base() ;
97
98 Tenseur msr(mp) ;
99 msr = ggrav * mass / r_bh ;
100 msr.set_std_base() ;
101
102 Tenseur tmp1(mp) ;
103 tmp1.set_etat_qcq() ;
104 tmp1.set() = 0 ;
105 tmp1.set_std_base() ;
106
107 for (int i=0; i<3; i++)
108 tmp1.set() += xsr_cov(i) % bsn(i) ;
109
110 tmp1.set_std_base() ;
111
112 Tenseur tmp2 = 2.*msr % tmp1 % tmp1 ;
113 tmp2.set_std_base() ;
114
115 for (int i=0; i<3; i++)
116 tmp2.set() += bsn(i) % bsn(i) ;
117
118 tmp2 = a_car % tmp2 ;
119
120 Tenseur gam0 = 1 / sqrt( 1 - unsurc2*tmp2 ) ;
121 gam0.set_std_base() ;
122
123 //--------------------------------------------
124 // Lorentz factor and 3-velocity of the fluid
125 // with respect to the Eulerian observer
126 //--------------------------------------------
127
128 if (irrotational) {
129
131
132 Tenseur xx_con(mp, 1, CON, ref_triad) ;
133 xx_con.set_etat_qcq() ;
134 xx_con.set(0) = xx + sepa ;
135 xx_con.set(1) = yy ;
136 xx_con.set(2) = zz ;
137 xx_con.set_std_base() ;
138
139 Tenseur xsr_con(mp, 1, CON, ref_triad) ;
140 xsr_con = xx_con / r_bh ;
141 xsr_con.set_std_base() ;
142
143 Tenseur tmp3(mp) ;
144 tmp3.set_etat_qcq() ;
145 tmp3.set() = 0 ;
146 tmp3.set_std_base() ;
147
148 for (int i=0; i<3; i++)
149 tmp3.set() += xsr_con(i) % d_psi(i) ;
150
151 tmp3.set_std_base() ;
152
153 Tenseur tmp4 = -2.*msr % tmp3 % tmp3 / (1.+2.*msr) ;
154 tmp4.set_std_base() ;
155
156 for (int i=0; i<3; i++)
157 tmp4.set() += d_psi(i) % d_psi(i) ;
158
159 tmp4 = tmp4 / a_car ;
160
161 gam_euler = sqrt( 1 + unsurc2 * tmp4 / (hhh%hhh) ) ;
162
163 gam_euler.set_std_base() ; // sets the standard spectral bases
164 // for a scalar field
165
166 Tenseur vtmp1 = d_psi / ( hhh % gam_euler % a_car ) ;
167 // COV (a_car correction) -> CON
168 Tenseur vtmp2 = -2.* msr % tmp3 % xsr_con / (1.+2.*msr)
169 / ( hhh % gam_euler % a_car ) ;
170 // CON
171
172 // The assignment of u_euler is performed component by component
173 // because u_euler is contravariant and d_psi is covariant
174 if (vtmp1.get_etat() == ETATZERO) {
176 }
177 else {
178 assert(vtmp1.get_etat() == ETATQCQ) ;
180 for (int i=0; i<3; i++) {
181 u_euler.set(i) = vtmp1(i) + vtmp2(i) ;
182 }
183 u_euler.set_triad( *(vtmp1.get_triad()) ) ;
184 }
185
187
188 }
189 else { // Rigid rotation
190 // --------------
191
192 gam_euler = gam0 ;
193 gam_euler.set_std_base() ; // sets the standard spectral bases
194 // for a scalar field
195
196 u_euler = - bsn ;
197
198 }
199
200 //--------------------------------------------------------
201 // Energy density E with respect to the Eulerian observer
202 //--------------------------------------------------------
203
205
206 //------------------------------------------------------------------
207 // Trace of the stress tensor with respect to the Eulerian observer
208 //------------------------------------------------------------------
209
210 Tenseur stmp1(mp) ;
211 stmp1.set_etat_qcq() ;
212 stmp1.set() = 0 ;
213 stmp1.set_std_base() ;
214
215 for (int i=0; i<3; i++)
216 stmp1.set() += xsr_cov(i) % u_euler(i) ;
217
218 stmp1.set_std_base() ;
219
220 Tenseur stmp2 = 2.*msr % stmp1 % stmp1 ;
221 stmp2.set_std_base() ;
222
223 for (int i=0; i<3; i++)
224 stmp2.set() += u_euler(i) % u_euler(i) ;
225
226 stmp2 = a_car % stmp2 ;
227
228 s_euler = 3 * press + ( ener_euler + press ) * stmp2 ;
230
231 //--------------------------------------
232 // Lorentz factor between the fluid and ---> gam
233 // co-orbiting observers
234 //--------------------------------------
235
236 Tenseur gtmp = 2.*msr % tmp1 % stmp1 ; //<- bsn^i = - U_0^i
237 gtmp.set_std_base() ;
238
239 for (int i=0; i<3; i++)
240 gtmp.set() += bsn(i) % u_euler(i) ; //<- bsn^i = - U_0^i
241
242 gtmp = a_car % gtmp ;
243
244 Tenseur tmp = ( 1+unsurc2*gtmp ) ; //<- (minus->plus) because of U_0^i
245 tmp.set_std_base() ;
246 Tenseur gam = gam0 % gam_euler % tmp ;
247
248 //--------------------------------------------
249 // Spatial projection of the fluid 3-velocity
250 // with respect to the co-orbiting observer
251 //--------------------------------------------
252
253 wit_w = gam_euler / gam % u_euler + gam0 % bsn ;
254
255 wit_w.set_std_base() ; // set the bases for spectral expansions
256
257 wit_w.annule(nzm1) ; // zero in the ZEC
258
259 //-----------------------------------------
260 // Logarithm of the Lorentz factor between
261 // the fluid and co-orbiting observers
262 //-----------------------------------------
263
264 if (relativistic) {
265
266 loggam = log( gam ) ;
267 loggam.set_std_base() ; // set the bases for spectral expansions
268
269 }
270 else {
271 cout << "BH-NS binary systems should be relativistic !!!" << endl ;
272 abort() ;
273 }
274
275 //-------------------------------------------------
276 // Velocity fields set to zero in external domains
277 //-------------------------------------------------
278
279 loggam.annule(nzm1) ; // zero in the ZEC only
280
281 wit_w.annule(nzm1) ; // zero outside the star
282
283 u_euler.annule(nzm1) ; // zero outside the star
284
285 if (loggam.get_etat() != ETATZERO) {
286 (loggam.set()).set_dzpuis(0) ;
287 }
288
289 if (!irrotational) {
290
291 gam = 1 ;
292 loggam = 0 ;
293 wit_w = 0 ;
294
295 }
296
297 // The derived quantities are obsolete
298 // -----------------------------------
299
301
302 }
303 else {
304
305 int nz = mp.get_mg()->get_nzone() ;
306 int nzm1 = nz - 1 ;
307
308 //--------------------------------
309 // Specific relativistic enthalpy ---> hhh
310 //--------------------------------
311
312 Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
313 hhh.set_std_base() ;
314
315 //----------------------------------------
316 // Lorentz factor between the co-orbiting ---> gam0
317 // observer and the Eulerian one
318 //----------------------------------------
319
320 Tenseur gam0 = 1 / sqrt( 1 - unsurc2*sprod(bsn, bsn) ) ;
321 gam0.set_std_base() ;
322
323 //--------------------------------------------
324 // Lorentz factor and 3-velocity of the fluid
325 // with respect to the Eulerian observer
326 //--------------------------------------------
327
328 if (irrotational) {
329
331
333 / (hhh%hhh) ) ;
334
335 gam_euler.set_std_base() ; // sets the standard spectral bases
336 // for a scalar field
337
338 Tenseur vtmp = d_psi / ( hhh % gam_euler % a_car ) ;
339 // COV (a_car correction) -> CON
340
341 // The assignment of u_euler is performed component by component
342 // because u_euler is contravariant and d_psi is covariant
343 if (vtmp.get_etat() == ETATZERO) {
345 }
346 else {
347 assert(vtmp.get_etat() == ETATQCQ) ;
349 for (int i=0; i<3; i++) {
350 u_euler.set(i) = vtmp(i) ;
351 }
352 u_euler.set_triad( *(vtmp.get_triad()) ) ;
353 }
354
356
357 }
358 else { // Rigid rotation
359 // --------------
360
361 gam_euler = gam0 ;
362 gam_euler.set_std_base() ; // sets the standard spectral bases
363 // for a scalar field
364
365 u_euler = - bsn ;
366
367 }
368
369 //--------------------------------------------------------
370 // Energy density E with respect to the Eulerian observer
371 //--------------------------------------------------------
372
374
375 //------------------------------------------------------------------
376 // Trace of the stress tensor with respect to the Eulerian observer
377 //------------------------------------------------------------------
378
379 s_euler = 3 * press + ( ener_euler + press )
380 * sprod(u_euler, u_euler) ;
382
383 //--------------------------------------
384 // Lorentz factor between the fluid and ---> gam
385 // co-orbiting observers
386 //--------------------------------------
387
388 Tenseur tmp = ( 1+unsurc2*sprod(bsn, u_euler) ) ;
389 //<- (minus->plus) because of U_0^i
390 tmp.set_std_base() ;
391 Tenseur gam = gam0 % gam_euler % tmp ;
392
393 //--------------------------------------------
394 // Spatial projection of the fluid 3-velocity
395 // with respect to the co-orbiting observer
396 //--------------------------------------------
397
398 wit_w = gam_euler / gam % u_euler + gam0 % bsn ;
399
400 wit_w.set_std_base() ; // set the bases for spectral expansions
401
402 wit_w.annule(nzm1) ; // zero in the ZEC
403
404 //-----------------------------------------
405 // Logarithm of the Lorentz factor between
406 // the fluid and co-orbiting observers
407 //-----------------------------------------
408
409 if (relativistic) {
410
411 loggam = log( gam ) ;
412 loggam.set_std_base() ; // set the bases for spectral expansions
413
414 }
415 else {
416 cout << "BH-NS binary systems should be relativistic !!!" << endl ;
417 abort() ;
418 }
419
420 //-------------------------------------------------
421 // Velocity fields set to zero in external domains
422 //-------------------------------------------------
423
424 loggam.annule(nzm1) ; // zero in the ZEC only
425
426 wit_w.annule(nzm1) ; // zero outside the star
427
428 u_euler.annule(nzm1) ; // zero outside the star
429
430 if (loggam.get_etat() != ETATZERO) {
431 (loggam.set()).set_dzpuis(0) ;
432 }
433
434 if (!irrotational) {
435
436 gam = 1 ;
437 loggam = 0 ;
438 wit_w = 0 ;
439
440 }
441
442 // The derived quantities are obsolete
443 // -----------------------------------
444
446
447 }
448
449}
450}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
void hydro_euler_extr(const double &mass, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tenseur bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition etoile.h:950
Tenseur d_psi
Gradient of (in the irrotational case) (Cartesian components with respect to ref_triad )
Definition etoile.h:838
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition etoile.h:828
bool irrotational
true for an irrotational star, false for a corotating one
Definition etoile.h:822
virtual void del_deriv() const
Deletes all the derived quantities.
Definition etoile_bin.C:447
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition etoile.h:849
Tenseur wit_w
Spatial projection of the fluid 3-velocity with respect to the co-orbiting observer.
Definition etoile.h:844
virtual Tenseur sprod(const Tenseur &t1, const Tenseur &t2) const
Performs the scalar product of two tensors by contracting the last index of t1 with the first index o...
Definition etoile_bin.C:748
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition etoile.h:471
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
Tenseur press
Fluid pressure.
Definition etoile.h:461
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition etoile.h:457
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition etoile.h:442
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord z
z coordinate centered on the grid
Definition map.h:728
Coord x
x coordinate centered on the grid
Definition map.h:726
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition tenseur.C:900
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tenseur.C:674
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.