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 
28 char 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 // --------------------
68 namespace Lorene {
69 AltBH_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 // -----------------------
180 AltBH_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 
207 void AltBH_QI::del_deriv() const {
208 
210 
211 
213 }
214 
215 
216 void AltBH_QI::set_der_0x0() const {
217 
218 }
219 
220  //--------------//
221  // Assignment //
222  //--------------//
223 
224 // Assignment to another AltBH_QI
225 // --------------------------------
226 void 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 // --------------
240 void AltBH_QI::sauve(FILE* ) const {
241 
242 
243 }
244 
245 // Printing
246 // --------
247 
248 ostream& AltBH_QI::operator>>(ostream& ost) const {
249 
250  using namespace Unites ;
251 
252  Compobj_QI::operator>>(ost) ;
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 }
Base class for axisymmetric stationary compact objects in Quasi-Isotropic coordinates (under developm...
Definition: compobj.h:280
AltBH_QI(Map &mp_i, const char *file_name, double a_spin_i)
Standard constructor.
Definition: altBH_QI.C:69
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition: altBH_QI.C:269
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: altBH_QI.C:207
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
double a_spin
Spin parameter of the model.
Definition: compobj.h:925
Scalar ak_car
Scalar .
Definition: compobj.h:315
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
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
Base class for coordinate mappings.
Definition: map.h:670
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: altBH_QI.C:248
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 allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:365
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
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
Scalar nphi
Metric coefficient .
Definition: compobj.h:296
Scalar b_car
Square of the metric factor B.
Definition: compobj.h:293
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.
virtual ~AltBH_QI()
Destructor.
Definition: altBH_QI.C:196
Scalar bbb
Metric factor B.
Definition: compobj.h:290
virtual void sauve(FILE *) const
Save in a file.
Definition: altBH_QI.C:240
virtual void homothetie(double lambda)=0
Sets a new radial scale.
void operator=(const AltBH_QI &)
Assignment to another AltBH_QI.
Definition: altBH_QI.C:226
void operator=(const Compobj_QI &)
Assignment to another Compobj_QI.
Definition: compobj_QI.C:195
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: altBH_QI.C:216
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
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
Alternative black hole spacetime in Quasi-Isotropic coordinates (under development).
Definition: compobj.h:917
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: compobj_QI.C:264
Scalar krphi
K_{(r)(phi)} read in the file.
Definition: compobj.h:927
Scalar nn
Lapse function N .
Definition: compobj.h:135
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
char description1[256]
String describing the model.
Definition: compobj.h:923
char description2[256]
String describing the model.
Definition: compobj.h:924
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: compobj_QI.C:163
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:132