LORENE
star_bhns_hydro.C
1 /*
2  * Method of class Star_bhns to compute hydro quantities
3  *
4  * (see file star_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Keisuke Taniguchi
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 version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: star_bhns_hydro.C,v 1.4 2016/12/05 16:18:16 j_novak Exp $
32  * $Log: star_bhns_hydro.C,v $
33  * Revision 1.4 2016/12/05 16:18:16 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.3 2014/10/13 08:53:41 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.2 2014/10/06 15:13:16 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.1 2007/06/22 01:31:42 k_taniguchi
43  * *** empty log message ***
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_hydro.C,v 1.4 2016/12/05 16:18:16 j_novak Exp $
47  *
48  */
49 
50 // C++ headers
51 //#include <>
52 
53 // C headers
54 #include <cmath>
55 
56 // Lorene headers
57 #include "star_bhns.h"
58 #include "utilitaires.h"
59 #include "unites.h"
60 
61 namespace Lorene {
62 void Star_bhns::hydro_euler_bhns(bool kerrschild, const double& mass_bh,
63  const double& sepa) {
64 
65  // Fundamental constants and units
66  // -------------------------------
67  using namespace Unites ;
68 
69  int nz = mp.get_mg()->get_nzone() ;
70  int nzm1 = nz - 1 ;
71 
72  //--------------------------------
73  // Specific relativistic enthalpy ---> hhh
74  //--------------------------------
75 
76  Scalar hhh = exp(ent) ; // = 1 at the Newtonian limit
77  hhh.std_spectral_base() ;
78 
79  if (kerrschild) {
80 
81  double mass = ggrav * mass_bh ;
82 
83  Scalar xx(mp) ;
84  xx = mp.x ;
85  xx.std_spectral_base() ;
86  Scalar yy(mp) ;
87  yy = mp.y ;
88  yy.std_spectral_base() ;
89  Scalar zz(mp) ;
90  zz = mp.z ;
91  zz.std_spectral_base() ;
92 
93  double yns = mp.get_ori_y() ;
94 
95  Scalar rbh(mp) ;
96  rbh = sqrt( (xx+sepa)*(xx+sepa) + (yy+yns)*(yy+yns) + zz*zz ) ;
97  rbh.std_spectral_base() ;
98 
99  Vector ll(mp, CON, mp.get_bvect_cart()) ;
100  ll.set_etat_qcq() ;
101  ll.set(1) = (xx+sepa) / rbh ;
102  ll.set(2) = (yy+yns) / rbh ;
103  ll.set(3) = zz / rbh ;
104  ll.std_spectral_base() ;
105 
106  Scalar msr(mp) ;
107  msr = mass / rbh ;
108  msr.std_spectral_base() ;
109 
110  //--------------------------------------------
111  // Lorentz factor and 3-velocity of the fluid
112  // with respect to the Eulerian observer
113  //--------------------------------------------
114 
115  if (irrotational) {
116 
118 
119  Scalar dpsi2(mp) ;
120  dpsi2 = d_psi(1)%d_psi(1) + d_psi(2)%d_psi(2)
121  + d_psi(3)%d_psi(3) ;
122  dpsi2.std_spectral_base() ;
123 
124  Scalar lldpsi(mp) ;
125  lldpsi = ll(1)%d_psi(1) + ll(2)%d_psi(2) + ll(3)%d_psi(3) ;
126  lldpsi.std_spectral_base() ;
127 
128  Scalar lap_bh2(mp) ;
129  lap_bh2 = 1. / (1.+2.*msr) ;
130  lap_bh2.std_spectral_base() ;
131 
132  Scalar tmp1(mp) ;
133  tmp1 = 2. * msr % lap_bh2 % lldpsi % lldpsi ;
134  tmp1.std_spectral_base() ;
135 
136  gam_euler = sqrt( 1.+(dpsi2-tmp1)/(hhh%hhh)/psi4 ) ;
138 
140  assert( u_euler.get_triad() == bsn.get_triad() ) ;
141 
142  for (int i=1; i<=3; i++)
143  u_euler.set(i) = (d_psi(i)-2.*msr%lap_bh2%lldpsi%ll(i))
144  / (hhh % gam_euler) / psi4 ;
145 
147 
148  }
149  else {
150 
151  // Rigid rotation
152  // --------------
153 
154  gam_euler = gam0 ;
156  u_euler = bsn ;
157 
158  }
159 
160  //--------------------------------------------------------
161  // Energy density E with respect to the Eulerian observer
162  // See Eq (53) from Gourgoulhon et al. (2001)
163  //--------------------------------------------------------
164 
167 
168  //---------------------------------------------------
169  // Trace of the stress-energy tensor with respect to
170  // the Eulerian observer
171  // See Eq (54) from Gourgoulhon et al. (2001)
172  //---------------------------------------------------
173 
174  Scalar ueuler2(mp) ;
175  ueuler2 = u_euler(1)%u_euler(1) + u_euler(2)%u_euler(2)
176  + u_euler(3)%u_euler(3) ;
177  ueuler2.std_spectral_base() ;
178 
179  Scalar llueuler(mp) ;
180  llueuler = ll(1)%u_euler(1) + ll(2)%u_euler(2) + ll(3)%u_euler(3) ;
181  llueuler.std_spectral_base() ;
182 
183  s_euler = 3. * press + (ener_euler + press) * psi4
184  * (ueuler2 + 2.*msr%llueuler%llueuler) ;
186 
187  //--------------------------------------------
188  // Lorentz factor between the fluid and ---> gam
189  // co-orbiting observers
190  // See Eq (58) from Gourgoulhon et al. (2001)
191  //--------------------------------------------
192 
193  if (irrotational) {
194 
195  Scalar bsnueuler(mp) ;
196  bsnueuler = bsn(1)%u_euler(1) + bsn(2)%u_euler(2)
197  + bsn(3)%u_euler(3) ;
198  bsnueuler.std_spectral_base() ;
199 
200  Scalar llbsn(mp) ;
201  llbsn = ll(1)%bsn(1) + ll(2)%bsn(2) + ll(3)%bsn(3) ;
202  llbsn.std_spectral_base() ;
203 
204  Scalar tmp2(mp) ;
205  tmp2 = 1. - psi4 * (bsnueuler + 2.*msr%llueuler%llbsn) ;
206  tmp2.std_spectral_base() ;
207 
208  gam = gam0 % gam_euler % tmp2 ;
210 
211  //--------------------------------------------
212  // Spetial projection of the fluid 3-velocity
213  // with respect to the co-orbiting observer
214  //--------------------------------------------
215 
216  wit_w = gam_euler / gam * u_euler - gam0 * bsn ;
218  wit_w.annule_domain(nzm1) ;
219 
220  //-----------------------------------------
221  // Logarithm of the Lorentz factor between
222  // the fluid and co-orbiting observer
223  //-----------------------------------------
224 
225  loggam = log( gam ) ;
227 
228  //-------------------------------------------------
229  // Velocity fields set to zero in external domains
230  //-------------------------------------------------
231 
232  loggam.annule_domain(nzm1) ;
233  wit_w.annule_domain(nzm1) ;
234  u_euler.annule_domain(nzm1) ;
235 
236  loggam.set_dzpuis(0) ;
237 
238  }
239  else { // Rigid rotation
240 
241  gam = 1. ;
243  loggam = 0. ;
244  wit_w.set_etat_zero() ;
245 
246  }
247 
248  } // End of Kerr-Schild
249  else { // Isotropic coordinates with the maximal slicing
250 
251  //--------------------------------------------
252  // Lorentz factor and 3-velocity of the fluid
253  // with respect to the Eulerian observer
254  //--------------------------------------------
255 
256  if (irrotational) {
257 
259 
260  Scalar dpsi2(mp) ;
261  dpsi2 = d_psi(1)%d_psi(1) + d_psi(2)%d_psi(2)
262  + d_psi(3)%d_psi(3) ;
263  dpsi2.std_spectral_base() ;
264 
265  gam_euler = sqrt( 1. + dpsi2/(hhh%hhh)/psi4 ) ;
267 
269  for (int i=1; i<=3; i++)
270  u_euler.set(i) = d_psi(i)/(hhh%gam_euler)/psi4 ;
271 
273 
274  }
275  else {
276 
277  // Rigid rotation
278  // --------------
279 
280  gam_euler = gam0 ;
282  u_euler = bsn ;
283 
284  }
285 
286  //--------------------------------------------------------
287  // Energy density E with respect to the Eulerian observer
288  // See Eq (53) from Gourgoulhon et al. (2001)
289  //--------------------------------------------------------
290 
293 
294  //---------------------------------------------------
295  // Trace of the stress-energy tensor with respect to
296  // the Eulerian observer
297  // See Eq (54) from Gourgoulhon et al. (2001)
298  //---------------------------------------------------
299 
300  Scalar ueuler2(mp) ;
301  ueuler2 = u_euler(1)%u_euler(1) + u_euler(2)%u_euler(2)
302  + u_euler(3)%u_euler(3) ;
303  ueuler2.std_spectral_base() ;
304 
305  s_euler = 3.*press + (ener_euler+press)*psi4*ueuler2 ;
307 
308 
309  //--------------------------------------------
310  // Lorentz factor between the fluid and ---> gam
311  // co-orbiting observers
312  // See Eq (58) from Gourgoulhon et al. (2001)
313  //--------------------------------------------
314 
315  if (irrotational) {
316 
317  Scalar bsnueuler(mp) ;
318  bsnueuler = bsn(1)%u_euler(1) + bsn(2)%u_euler(2)
319  + bsn(3)%u_euler(3) ;
320  bsnueuler.std_spectral_base() ;
321 
322  Scalar tmp3(mp) ;
323  tmp3 = 1. - psi4 % bsnueuler ;
324  tmp3.std_spectral_base() ;
325 
326  gam = gam0 % gam_euler % tmp3 ;
328 
329  //--------------------------------------------
330  // Spetial projection of the fluid 3-velocity
331  // with respect to the co-orbiting observer
332  //--------------------------------------------
333 
334  wit_w = gam_euler / gam * u_euler - gam0 * bsn ;
336  wit_w.annule_domain(nzm1) ;
337 
338  //-----------------------------------------
339  // Logarithm of the Lorentz factor between
340  // the fluid and co-orbiting observer
341  //-----------------------------------------
342 
343  loggam = log( gam ) ;
345 
346  //-------------------------------------------------
347  // Velocity fields set to zero in external domains
348  //-------------------------------------------------
349 
350  loggam.annule_domain(nzm1) ;
351  wit_w.annule_domain(nzm1) ;
352  u_euler.annule_domain(nzm1) ;
353 
354  loggam.set_dzpuis(0) ;
355 
356  }
357  else { // Rigid rotation
358 
359  gam = 1. ;
361  loggam = 0. ;
362  wit_w.set_etat_zero() ;
363 
364  }
365 
366  } // End of Isotropic coordinates
367 
368  // The derived quantities are obsolete
369  // -----------------------------------
370 
371  del_deriv() ;
372 
373 }
374 }
Scalar psi4
Fourth power of the total conformal factor.
Definition: star_bhns.h:176
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
void hydro_euler_bhns(bool kerrschild, const double &mass_bh, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Map & mp
Mapping associated with the star.
Definition: star.h:180
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star_bhns.C:344
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Vector bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition: star_bhns.h:99
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:788
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Vector wit_w
Spatial projection of the fluid 3-velocity with respect to the co-orbiting observer.
Definition: star_bhns.h:88
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: star_bhns.h:93
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
Scalar ent
Log-enthalpy.
Definition: star.h:190
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Tensor field of valence 1.
Definition: vector.h:188
Scalar gam
Lorentz factor between the fluid and the co-orbiting observer.
Definition: star_bhns.h:102
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition: star_bhns.h:82
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
Scalar ener
Total energy density in the fluid frame.
Definition: star.h:193
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
Scalar press
Fluid pressure.
Definition: star.h:194
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
Definition: star_bhns.h:107
bool irrotational
true for an irrotational star, false for a corotating one
Definition: star_bhns.h:72
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord x
x coordinate centered on the grid
Definition: map.h:744
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
Coord z
z coordinate centered on the grid
Definition: map.h:746
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198