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 
29 char 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 
105 namespace 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 
138  d_psi.set_std_base() ;
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) ;
158  u_euler.set_etat_qcq() ;
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 
165  u_euler.set_std_base() ;
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 
184  ener_euler = gam_euler % gam_euler % ( ener + press ) - press ;
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 ) *
198  sprod(u_euler, u_euler) ;
199  s_euler.set_std_base() ;
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...
Definition: et_bin_hydro.C:106
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
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 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