LORENE
scalarBH.C
1 /*
2  * Methods of the class ScalarBH
3  *
4  * (see file compobj.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2015 Frederic Vincent, 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: scalarBH.C,v 1.7 2016/12/05 16:17:49 j_novak Exp $
32  * $Log: scalarBH.C,v $
33  * Revision 1.7 2016/12/05 16:17:49 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.6 2016/05/10 12:52:32 f_vincent
37  * scalarBH: adding a flag to treat both boson stars and scalar BH
38  *
39  * Revision 1.5 2015/12/15 06:45:47 f_vincent
40  * Few modifs to scalarBH.C to handle spacetime with horizon
41  *
42  * Revision 1.4 2015/11/09 16:00:57 f_vincent
43  * Updated ScalarBH class
44  *
45  * Revision 1.3 2015/11/05 17:30:46 f_vincent
46  * Updated class scalarBH.
47  *
48  * Revision 1.2 2015/10/27 10:53:23 f_vincent
49  * Updated class scalarBH
50  *
51  * Revision 1.1 2015/10/22 09:18:36 f_vincent
52  * New class ScalarBH
53  *
54  *
55  * $Header: /cvsroot/Lorene/C++/Source/Compobj/scalarBH.C,v 1.7 2016/12/05 16:17:49 j_novak Exp $
56  *
57  */
58 
59 
60 // C headers
61 #include <cassert>
62 
63 // Lorene headers
64 #include "compobj.h"
65 #include "unites.h"
66 #include "nbr_spx.h"
67 
68 //--------------//
69 // Constructors //
70 //--------------//
71 
72 // Standard constructor
73 // --------------------
74 namespace Lorene {
75  ScalarBH::ScalarBH(Map& mpi, const char* file_name) :
76  Compobj(mpi),
77  ff0(mpi),
78  ff1(mpi),
79  ff2(mpi),
80  ww(mpi),
81  sfield(mpi),
82  rHor(0.)
83  {
84 
85  ifstream file(file_name) ;
86  if ( !file.good() ) {
87  cerr << "Problem in opening the file " << file_name << endl ;
88  abort() ;
89  }
90 
91  const Mg3d* mg = mp.get_mg() ;
92  double rH2 ;
93  int nrfile, nthetafile;
94  file >> nrfile >> nthetafile ;
95  file >> rHor ;
96  rH2 = rHor*rHor ;
97 
98  if (rHor<0. || nrfile<0. || nthetafile<0.){
99  cerr << "In scalarBH::scalarBH(Map,char*): "
100  << "Bad parameters rHor, nrfile, nthetafile" << endl;
101  abort();
102  }
103 
104  int isphi;
105  file >> isphi ;
106 
107  cout << "nr, ntheta from file = " << nrfile << " " << nthetafile << endl;
108  cout << "rHor from file = " << rHor << endl;
109  if (isphi==1) {
110  cout << "Scalar field values provided" << endl;
111  }else if (isphi==0){
112  cout << "Scalar field values not provided, put to zero" << endl;
113  }else{
114  cerr << "In scalarBH::scalarBH(Map,char*): "
115  << "Bad parameter isphi" << endl;
116  abort();
117  }
118 
119  double* Xfile = new double[nrfile * nthetafile] ;
120  double* thetafile = new double[nrfile * nthetafile] ;
121  double* f0file = new double[nrfile * nthetafile] ;
122  double* f1file = new double[nrfile * nthetafile] ;
123  double* f2file = new double[nrfile * nthetafile] ;
124  double* wwfile = new double[nrfile * nthetafile] ;
125  double* sfieldfile = new double[nrfile * nthetafile] ;
126 
127  cout << "Reading metric data... ";
128  for (int ii=0;ii<nrfile*nthetafile;ii++){
129  // there are empty lines in Carlos file, but it doesn't seem to be a pb
130  file >> Xfile[ii] ;
131  file >> thetafile[ii] ;
132  file >> f1file[ii] ;
133  file >> f2file[ii] ;
134  file >> f0file[ii] ;
135  if (isphi==1) {
136  file >> sfieldfile[ii] ;
137  }else{
138  sfieldfile[ii] = 0.;
139  }
140  file >> wwfile[ii] ;
141  //cout << ii << " " << Xfile[ii] << " " << thetafile[ii] << " " << f0file[ii] << " " << f1file[ii] << " " << f2file[ii] << " " << wwfile[ii] << " " << sfieldfile[ii] << endl;
142  }
143  cout << "done" << endl;
144  file.close() ;
145 
146  double Xbefmax = Xfile[nrfile-2];
147 
148  int nz = mg->get_nzone() ;
149  //cout << "nz : " << nz << endl ;
150  ff0.allocate_all() ; // Memory allocation for F_0
151  ff1.allocate_all() ; // Memory allocation for F_1
152  ff2.allocate_all() ; // Memory allocation for F_2
153  ww.allocate_all() ; // Memory allocation for W
154  sfield.allocate_all() ; // Memory allocation for scalar field phi
155 
156  double delta_theta = 0.1*fabs(thetafile[nrfile] - thetafile[0]); // small wrt the theta step in Carlos grid
157 
158  cout << "Starting interpolating Lorene grid... " ;
159  Mtbl rr(mp.r);
160  Mtbl theta(mp.tet);
161  int l_min; // this should be 0 for a spacetime without horizon, 1 with
162  if (rHor>0.){
163  l_min = 1;
164  }else{
165  l_min = 0;
166  }
167  for (int l=l_min; l<nz; l++) {
168  int nr = mg->get_nr(l) ;
169  int nt = mg->get_nt(l) ;
170  int np = mg->get_np(l) ;
171  //cout << "Starting loop k j i" << endl;
172  for (int k=0; k<np; k++){
173  for (int j=0; j<nt; j++){
174  double th0 = theta(l, k, j, 0);
175  for (int i=0; i<nr; i++){
176  double r0 = rr(l, k, j, i);
177  double x0, xx0;
178  if (r0 < __infinity){
179  x0 = sqrt(r0*r0 - rH2);
180  xx0 = x0 / (x0+1);
181  }else{
182  xx0 = 1.;
183  }
184 
185  //cout << "Loene radial stuff= " << r0 << " " << x0 << " " << xx0 << endl;
186  //cout << "Lorene points= " << xx0 << " " << th0 << endl;
187 
188  int ith=0;
189  int ithbis=0;
190  double thc = thetafile[ith];
191  while (fabs(th0 - thc) > delta_theta){
192  ith += nrfile;
193  ithbis++;
194  thc = thetafile[ith];
195  }
196  int ir=ith;
197  double xc = Xfile[ir];
198  if (xc > 0.){
199  cerr << "In scalarBH::ScalarBH(): "
200  "r should be zero here" << endl;
201  abort();
202  }
203  int doradinterp=1;
204  if (xx0 == 0.){
205  doradinterp=0;
206  }else if(xx0 == 1.){
207  ir += nrfile-1;
208  xc = Xfile[ir];
209  doradinterp=0;
210  }else if(xx0 > Xbefmax && xx0< 1.){
211  ir+=nrfile-2;
212  xc = Xfile[ir];
213  }else{
214  //cout << "Indice= " << ithbis << " " << ir << " " << xc << endl;
215 
216  while (xx0 > xc){
217  ir ++;
218  xc = Xfile[ir];
219  //cout << "xc, xx0 in loop= " << xc << " " << xx0 << endl;
220  }
221  }
222 
223  double f0interp, f1interp, f2interp, winterp, sfieldinterp;
224  if (doradinterp){
225  double xcinf = Xfile[ir-1], xcsup = Xfile[ir+1];
226  int irext1, irext2;
227  if (ir-3>0) {irext1=ir-2; irext2=ir-3;}
228  else if (ir+3<nrfile) {irext1=ir+2; irext2=ir+3;}
229  else{
230  cerr << "scalarBH::scalarBH(): bad radial indice" << endl;
231  abort();
232  }
233  double xcext1 = Xfile[irext1], xcext2 = Xfile[irext2];
234 
235  // At this stage we have either xcext2<xcext1<xcinf<xx0<xc<xcsup
236  // or xcinf<xx0<xc<xcsup<xcext1<xcext2
237 
238  //cout << "index, X= " << irext1 << " " << xcext1 <<" " << irext2 << " " << xcext2 << endl;
239  //cout << "X stuff= " << xcext2 << " " << xcext1 << " " << xcinf << " " << xx0 << " " << xc << " " << xcsup << endl;
240 
241  if (fabs(thc-th0)>delta_theta){
242  cerr << "scalarBH::ScalarBH(): theta problem in grid" << endl;
243  cerr << "Theta info: " << thc << " " << th0 << endl;
244  abort();
245  }
246  if (xx0 <= Xbefmax &&
247  (xx0 < xcinf || xx0 > xc || xx0 > xcsup)){
248  cerr << "scalarBH::ScalarBH(): rad problem in grid" << endl;
249  cerr << "Radial info: " << xcinf << " " << xx0 << " "
250  << xc << " " << xcsup << endl;
251  abort();
252  }else if (xx0 > Xbefmax &&
253  (xx0 < xcinf || xx0 < xc || xx0 > xcsup)){
254  cerr << "scalarBH::ScalarBH(): special rad "
255  "problem in grid" << endl;
256  cerr << "Radial info: " << xcinf << " " << xx0 << " "
257  << xc << " " << xcsup << endl;
258  abort();
259  }
260  //Radial polynomials
261  double polyrinf = (xx0-xc)*(xx0-xcsup)*(xx0-xcext1)*(xx0-xcext2)/((xcinf-xc)*(xcinf-xcsup)*(xcinf-xcext1)*(xcinf-xcext2));
262  double polyrmid = (xx0-xcinf)*(xx0-xcsup)*(xx0-xcext1)*(xx0-xcext2)/((xc-xcinf)*(xc-xcsup)*(xc-xcext1)*(xc-xcext2));
263  double polyrsup = (xx0-xcinf)*(xx0-xc)*(xx0-xcext1)*(xx0-xcext2)/((xcsup-xcinf)*(xcsup-xc)*(xcsup-xcext1)*(xcsup-xcext2));
264  double polyrext1 = (xx0-xcinf)*(xx0-xc)*(xx0-xcsup)*(xx0-xcext2)/((xcext1-xcinf)*(xcext1-xc)*(xcext1-xcsup)*(xcext1-xcext2));
265  double polyrext2 = (xx0-xcinf)*(xx0-xc)*(xx0-xcsup)*(xx0-xcext1)/((xcext2-xcinf)*(xcext2-xc)*(xcext2-xcsup)*(xcext2-xcext1));
266 
267  // Grid values of all Scalars
268  double f0ext1 = f0file[irext1], f0ext2 = f0file[irext2],
269  f0inf = f0file[ir-1],
270  f0mid = f0file[ir],
271  f0sup = f0file[ir+1], f1ext1=f1file[irext1],
272  f1ext2=f1file[irext2],
273  f1inf = f1file[ir-1], f1mid = f1file[ir],
274  f1sup = f1file[ir+1], f2ext1=f2file[irext1],
275  f2ext2=f2file[irext2],
276  f2inf = f2file[ir-1], f2mid = f2file[ir],
277  f2sup = f2file[ir+1], wext1=wwfile[irext1],
278  wext2=wwfile[irext2],
279  winf = wwfile[ir-1], wmid = wwfile[ir],
280  wsup = wwfile[ir+1], sfext1=sfieldfile[irext1],
281  sfext2=sfieldfile[irext2],
282  sfinf = sfieldfile[ir-1],
283  sfmid = sfieldfile[ir], sfsup = sfieldfile[ir+1];
284 
285  /*cout << "Interpolating" << endl;
286  cout << "Carlos points= " << xcinf << " " << xx0 << " " << xc << " " << xcsup << endl;
287  cout << "f0= " << f0inf << " " << f0mid << " " << f0sup << endl;*/
288 
289  // Interpolate Scalars
290  f0interp = f0ext1*polyrext1 + f0ext2*polyrext2 + f0inf*polyrinf
291  + f0mid*polyrmid + f0sup*polyrsup;
292  f1interp = f1ext1*polyrext1 + f1ext2*polyrext2 + f1inf*polyrinf
293  + f1mid*polyrmid + f1sup*polyrsup;
294  f2interp = f2ext1*polyrext1 + f2ext2*polyrext2 + f2inf*polyrinf
295  + f2mid*polyrmid + f2sup*polyrsup;
296  winterp = wext1*polyrext1 + wext2*polyrext2 + winf*polyrinf
297  + wmid*polyrmid + wsup*polyrsup;
298  sfieldinterp = sfext1*polyrext1 + sfext2*polyrext2
299  + sfinf*polyrinf
300  + sfmid*polyrmid + sfsup*polyrsup;
301  }else{
302 
303  /*cout << "Not interpolating" << endl;
304  cout << "Carlos point= " << xc << endl;
305  cout << "W= " << wwfile[ir] << endl;*/
306  // No interpolation at grid ends
307  f0interp = f0file[ir] ;
308  f1interp = f1file[ir] ;
309  f2interp = f2file[ir] ;
310  winterp = wwfile[ir] ;
311  sfieldinterp = sfieldfile[ir] ;
312  }
313 
314  ff0.set_grid_point(l,k,j,i) = f0interp ;
315  ff1.set_grid_point(l,k,j,i) = f1interp ;
316  ff2.set_grid_point(l,k,j,i) = f2interp ;
317  ww.set_grid_point(l,k,j,i) = winterp ;
318  sfield.set_grid_point(l,k,j,i) = sfieldinterp ;
319  }
320  }
321  }
322  }
323 
324  // Deleting arrays useless now
325  delete[] Xfile;
326  delete[] thetafile;
327  delete[] f0file;
328  delete[] f1file;
329  delete[] f2file;
330  delete[] wwfile;
331  delete[] sfieldfile;
332 
337  sfield.std_spectral_base() ; // to be modified: parity of |Phi|
338 
339 
340  cout << "done." << endl;
341 
342  // At this point the Scalar ff0, ff1, ff2, ww, sfield
343  // are initialized on the Lorene grid to proper interpolated values
344 
345  cout << "Starting updating metric... " ;
346  update_metric();
347  cout << "done." << endl;
348 
349  // Pointers of derived quantities initialized to zero :
350  set_der_0x0() ;
351  }
352 
353  // Copy constructor
354  // --------------------
355  ScalarBH::ScalarBH(const ScalarBH& other) :
356  Compobj(other),
357  ff0(other.ff0),
358  ff1(other.ff0),
359  ff2(other.ff0),
360  ww(other.ff0),
361  sfield(other.ff0),
362  rHor(other.rHor)
363  {
364  // Pointers of derived quantities initialized to zero :
365  set_der_0x0() ;
366  }
367 
368 
369  // Constructor from a file
370  // -----------------------
371  ScalarBH::ScalarBH(Map& mpi, FILE* ) :
372  Compobj(mpi),
373  ff0(mpi),
374  ff1(mpi),
375  ff2(mpi),
376  ww(mpi),
377  sfield(mpi),
378  rHor(0.)
379  {
380  // Pointers of derived quantities initialized to zero :
381  set_der_0x0() ;
382 
383  // Read of the saved fields:
384  // ------------------------
385 
386  }
387 
388  //------------//
389  // Destructor //
390  //------------//
391 
393 
394  del_deriv() ;
395 
396  }
397 
398 
399  //----------------------------------//
400  // Management of derived quantities //
401  //----------------------------------//
402 
403  void ScalarBH::del_deriv() const {
404 
405  Compobj::del_deriv() ;
406 
407 
409  }
410 
411 
412  void ScalarBH::set_der_0x0() const {
413 
414  }
415 
416  //--------------//
417  // Assignment //
418  //--------------//
419 
420  // Assignment to another ScalarBH
421  // --------------------------------
422  void ScalarBH::operator=(const ScalarBH& other) {
423 
424  // Assignment of quantities common to all the derived classes of Compobj
425  Compobj::operator=(other) ;
426 
427  del_deriv() ; // Deletes all derived quantities
428  }
429 
430  //--------------//
431  // Outputs //
432  //--------------//
433 
434  // Save in a file
435  // --------------
436  void ScalarBH::sauve(FILE* ) const {
437 
438 
439  }
440 
441  // Printing
442  // --------
443 
444  ostream& ScalarBH::operator>>(ostream& ost) const {
445 
446  using namespace Unites ;
447 
448  Compobj::operator>>(ost) ;
449 
450  ost << endl << "Black hole with scalar hair (class ScalarBH) " << endl ;
451  // ost << description1 << endl ;
452  // ost << description2 << endl ;
453 
454  return ost ;
455 
456  }
457 
458  //-------------------------//
459  // Computational methods //
460  //-------------------------//
461 
462  // Updates the extrinsic curvature
463  // -------------------------------
464 
465  //void ScalarBH::extrinsic_curvature() {
466 
467  // FV: commenting out October 2015 to compile
468 
469  // // Special treatment for axisymmetric case:
470 
471  // if ( (mp.get_mg())->get_np(0) == 1) {
472 
473  // // What follows is valid only for a mapping of class Map_radial :
474  // assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
475 
476  // Scalar tmp = krphi ;
477  // tmp.mult_sint() ; // multiplication by sin(theta)
478  // kk.set(1,3) = tmp ;
479 
480  // kk.set(2,3) = 0 ;
481 
482  // kk.set(1,1) = 0 ;
483  // kk.set(1,2) = 0 ;
484  // kk.set(2,2) = 0 ;
485  // kk.set(3,3) = 0 ;
486  // }
487  // else {
488 
489  // // General case:
490 
491  // Compobj::extrinsic_curvature() ;
492  // }
493 
494  // // Computation of A^2 K_{ij} K^{ij}
495  // // --------------------------------
496 
497  // ak_car = 2 * ( kk(1,3)*kk(1,3) + kk(2,3)*kk(2,3) ) / b_car ;
498 
499  // del_deriv() ;
500 
501  //}
502 
503 
504 void ScalarBH::update_metric() {
505  Mtbl rr(mp.r);
506  Scalar NN(mp);
507  NN = 1 - rHor/rr;
508  if (rHor>0.){
509  NN.set_domain(0) = 1;
510  }
511 
512  nn = exp(ff0)*sqrt(NN);
514 
515  Sym_tensor gam(mp, COV, mp.get_bvect_spher()) ;
516  // Component in an orthonormal basis, thus, no r^2, r^2sin2theta terms
517  gam.set(1,1) = exp(2*ff1)/NN ;
518  gam.set(1,1).std_spectral_base() ;
519  gam.set(1,2) = 0 ;
520  gam.set(1,3) = 0 ;
521  gam.set(2,2) = exp(2*ff1); //gam(1,1) ;
522  gam.set(2,2).std_spectral_base() ;
523  gam.set(2,3) = 0 ;
524  gam.set(3,3) = exp(2*ff2) ;
525  gam.set(3,3).std_spectral_base() ;
526 
527  gamma = gam ;
528 
529  assert(*(beta.get_triad()) == mp.get_bvect_spher()) ;
530 
531  beta.set(1) = 0 ;
532  beta.set(2) = 0 ;
533  Scalar nphi_ortho(ww) ;
534  nphi_ortho.mult_rsint() ;
535  beta.set(3) = - nphi_ortho ;
536 
537  // Tensor B^{-2} K_{ij} and Scalar A^2 K_{ij} K^{ij}
538  // -------------------------------------------------
539 
541 
542 
543  // The derived quantities are no longer up to date :
544  // -----------------------------------------------
545 
546  del_deriv() ;
547 
548 }
549 
550 
551 } // End nammespace Lorene
virtual ~ScalarBH()
Destructor.
Definition: scalarBH.C:392
Scalar ff2
Metric field F_2 of Herdeiro & Radu (2015)
Definition: compobj.h:1033
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
virtual void sauve(FILE *) const
Save in a file.
Definition: scalarBH.C:436
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Multi-domain array.
Definition: mtbl.h:118
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
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
Scalar ww
Metric field W of Herdeiro & Radu (2015)
Definition: compobj.h:1034
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 allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:371
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
Definition: compobj.C:293
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: scalarBH.C:403
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: scalarBH.C:412
Coord tet
coordinate centered on the grid
Definition: map.h:737
void operator=(const ScalarBH &)
Assignment to another ScalarBH.
Definition: scalarBH.C:422
Scalar ff1
Metric field F_1 of Herdeiro & Radu (2015)
Definition: compobj.h:1032
Base class for stationary compact objects (under development).
Definition: compobj.h:129
void operator=(const Compobj &)
Assignment to another Compobj.
Definition: compobj.C:178
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
Vector beta
Shift vector .
Definition: compobj.h:141
ScalarBH(Map &mp_i, const char *file_name)
Standard constructor.
Definition: scalarBH.C:75
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Scalar ff0
Metric field F_0 of Herdeiro & Radu (2015)
Definition: compobj.h:1031
double rHor
Event horizon coordinate radius.
Definition: compobj.h:1036
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: compobj.C:158
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: compobj.C:242
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: scalarBH.C:444
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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:690
Black hole with scalar hair spacetime (under development).
Definition: compobj.h:1023
Scalar nn
Lapse function N .
Definition: compobj.h:138
Scalar sfield
Scalar field (modulus of Phi)
Definition: compobj.h:1035
Metric gamma
3-metric
Definition: compobj.h:144
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:135
Coord r
r coordinate centered on the grid
Definition: map.h:736