LORENE
altBH_QI.C
1/*
2 * Methods of the class AltBH_QI
3 *
4 * (see file compobj.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2013 Odele Straub, 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 altBH_QI_C[] = "$Header: /cvsroot/Lorene/C++/Source/Compobj/altBH_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $" ;
29
30/*
31 * $Id: altBH_QI.C,v 1.5 2014/10/13 08:52:49 j_novak Exp $
32 * $Log: altBH_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 2014/01/17 07:33:13 o_straub
37 * adjustment to read in files with 4 lines in the header
38 *
39 * Revision 1.3 2014/01/14 16:36:11 e_gourgoulhon
40 * Corrected initialization of bbb
41 *
42 * Revision 1.2 2013/04/17 13:02:47 e_gourgoulhon
43 * Added member krphi + method extrinsic_curvature
44 *
45 * Revision 1.1 2013/04/16 15:27:27 e_gourgoulhon
46 * New class AltBH_QI
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Compobj/altBH_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#include "nbr_spx.h"
61
62 //--------------//
63 // Constructors //
64 //--------------//
65
66// Standard constructor
67// --------------------
68namespace Lorene {
69AltBH_QI::AltBH_QI(Map& mpi, const char* file_name, double a_spin_i) :
70 Compobj_QI(mpi) ,
71 a_spin(a_spin_i),
72 krphi(mpi)
73{
74
75 ifstream file(file_name) ;
76 if ( !file.good() ) {
77 cerr << "Problem in opening the file " << file_name << endl ;
78 abort() ;
79 }
80
81 file.getline(description1, 256) ;
82 file.getline(description2, 256) ;
83 cout << "description1 : " << description1 << endl ;
84 cout << "description2 : " << description2 << endl ;
85// description1[0] = " " ;
86// description2[0] = " " ;
87// cout << "description1 : " << description1 << endl ;
88// cout << "description2 : " << description2 << endl ;
89
90 const Mg3d* mg = mp.get_mg() ;
91 double tmp ;
92 file.ignore(1000,'\n') ;
93 file.ignore(1000,'\n') ;
94 int nz = mg->get_nzone() ;
95 cout << "nz : " << nz << endl ;
96 a_car.set_etat_qcq() ; // Memory allocation for A^2
97 nn.set_etat_qcq() ; // Memory allocation for N
98 nphi.allocate_all() ; // Memory allocation for N^phi
99 krphi.allocate_all() ; // Memory allocation for K_rphi
100 double r_inner = 0 ;
101 for (int l=1; l<nz; l++) {
102 cout << "l = " << l << endl ;
103 int nr = mg->get_nr(l) ;
104 int nt = mg->get_nt(l) ;
105 int np = mg->get_np(l) ;
106 double* r_iso = new double[nr] ;
107 double* r_areal = new double[nr] ;
108 double* psi4 = new double[nr] ;
109 double* alpha = new double[nr] ;
110 double* Krphi = new double[nr] ;
111 double* beta_phi = new double[nr] ;
112 int nr_max = nr ;
113 if (l==nz-1) {
114 nr_max = nr - 1 ;
115 r_areal[nr-1] = __infinity ;
116 r_iso[nr-1] = __infinity ;
117 psi4[nr-1] = 1 ;
118 alpha[nr-1] = 1 ;
119 Krphi[nr-1] = 0 ;
120 beta_phi[nr-1] = 0 ;
121 }
122 for (int i=0; i<nr_max; i++) {
123 file >> r_areal[i] ;
124 file >> r_iso[i] ;
125 file >> tmp ;
126 file >> psi4[i] ;
127 file >> alpha[i] ;
128 file >> tmp ;
129 file >> Krphi[i] ;
130 file >> beta_phi[i] ;
131 cout << "r_iso, psi4, beta_phi : " << r_iso[i] << " " << psi4[i] << " " << beta_phi[i] << endl ;
132 }
133
134 if (l==1) r_inner = r_iso[0] ;
135
136 for (int k=0; k<np; k++) {
137 for (int j=0; j<nt; j++) {
138 for (int i=0; i<nr; i++) {
139 a_car.set_grid_point(l,k,j,i) = psi4[i] ;
140 nn.set_grid_point(l,k,j,i) = alpha[i] ;
141 nphi.set_grid_point(l,k,j,i) = - a_spin * beta_phi[i] / psi4[i] / (r_iso[i]*r_iso[i]) ;
142 krphi.set_grid_point(l,k,j,i) = a_spin * Krphi[i] / r_iso[i] ;
143 }
144 }
145 }
146 }
147 file.close() ;
148
149 mp.homothetie(r_inner / mp.val_r(1,-1.,0.,0.)) ;
150
151 // Set the shift to zero in domain l=0:
152 nphi.annule(0,0) ;
153
158
159 b_car = a_car ; // slow rotation limit: B^2 = A^2
160 bbb = sqrt(b_car) ;
162
163 // Pointers of derived quantities initialized to zero :
164 set_der_0x0() ;
165}
166
167// Copy constructor
168// --------------------
170 Compobj_QI(other),
171 krphi(other.krphi)
172{
173 // Pointers of derived quantities initialized to zero :
174 set_der_0x0() ;
175}
176
177
178// Constructor from a file
179// -----------------------
180AltBH_QI::AltBH_QI(Map& mpi, FILE* ) :
181 Compobj_QI(mpi),
182 krphi(mpi)
183{
184 // Pointers of derived quantities initialized to zero :
185 set_der_0x0() ;
186
187 // Read of the saved fields:
188 // ------------------------
189
190}
191
192 //------------//
193 // Destructor //
194 //------------//
195
197
198 del_deriv() ;
199
200}
201
202
203 //----------------------------------//
204 // Management of derived quantities //
205 //----------------------------------//
206
208
210
211
213}
214
215
217
218}
219
220 //--------------//
221 // Assignment //
222 //--------------//
223
224// Assignment to another AltBH_QI
225// --------------------------------
226void AltBH_QI::operator=(const AltBH_QI& other) {
227
228 // Assignment of quantities common to all the derived classes of Compobj_QI
229 Compobj_QI::operator=(other) ;
230
231 del_deriv() ; // Deletes all derived quantities
232}
233
234 //--------------//
235 // Outputs //
236 //--------------//
237
238// Save in a file
239// --------------
240void AltBH_QI::sauve(FILE* ) const {
241
242
243}
244
245// Printing
246// --------
247
248ostream& AltBH_QI::operator>>(ostream& ost) const {
249
250 using namespace Unites ;
251
253
254 ost << endl << "Alternative black hole spacetime in quasi-isotropic coordinates (class AltBH_QI) " << endl ;
255 ost << description1 << endl ;
256 ost << description2 << endl ;
257
258 return ost ;
259
260}
261
262 //-------------------------//
263 // Computational methods //
264 //-------------------------//
265
266// Updates the extrinsic curvature
267// -------------------------------
268
270
271 // Special treatment for axisymmetric case:
272
273 if ( (mp.get_mg())->get_np(0) == 1) {
274
275 // What follows is valid only for a mapping of class Map_radial :
276 assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
277
278 Scalar tmp = krphi ;
279 tmp.mult_sint() ; // multiplication by sin(theta)
280 kk.set(1,3) = tmp ;
281
282 kk.set(2,3) = 0 ;
283
284 kk.set(1,1) = 0 ;
285 kk.set(1,2) = 0 ;
286 kk.set(2,2) = 0 ;
287 kk.set(3,3) = 0 ;
288 }
289 else {
290
291 // General case:
292
294 }
295
296 // Computation of A^2 K_{ij} K^{ij}
297 // --------------------------------
298
299 ak_car = 2 * ( kk(1,3)*kk(1,3) + kk(2,3)*kk(2,3) ) / b_car ;
300
301 del_deriv() ;
302
303}
304}
Alternative black hole spacetime in Quasi-Isotropic coordinates (under development).
Definition compobj.h:917
char description2[256]
String describing the model.
Definition compobj.h:924
void operator=(const AltBH_QI &)
Assignment to another AltBH_QI.
Definition altBH_QI.C:226
AltBH_QI(Map &mp_i, const char *file_name, double a_spin_i)
Standard constructor.
Definition altBH_QI.C:69
virtual ~AltBH_QI()
Destructor.
Definition altBH_QI.C:196
char description1[256]
String describing the model.
Definition compobj.h:923
virtual void del_deriv() const
Deletes all the derived quantities.
Definition altBH_QI.C:207
Scalar krphi
K_{(r)(phi)} read in the file.
Definition compobj.h:927
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition altBH_QI.C:248
double a_spin
Spin parameter of the model.
Definition compobj.h:925
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition altBH_QI.C:216
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition altBH_QI.C:269
virtual void sauve(FILE *) const
Save in a file.
Definition altBH_QI.C:240
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
Scalar ak_car
Scalar .
Definition compobj.h:315
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
Sym_tensor kk
Extrinsic curvature tensor
Definition compobj.h:153
Scalar nn
Lapse function N .
Definition compobj.h:135
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition compobj.C:290
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition compobj.h:132
Base class for pure radial mappings.
Definition map.h:1536
Base class for coordinate mappings.
Definition map.h:670
virtual void homothetie(double lambda)=0
Sets a new radial scale.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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 annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
void mult_sint()
Multiplication by .
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.