blackhole_bc.C
/*
* Methods of class Black_hole to compute the inner boundary condition
* at the excised surface
*
* (see file blackhole.h for documentation).
*
*/
/*
* Copyright (c) 2005-2007 Keisuk Taniguchi
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 2
* as published by the Free Software Foundation.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char 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 $" ;
/*
* $Id: blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $
* $Log: blackhole_bc.C,v $
* Revision 1.5 2014/10/13 08:52:45 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.4 2014/10/06 15:13:02 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.3 2008/07/03 14:53:47 k_taniguchi
* Modification of a signature in bc_shift_x and bc_shift_y.
*
* Revision 1.2 2008/05/15 19:25:43 k_taniguchi
* Change of some parameters.
*
* Revision 1.1 2007/06/22 01:18:23 k_taniguchi
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_bc.C,v 1.5 2014/10/13 08:52:45 j_novak Exp $
*
*/
// C++ headers
//#include <>
// C headers
#include <cmath>
// Lorene headers
#include "blackhole.h"
#include "valeur.h"
#include "grilles.h"
#include "unites.h"
//----------------------------------//
// Inner boundary condition //
//----------------------------------//
namespace Lorene {
const Valeur Black_hole::bc_lapconf(bool neumann, bool first) const {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
const Mg3d* mg_angu = mg->get_angu() ;
Valeur bc(mg_angu) ;
if (kerrschild) {
cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
abort() ;
/*
if (neumann) {
if (first) {
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
Scalar tmp(mp) ;
tmp = - pow(lapse_bh,3.) * ggrav * mass_bh / rr / rr ;
// dlapse/dr = 0
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
}
else {
bc = 0. ; // dlapse/dr = 0.25*lapse/rr
}
}
else {
if (first) { // The poisson solver in LORENE assumes the
// asymptotic behavior of the function -> 0
bc = 0.5 - 1./sqrt(2.) ; // <- bc of the real function = 0.5
}
else {
bc = 0. ; // <- bc of the real function = 1./sqrt(2.)
}
}
*/
}
else { // Isotropic coordinates with the maximal slicing
if (neumann) {
if (first) {
bc = 0. ; // d(lapconf)/dr = 0
}
else {
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
Scalar tmp(mp) ;
tmp = 0.5 * lapconf / rr ;
// d(lapconf)/dr = lapconf/2/rah
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
}
}
else {
if (first) {
cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
abort() ;
// bc = - 0.5 ; // lapconf = 0.5
}
else {
cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
abort() ;
/*
Scalar r_are(mp) ;
r_are = r_coord(neumann, first) ;
r_are.std_spectral_base() ;
r_are.annule_domain(0) ;
r_are.raccord(1) ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
Scalar tmp(mp) ;
tmp = sqrt(0.5*r_are) - 1. ; // lapse = 1./sqrt(2.)
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
*/
}
}
}
bc.std_base_scal() ;
return bc ;
}
const Valeur Black_hole::bc_shift_x(double omega_r) const {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
const Mg3d* mg_angu = mg->get_angu() ;
Valeur bc(mg_angu) ;
Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
Scalar st(mp) ;
st = mp.sint ;
st.std_spectral_base() ;
Scalar cp(mp) ;
cp = mp.cosp ;
cp.std_spectral_base() ;
Scalar yy(mp) ;
yy = mp.y ;
yy.std_spectral_base() ;
Scalar tmp(mp) ;
if (kerrschild) {
cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
abort() ;
/*
tmp = lapse_bh * (lapse / confo / confo) * st * cp
- omega_r * yy - shift_bh(1) ;
*/
// tmp = lap_bh * lap_bh * st * cp - omega_r * yy ;
}
else { // Isotropic coordinates with the maximal slicing
// Note: the signature of omega_r is opposite to that in the
// binary case because of the direction of the spin
tmp = lapconf / pow(confo, 3.) * st * cp + omega_r * yy ;
// tmp = lapconf / pow(confo, 3.) * st * cp - omega_r * yy ;
}
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
bc.base = *bases[0] ;
// bc.std_base_scal() ;
for (int i=0; i<3; i++)
delete bases[i] ;
delete [] bases ;
return bc ;
}
const Valeur Black_hole::bc_shift_y(double omega_r) const {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
const Mg3d* mg_angu = mg->get_angu() ;
Valeur bc(mg_angu) ;
Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
Scalar st(mp) ;
st = mp.sint ;
st.std_spectral_base() ;
Scalar sp(mp) ;
sp = mp.sinp ;
sp.std_spectral_base() ;
Scalar xx(mp) ;
xx = mp.x ;
xx.std_spectral_base() ;
Scalar tmp(mp) ;
if (kerrschild) {
cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
abort() ;
/*
tmp = lapse_bh * (lapse / confo / confo) * st * sp
+ omega_r * xx - shift_bh(2) ;
*/
// tmp = lap_bh * lap_bh * st * sp + omega_r * xx ;
}
else {
// Note: the signature of omega_r is opposite to that in the
// binary case because of the direction of the spin
tmp = lapconf / pow(confo, 3.) * st * sp - omega_r * xx ;
// tmp = lapconf / pow(confo, 3.) * st * sp + omega_r * xx ;
}
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
bc.base = *bases[1] ;
// bc.std_base_scal() ;
for (int i=0; i<3; i++)
delete bases[i] ;
delete [] bases ;
return bc ;
}
const Valeur Black_hole::bc_shift_z() const {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
const Mg3d* mg_angu = mg->get_angu() ;
Valeur bc(mg_angu) ;
Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
Scalar ct(mp) ;
ct = mp.cost ;
ct.std_spectral_base() ;
Scalar tmp(mp) ;
if (kerrschild) {
cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
abort() ;
/*
tmp = lapse_bh * (lapse / confo / confo) * ct
- shift_bh(3) ;
*/
// tmp = lap_bh * lap_bh * ct ;
}
else {
tmp = lapconf / pow(confo, 3.) * ct ;
}
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
bc.base = *bases[2] ;
// bc.std_base_scal() ;
for (int i=0; i<3; i++)
delete bases[i] ;
delete [] bases ;
return bc ;
}
const Valeur Black_hole::bc_confo() const {
// Fundamental constants and units
// -------------------------------
using namespace Unites ;
const Mg3d* mg = mp.get_mg() ;
const Mg3d* mg_angu = mg->get_angu() ;
Valeur bc(mg_angu) ;
double mass = ggrav * mass_bh ;
Scalar rr(mp) ;
rr = mp.r ;
rr.std_spectral_base() ;
Scalar st(mp) ;
st = mp.sint ;
st.std_spectral_base() ;
Scalar ct(mp) ;
ct = mp.cost ;
ct.std_spectral_base() ;
Scalar sp(mp) ;
sp = mp.sinp ;
sp.std_spectral_base() ;
Scalar cp(mp) ;
cp = mp.cosp ;
cp.std_spectral_base() ;
int nt = mg->get_nt(0) ;
int np = mg->get_np(0) ;
Scalar tmp(mp) ;
if (kerrschild) { // Assumes that r_BH = 1.
cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
abort() ;
/*
Scalar divshift(mp) ;
divshift = shift_rs(1).deriv(1) + shift_rs(2).deriv(2)
+ shift_rs(3).deriv(3) ;
divshift.std_spectral_base() ;
Scalar llshift(mp) ;
llshift = st*cp*shift_rs(1) + st*sp*shift_rs(2) + ct*shift_rs(3) ;
llshift.std_spectral_base() ;
Scalar lldllsh = llshift.dsdr() ;
lldllsh.std_spectral_base() ;
Scalar tmp1 = divshift ;
Scalar tmp2 = -3.*lldllsh ;
Scalar tmp5 = 0.5*confo*(lapse_bh*confo*confo/lapse - 1.)/rr ;
tmp1.set_dzpuis(tmp5.get_dzpuis()) ;
tmp2.set_dzpuis(tmp5.get_dzpuis()) ;
Scalar tmp3 = 2. * lapse_bh * lapse_bh * mass * llshift / rr / rr ;
Scalar tmp4 = 4. * pow(lapse_bh,3.) * mass * (1.+3.*mass/rr)
* lapse_rs / rr / rr ;
tmp3.set_dzpuis(tmp5.get_dzpuis()) ;
tmp4.set_dzpuis(tmp5.get_dzpuis()) ;
tmp = tmp5 + pow(confo,3.)*(tmp1+tmp2+tmp3+tmp4)/12./lapse/lapse_bh ;
*/
// tmp = -0.5 * (1. - 2. * mass / rr) / rr ;
}
else { // Isotropic coordinates with the maximal slicing
Scalar divshift(mp) ;
divshift = shift(1).deriv(1) + shift(2).deriv(2)
+ shift(3).deriv(3) ;
divshift.std_spectral_base() ;
Scalar llshift(mp) ;
llshift = st*cp*shift(1) + st*sp*shift(2) + ct*shift(3) ;
llshift.std_spectral_base() ;
Scalar lldllsh = llshift.dsdr() ;
lldllsh.std_spectral_base() ;
Scalar tmp1 = divshift ;
Scalar tmp2 = -3.*lldllsh ;
Scalar tmp5 = - 0.5 * confo / rr ;
tmp1.set_dzpuis(tmp5.get_dzpuis()) ;
tmp2.set_dzpuis(tmp5.get_dzpuis()) ;
tmp = tmp5 + pow(confo, 4.) * (tmp1 + tmp2) / 12. / lapconf ;
}
bc = 1. ;
for (int j=0; j<nt; j++) {
for (int k=0; k<np; k++) {
bc.set(0,k,j,0) = tmp.val_grid_point(1,k,j,0) ;
}
}
bc.std_base_scal() ;
return bc ;
}
}
Lorene::Base_val
Bases of the spectral expansions.
Definition base_val.h:322
Lorene::Black_hole::bc_confo
const Valeur bc_confo() const
Boundary condition on the apparent horizon of the black hole for the conformal factor: 2-D Valeur.
Definition blackhole_bc.C:410
Lorene::Black_hole::bc_shift_x
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...
Definition blackhole_bc.C:211
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Lorene::Black_hole::lapconf
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition blackhole.h:97
Lorene::Black_hole::bc_lapconf
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.
Definition blackhole_bc.C:71
Lorene::Black_hole::shift
Vector shift
Shift vector generated by the black hole.
Definition blackhole.h:109
Lorene::Black_hole::bc_shift_z
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...
Definition blackhole_bc.C:350
Lorene::Black_hole::bc_shift_y
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...
Definition blackhole_bc.C:281
Lorene::Black_hole::kerrschild
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Lorene::Black_hole::confo
Scalar confo
Conformal factor generated by the black hole.
Definition blackhole.h:118
Lorene::Black_hole::mass_bh
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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::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::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:557
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
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition scalar_deriv.C:113
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::set
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
Lorene::Valeur::base
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
Lorene::Valeur::std_base_scal
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.