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