LORENE
star_QI.C
1/*
2 * Methods of the class Star_QI
3 *
4 * (see file compobj.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2012 Claire Some, Eric Gourgoulhon
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_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $" ;
29
30/*
31 * $Id: star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
32 * $Log: star_QI.C,v $
33 * Revision 1.5 2014/10/13 08:52:49 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2012/12/03 15:27:30 c_some
37 * Small changes
38 *
39 * Revision 1.3 2012/11/22 16:04:51 c_some
40 * Minor modifications
41 *
42 * Revision 1.2 2012/11/21 14:55:27 c_some
43 * Added methods fait_shift and fait_nphi
44 *
45 * Revision 1.1 2012/11/20 16:29:50 c_some
46 * New class Star_QI
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
50 *
51 */
52
53
54// C headers
55#include <cassert>
56
57// Lorene headers
58#include "compobj.h"
59#include "unites.h"
60
61 //--------------//
62 // Constructors //
63 //--------------//
64
65// Standard constructor
66// --------------------
67namespace Lorene {
69 Compobj_QI(mpi) ,
70 logn(mpi),
71 tnphi(mpi),
72 nuf(mpi),
73 nuq(mpi),
74 dzeta(mpi),
75 tggg(mpi),
76 w_shift(mpi, CON, mp.get_bvect_cart()),
77 khi_shift(mpi),
78 ssjm1_nuf(mpi),
79 ssjm1_nuq(mpi),
80 ssjm1_dzeta(mpi),
81 ssjm1_tggg(mpi),
82 ssjm1_khi(mpi),
83 ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
84{
85 // Pointers of derived quantities initialized to zero :
86 set_der_0x0() ;
87
88 // Initialization to a flat metric :
89 logn = 0 ;
90 tnphi = 0 ;
91 nuf = 0 ;
92 nuq = 0 ;
93 dzeta = 0 ;
94 tggg = 0 ;
95
97 khi_shift = 0 ;
98
100
101 ssjm1_nuf = 0 ;
102 ssjm1_nuq = 0 ;
103 ssjm1_dzeta = 0 ;
104 ssjm1_tggg = 0 ;
105 ssjm1_khi = 0 ;
106
108
109}
110
111// Copy constructor
112// --------------------
114 Compobj_QI(st),
115 logn(st.logn),
116 tnphi(st.tnphi),
117 nuf(st.nuf),
118 nuq(st.nuq),
119 dzeta(st.dzeta),
120 tggg(st.tggg),
121 w_shift(st.w_shift),
122 khi_shift(st.khi_shift),
123 ssjm1_nuf(st.ssjm1_nuf),
124 ssjm1_nuq(st.ssjm1_nuq),
125 ssjm1_dzeta(st.ssjm1_dzeta),
126 ssjm1_tggg(st.ssjm1_tggg),
127 ssjm1_khi(st.ssjm1_khi),
128 ssjm1_wshift(st.ssjm1_wshift)
129{
130 // Pointers of derived quantities initialized to zero :
131 set_der_0x0() ;
132}
133
134
135// Constructor from a file
136// -----------------------
137Star_QI::Star_QI(Map& mpi, FILE* fich) :
138 Compobj_QI(mpi) ,
139 logn(mpi),
140 tnphi(mpi),
141 nuf(mpi),
142 nuq(mpi),
143 dzeta(mpi),
144 tggg(mpi),
145 w_shift(mpi, CON, mp.get_bvect_cart()),
146 khi_shift(mpi),
147 ssjm1_nuf(mpi),
148 ssjm1_nuq(mpi),
149 ssjm1_dzeta(mpi),
150 ssjm1_tggg(mpi),
151 ssjm1_khi(mpi),
152 ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
153{
154 // Pointers of derived quantities initialized to zero :
155 set_der_0x0() ;
156
157 // Read of the saved fields:
158 // ------------------------
159
160 Scalar nuf_file(mp, *(mp.get_mg()), fich) ;
161 nuf = nuf_file ;
162
163 Scalar nuq_file(mp, *(mp.get_mg()), fich) ;
164 nuq = nuq_file ;
165
166 logn = nuf + nuq ; //## to be checked !
167
168 Scalar dzeta_file(mp, *(mp.get_mg()), fich) ;
169 dzeta = dzeta_file ;
170
171 Scalar tggg_file(mp, *(mp.get_mg()), fich) ;
172 tggg = tggg_file ;
173
174 Vector w_shift_file(mp, mp.get_bvect_cart(), fich) ;
175 w_shift = w_shift_file ;
176
177 Scalar khi_shift_file(mp, *(mp.get_mg()), fich) ;
178 khi_shift = khi_shift_file ;
179
180 fait_shift() ; // constructs shift from w_shift and khi_shift
181 fait_nphi() ; // constructs N^phi from (N^x,N^y,N^z)
182
183 Scalar ssjm1_nuf_file(mp, *(mp.get_mg()), fich) ;
184 ssjm1_nuf = ssjm1_nuf_file ;
185
186 Scalar ssjm1_nuq_file(mp, *(mp.get_mg()), fich) ;
187 ssjm1_nuq = ssjm1_nuq_file ;
188
189 Scalar ssjm1_dzeta_file(mp, *(mp.get_mg()), fich) ;
190 ssjm1_dzeta = ssjm1_dzeta_file ;
191
192 Scalar ssjm1_tggg_file(mp, *(mp.get_mg()), fich) ;
193 ssjm1_tggg = ssjm1_tggg_file ;
194
195 Scalar ssjm1_khi_file(mp, *(mp.get_mg()), fich) ;
196 ssjm1_khi = ssjm1_khi_file ;
197
198 Vector ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ;
199 ssjm1_wshift = ssjm1_wshift_file ;
200
201
202 // Initialization of N, A^2, B^2, gamma_ij, tkij and ak_car
203 update_metric() ;
204
205}
206
207 //------------//
208 // Destructor //
209 //------------//
210
212
213 del_deriv() ;
214
215}
216
217
218 //----------------------------------//
219 // Management of derived quantities //
220 //----------------------------------//
221
222void Star_QI::del_deriv() const {
223
225
226 if (p_grv2 != 0x0) delete p_grv2 ;
227 if (p_grv3 != 0x0) delete p_grv3 ;
228 if (p_mom_quad != 0x0) delete p_mom_quad ;
229 if (p_mass_g != 0x0) delete p_mass_g ;
230
232}
233
234
236
237 p_grv2 = 0x0 ;
238 p_grv3 = 0x0 ;
239 p_mom_quad = 0x0 ;
240 p_mass_g = 0x0 ;
241
242}
243
244 //--------------//
245 // Assignment //
246 //--------------//
247
248// Assignment to another Star_QI
249// --------------------------------
251
252 // Assignment of quantities common to all the derived classes of Compobj_QI
254
255 logn = st.logn ;
256 tnphi = st.tnphi ;
257 nuf = st.nuf ;
258 nuq = st.nuq ;
259 dzeta = st.dzeta ;
260 tggg = st.tggg ;
261 w_shift = st.w_shift ;
262 khi_shift = st.khi_shift ;
263 ssjm1_nuf = st.ssjm1_nuf ;
264 ssjm1_nuq = st.ssjm1_nuq ;
267 ssjm1_khi = st.ssjm1_khi ;
269
270 del_deriv() ; // Deletes all derived quantities
271}
272
273 //--------------//
274 // Outputs //
275 //--------------//
276
277// Save in a file
278// --------------
279void Star_QI::sauve(FILE* fich) const {
280
281 nuf.sauve(fich) ;
282 nuq.sauve(fich) ;
283 dzeta.sauve(fich) ;
284 tggg.sauve(fich) ;
285 w_shift.sauve(fich) ;
286 khi_shift.sauve(fich) ;
287
288 ssjm1_nuf.sauve(fich) ;
289 ssjm1_nuq.sauve(fich) ;
290 ssjm1_dzeta.sauve(fich) ;
291 ssjm1_tggg.sauve(fich) ;
292 ssjm1_khi.sauve(fich) ;
293 ssjm1_wshift.sauve(fich) ;
294
295}
296
297// Printing
298// --------
299
300ostream& Star_QI::operator>>(ostream& ost) const {
301
302 using namespace Unites ;
303
305
306 ost << endl << "Axisymmetric stationary compact star in quasi-isotropic coordinates (class Star_QI) " << endl ;
307
308 ost << "Central values of various fields : " << endl ;
309 ost << "-------------------------------- " << endl ;
310 ost << " ln(N) : " << logn.val_grid_point(0,0,0,0) << endl ;
311 ost << " nuf : " << nuf.val_grid_point(0,0,0,0) << endl ;
312 ost << " nuq : " << nuq.val_grid_point(0,0,0,0) << endl ;
313 ost << " zeta = ln(AN): " << dzeta.val_grid_point(0,0,0,0) << endl << endl ;
314
315 ost << "Error on the virial identity GRV2 : " << grv2() << endl ;
316 ost << "Error on the virial identity GRV3 : " << grv3(&ost) << endl ;
317
318 double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
319 / double(1.e38) ) ;
320 ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
321 << endl ;
322 ost << "c^4 Q / (G^2 M^3) : "
323 << mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) )
324 << endl ;
325
326 ost << "Total angular momentum J : "
327 << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
328 << endl ;
329 ost << "c J / (G M^2) : "
330 << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ;
331
332 return ost ;
333
334}
335
336 //-------------------------//
337 // Computational methods //
338 //-------------------------//
339
341
343
344 d_khi.dec_dzpuis(2) ; // divide by r^2 in the external compactified domain
345
346 // x_k dW^k/dx_i
347 Scalar xx(mp) ;
348 Scalar yy(mp) ;
349 Scalar zz(mp) ;
350 Scalar sintcosp(mp) ;
351 Scalar sintsinp(mp) ;
352 Scalar cost(mp) ;
353 xx = mp.x ;
354 yy = mp.y ;
355 zz = mp.z ;
356 sintcosp = mp.sint * mp.cosp ;
357 sintsinp = mp.sint * mp.sinp ;
358 cost = mp.cost ;
359
360 int nz = mp.get_mg()->get_nzone() ;
361 Vector xk(mp, COV, mp.get_bvect_cart()) ;
362 xk.set(1) = xx ;
363 xk.set(1).set_domain(nz-1) = sintcosp.domain(nz-1) ;
364 xk.set(1).set_dzpuis(-1) ;
365 xk.set(2) = yy ;
366 xk.set(2).set_domain(nz-1) = sintsinp.domain(nz-1) ;
367 xk.set(2).set_dzpuis(-1) ;
368 xk.set(3) = zz ;
369 xk.set(3).set_domain(nz-1) = cost.domain(nz-1) ;
370 xk.set(3).set_dzpuis(-1) ;
371 xk.std_spectral_base() ;
372
374
375 Vector x_d_w = contract(xk, 0, d_w, 0) ;
376 x_d_w.dec_dzpuis() ;
377
378 double lambda = double(1) / double(3) ;
379
380 beta = - (lambda+2)/2./(lambda+1) * w_shift
381 + (lambda/2./(lambda+1)) * (d_khi + x_d_w) ;
382
383}
384
385
387
388 if ( (beta(1).get_etat() == ETATZERO) && (beta(2).get_etat() == ETATZERO) ) {
389 tnphi = 0 ;
390 nphi = 0 ;
391 return ;
392 }
393
394 // Computation of tnphi
395 // --------------------
396
397 mp.comp_p_from_cartesian( -beta(1), -beta(2), tnphi ) ;
398
399 // Computation of nphi
400 // -------------------
401
402 nphi = tnphi ;
403 nphi.div_rsint() ;
404
405}
406
407
408// Updates the 3-metric and the shift
409
411
412 // Lapse function N
413 // ----------------
414
415 nn = exp( logn ) ;
416
417 nn.std_spectral_base() ; // set the bases for spectral expansions
418
419
420 // Metric factor A^2
421 // -----------------
422
423 a_car = exp( 2*( dzeta - logn ) ) ;
424
425 a_car.std_spectral_base() ; // set the bases for spectral expansions
426
427 // Metric factor B
428 // ---------------
429
430 Scalar tmp = tggg ;
431 tmp.div_rsint() ; //... Division of tG by r sin(theta)
432
433 bbb = (1 + tmp) / nn ;
434
435 bbb.std_spectral_base() ; // set the bases for spectral expansions
436
437 b_car = bbb * bbb ;
438
439
440 Compobj_QI::update_metric() ; // updates gamma_{ij}
441
442
443 // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
444 // -------------------------------------------------
445
446 //## extrinsic_curvature() ; // should be done by Compobj_QI::update_metric()
447
448
449 // The derived quantities are no longer up to date :
450 // -----------------------------------------------
451
452 del_deriv() ;
453
454}
455
456
457
458}
Base class for axisymmetric stationary compact objects in Quasi-Isotropic coordinates (under developm...
Definition compobj.h:280
void operator=(const Compobj_QI &)
Assignment to another Compobj_QI.
Definition compobj_QI.C:195
Scalar nphi
Metric coefficient .
Definition compobj.h:296
virtual void del_deriv() const
Deletes all the derived quantities.
Definition compobj_QI.C:163
virtual void update_metric()
Updates the 3-metric from A and B and the shift vector from .
Definition compobj_QI.C:302
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition compobj_QI.C:264
Scalar b_car
Square of the metric factor B.
Definition compobj.h:293
Scalar bbb
Metric factor B.
Definition compobj.h:290
Scalar a_car
Square of the metric factor A.
Definition compobj.h:287
Scalar nn
Lapse function N .
Definition compobj.h:135
Vector beta
Shift vector .
Definition compobj.h:138
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition compobj.h:132
Base class for coordinate mappings.
Definition map.h:670
Coord cosp
Definition map.h:724
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord sint
Definition map.h:721
virtual void comp_p_from_cartesian(const Scalar &v_x, const Scalar &v_y, Scalar &v_p) const =0
Computes the Spherical component (with respect to bvect_spher ) of a vector given by its cartesian c...
Coord sinp
Definition map.h:723
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition map.C:331
Coord z
z coordinate centered on the grid
Definition map.h:728
Coord x
x coordinate centered on the grid
Definition map.h:726
Coord cost
Definition map.h:722
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 field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
void div_rsint()
Division by everywhere; dzpuis is not changed.
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition scalar.h:625
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Base class for axisymmetric stationary compact stars in Quasi-Isotropic coordinates (under developmen...
Definition compobj.h:487
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata's prescription [Prog.
Definition star_QI.C:340
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition compobj.h:529
virtual double grv2() const
Error on the virial identity GRV2.
double * p_grv2
Error on the virial identity GRV2.
Definition compobj.h:585
Scalar logn
Logarithm of the lapse N .
Definition compobj.h:495
double * p_mass_g
Gravitational mass (ADM mass as a volume integral)
Definition compobj.h:588
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition compobj.h:510
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition compobj.h:569
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Star_QI(Map &mp_i)
Standard constructor.
Definition star_QI.C:68
virtual void sauve(FILE *) const
Save in a file.
Definition star_QI.C:279
virtual double mom_quad() const
Quadrupole moment.
virtual ~Star_QI()
Destructor.
Definition star_QI.C:211
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Definition compobj.h:505
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition compobj.h:551
void update_metric()
Computes metric coefficients from known potentials.
Definition star_QI.C:410
double * p_grv3
Error on the virial identity GRV3.
Definition compobj.h:586
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition star_QI.C:300
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition compobj.h:578
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition star_QI.C:386
virtual double angu_mom() const
Angular momentum.
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition compobj.h:539
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition compobj.h:545
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star_QI.C:222
virtual double mass_g() const
Gravitational mass.
void operator=(const Star_QI &)
Assignment to another Star_QI.
Definition star_QI.C:250
Scalar tggg
Metric potential .
Definition compobj.h:516
Scalar tnphi
Component of the shift vector.
Definition compobj.h:500
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition compobj.h:561
Scalar ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta .
Definition compobj.h:556
double * p_mom_quad
Quadrupole moment
Definition compobj.h:587
Scalar dzeta
Metric potential .
Definition compobj.h:513
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition star_QI.C:235
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.