LORENE
compobj_QI.C
1 /*
2  * Methods of the class Compobj_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 
28 char compobj_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $" ;
29 
30 /*
31  * $Id: compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
32  * $Log: compobj_QI.C,v $
33  * Revision 1.9 2014/10/13 08:52:49 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.8 2014/05/16 11:55:18 o_straub
37  * fixed: GYOTO output from compobj & compobj_QI
38  *
39  * Revision 1.7 2013/07/25 19:44:11 o_straub
40  * calculation of the marginally bound radius
41  *
42  * Revision 1.6 2013/04/04 15:32:32 e_gourgoulhon
43  * Better computation of the ISCO
44  *
45  * Revision 1.5 2013/04/04 08:53:47 e_gourgoulhon
46  * Minor improvements
47  *
48  * Revision 1.4 2013/04/03 12:10:13 e_gourgoulhon
49  * Added member kk to Compobj; suppressed tkij
50  *
51  * Revision 1.3 2012/11/22 16:04:51 c_some
52  * Minor modifications
53  *
54  * Revision 1.2 2012/11/20 16:28:48 c_some
55  * -- tkij is created on the Cartesian triad.
56  * -- implemented method extrinsic_curvature()
57  *
58  * Revision 1.1 2012/11/16 16:14:11 c_some
59  * New class Compobj_QI
60  *
61  *
62  * $Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI.C,v 1.9 2014/10/13 08:52:49 j_novak Exp $
63  *
64  */
65 
66 
67 // C headers
68 #include <cassert>
69 #include <cmath>
70 #include <cstdio>
71 
72 // Lorene headers
73 #include "compobj.h"
74 #include "nbr_spx.h"
75 #include "utilitaires.h"
76 
77 
78 
79  //--------------//
80  // Constructors //
81  //--------------//
82 
83 // Standard constructor
84 // --------------------
85 namespace Lorene {
87  Compobj(map_i) ,
88  a_car(map_i) ,
89  bbb(map_i) ,
90  b_car(map_i) ,
91  nphi(map_i) ,
92  ak_car(map_i)
93 {
94  // Pointers of derived quantities initialized to zero :
95  set_der_0x0() ;
96 
97  // Initialization to a flat metric :
98  a_car = 1 ;
100  bbb = 1 ;
102  b_car = bbb*bbb ;
103  nphi = 0 ;
104  ak_car = 0 ;
105 
106 }
107 
108 // Copy constructor
109 // --------------------
111  Compobj(co),
112  a_car(co.a_car) ,
113  bbb(co.bbb) ,
114  b_car(co.b_car) ,
115  nphi(co.nphi) ,
116  ak_car(co.ak_car)
117 {
118  // Pointers of derived quantities initialized to zero :
119  set_der_0x0() ;
120 }
121 
122 
123 // Constructor from a file
124 // -----------------------
125 Compobj_QI::Compobj_QI(Map& map_i, FILE* fich) :
126  Compobj(map_i) ,
127  a_car(map_i, *(map_i.get_mg()), fich) ,
128  bbb(map_i, *(map_i.get_mg()), fich) ,
129  b_car(map_i) ,
130  nphi(map_i, *(map_i.get_mg()), fich) ,
131  ak_car(map_i)
132 {
133  // Pointers of derived quantities initialized to zero :
134  set_der_0x0() ;
135 
136  Scalar nn_file(mp, *(mp.get_mg()), fich) ;
137  nn = nn_file ;
138 
139  b_car = bbb*bbb ;
140 
141  // Initialization of gamma_ij:
142  update_metric() ;
143 
144  // Computation of K_ij and ak_car:
146 }
147 
148  //------------//
149  // Destructor //
150  //------------//
151 
153 
154  del_deriv() ;
155 
156 }
157 
158 
159  //----------------------------------//
160  // Management of derived quantities //
161  //----------------------------------//
162 
163 void Compobj_QI::del_deriv() const {
164 
165  Compobj::del_deriv() ;
166 
167  if (p_angu_mom != 0x0) delete p_angu_mom ;
168  if (p_r_mb != 0x0) delete p_r_mb ;
169  if (p_r_isco != 0x0) delete p_r_isco ;
170  if (p_f_isco != 0x0) delete p_f_isco ;
171  if (p_lspec_isco != 0x0) delete p_lspec_isco ;
172  if (p_espec_isco != 0x0) delete p_espec_isco ;
173 
175 }
176 
177 
179 
180  p_angu_mom = 0x0 ;
181  p_r_mb = 0x0 ;
182  p_r_isco = 0x0 ;
183  p_f_isco = 0x0 ;
184  p_lspec_isco = 0x0 ;
185  p_espec_isco = 0x0 ;
186 
187 }
188 
189  //--------------//
190  // Assignment //
191  //--------------//
192 
193 // Assignment to another Compobj_QI
194 // --------------------------------
196 
197  // Assignment of quantities common to all the derived classes of Compobj
198  Compobj::operator=(co) ;
199 
200  a_car = co.a_car ;
201  bbb = co.bbb ;
202  b_car = co.b_car ;
203  nphi = co.nphi ;
204  ak_car = co.ak_car ;
205 
206  del_deriv() ; // Deletes all derived quantities
207 }
208 
209  //--------------//
210  // Outputs //
211  //--------------//
212 
213 
214 // Save in a file
215 // --------------
216 void Compobj_QI::sauve(FILE* fich) const {
217 
218  a_car.sauve(fich) ;
219  bbb.sauve(fich) ;
220  nphi.sauve(fich) ;
221 
222  nn.sauve(fich) ;
223 
224 }
225 
226 
227 // Save in a file for GYOTO input
228 // -------------------------------
229 
230 // Redefinition of Compobj::gyoto_data
231 
232 void Compobj_QI::gyoto_data(const char* file_name) const {
233 
234  FILE* file_out = fopen(file_name, "w") ;
235  double total_time = 0. ; // for compatibility
236  double RISCO=r_isco(0) ;
237  double RMB=r_mb(0);
238 
239  fwrite_be(&total_time, sizeof(double), 1, file_out) ;
240  mp.get_mg()->sauve(file_out) ;
241  mp.sauve(file_out) ;
242  nn.sauve(file_out) ;
243  beta.sauve(file_out) ;
244  gamma.cov().sauve(file_out) ;
245  gamma.con().sauve(file_out) ;
246  kk.sauve(file_out) ;
247  fwrite_be(&RISCO, sizeof(double), 1, file_out) ;
248  fwrite_be(&RMB, sizeof(double), 1, file_out) ;
249  fclose(file_out) ;
250 
251  cout << "WRITING RISCO TO GYOTO FILE : " << RISCO << endl ;
252  cout << "WRITING RMB TO GYOTO FILE : " << RMB << endl ;
253  cout << "WRITING TO GYOTO FILE - end of part " << endl ;
254 }
255 
256 
257 
258 
259 
260 
261 // Printing
262 // --------
263 
264 ostream& Compobj_QI::operator>>(ostream& ost) const {
265 
266  Compobj::operator>>(ost) ;
267 
268  ost << endl << "Axisymmetric stationary compact object in quasi-isotropic coordinates (class Compobj_QI) " << endl ;
269 
270  ost << "Central values of various fields : " << endl ;
271  ost << "-------------------------------- " << endl ;
272  ost << " metric coefficient A^2 : " << a_car.val_grid_point(0,0,0,0) << endl ;
273  ost << " metric coefficient B^2 : " << b_car.val_grid_point(0,0,0,0) << endl ;
274  ost << " metric coefficient N^phi : " << nphi.val_grid_point(0,0,0,0) << endl ;
275  ost << " A^2 K_{ij} K^{ij} = " << ak_car.val_grid_point(0,0,0,0) << endl << endl ;
276 
277 
278  double RISCO = r_isco(0, &ost) ;
279  ost << "Coordinate r at the innermost stable circular orbit (ISCO) : " <<
280  RISCO << endl ;
281  // ost << "Circumferential radius of the innermost stable circular orbit (ISCO) : " <<
282  // bbb.val_point(risco, M_PI/2, 0)*risco << endl ;
283  // ost << "Orbital frequency at the ISCO : " << f_isco(0) << endl ;
284  // ost << "Specific energy of a particle on the ISCO : " << espec_isco(0) << endl ;
285  // ost << "Specific angular momentum of a particle on the ISCO : " << lspec_isco(0) << endl ;
286 
287 
288  double RMB = r_mb(0, &ost) ;
289  ost << "Coordinate r at the marginally bound circular orbit (R_mb) : " << RMB << endl ;
290 
291 
292  // ost << "A^2 : " << a_car << endl ;
293  // ost << "B^2 : " << b_car << endl ;
294  // ost << "nphi : " << nphi << endl ;
295 
296  return ost ;
297 
298 }
299 
300 // Updates the 3-metric and the shift
301 
303 
304  Sym_tensor gam(mp, COV, mp.get_bvect_spher()) ;
305  gam.set(1,1) = a_car ;
306  gam.set(1,2) = 0 ;
307  gam.set(1,3) = 0 ;
308  gam.set(2,2) = a_car ;
309  gam.set(2,3) = 0 ;
310  gam.set(3,3) = b_car ;
311 
312  gamma = gam ;
313 
314  assert(*(beta.get_triad()) == mp.get_bvect_spher()) ;
315 
316  beta.set(1) = 0 ;
317  beta.set(2) = 0 ;
318  Scalar nphi_ortho(nphi) ;
319  nphi_ortho.mult_rsint() ;
320  beta.set(3) = - nphi_ortho ;
321 
322  // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
323  // -------------------------------------------------
324 
326 
327 
328  // The derived quantities are no longer up to date :
329  // -----------------------------------------------
330 
331  del_deriv() ;
332 
333 }
334 
335 
336 // Updates the extrinsic curvature
337 // -------------------------------
338 
340 
341  // Special treatment for axisymmetric case:
342 
343  if ( (mp.get_mg())->get_np(0) == 1) {
344 
345  Scalar dnpdr = nphi.dsdr() ; // d/dr (N^phi)
346  Scalar dnpdt = nphi.dsdt() ; // d/dtheta (N^phi)
347 
348  // What follows is valid only for a mapping of class Map_radial :
349  assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
350 
351  dnpdr.mult_rsint() ; // multiplication by r sin(theta)
352  kk.set(1,3) = - b_car * dnpdr / (2*nn) ;
353 
354  dnpdt.mult_sint() ; // multiplication by sin(theta)
355  kk.set(2,3) = - b_car * dnpdt / (2*nn) ;
356  kk.set(2,3).inc_dzpuis(2) ; // to have the same dzpuis as kk(1,3)
357 
358  kk.set(1,1) = 0 ;
359  kk.set(1,2) = 0 ;
360  kk.set(2,2) = 0 ;
361  kk.set(3,3) = 0 ;
362  }
363  else {
364 
365  // General case:
366 
368  }
369 
370  // Computation of A^2 K_{ij} K^{ij}
371  // --------------------------------
372 
373  ak_car = 2 * ( kk(1,3)*kk(1,3) + kk(2,3)*kk(2,3) ) / b_car ;
374 
375  del_deriv() ;
376 
377 }
378 }
Base class for axisymmetric stationary compact objects in Quasi-Isotropic coordinates (***under devel...
Definition: compobj.h:280
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:372
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
virtual void sauve(FILE *) const
Save in a file.
Definition: map.C:224
Compobj_QI(Map &map_i)
Standard constructor.
Definition: compobj_QI.C:86
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: compobj_QI.C:178
Scalar ak_car
Scalar .
Definition: compobj.h:315
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene prototypes.
Definition: app_hor.h:64
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
virtual void extrinsic_curvature()
Computes the extrinsic curvature and ak_car from nphi , nn and b_car .
Definition: compobj_QI.C:339
Base class for coordinate mappings.
Definition: map.h:670
double * p_lspec_isco
Specific angular momentum of a particle at the ISCO.
Definition: compobj.h:327
void mult_sint()
Multiplication by .
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
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:906
void gyoto_data(const char *file_name) const
Save in a file for GYOTO.
Definition: compobj_QI.C:232
Sym_tensor kk
Extrinsic curvature tensor .
Definition: compobj.h:153
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition: compobj.C:290
Scalar a_car
Square of the metric factor A.
Definition: compobj.h:287
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
Scalar nphi
Metric coefficient .
Definition: compobj.h:296
Scalar b_car
Square of the metric factor B.
Definition: compobj.h:293
Base class for stationary compact objects (***under development***).
Definition: compobj.h:126
double * p_angu_mom
Angular momentum.
Definition: compobj.h:321
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:686
void operator=(const Compobj &)
Assignment to another Compobj.
Definition: compobj.C:175
virtual double r_isco(int lmin, ostream *ost=0x0) const
Coordinate r of the innermost stable circular orbit (ISCO).
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
void sauve(FILE *fd, bool save_base=false) const
Saves into a file.
Definition: mg3d.C:371
Scalar bbb
Metric factor B.
Definition: compobj.h:290
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
double * p_espec_isco
Specific energy of a particle at the ISCO.
Definition: compobj.h:325
Vector beta
Shift vector .
Definition: compobj.h:138
void operator=(const Compobj_QI &)
Assignment to another Compobj_QI.
Definition: compobj_QI.C:195
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: compobj.C:155
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
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
virtual ~Compobj_QI()
Destructor.
Definition: compobj_QI.C:152
double * p_r_isco
Coordinate r of the ISCO.
Definition: compobj.h:322
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: compobj.C:239
double * p_r_mb
Coordinate r of the marginally bound orbit.
Definition: compobj.h:328
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
double * p_f_isco
Orbital frequency of the ISCO.
Definition: compobj.h:323
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
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 nn
Lapse function N .
Definition: compobj.h:135
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
virtual double r_mb(int lmin, ostream *ost=0x0) const
Coordinate r of the marginally bound circular orbit (R_mb).
Metric gamma
3-metric
Definition: compobj.h:141
virtual void sauve(FILE *) const
Save in a file.
Definition: compobj_QI.C:216
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: compobj_QI.C:163
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:132