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