LORENE
star_rot_dirac.C
1 /*
2  * Methods of class Star_rot_Dirac
3  *
4  * (see file star_rot_dirac.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
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 char star_rot_dirac_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_rot_dirac.C,v 1.10 2014/10/13 08:53:39 j_novak Exp $" ;
29 
30 /*
31  * $Id: star_rot_dirac.C,v 1.10 2014/10/13 08:53:39 j_novak Exp $
32  * $Log: star_rot_dirac.C,v $
33  * Revision 1.10 2014/10/13 08:53:39 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.9 2014/10/06 15:13:17 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.8 2013/04/25 15:46:06 j_novak
40  * Added special treatment in the case np = 1, for type_p = NONSYM.
41  *
42  * Revision 1.7 2008/05/30 08:27:38 j_novak
43  * New global quantities rp_circ and ellipt (circumferential polar coordinate and
44  * ellipticity).
45  *
46  * Revision 1.6 2007/11/06 16:23:59 j_novak
47  * Added the flag spectral_filter giving the order of possible spectral filtering
48  * of the hydro sources of metric equations (some members *_euler). The filtering
49  * is done in strot_dirac_hydro, if this flag is non-zero.
50  *
51  * Revision 1.5 2007/11/06 10:15:19 j_novak
52  * Change the order of updates in the constructor from a file, to avoid
53  * inconsistencies.
54  *
55  * Revision 1.4 2007/10/30 16:55:23 j_novak
56  * Completed the input/ouput to a file
57  *
58  * Revision 1.3 2005/02/09 13:37:37 lm_lin
59  *
60  * Add pointers p_tsw, p_aplat, and p_r_circ; add more screen output
61  * information.
62  *
63  * Revision 1.2 2005/02/02 09:22:29 lm_lin
64  *
65  * Add the GRV3 error to screen output
66  *
67  * Revision 1.1 2005/01/31 08:51:48 j_novak
68  * New files for rotating stars in Dirac gauge (still under developement).
69  *
70  *
71  * $Header: /cvsroot/Lorene/C++/Source/Star/star_rot_dirac.C,v 1.10 2014/10/13 08:53:39 j_novak Exp $
72  *
73  */
74 
75 
76 // C headers
77 #include <cmath>
78 #include <cassert>
79 
80 // Lorene headers
81 #include "star_rot_dirac.h"
82 #include "unites.h"
83 #include "utilitaires.h"
84 
85 
86  //--------------//
87  // Constructors //
88  //--------------//
89 
90 // Standard constructor
91 //-------------------------
92 namespace Lorene {
93 Star_rot_Dirac::Star_rot_Dirac(Map& mpi, int nzet_i, const Eos& eos_i, int filter)
94  : Star(mpi, nzet_i, eos_i),
95  spectral_filter(filter),
96  psi4(mpi),
97  psi2(mpi),
98  qqq(mpi),
99  ln_psi(mpi),
100  j_euler(mpi, CON, mpi.get_bvect_spher()),
101  v2(mpi),
102  flat(mpi.flat_met_spher()),
103  tgamma(flat),
104  aa(mpi, CON, mpi.get_bvect_spher()),
105  taa(mpi, COV, mpi.get_bvect_spher()),
106  aa_quad(mpi),
107  hh(mpi, mpi.get_bvect_spher(), flat)
108 {
109  assert (spectral_filter>=0) ;
110  assert (spectral_filter<1000) ;
111 
112  // Initialization to a static state
113  omega = 0 ;
114  v2 = 0 ;
115 
116  // All the matter quantities are initialized to zero
118 
119  // Initialization to a flat case
120  psi4 = 1 ;
121  psi2 = 1 ;
122  qqq = 1 ;
123  ln_psi = 0 ;
124  aa.set_etat_zero() ;
125  taa.set_etat_zero() ;
126  aa_quad = 0 ;
127  hh.set_etat_zero() ;
128 
129  // Pointers of derived quantities initialized to zero :
130  set_der_0x0() ;
131 
132 }
133 
134 
135 // Copy constructor
136 //-----------------
138  : Star(star),
139  spectral_filter(star.spectral_filter),
140  psi4(star.psi4),
141  psi2(star.psi2),
142  qqq(star.qqq),
143  ln_psi(star.ln_psi),
144  j_euler(star.j_euler),
145  v2(star.v2),
146  flat(star.flat),
147  tgamma(star.tgamma),
148  aa(star.aa),
149  taa(star.taa),
150  aa_quad(star.aa_quad),
151  hh(star.hh)
152 {
153 
154  omega = star.omega ;
155 
156  // Pointers of derived quantities initialized to zero :
157  set_der_0x0() ;
158 
159 }
160 
161 
162 //Constructor from a file
163 //------------------------
164 Star_rot_Dirac::Star_rot_Dirac(Map& mpi, const Eos& eos_i, FILE* fich)
165  : Star(mpi, eos_i, fich),
166  psi4(mpi),
167  psi2(mpi),
168  qqq(mpi, *(mpi.get_mg()), fich),
169  ln_psi(mpi),
170  j_euler(mpi, CON, mpi.get_bvect_spher()),
171  v2(mpi),
172  flat(mpi.flat_met_spher()),
173  tgamma(flat),
174  aa(mpi, CON, mpi.get_bvect_spher()),
175  taa(mpi, COV, mpi.get_bvect_spher()),
176  aa_quad(mpi),
177  hh(mpi, mpi.get_bvect_spher(), flat, fich)
178 {
179 
180  // Pointers of derived quantities initialized to zero
181  //----------------------------------------------------
182  set_der_0x0() ;
183 
184  fread_be(&spectral_filter, sizeof(int), 1, fich) ;
185 
186  // Metric fields are read in the file:
187  fread_be(&omega, sizeof(double), 1, fich) ;
188  Vector shift_tmp(mpi, mpi.get_bvect_spher(), fich) ;
189  beta = shift_tmp ;
190 
191  update_metric() ;
192 
194 
195  hydro_euler() ;
196 
197 
198 
199 
200 }
201 
202 
203  //------------//
204  // Destructor //
205  //------------//
206 
208 
210 
211 }
212 
213 
214  //----------------------------------//
215  // Management of derived quantities //
216  //----------------------------------//
217 
219 
220  if (p_angu_mom != 0x0) delete p_angu_mom ;
221  if (p_grv2 != 0x0) delete p_grv2 ;
222  if (p_grv3 != 0x0) delete p_grv3 ;
223  if (p_tsw != 0x0) delete p_tsw ;
224  if (p_r_circ != 0x0) delete p_r_circ ;
225  if (p_rp_circ != 0x0) delete p_rp_circ ;
226 
227  set_der_0x0() ;
228 
229  Star::del_deriv() ;
230 
231 }
232 
233 
235 
236  p_angu_mom = 0x0 ;
237  p_grv2 = 0x0 ;
238  p_grv3 = 0x0 ;
239  p_tsw = 0x0 ;
240  p_r_circ = 0x0 ;
241  p_rp_circ = 0x0 ;
242 
243 }
244 
245 
247 
249  v2.set_etat_nondef() ;
250 
251  del_deriv() ;
252 
254 
255 }
256 
257 
258 
259  //---------------//
260  // Assignment //
261  //---------------//
262 
263 // Assignment to another Star_rot_Dirac
264 // ------------------------------------
265 
267 
268  // Assignment of quantities common to all the derived classes of Star
269  Star::operator=(star) ;
270 
271  // Assignment of proper quantities of class Star_rot_Dirac
273  omega = star.omega ;
274  psi4 = star.psi4 ;
275  psi2 = star.psi2 ;
276  qqq = star.qqq ;
277  ln_psi = star.ln_psi ;
278  j_euler = star.j_euler ;
279  v2 = star.v2 ;
280  tgamma = star.tgamma ;
281  aa = star.aa ;
282  aa_quad = star.aa_quad ;
283  hh = star.hh ;
284 
285  assert(&flat == &star.flat) ;
286 
287  del_deriv() ; // Deletes all derived quantities
288 
289 }
290 
291 
292  //-----------//
293  // Outputs //
294  //-----------//
295 
296 // Save in a file
297 // --------------
298 
299 void Star_rot_Dirac::sauve(FILE* fich) const {
300 
301  Star::sauve(fich) ;
302 
303  qqq.sauve(fich) ;
304  hh.sauve(fich) ;
305  fwrite_be(&spectral_filter, sizeof(int), 1, fich) ;
306  fwrite_be(&omega, sizeof(double), 1, fich) ;
307  beta.sauve(fich) ;
308 
309 }
310 
311 
312 // Printing
313 // ---------
314 
315 ostream& Star_rot_Dirac::operator>>(ostream& ost) const {
316 
317  using namespace Unites ;
318 
319  Star::operator>>(ost) ;
320 
321  ost << "Rotating star in Dirac gauge" << endl ;
322 
323  // Only uniformly rotating star for the moment....
324  ost << endl ;
325  ost << "Uniformly rotating star" << endl ;
326  ost << "-----------------------" << endl ;
327  if (spectral_filter > 0)
328  ost << "hydro sources of equations are filtered\n"
329  << "with " << spectral_filter << "-order exponential filter" << endl ;
330 
331  double freq = omega/ (2.*M_PI) ;
332  ost << "Omega : " << omega * f_unit
333  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
334  ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
335  << endl ;
336 
337  ost << "Error on the virial identity GRV2 : " << endl ;
338  ost << "GRV2 = " << grv2() << endl ;
339  ost << "Error on the virial identity GRV3 : " << endl ;
340  ost << "GRV3 = " << grv3() << endl ;
341 
342  ost << "Angular momentum J : "
343  << angu_mom()/( qpig / (4*M_PI) *msol*msol) << " G M_sol^2 / c"
344  << endl ;
345  ost << "c J / (G M^2) : "
346  << angu_mom()/( qpig / (4*M_PI) * pow(mass_g(), 2.) ) << endl ;
347 
348  if (omega != 0.) {
349  double mom_iner = angu_mom() / omega ;
350  double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.))
351  / double(1.e38) ) ;
352  ost << "Moment of inertia: " << mom_iner_38si << " 10^38 kg m^2"
353  << endl ;
354  }
355 
356  ost << "Ratio T/W : " << tsw() << endl ;
357  ost << "Circumferential equatorial radius R_circ : "
358  << r_circ()/km << " km" << endl ;
359  if (mp.get_mg()->get_np(0) == 1)
360  ost << "Circumferential polar radius Rp_circ : "
361  << rp_circ()/km << " km" << endl ;
362  ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
363  << endl ;
364  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
365  if (mp.get_mg()->get_np(0) == 1)
366  ost << "Ellipticity sqrt(1-(Rp_circ/R_circ)^2) : " << ellipt() << endl ;
367 
368  double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
369  ost << "Compaction parameter M_g / R_circ : " << compact << endl ;
370 
371 
372  // More to come here.....
373 
374  return ost ;
375 
376 }
377 }
Equation of state base class.
Definition: eos.h:190
Base class for coordinate mappings.
Definition: map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:686
virtual void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
Definition: scalar.C:344
Class for relativistic rotating stars in Dirac gauge and maximal slicing.
virtual double mass_g() const
Gravitational mass.
virtual double ellipt() const
Ellipticity e.
virtual double grv3() const
Error on the virial identity GRV3.
Star_rot_Dirac(Map &mp_i, int nzet_i, const Eos &eos_i, int filter=0)
Standard constructor.
Sym_tensor_trans hh
is defined by .
virtual double angu_mom() const
Angular momentum.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
virtual void del_deriv() const
Deletes all the derived quantities.
double * p_grv3
Error on the virial identity GRV3.
virtual double r_circ() const
Circumferential equatorial radius.
int spectral_filter
Spectral exponential filtering order.
double omega
Rotation angular velocity ([f_unit] )
double * p_tsw
Ratio T/W.
double * p_r_circ
Circumferential equatorial radius.
virtual double tsw() const
Ratio T/W.
Vector j_euler
Momentum density 3-vector with respect to the Eulerian observer.
const Metric_flat & flat
flat metric (spherical components)
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
virtual void sauve(FILE *) const
Save in a file.
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
virtual double rp_circ() const
Circumferential polar radius.
void update_metric()
Computes metric quantities from known potentials.
Scalar psi4
Conformal factor .
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
double * p_grv2
Error on the virial identity GRV2.
virtual ~Star_rot_Dirac()
Destructor.
double * p_angu_mom
Angular momentum.
double * p_rp_circ
Circumferential polar radius.
virtual double aplat() const
Flattening r_pole/r_eq.
void operator=(const Star_rot_Dirac &)
Assignment to another Star_rot_Dirac.
virtual double grv2() const
Error on the virial identity GRV2.
Base class for stars.
Definition: star.h:175
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition: star.C:330
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star.C:298
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: star.C:417
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
virtual void sauve(FILE *) const
Save in a file.
Definition: star.C:397
Map & mp
Mapping associated with the star.
Definition: star.h:180
void operator=(const Star &)
Assignment to another Star.
Definition: star.C:351
Vector beta
Shift vector.
Definition: star.h:228
Tensor field of valence 1.
Definition: vector.h:188
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:70
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:906
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:372
virtual void set_etat_nondef()
Sets the logical state of all components to ETATNONDEF (undefined state).
Definition: tensor.C:489
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.