LORENE
blackhole_bc.C
1/*
2 * Methods of class Black_hole to compute the inner boundary condition
3 * at the excised surface
4 *
5 * (see file blackhole.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2005-2007 Keisuk 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_bc_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $" ;
30
31/*
32 * $Id: blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $
33 * $Log: blackhole_bc.C,v $
34 * Revision 1.5 2014/10/13 08:52:45 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.4 2014/10/06 15:13:02 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.3 2008/07/03 14:53:47 k_taniguchi
41 * Modification of a signature in bc_shift_x and bc_shift_y.
42 *
43 * Revision 1.2 2008/05/15 19:25:43 k_taniguchi
44 * Change of some parameters.
45 *
46 * Revision 1.1 2007/06/22 01:18:23 k_taniguchi
47 * *** empty log message ***
48 *
49 *
50 * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_bc.C,v 1.5 2014/10/13 08:52:45 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 "valeur.h"
63#include "grilles.h"
64#include "unites.h"
65
66 //----------------------------------//
67 // Inner boundary condition //
68 //----------------------------------//
69
70namespace Lorene {
71const Valeur Black_hole::bc_lapconf(bool neumann, bool first) const {
72
73 // Fundamental constants and units
74 // -------------------------------
75 using namespace Unites ;
76
77 const Mg3d* mg = mp.get_mg() ;
78 const Mg3d* mg_angu = mg->get_angu() ;
79 Valeur bc(mg_angu) ;
80
81 if (kerrschild) {
82
83 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
84 abort() ;
85
86 /*
87 if (neumann) {
88
89 if (first) {
90
91 Scalar rr(mp) ;
92 rr = mp.r ;
93 rr.std_spectral_base() ;
94
95 int nt = mg->get_nt(0) ;
96 int np = mg->get_np(0) ;
97
98 Scalar tmp(mp) ;
99
100 tmp = - pow(lapse_bh,3.) * ggrav * mass_bh / rr / rr ;
101 // dlapse/dr = 0
102
103 bc = 1. ;
104 for (int j=0; j<nt; j++) {
105 for (int k=0; k<np; k++) {
106 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
107 }
108 }
109 }
110 else {
111
112 bc = 0. ; // dlapse/dr = 0.25*lapse/rr
113
114 }
115 }
116 else {
117 if (first) { // The poisson solver in LORENE assumes the
118 // asymptotic behavior of the function -> 0
119 bc = 0.5 - 1./sqrt(2.) ; // <- bc of the real function = 0.5
120
121 }
122 else {
123
124 bc = 0. ; // <- bc of the real function = 1./sqrt(2.)
125
126 }
127 }
128 */
129 }
130 else { // Isotropic coordinates with the maximal slicing
131
132 if (neumann) {
133
134 if (first) {
135
136 bc = 0. ; // d(lapconf)/dr = 0
137
138 }
139 else {
140
141 Scalar rr(mp) ;
142 rr = mp.r ;
143 rr.std_spectral_base() ;
144
145 int nt = mg->get_nt(0) ;
146 int np = mg->get_np(0) ;
147
148 Scalar tmp(mp) ;
149
150 tmp = 0.5 * lapconf / rr ;
151 // d(lapconf)/dr = lapconf/2/rah
152
153 bc = 1. ;
154 for (int j=0; j<nt; j++) {
155 for (int k=0; k<np; k++) {
156 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
157 }
158 }
159
160
161 }
162
163 }
164 else {
165
166 if (first) {
167
168 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
169 abort() ;
170 // bc = - 0.5 ; // lapconf = 0.5
171
172 }
173 else {
174
175 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
176 abort() ;
177 /*
178 Scalar r_are(mp) ;
179 r_are = r_coord(neumann, first) ;
180 r_are.std_spectral_base() ;
181 r_are.annule_domain(0) ;
182 r_are.raccord(1) ;
183
184 int nt = mg->get_nt(0) ;
185 int np = mg->get_np(0) ;
186
187 Scalar tmp(mp) ;
188
189 tmp = sqrt(0.5*r_are) - 1. ; // lapse = 1./sqrt(2.)
190
191 bc = 1. ;
192 for (int j=0; j<nt; j++) {
193 for (int k=0; k<np; k++) {
194 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
195 }
196 }
197 */
198
199 }
200
201 }
202
203 }
204
205 bc.std_base_scal() ;
206 return bc ;
207
208}
209
210
211const Valeur Black_hole::bc_shift_x(double omega_r) const {
212
213 // Fundamental constants and units
214 // -------------------------------
215 using namespace Unites ;
216
217 const Mg3d* mg = mp.get_mg() ;
218 const Mg3d* mg_angu = mg->get_angu() ;
219 Valeur bc(mg_angu) ;
220
221 Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
222
223 Scalar rr(mp) ;
224 rr = mp.r ;
225 rr.std_spectral_base() ;
226 Scalar st(mp) ;
227 st = mp.sint ;
228 st.std_spectral_base() ;
229 Scalar cp(mp) ;
230 cp = mp.cosp ;
231 cp.std_spectral_base() ;
232 Scalar yy(mp) ;
233 yy = mp.y ;
234 yy.std_spectral_base() ;
235
236 Scalar tmp(mp) ;
237
238 if (kerrschild) {
239
240 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
241 abort() ;
242 /*
243 tmp = lapse_bh * (lapse / confo / confo) * st * cp
244 - omega_r * yy - shift_bh(1) ;
245 */
246
247 // tmp = lap_bh * lap_bh * st * cp - omega_r * yy ;
248 }
249 else { // Isotropic coordinates with the maximal slicing
250
251 // Note: the signature of omega_r is opposite to that in the
252 // binary case because of the direction of the spin
253 tmp = lapconf / pow(confo, 3.) * st * cp + omega_r * yy ;
254 // tmp = lapconf / pow(confo, 3.) * st * cp - omega_r * yy ;
255
256 }
257
258 int nt = mg->get_nt(0) ;
259 int np = mg->get_np(0) ;
260
261 bc = 1. ;
262 for (int j=0; j<nt; j++) {
263 for (int k=0; k<np; k++) {
264 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
265 }
266 }
267
268 bc.base = *bases[0] ;
269 // bc.std_base_scal() ;
270
271 for (int i=0; i<3; i++)
272 delete bases[i] ;
273
274 delete [] bases ;
275
276 return bc ;
277
278}
279
280
281const Valeur Black_hole::bc_shift_y(double omega_r) const {
282
283 // Fundamental constants and units
284 // -------------------------------
285 using namespace Unites ;
286
287 const Mg3d* mg = mp.get_mg() ;
288 const Mg3d* mg_angu = mg->get_angu() ;
289 Valeur bc(mg_angu) ;
290
291 Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
292
293 Scalar rr(mp) ;
294 rr = mp.r ;
295 rr.std_spectral_base() ;
296 Scalar st(mp) ;
297 st = mp.sint ;
298 st.std_spectral_base() ;
299 Scalar sp(mp) ;
300 sp = mp.sinp ;
301 sp.std_spectral_base() ;
302 Scalar xx(mp) ;
303 xx = mp.x ;
304 xx.std_spectral_base() ;
305
306 Scalar tmp(mp) ;
307
308 if (kerrschild) {
309
310 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
311 abort() ;
312 /*
313 tmp = lapse_bh * (lapse / confo / confo) * st * sp
314 + omega_r * xx - shift_bh(2) ;
315 */
316 // tmp = lap_bh * lap_bh * st * sp + omega_r * xx ;
317 }
318 else {
319
320 // Note: the signature of omega_r is opposite to that in the
321 // binary case because of the direction of the spin
322 tmp = lapconf / pow(confo, 3.) * st * sp - omega_r * xx ;
323 // tmp = lapconf / pow(confo, 3.) * st * sp + omega_r * xx ;
324
325 }
326
327 int nt = mg->get_nt(0) ;
328 int np = mg->get_np(0) ;
329
330 bc = 1. ;
331 for (int j=0; j<nt; j++) {
332 for (int k=0; k<np; k++) {
333 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
334 }
335 }
336
337 bc.base = *bases[1] ;
338 // bc.std_base_scal() ;
339
340 for (int i=0; i<3; i++)
341 delete bases[i] ;
342
343 delete [] bases ;
344
345 return bc ;
346
347}
348
349
351
352 // Fundamental constants and units
353 // -------------------------------
354 using namespace Unites ;
355
356 const Mg3d* mg = mp.get_mg() ;
357 const Mg3d* mg_angu = mg->get_angu() ;
358 Valeur bc(mg_angu) ;
359
360 Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
361
362 Scalar rr(mp) ;
363 rr = mp.r ;
364 rr.std_spectral_base() ;
365 Scalar ct(mp) ;
366 ct = mp.cost ;
367 ct.std_spectral_base() ;
368
369 Scalar tmp(mp) ;
370
371 if (kerrschild) {
372
373 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
374 abort() ;
375 /*
376 tmp = lapse_bh * (lapse / confo / confo) * ct
377 - shift_bh(3) ;
378 */
379 // tmp = lap_bh * lap_bh * ct ;
380 }
381 else {
382
383 tmp = lapconf / pow(confo, 3.) * ct ;
384
385 }
386
387 int nt = mg->get_nt(0) ;
388 int np = mg->get_np(0) ;
389
390 bc = 1. ;
391 for (int j=0; j<nt; j++) {
392 for (int k=0; k<np; k++) {
393 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
394 }
395 }
396
397 bc.base = *bases[2] ;
398 // bc.std_base_scal() ;
399
400 for (int i=0; i<3; i++)
401 delete bases[i] ;
402
403 delete [] bases ;
404
405 return bc ;
406
407}
408
409
411
412 // Fundamental constants and units
413 // -------------------------------
414 using namespace Unites ;
415
416 const Mg3d* mg = mp.get_mg() ;
417 const Mg3d* mg_angu = mg->get_angu() ;
418 Valeur bc(mg_angu) ;
419
420 double mass = ggrav * mass_bh ;
421
422 Scalar rr(mp) ;
423 rr = mp.r ;
424 rr.std_spectral_base() ;
425 Scalar st(mp) ;
426 st = mp.sint ;
427 st.std_spectral_base() ;
428 Scalar ct(mp) ;
429 ct = mp.cost ;
430 ct.std_spectral_base() ;
431 Scalar sp(mp) ;
432 sp = mp.sinp ;
433 sp.std_spectral_base() ;
434 Scalar cp(mp) ;
435 cp = mp.cosp ;
436 cp.std_spectral_base() ;
437
438 int nt = mg->get_nt(0) ;
439 int np = mg->get_np(0) ;
440
441 Scalar tmp(mp) ;
442
443 if (kerrschild) { // Assumes that r_BH = 1.
444
445 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
446 abort() ;
447 /*
448 Scalar divshift(mp) ;
449 divshift = shift_rs(1).deriv(1) + shift_rs(2).deriv(2)
450 + shift_rs(3).deriv(3) ;
451 divshift.std_spectral_base() ;
452
453 Scalar llshift(mp) ;
454 llshift = st*cp*shift_rs(1) + st*sp*shift_rs(2) + ct*shift_rs(3) ;
455 llshift.std_spectral_base() ;
456
457 Scalar lldllsh = llshift.dsdr() ;
458 lldllsh.std_spectral_base() ;
459
460 Scalar tmp1 = divshift ;
461 Scalar tmp2 = -3.*lldllsh ;
462
463 Scalar tmp5 = 0.5*confo*(lapse_bh*confo*confo/lapse - 1.)/rr ;
464 tmp1.set_dzpuis(tmp5.get_dzpuis()) ;
465 tmp2.set_dzpuis(tmp5.get_dzpuis()) ;
466
467 Scalar tmp3 = 2. * lapse_bh * lapse_bh * mass * llshift / rr / rr ;
468 Scalar tmp4 = 4. * pow(lapse_bh,3.) * mass * (1.+3.*mass/rr)
469 * lapse_rs / rr / rr ;
470
471 tmp3.set_dzpuis(tmp5.get_dzpuis()) ;
472 tmp4.set_dzpuis(tmp5.get_dzpuis()) ;
473
474 tmp = tmp5 + pow(confo,3.)*(tmp1+tmp2+tmp3+tmp4)/12./lapse/lapse_bh ;
475 */
476 // tmp = -0.5 * (1. - 2. * mass / rr) / rr ;
477
478 }
479 else { // Isotropic coordinates with the maximal slicing
480
481 Scalar divshift(mp) ;
482 divshift = shift(1).deriv(1) + shift(2).deriv(2)
483 + shift(3).deriv(3) ;
484 divshift.std_spectral_base() ;
485
486 Scalar llshift(mp) ;
487 llshift = st*cp*shift(1) + st*sp*shift(2) + ct*shift(3) ;
488 llshift.std_spectral_base() ;
489
490 Scalar lldllsh = llshift.dsdr() ;
491 lldllsh.std_spectral_base() ;
492
493 Scalar tmp1 = divshift ;
494 Scalar tmp2 = -3.*lldllsh ;
495
496 Scalar tmp5 = - 0.5 * confo / rr ;
497
498 tmp1.set_dzpuis(tmp5.get_dzpuis()) ;
499 tmp2.set_dzpuis(tmp5.get_dzpuis()) ;
500
501 tmp = tmp5 + pow(confo, 4.) * (tmp1 + tmp2) / 12. / lapconf ;
502
503 }
504
505 bc = 1. ;
506 for (int j=0; j<nt; j++) {
507 for (int k=0; k<np; k++) {
508 bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
509 }
510 }
511
512 bc.std_base_scal() ;
513 return bc ;
514
515}
516}
Bases of the spectral expansions.
Definition base_val.h:322
const Valeur bc_confo() const
Boundary condition on the apparent horizon of the black hole for the conformal factor: 2-D Valeur.
const Valeur bc_shift_x(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fx\f direct...
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition blackhole.h:97
const Valeur bc_lapconf(bool neumann, bool first) const
Boundary condition on the apparent horizon of the black hole for the lapse function: 2-D Valeur.
Vector shift
Shift vector generated by the black hole.
Definition blackhole.h:109
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fz\f direct...
const Valeur bc_shift_y(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fy\f direct...
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Scalar confo
Conformal factor generated by the black hole.
Definition blackhole.h:118
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Coord cosp
Definition map.h:724
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
Coord sinp
Definition map.h:723
Coord x
x coordinate centered on the grid
Definition map.h:726
Coord cost
Definition map.h:722
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
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
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
const Scalar & dsdr() const
Returns of *this .
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.