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 
31 
32 /*
33  * $Id: et_bin_bhns_extr_hydro.C,v 1.4 2016/12/05 16:17:52 j_novak Exp $
34  * $Log: et_bin_bhns_extr_hydro.C,v $
35  * Revision 1.4 2016/12/05 16:17:52 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.3 2014/10/13 08:52:55 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.2 2005/02/28 23:14:16 k_taniguchi
42  * Modification to include the case of the conformally flat background metric
43  *
44  * Revision 1.1 2004/11/30 20:49:34 k_taniguchi
45  * *** empty log message ***
46  *
47  *
48  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_hydro.C,v 1.4 2016/12/05 16:17:52 j_novak Exp $
49  *
50  */
51 
52 // Lorene headers
53 #include "et_bin_bhns_extr.h"
54 #include "etoile.h"
55 #include "coord.h"
56 #include "unites.h"
57 
58 namespace Lorene {
59 void Et_bin_bhns_extr::hydro_euler_extr(const double& mass,
60  const double& sepa) {
61 
62  using namespace Unites ;
63 
64  if (kerrschild) {
65 
66  int nz = mp.get_mg()->get_nzone() ;
67  int nzm1 = nz - 1 ;
68 
69  //--------------------------------
70  // Specific relativistic enthalpy ---> hhh
71  //--------------------------------
72 
73  Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
74  hhh.set_std_base() ;
75 
76  //----------------------------------------
77  // Lorentz factor between the co-orbiting ---> gam0
78  // observer and the Eulerian one
79  //----------------------------------------
80 
81  const Coord& xx = mp.x ;
82  const Coord& yy = mp.y ;
83  const Coord& zz = mp.z ;
84 
85  Tenseur r_bh(mp) ;
86  r_bh.set_etat_qcq() ;
87  r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
88  r_bh.set_std_base() ;
89 
90  Tenseur xx_cov(mp, 1, COV, ref_triad) ;
91  xx_cov.set_etat_qcq() ;
92  xx_cov.set(0) = xx + sepa ;
93  xx_cov.set(1) = yy ;
94  xx_cov.set(2) = zz ;
95  xx_cov.set_std_base() ;
96 
97  Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
98  xsr_cov = xx_cov / r_bh ;
99  xsr_cov.set_std_base() ;
100 
101  Tenseur msr(mp) ;
102  msr = ggrav * mass / r_bh ;
103  msr.set_std_base() ;
104 
105  Tenseur tmp1(mp) ;
106  tmp1.set_etat_qcq() ;
107  tmp1.set() = 0 ;
108  tmp1.set_std_base() ;
109 
110  for (int i=0; i<3; i++)
111  tmp1.set() += xsr_cov(i) % bsn(i) ;
112 
113  tmp1.set_std_base() ;
114 
115  Tenseur tmp2 = 2.*msr % tmp1 % tmp1 ;
116  tmp2.set_std_base() ;
117 
118  for (int i=0; i<3; i++)
119  tmp2.set() += bsn(i) % bsn(i) ;
120 
121  tmp2 = a_car % tmp2 ;
122 
123  Tenseur gam0 = 1 / sqrt( 1 - unsurc2*tmp2 ) ;
124  gam0.set_std_base() ;
125 
126  //--------------------------------------------
127  // Lorentz factor and 3-velocity of the fluid
128  // with respect to the Eulerian observer
129  //--------------------------------------------
130 
131  if (irrotational) {
132 
133  d_psi.set_std_base() ;
134 
135  Tenseur xx_con(mp, 1, CON, ref_triad) ;
136  xx_con.set_etat_qcq() ;
137  xx_con.set(0) = xx + sepa ;
138  xx_con.set(1) = yy ;
139  xx_con.set(2) = zz ;
140  xx_con.set_std_base() ;
141 
142  Tenseur xsr_con(mp, 1, CON, ref_triad) ;
143  xsr_con = xx_con / r_bh ;
144  xsr_con.set_std_base() ;
145 
146  Tenseur tmp3(mp) ;
147  tmp3.set_etat_qcq() ;
148  tmp3.set() = 0 ;
149  tmp3.set_std_base() ;
150 
151  for (int i=0; i<3; i++)
152  tmp3.set() += xsr_con(i) % d_psi(i) ;
153 
154  tmp3.set_std_base() ;
155 
156  Tenseur tmp4 = -2.*msr % tmp3 % tmp3 / (1.+2.*msr) ;
157  tmp4.set_std_base() ;
158 
159  for (int i=0; i<3; i++)
160  tmp4.set() += d_psi(i) % d_psi(i) ;
161 
162  tmp4 = tmp4 / a_car ;
163 
164  gam_euler = sqrt( 1 + unsurc2 * tmp4 / (hhh%hhh) ) ;
165 
166  gam_euler.set_std_base() ; // sets the standard spectral bases
167  // for a scalar field
168 
169  Tenseur vtmp1 = d_psi / ( hhh % gam_euler % a_car ) ;
170  // COV (a_car correction) -> CON
171  Tenseur vtmp2 = -2.* msr % tmp3 % xsr_con / (1.+2.*msr)
172  / ( hhh % gam_euler % a_car ) ;
173  // CON
174 
175  // The assignment of u_euler is performed component by component
176  // because u_euler is contravariant and d_psi is covariant
177  if (vtmp1.get_etat() == ETATZERO) {
179  }
180  else {
181  assert(vtmp1.get_etat() == ETATQCQ) ;
183  for (int i=0; i<3; i++) {
184  u_euler.set(i) = vtmp1(i) + vtmp2(i) ;
185  }
186  u_euler.set_triad( *(vtmp1.get_triad()) ) ;
187  }
188 
190 
191  }
192  else { // Rigid rotation
193  // --------------
194 
195  gam_euler = gam0 ;
196  gam_euler.set_std_base() ; // sets the standard spectral bases
197  // for a scalar field
198 
199  u_euler = - bsn ;
200 
201  }
202 
203  //--------------------------------------------------------
204  // Energy density E with respect to the Eulerian observer
205  //--------------------------------------------------------
206 
207  ener_euler = gam_euler % gam_euler % ( ener + press ) - press ;
208 
209  //------------------------------------------------------------------
210  // Trace of the stress tensor with respect to the Eulerian observer
211  //------------------------------------------------------------------
212 
213  Tenseur stmp1(mp) ;
214  stmp1.set_etat_qcq() ;
215  stmp1.set() = 0 ;
216  stmp1.set_std_base() ;
217 
218  for (int i=0; i<3; i++)
219  stmp1.set() += xsr_cov(i) % u_euler(i) ;
220 
221  stmp1.set_std_base() ;
222 
223  Tenseur stmp2 = 2.*msr % stmp1 % stmp1 ;
224  stmp2.set_std_base() ;
225 
226  for (int i=0; i<3; i++)
227  stmp2.set() += u_euler(i) % u_euler(i) ;
228 
229  stmp2 = a_car % stmp2 ;
230 
231  s_euler = 3 * press + ( ener_euler + press ) * stmp2 ;
233 
234  //--------------------------------------
235  // Lorentz factor between the fluid and ---> gam
236  // co-orbiting observers
237  //--------------------------------------
238 
239  Tenseur gtmp = 2.*msr % tmp1 % stmp1 ; //<- bsn^i = - U_0^i
240  gtmp.set_std_base() ;
241 
242  for (int i=0; i<3; i++)
243  gtmp.set() += bsn(i) % u_euler(i) ; //<- bsn^i = - U_0^i
244 
245  gtmp = a_car % gtmp ;
246 
247  Tenseur tmp = ( 1+unsurc2*gtmp ) ; //<- (minus->plus) because of U_0^i
248  tmp.set_std_base() ;
249  Tenseur gam = gam0 % gam_euler % tmp ;
250 
251  //--------------------------------------------
252  // Spatial projection of the fluid 3-velocity
253  // with respect to the co-orbiting observer
254  //--------------------------------------------
255 
256  wit_w = gam_euler / gam % u_euler + gam0 % bsn ;
257 
258  wit_w.set_std_base() ; // set the bases for spectral expansions
259 
260  wit_w.annule(nzm1) ; // zero in the ZEC
261 
262  //-----------------------------------------
263  // Logarithm of the Lorentz factor between
264  // the fluid and co-orbiting observers
265  //-----------------------------------------
266 
267  if (relativistic) {
268 
269  loggam = log( gam ) ;
270  loggam.set_std_base() ; // set the bases for spectral expansions
271 
272  }
273  else {
274  cout << "BH-NS binary systems should be relativistic !!!" << endl ;
275  abort() ;
276  }
277 
278  //-------------------------------------------------
279  // Velocity fields set to zero in external domains
280  //-------------------------------------------------
281 
282  loggam.annule(nzm1) ; // zero in the ZEC only
283 
284  wit_w.annule(nzm1) ; // zero outside the star
285 
286  u_euler.annule(nzm1) ; // zero outside the star
287 
288  if (loggam.get_etat() != ETATZERO) {
289  (loggam.set()).set_dzpuis(0) ;
290  }
291 
292  if (!irrotational) {
293 
294  gam = 1 ;
295  loggam = 0 ;
296  wit_w = 0 ;
297 
298  }
299 
300  // The derived quantities are obsolete
301  // -----------------------------------
302 
304 
305  }
306  else {
307 
308  int nz = mp.get_mg()->get_nzone() ;
309  int nzm1 = nz - 1 ;
310 
311  //--------------------------------
312  // Specific relativistic enthalpy ---> hhh
313  //--------------------------------
314 
315  Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
316  hhh.set_std_base() ;
317 
318  //----------------------------------------
319  // Lorentz factor between the co-orbiting ---> gam0
320  // observer and the Eulerian one
321  //----------------------------------------
322 
323  Tenseur gam0 = 1 / sqrt( 1 - unsurc2*sprod(bsn, bsn) ) ;
324  gam0.set_std_base() ;
325 
326  //--------------------------------------------
327  // Lorentz factor and 3-velocity of the fluid
328  // with respect to the Eulerian observer
329  //--------------------------------------------
330 
331  if (irrotational) {
332 
333  d_psi.set_std_base() ;
334 
336  / (hhh%hhh) ) ;
337 
338  gam_euler.set_std_base() ; // sets the standard spectral bases
339  // for a scalar field
340 
341  Tenseur vtmp = d_psi / ( hhh % gam_euler % a_car ) ;
342  // COV (a_car correction) -> CON
343 
344  // The assignment of u_euler is performed component by component
345  // because u_euler is contravariant and d_psi is covariant
346  if (vtmp.get_etat() == ETATZERO) {
348  }
349  else {
350  assert(vtmp.get_etat() == ETATQCQ) ;
352  for (int i=0; i<3; i++) {
353  u_euler.set(i) = vtmp(i) ;
354  }
355  u_euler.set_triad( *(vtmp.get_triad()) ) ;
356  }
357 
359 
360  }
361  else { // Rigid rotation
362  // --------------
363 
364  gam_euler = gam0 ;
365  gam_euler.set_std_base() ; // sets the standard spectral bases
366  // for a scalar field
367 
368  u_euler = - bsn ;
369 
370  }
371 
372  //--------------------------------------------------------
373  // Energy density E with respect to the Eulerian observer
374  //--------------------------------------------------------
375 
376  ener_euler = gam_euler % gam_euler % ( ener + press ) - press ;
377 
378  //------------------------------------------------------------------
379  // Trace of the stress tensor with respect to the Eulerian observer
380  //------------------------------------------------------------------
381 
382  s_euler = 3 * press + ( ener_euler + press )
383  * sprod(u_euler, u_euler) ;
385 
386  //--------------------------------------
387  // Lorentz factor between the fluid and ---> gam
388  // co-orbiting observers
389  //--------------------------------------
390 
391  Tenseur tmp = ( 1+unsurc2*sprod(bsn, u_euler) ) ;
392  //<- (minus->plus) because of U_0^i
393  tmp.set_std_base() ;
394  Tenseur gam = gam0 % gam_euler % tmp ;
395 
396  //--------------------------------------------
397  // Spatial projection of the fluid 3-velocity
398  // with respect to the co-orbiting observer
399  //--------------------------------------------
400 
401  wit_w = gam_euler / gam % u_euler + gam0 % bsn ;
402 
403  wit_w.set_std_base() ; // set the bases for spectral expansions
404 
405  wit_w.annule(nzm1) ; // zero in the ZEC
406 
407  //-----------------------------------------
408  // Logarithm of the Lorentz factor between
409  // the fluid and co-orbiting observers
410  //-----------------------------------------
411 
412  if (relativistic) {
413 
414  loggam = log( gam ) ;
415  loggam.set_std_base() ; // set the bases for spectral expansions
416 
417  }
418  else {
419  cout << "BH-NS binary systems should be relativistic !!!" << endl ;
420  abort() ;
421  }
422 
423  //-------------------------------------------------
424  // Velocity fields set to zero in external domains
425  //-------------------------------------------------
426 
427  loggam.annule(nzm1) ; // zero in the ZEC only
428 
429  wit_w.annule(nzm1) ; // zero outside the star
430 
431  u_euler.annule(nzm1) ; // zero outside the star
432 
433  if (loggam.get_etat() != ETATZERO) {
434  (loggam.set()).set_dzpuis(0) ;
435  }
436 
437  if (!irrotational) {
438 
439  gam = 1 ;
440  loggam = 0 ;
441  wit_w = 0 ;
442 
443  }
444 
445  // The derived quantities are obsolete
446  // -----------------------------------
447 
449 
450  }
451 
452 }
453 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:916
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur &#39;s...
Definition: etoile.h:831
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:690
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1186
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:471
Tenseur wit_w
Spatial projection of the fluid 3-velocity with respect to the co-orbiting observer.
Definition: etoile.h:847
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition: etoile.h:445
Tenseur press
Fluid pressure.
Definition: etoile.h:464
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...
bool irrotational
true for an irrotational star, false for a corotating one
Definition: etoile.h:825
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:477
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: etoile_bin.C:450
Tenseur bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition: etoile.h:953
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:707
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:474
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int get_etat() const
Returns the logical state.
Definition: tenseur.h:710
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
Tenseur a_car
Total conformal factor .
Definition: etoile.h:518
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:440
Tenseur ener
Total energy density in the fluid frame.
Definition: etoile.h:463
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: etoile.h:852
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:460
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord x
x coordinate centered on the grid
Definition: map.h:744
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:468
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:661
Coord z
z coordinate centered on the grid
Definition: map.h:746
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:751
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Tenseur d_psi
Gradient of (in the irrotational case) (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:841