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 
30 char 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 
55 namespace Lorene {
56 void 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 
130  d_psi.set_std_base() ;
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 
204  ener_euler = gam_euler % gam_euler % ( ener + press ) - press ;
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 
330  d_psi.set_std_base() ;
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 
373  ener_euler = gam_euler % gam_euler % ( ener + press ) - press ;
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
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Coord z
z coordinate centered on the grid
Definition: map.h:728
Coord x
x coordinate centered on the grid
Definition: map.h:726
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:704
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.