LORENE
blackhole_r_coord.C
1 /*
2  * Method of class Black_hole to express the radial coordinate
3  * in Kerr-Schild coordinates by that in spatially isotropic coordinates
4  *
5  * (see file blackhole.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2006-2007 Keisuke Taniguchi
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: blackhole_r_coord.C,v 1.5 2016/12/05 16:17:48 j_novak Exp $
33  * $Log: blackhole_r_coord.C,v $
34  * Revision 1.5 2016/12/05 16:17:48 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.4 2014/10/13 08:52:46 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.3 2014/10/06 15:13:02 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.2 2008/05/15 19:29:58 k_taniguchi
44  * Change of some parameters.
45  *
46  * Revision 1.1 2007/06/22 01:20:33 k_taniguchi
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_r_coord.C,v 1.5 2016/12/05 16:17:48 j_novak Exp $
51  *
52  */
53 
54 // C++ headers
55 //#include <>
56 
57 // C headers
58 #include <cmath>
59 
60 // Lorene headers
61 #include "blackhole.h"
62 #include "unites.h"
63 #include "utilitaires.h"
64 
65 // Local function
66 namespace Lorene {
67 double gg(double, const double) ;
68 
69 const Scalar Black_hole::r_coord(bool neumann, bool first) const {
70 
71  // Fundamental constants and units
72  // -------------------------------
73  using namespace Unites ;
74 
75  const Mg3d* mg = mp.get_mg() ;
76  int nz = mg->get_nzone() ; // total number of domains
77  int nr = mg->get_nr(0) ;
78  int nt = mg->get_nt(0) ;
79  int np = mg->get_np(0) ;
80 
81  double mass = ggrav * mass_bh ;
82 
83  Scalar r_iso(mp) ;
84  r_iso = mp.r ;
85  r_iso.std_spectral_base() ;
86 
87  Scalar r_are(mp) ;
88  r_are = r_iso ; // Initialization
89  r_are.std_spectral_base() ;
90 
91  // Sets C/M^2 for each case of the lapse boundary condition
92  // --------------------------------------------------------
93  double cc ;
94 
95  if (neumann) { // Neumann boundary condition
96  if (first) { // First condition
97  // d(\alpha \psi)/dr = 0
98  // ---------------------
99  cc = 2. * (sqrt(13.) - 1.) / 3. ;
100  }
101  else { // Second condition
102  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
103  // -----------------------------------------
104  cc = 4. / 3. ;
105  }
106  }
107  else { // Dirichlet boundary condition
108  if (first) { // First condition
109  // (\alpha \psi) = 1/2
110  // -------------------
111  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
112  abort() ;
113  }
114  else { // Second condition
115  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
116  // ----------------------------------
117  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
118  abort() ;
119  }
120  }
121 
122  int ll ;
123  double diff ;
124  double ratio ;
125  double precis = 1.e-15 ;
126  double dp ;
127  double tmp ;
128  double tr ;
129 
130  int nn = 1000 ;
131  assert(nn%4 == 0) ;
132  int mm = nn/4 ;
133  double x1, x2, x3, x4, x5 ;
134  double hh, integ ;
135 
136  // Boole's Rule (Newton-Cotes Integral) for integration
137  // ----------------------------------------------------
138 
139  for (int l=1; l<nz; l++) {
140 
141  for (int i=0; i<nr; i++) {
142 
143  ratio = 1. ;
144  dp = 10. ;
145  tr = r_iso.val_grid_point(l,0,0,i) ;
146 
147  while ( dp > precis ) {
148 
149  diff = 1. ; // Initialization
150  ll = 0 ;
151  dp = 0.1 * dp ;
152 
153  while ( diff > precis ) {
154 
155  ll++ ;
156  tmp = ratio + ll * dp ;
157 
158  double r_max = 2.*mass/tmp/tr ;
159 
160  hh = r_max / double(nn) ;
161  integ = 0. ;
162 
163  for (int n=0; n<mm; n++) {
164 
165  x1 = hh * double(4*n) ;
166  x2 = hh * double(4*n+1) ;
167  x3 = hh * double(4*n+2) ;
168  x4 = hh * double(4*n+3) ;
169  x5 = hh * double(4*n+4) ;
170 
171  integ += (hh/45.) * (14.*gg(x1,cc) + 64.*gg(x2,cc)
172  + 24.*gg(x3,cc) + 64.*gg(x4,cc)
173  + 14.*gg(x5,cc)) ;
174 
175  }
176 
177  diff = -log( tmp ) - integ ;
178 
179  // cout << "diff: " << diff << " x: " << tmp << endl ;
180 
181  }
182 
183  ratio += (ll - 1) * dp ;
184 
185  }
186 
187  for (int j=0; j<nt; j++) {
188  for (int k=0; k<np; k++) {
189 
190  r_are.set_grid_point(l,k,j,i) = ratio ;
191 
192  }
193  }
194 
195  // arrete() ;
196 
197  }
198  }
199 
200  r_are.std_spectral_base() ;
201  r_are.annule_domain(0) ;
202  r_are.raccord(1) ;
203 
204  /*
205  cout << "r_are:" << endl ;
206  for (int l=0; l<nz; l++) {
207  cout << r_are.val_grid_point(l,0,0,0) << " "
208  << r_are.val_grid_point(l,0,0,nr-1) << endl ;
209  }
210  */
211 
212  return r_are ;
213 
214 }
215 
216 //*****************************************************************
217 
218 double gg(double xx, const double cc) {
219 
220  double tcc2 = cc*cc/16. ;
221  double tmp = sqrt(1. - xx + tcc2*pow(xx, 4.)) ;
222 
223  double resu = (-1. + tcc2 * pow(xx, 3.)) / tmp / (1. + tmp) ;
224 
225  return resu ;
226 
227 }
228 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
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
Lorene prototypes.
Definition: app_hor.h:67
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
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
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
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
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
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
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
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
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
Coord r
r coordinate centered on the grid
Definition: map.h:736