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 char blackhole_r_coord_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_r_coord.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
30 
31 /*
32  * $Id: blackhole_r_coord.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
33  * $Log: blackhole_r_coord.C,v $
34  * Revision 1.4 2014/10/13 08:52:46 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.3 2014/10/06 15:13:02 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.2 2008/05/15 19:29:58 k_taniguchi
41  * Change of some parameters.
42  *
43  * Revision 1.1 2007/06/22 01:20:33 k_taniguchi
44  * *** empty log message ***
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_r_coord.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
48  *
49  */
50 
51 // C++ headers
52 //#include <>
53 
54 // C headers
55 #include <cmath>
56 
57 // Lorene headers
58 #include "blackhole.h"
59 #include "unites.h"
60 #include "utilitaires.h"
61 
62 // Local function
63 namespace Lorene {
64 double gg(double, const double) ;
65 
66 const Scalar Black_hole::r_coord(bool neumann, bool first) const {
67 
68  // Fundamental constants and units
69  // -------------------------------
70  using namespace Unites ;
71 
72  const Mg3d* mg = mp.get_mg() ;
73  int nz = mg->get_nzone() ; // total number of domains
74  int nr = mg->get_nr(0) ;
75  int nt = mg->get_nt(0) ;
76  int np = mg->get_np(0) ;
77 
78  double mass = ggrav * mass_bh ;
79 
80  Scalar r_iso(mp) ;
81  r_iso = mp.r ;
82  r_iso.std_spectral_base() ;
83 
84  Scalar r_are(mp) ;
85  r_are = r_iso ; // Initialization
86  r_are.std_spectral_base() ;
87 
88  // Sets C/M^2 for each case of the lapse boundary condition
89  // --------------------------------------------------------
90  double cc ;
91 
92  if (neumann) { // Neumann boundary condition
93  if (first) { // First condition
94  // d(\alpha \psi)/dr = 0
95  // ---------------------
96  cc = 2. * (sqrt(13.) - 1.) / 3. ;
97  }
98  else { // Second condition
99  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
100  // -----------------------------------------
101  cc = 4. / 3. ;
102  }
103  }
104  else { // Dirichlet boundary condition
105  if (first) { // First condition
106  // (\alpha \psi) = 1/2
107  // -------------------
108  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
109  abort() ;
110  }
111  else { // Second condition
112  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
113  // ----------------------------------
114  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
115  abort() ;
116  }
117  }
118 
119  int ll ;
120  double diff ;
121  double ratio ;
122  double precis = 1.e-15 ;
123  double dp ;
124  double tmp ;
125  double tr ;
126 
127  int nn = 1000 ;
128  assert(nn%4 == 0) ;
129  int mm = nn/4 ;
130  double x1, x2, x3, x4, x5 ;
131  double hh, integ ;
132 
133  // Boole's Rule (Newton-Cotes Integral) for integration
134  // ----------------------------------------------------
135 
136  for (int l=1; l<nz; l++) {
137 
138  for (int i=0; i<nr; i++) {
139 
140  ratio = 1. ;
141  dp = 10. ;
142  tr = r_iso.val_grid_point(l,0,0,i) ;
143 
144  while ( dp > precis ) {
145 
146  diff = 1. ; // Initialization
147  ll = 0 ;
148  dp = 0.1 * dp ;
149 
150  while ( diff > precis ) {
151 
152  ll++ ;
153  tmp = ratio + ll * dp ;
154 
155  double r_max = 2.*mass/tmp/tr ;
156 
157  hh = r_max / double(nn) ;
158  integ = 0. ;
159 
160  for (int n=0; n<mm; n++) {
161 
162  x1 = hh * double(4*n) ;
163  x2 = hh * double(4*n+1) ;
164  x3 = hh * double(4*n+2) ;
165  x4 = hh * double(4*n+3) ;
166  x5 = hh * double(4*n+4) ;
167 
168  integ += (hh/45.) * (14.*gg(x1,cc) + 64.*gg(x2,cc)
169  + 24.*gg(x3,cc) + 64.*gg(x4,cc)
170  + 14.*gg(x5,cc)) ;
171 
172  }
173 
174  diff = -log( tmp ) - integ ;
175 
176  // cout << "diff: " << diff << " x: " << tmp << endl ;
177 
178  }
179 
180  ratio += (ll - 1) * dp ;
181 
182  }
183 
184  for (int j=0; j<nt; j++) {
185  for (int k=0; k<np; k++) {
186 
187  r_are.set_grid_point(l,k,j,i) = ratio ;
188 
189  }
190  }
191 
192  // arrete() ;
193 
194  }
195  }
196 
197  r_are.std_spectral_base() ;
198  r_are.annule_domain(0) ;
199  r_are.raccord(1) ;
200 
201  /*
202  cout << "r_are:" << endl ;
203  for (int l=0; l<nz; l++) {
204  cout << r_are.val_grid_point(l,0,0,0) << " "
205  << r_are.val_grid_point(l,0,0,nr-1) << endl ;
206  }
207  */
208 
209  return r_are ;
210 
211 }
212 
213 //*****************************************************************
214 
215 double gg(double xx, const double cc) {
216 
217  double tcc2 = cc*cc/16. ;
218  double tmp = sqrt(1. - xx + tcc2*pow(xx, 4.)) ;
219 
220  double resu = (-1. + tcc2 * pow(xx, 3.)) / tmp / (1. + tmp) ;
221 
222  return resu ;
223 
224 }
225 }
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: scalar_raccord.C:62
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:273
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Black_hole::mass_bh
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Scalar::set_grid_point
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
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Unites
Standard units of space, time and mass.
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Black_hole::r_coord
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
Definition: blackhole_r_coord.C:66
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Scalar::std_spectral_base
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
Lorene::Scalar::val_grid_point
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:637