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