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
28char 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//-------------------------
92namespace Lorene {
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{
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() ;
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//------------------------
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
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
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
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
300
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
316
317 using namespace Unites ;
318
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
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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
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].
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.