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
28char star_bhns_hydro_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_hydro.C,v 1.3 2014/10/13 08:53:41 j_novak Exp $" ;
29
30/*
31 * $Id: star_bhns_hydro.C,v 1.3 2014/10/13 08:53:41 j_novak Exp $
32 * $Log: star_bhns_hydro.C,v $
33 * Revision 1.3 2014/10/13 08:53:41 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.2 2014/10/06 15:13:16 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.1 2007/06/22 01:31:42 k_taniguchi
40 * *** empty log message ***
41 *
42 *
43 * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_hydro.C,v 1.3 2014/10/13 08:53:41 j_novak Exp $
44 *
45 */
46
47// C++ headers
48//#include <>
49
50// C headers
51#include <cmath>
52
53// Lorene headers
54#include "star_bhns.h"
55#include "utilitaires.h"
56#include "unites.h"
57
58namespace Lorene {
59void Star_bhns::hydro_euler_bhns(bool kerrschild, const double& mass_bh,
60 const double& sepa) {
61
62 // Fundamental constants and units
63 // -------------------------------
64 using namespace Unites ;
65
66 int nz = mp.get_mg()->get_nzone() ;
67 int nzm1 = nz - 1 ;
68
69 //--------------------------------
70 // Specific relativistic enthalpy ---> hhh
71 //--------------------------------
72
73 Scalar hhh = exp(ent) ; // = 1 at the Newtonian limit
74 hhh.std_spectral_base() ;
75
76 if (kerrschild) {
77
78 double mass = ggrav * mass_bh ;
79
80 Scalar xx(mp) ;
81 xx = mp.x ;
82 xx.std_spectral_base() ;
83 Scalar yy(mp) ;
84 yy = mp.y ;
85 yy.std_spectral_base() ;
86 Scalar zz(mp) ;
87 zz = mp.z ;
88 zz.std_spectral_base() ;
89
90 double yns = mp.get_ori_y() ;
91
92 Scalar rbh(mp) ;
93 rbh = sqrt( (xx+sepa)*(xx+sepa) + (yy+yns)*(yy+yns) + zz*zz ) ;
94 rbh.std_spectral_base() ;
95
96 Vector ll(mp, CON, mp.get_bvect_cart()) ;
97 ll.set_etat_qcq() ;
98 ll.set(1) = (xx+sepa) / rbh ;
99 ll.set(2) = (yy+yns) / rbh ;
100 ll.set(3) = zz / rbh ;
101 ll.std_spectral_base() ;
102
103 Scalar msr(mp) ;
104 msr = mass / rbh ;
105 msr.std_spectral_base() ;
106
107 //--------------------------------------------
108 // Lorentz factor and 3-velocity of the fluid
109 // with respect to the Eulerian observer
110 //--------------------------------------------
111
112 if (irrotational) {
113
115
116 Scalar dpsi2(mp) ;
117 dpsi2 = d_psi(1)%d_psi(1) + d_psi(2)%d_psi(2)
118 + d_psi(3)%d_psi(3) ;
119 dpsi2.std_spectral_base() ;
120
121 Scalar lldpsi(mp) ;
122 lldpsi = ll(1)%d_psi(1) + ll(2)%d_psi(2) + ll(3)%d_psi(3) ;
123 lldpsi.std_spectral_base() ;
124
125 Scalar lap_bh2(mp) ;
126 lap_bh2 = 1. / (1.+2.*msr) ;
127 lap_bh2.std_spectral_base() ;
128
129 Scalar tmp1(mp) ;
130 tmp1 = 2. * msr % lap_bh2 % lldpsi % lldpsi ;
131 tmp1.std_spectral_base() ;
132
133 gam_euler = sqrt( 1.+(dpsi2-tmp1)/(hhh%hhh)/psi4 ) ;
135
138
139 for (int i=1; i<=3; i++)
140 u_euler.set(i) = (d_psi(i)-2.*msr%lap_bh2%lldpsi%ll(i))
141 / (hhh % gam_euler) / psi4 ;
142
144
145 }
146 else {
147
148 // Rigid rotation
149 // --------------
150
151 gam_euler = gam0 ;
153 u_euler = bsn ;
154
155 }
156
157 //--------------------------------------------------------
158 // Energy density E with respect to the Eulerian observer
159 // See Eq (53) from Gourgoulhon et al. (2001)
160 //--------------------------------------------------------
161
164
165 //---------------------------------------------------
166 // Trace of the stress-energy tensor with respect to
167 // the Eulerian observer
168 // See Eq (54) from Gourgoulhon et al. (2001)
169 //---------------------------------------------------
170
171 Scalar ueuler2(mp) ;
172 ueuler2 = u_euler(1)%u_euler(1) + u_euler(2)%u_euler(2)
173 + u_euler(3)%u_euler(3) ;
174 ueuler2.std_spectral_base() ;
175
177 llueuler = ll(1)%u_euler(1) + ll(2)%u_euler(2) + ll(3)%u_euler(3) ;
178 llueuler.std_spectral_base() ;
179
180 s_euler = 3. * press + (ener_euler + press) * psi4
181 * (ueuler2 + 2.*msr%llueuler%llueuler) ;
183
184 //--------------------------------------------
185 // Lorentz factor between the fluid and ---> gam
186 // co-orbiting observers
187 // See Eq (58) from Gourgoulhon et al. (2001)
188 //--------------------------------------------
189
190 if (irrotational) {
191
193 bsnueuler = bsn(1)%u_euler(1) + bsn(2)%u_euler(2)
194 + bsn(3)%u_euler(3) ;
195 bsnueuler.std_spectral_base() ;
196
197 Scalar llbsn(mp) ;
198 llbsn = ll(1)%bsn(1) + ll(2)%bsn(2) + ll(3)%bsn(3) ;
199 llbsn.std_spectral_base() ;
200
201 Scalar tmp2(mp) ;
202 tmp2 = 1. - psi4 * (bsnueuler + 2.*msr%llueuler%llbsn) ;
203 tmp2.std_spectral_base() ;
204
205 gam = gam0 % gam_euler % tmp2 ;
207
208 //--------------------------------------------
209 // Spetial projection of the fluid 3-velocity
210 // with respect to the co-orbiting observer
211 //--------------------------------------------
212
213 wit_w = gam_euler / gam * u_euler - gam0 * bsn ;
216
217 //-----------------------------------------
218 // Logarithm of the Lorentz factor between
219 // the fluid and co-orbiting observer
220 //-----------------------------------------
221
222 loggam = log( gam ) ;
224
225 //-------------------------------------------------
226 // Velocity fields set to zero in external domains
227 //-------------------------------------------------
228
232
233 loggam.set_dzpuis(0) ;
234
235 }
236 else { // Rigid rotation
237
238 gam = 1. ;
240 loggam = 0. ;
242
243 }
244
245 } // End of Kerr-Schild
246 else { // Isotropic coordinates with the maximal slicing
247
248 //--------------------------------------------
249 // Lorentz factor and 3-velocity of the fluid
250 // with respect to the Eulerian observer
251 //--------------------------------------------
252
253 if (irrotational) {
254
256
257 Scalar dpsi2(mp) ;
258 dpsi2 = d_psi(1)%d_psi(1) + d_psi(2)%d_psi(2)
259 + d_psi(3)%d_psi(3) ;
260 dpsi2.std_spectral_base() ;
261
262 gam_euler = sqrt( 1. + dpsi2/(hhh%hhh)/psi4 ) ;
264
266 for (int i=1; i<=3; i++)
267 u_euler.set(i) = d_psi(i)/(hhh%gam_euler)/psi4 ;
268
270
271 }
272 else {
273
274 // Rigid rotation
275 // --------------
276
277 gam_euler = gam0 ;
279 u_euler = bsn ;
280
281 }
282
283 //--------------------------------------------------------
284 // Energy density E with respect to the Eulerian observer
285 // See Eq (53) from Gourgoulhon et al. (2001)
286 //--------------------------------------------------------
287
290
291 //---------------------------------------------------
292 // Trace of the stress-energy tensor with respect to
293 // the Eulerian observer
294 // See Eq (54) from Gourgoulhon et al. (2001)
295 //---------------------------------------------------
296
297 Scalar ueuler2(mp) ;
298 ueuler2 = u_euler(1)%u_euler(1) + u_euler(2)%u_euler(2)
299 + u_euler(3)%u_euler(3) ;
300 ueuler2.std_spectral_base() ;
301
304
305
306 //--------------------------------------------
307 // Lorentz factor between the fluid and ---> gam
308 // co-orbiting observers
309 // See Eq (58) from Gourgoulhon et al. (2001)
310 //--------------------------------------------
311
312 if (irrotational) {
313
315 bsnueuler = bsn(1)%u_euler(1) + bsn(2)%u_euler(2)
316 + bsn(3)%u_euler(3) ;
317 bsnueuler.std_spectral_base() ;
318
319 Scalar tmp3(mp) ;
320 tmp3 = 1. - psi4 % bsnueuler ;
321 tmp3.std_spectral_base() ;
322
323 gam = gam0 % gam_euler % tmp3 ;
325
326 //--------------------------------------------
327 // Spetial projection of the fluid 3-velocity
328 // with respect to the co-orbiting observer
329 //--------------------------------------------
330
331 wit_w = gam_euler / gam * u_euler - gam0 * bsn ;
334
335 //-----------------------------------------
336 // Logarithm of the Lorentz factor between
337 // the fluid and co-orbiting observer
338 //-----------------------------------------
339
340 loggam = log( gam ) ;
342
343 //-------------------------------------------------
344 // Velocity fields set to zero in external domains
345 //-------------------------------------------------
346
350
351 loggam.set_dzpuis(0) ;
352
353 }
354 else { // Rigid rotation
355
356 gam = 1. ;
358 loggam = 0. ;
360
361 }
362
363 } // End of Isotropic coordinates
364
365 // The derived quantities are obsolete
366 // -----------------------------------
367
368 del_deriv() ;
369
370}
371}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord y
y coordinate centered on the grid
Definition map.h:727
double get_ori_y() const
Returns the y coordinate of the origin.
Definition map.h:770
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Vector wit_w
Spatial projection of the fluid 3-velocity with respect to the co-orbiting observer.
Definition star_bhns.h:88
Vector bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition star_bhns.h:99
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star_bhns.C:341
Scalar psi4
Fourth power of the total conformal factor.
Definition star_bhns.h:176
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition star_bhns.h:82
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
Definition star_bhns.h:107
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...
Scalar gam
Lorentz factor between the fluid and the co-orbiting observer.
Definition star_bhns.h:102
bool irrotational
true for an irrotational star, false for a corotating one
Definition star_bhns.h:72
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star_bhns.h:93
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition star.h:204
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Scalar press
Fluid pressure.
Definition star.h:194
Scalar ent
Log-enthalpy.
Definition star.h:190
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition star.h:207
Map & mp
Mapping associated with the star.
Definition star.h:180
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.