LORENE
et_bin_hydro.C
1/*
2 * Methods of the class Etoile_bin for computing hydro quantities
3 *
4 * (see file etoile.h for documentation)
5 */
6
7/*
8 * Copyright (c) 2000-2001 Eric Gourgoulhon
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
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
29char et_bin_hydro_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_hydro.C,v 1.6 2014/10/13 08:52:55 j_novak Exp $" ;
30
31/*
32 * $Id: et_bin_hydro.C,v 1.6 2014/10/13 08:52:55 j_novak Exp $
33 * $Log: et_bin_hydro.C,v $
34 * Revision 1.6 2014/10/13 08:52:55 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.5 2005/08/29 15:10:16 p_grandclement
38 * Addition of things needed :
39 * 1) For BBH with different masses
40 * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
41 * WORKING YET !!!
42 *
43 * Revision 1.4 2003/03/03 19:46:09 f_limousin
44 * Set standard bases for s_euler.
45 *
46 * Revision 1.3 2003/01/17 13:34:56 f_limousin
47 * Replace A^2*flat_scalar_prod by sprod
48 *
49 * Revision 1.2 2002/12/10 15:29:29 k_taniguchi
50 * Change the multiplication "*" to "%"
51 * and flat_scalar_prod to flat_scalar_prod_desal.
52 *
53 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
54 * LORENE
55 *
56 * Revision 2.11 2000/12/22 13:10:32 eric
57 * Modif des annulations en dehors de l'etoile.
58 *
59 * Revision 2.10 2000/03/29 11:47:53 eric
60 * Suppression affichage.
61 *
62 * Revision 2.9 2000/03/01 16:12:09 eric
63 * Appel de set_std_base sur u_euler dans le cas irrotationnel.
64 *
65 * Revision 2.8 2000/02/24 09:34:10 keisuke
66 * Ajout de l'appel a set_std_base() sur wit_w.
67 *
68 * Revision 2.7 2000/02/21 13:59:09 eric
69 * Appel de set_dzpuis(0) sur loggam.
70 *
71 * Revision 2.6 2000/02/21 08:43:17 keisuke
72 * Ajout de l'appel a set_std_base() sur loggam.
73 *
74 * Revision 2.5 2000/02/17 14:42:45 eric
75 * Ajout de l'appel a set_std_base() sur gam_euler.
76 *
77 * Revision 2.4 2000/02/14 11:06:15 eric
78 * Suppression de l'appel a annule(nzet.nzm1) sur gam_euler dans le cas en
79 * corotation.
80 * Ajout de l'appel a annule(nzet,nzm1) sur wit_w.
81 *
82 * Revision 2.3 2000/02/12 18:36:23 eric
83 * Appel de set_std_base sur loggam.
84 *
85 * Revision 2.2 2000/02/08 19:29:03 eric
86 * Calcul sur les tenseurs.
87 * wit_w est calcule.
88 *
89 * Revision 2.1 2000/02/04 16:37:28 eric
90 * Premiere version implementee (non testee !)/
91 *
92 * Revision 2.0 2000/01/31 15:57:30 eric
93 * *** empty log message ***
94 *
95 *
96 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_hydro.C,v 1.6 2014/10/13 08:52:55 j_novak Exp $
97 *
98 */
99
100// Headers C
101
102// Headers Lorene
103#include "etoile.h"
104
105namespace Lorene {
107
108 int nz = mp.get_mg()->get_nzone() ;
109 int nzm1 = nz - 1 ;
110
111 //----------------------------------
112 // Specific relativistic enthalpy ---> hhh
113 //----------------------------------
114
115 Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
116 hhh.set_std_base() ;
117
118 //---------------------------------------------------
119 // Lorentz factor between the co-orbiting ---> gam0
120 // observer and the Eulerian one
121 // See Eq (23) and (24) from Gourgoulhon et al. (2001)
122 //---------------------------------------------------
123
124 Tenseur gam0 = 1/sqrt(1-unsurc2*sprod(bsn,bsn)) ;
125 gam0.set_std_base() ;
126
127 //------------------------------------------
128 // Lorentz factor and 3-velocity of the fluid
129 // with respect to the Eulerian observer
130 //------------------------------------------
131
132 if (irrotational) {
133
134//## cout << "Etoile_bin::hydro_euler : ### warning : " << endl ;
135// cout << " d_psi.d_psi would be better computed in spher. coord. !"
136//## << endl ;
137
139
140 // See Eq (32) from Gourgoulhon et al. (2001)
142 / (hhh%hhh) ) ;
143
144 gam_euler.set_std_base() ; // sets the standard spectral bases for
145 // a scalar field
146
147
148 // See Eq (31) from Gourgoulhon et al. (2001)
149 Tenseur vtmp = d_psi / ( hhh % gam_euler % a_car ) ;
150
151 // The assignment of u_euler is performed component by component
152 // because u_euler is contravariant and d_psi is covariant
153 if (vtmp.get_etat() == ETATZERO) {
155 }
156 else{
157 assert(vtmp.get_etat() == ETATQCQ) ;
159 for (int i=0; i<3; i++) {
160 u_euler.set(i) = vtmp(i) ;
161 }
162 u_euler.set_triad( *(vtmp.get_triad()) ) ;
163 }
164
166
167 }
168 else { // Rigid rotation
169 // --------------
170
171 gam_euler = gam0 ;
172 gam_euler.set_std_base() ; // sets the standard spectral bases for
173 // a scalar field
174
175 u_euler = -bsn ;
176
177 }
178
179 //------------------------------------
180 // Energy density E with respect to the Eulerian observer
181 // See Eq (53) from Gourgoulhon et al. (2001)
182 //------------------------------------
183
185
186 //------------------------------------
187 // Trace of the stress tensor with respect to the Eulerian observer
188 // See Eq (54) from Gourgoulhon et al. (2001)
189 //------------------------------------
190
191 //Tenseur tmp00 = sprod(u_euler, u_euler) ;
192 //cout << hex << u_euler(0).va.base.b[0] << endl ;
193 //cout << hex << u_euler(1).va.base.b[0] << endl ;
194 //cout << hex << u_euler(2).va.base.b[0] << endl ;
195 //cout << hex << tmp00().va.base.b[0] << endl ;
196
197 s_euler = 3 * press + ( ener_euler + press ) *
200
201
202 //-------------------------------------------
203 // Lorentz factor between the fluid and ---> gam
204 // co-orbiting observers
205 // See Eq (58) from Gourgoulhon et al. (2001)
206 //--------------------------------------------
207
208 Tenseur tmp = ( 1+unsurc2*sprod(bsn,u_euler) ) ;
209 tmp.set_std_base() ;
210 Tenseur gam = gam0 % gam_euler % tmp ;
211
212 //-------------------------------------------
213 // Spatial projection of the fluid 3-velocity
214 // with respect to the co-orbiting observer
215 //--------------------------------------------
216
217 wit_w = gam_euler / gam % u_euler + gam0 % bsn ;
218
219 wit_w.set_std_base() ; // set the bases for spectral expansions
220
221 wit_w.annule(nzm1) ; // zero in the ZEC
222
223
224 //-------------------------------------------
225 // Logarithm of the Lorentz factor between
226 // the fluid and co-orbiting observers
227 //--------------------------------------------
228
229 if (relativistic) {
230
231 loggam = log( gam ) ;
232
233 loggam.set_std_base() ; // set the bases for spectral expansions
234 }
235 else {
236
238
239 loggam.set_std_base() ; // set the bases for spectral expansions
240
241//### Forcage a zero des termes en sin(m*phi) :
242// loggam.coef() ;
243// int np = mgrille->np[0] ;
244// int nt = mgrille->nt[0] ;
245// int nr = mgrille->nr[0] ;
246// int ntnr = nt * nr ;
247// for (int k = 1; k < np+2; k+=2) {
248// for (int j = 0; j < nt; j++) {
249// for (int i = 0; i<nr; i++) {
250// loggam.c_cf[0]->t[0]->t[ntnr*k + nr*j + i] = 0 ;
251// }
252// }
253// }
254// loggam.c_ajx[0] = 0 ;
255// loggam.c_aj = 0 ;
256// loggam.coef_i() ;
257//###
258 }
259
260
261 //-------------------------------------------------
262 // Velocity fields set to zero in external domains
263 //-------------------------------------------------
264
265 loggam.annule(nzm1) ; // zero in the ZEC only
266
267 wit_w.annule(nzm1) ; // zero outside the star
268
269 u_euler.annule(nzm1) ; // zero outside the star
270
271
272 if (loggam.get_etat() != ETATZERO) {
273 (loggam.set()).set_dzpuis(0) ;
274 }
275
276 //################
277 // verification: test on gam_euler
278
279 // if (irrotational) {
280
281 // Tenseur gam_test = 1. / sqrt( 1 - unsurc2 * sprod(u_euler, u_euler) ) ;
282
283 // cout << "Etoile_bin::hydro_euler: test of gam_euler : " << endl ;
284 // cout << " maximum error : " << endl ;
285 // cout << max(gam_test() - gam_euler()) << endl ;
286 //cout << " relative error : " << endl ;
287 // cout << diffrel(gam_test(), gam_euler()) << endl ;
288
289 // }
290
291 //##################
292
293
294 //### Test
295
296 if (!irrotational) {
297
298 // Tenseur diff = gam - 1 ;
299 // cout << "Etoile_bin::hydro_euler: rigid motion: difference gam <-> 1 : "
300 // << endl ;
301 // cout << norme(diff()) / norme(gam()) << endl ;
302 //
303 // cout << "Etoile_bin::hydro_euler: rigid motion: difference wit_w <-> 0 : "
304 // << endl ;
305 // cout << "Component x : " << endl << norme(wit_w(0)) << endl ;
306 // cout << "Component y : " << endl << norme(wit_w(1)) << endl ;
307 // cout << "Component z : " << endl << norme(wit_w(2)) << endl ;
308
309//####
310 gam = 1 ;
311 loggam = 0 ;
312 wit_w = 0 ;
313 }
314
315 // The derived quantities are obsolete
316 // -----------------------------------
317
318 del_deriv() ;
319
320
321}
322}
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
bool irrotational
true for an irrotational star, false for a corotating one
Definition etoile.h:822
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
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
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 log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Lorene prototypes.
Definition app_hor.h:64