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
29char 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
63namespace Lorene {
64double gg(double, const double) ;
65
66const 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
215double 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}
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Coord r
r coordinate centered on the grid
Definition map.h:718
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
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
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
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
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.