LORENE
hole_bhns_rk_phi.C
1/*
2 * Methods of class Hole_bhns to compute a forth-order Runge-Kutta
3 * integration to the phi direction for the solution of the Killing vectors
4 *
5 * (see file hole_bhns.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 hole_bhns_rk_phi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_rk_phi.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $" ;
30
31/*
32 * $Id: hole_bhns_rk_phi.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
33 * $Log: hole_bhns_rk_phi.C,v $
34 * Revision 1.4 2014/10/13 08:53:00 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.3 2014/10/06 15:13:10 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.2 2008/07/02 20:47:55 k_taniguchi
41 * Typos removed.
42 *
43 * Revision 1.1 2008/05/15 19:10:31 k_taniguchi
44 * *** empty log message ***
45 *
46 *
47 * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_rk_phi.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
48 *
49 */
50
51// C++ headers
52//#include <>
53
54// C headers
55#include <cmath>
56
57// Lorene headers
58#include "hole_bhns.h"
59#include "unites.h"
60#include "utilitaires.h"
61
62 //--------------------------------------------------//
63 // Forth-order Runge-Kutta on the equator //
64 //--------------------------------------------------//
65
66namespace Lorene {
67Tbl Hole_bhns::runge_kutta_phi(const Tbl& xi_i, const double& phi_i,
68 const int& nrk_phi) const {
69
70 using namespace Unites ;
71
72 const Mg3d* mg = mp.get_mg() ;
73 int np = mg->get_np(1) ;
74
75 Tbl xi_f(3) ; // xi_f(0)=xi_hat{theta}, xi_f(1)=xi_hat{phi}, xi_f(2)=L
76 xi_f.set_etat_qcq() ;
77
78 if (kerrschild) {
79
80 cout << "Not yet prepared!!!" << endl ;
81 abort() ;
82
83 }
84 else { // Isotropic coordinates
85
86 // Initial data at phi=0 on the equator
87 double xi_t0 = xi_i(0) ; // xi_hat{theta}
88 double xi_p0 = xi_i(1) ; // xi_hat{phi}
89 double xi_l0 = xi_i(2) ; // L
90 double phi0 = phi_i ;
91
92 double dp = 2. * M_PI / double(np) / double(nrk_phi) ;
93
94 double rah = rad_ah() ;
95
96 Scalar dlnconfo(mp) ;
97 dlnconfo = confo_tot.dsdt() / confo_tot ;
98 dlnconfo.std_spectral_base() ;
99
100 Scalar laplnconfo(mp) ;
101 laplnconfo = confo_tot.lapang() / confo_tot ;
102 laplnconfo.std_spectral_base() ;
103
104 Scalar confo2(mp) ;
105 confo2 = confo_tot * confo_tot ;
106 confo2.std_spectral_base() ;
107
108 double xi_t1, xi_t2, xi_t3, xi_t4, xi_tf ;
109 double xi_p1, xi_p2, xi_p3, xi_p4, xi_pf ;
110 double xi_l1, xi_l2, xi_l3, xi_l4, xi_lf ;
111 double f1, f2, f3, f4 ;
112 double g1, g2, g3, g4 ;
113 double h1, h2, h3, h4 ;
114
115 // Forth-order Runge-Kutta
116 // (nrk_phi times steps between two collocation points)
117 // ----------------------------------------------------
118
119 for (int i=0; i<nrk_phi; i++) {
120
121 // First
122 f1 = - xi_l0 * rah * confo2.val_point(rah, M_PI/2., phi0)
123 + 2. * xi_p0 * dlnconfo.val_point(rah, M_PI/2., phi0) ;
124 g1 = -2. * xi_t0 * dlnconfo.val_point(rah, M_PI/2., phi0) ;
125 h1 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0)) * xi_t0
126 / rah / confo2.val_point(rah, M_PI/2., phi0) ;
127
128 xi_t1 = dp * f1 ;
129 xi_p1 = dp * g1 ;
130 xi_l1 = dp * h1 ;
131
132 // Second
133 f2 = - (xi_l0+0.5*xi_l1) * rah
134 * confo2.val_point(rah, M_PI/2., phi0+0.5*dp)
135 + 2. * (xi_p0+0.5*xi_p1)
136 * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
137 g2 = -2. * (xi_t0+0.5*xi_t1)
138 * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
139 h2 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0+0.5*dp))
140 * (xi_t0+0.5*xi_t1) / rah
141 / confo2.val_point(rah, M_PI/2., phi0+0.5*dp) ;
142
143 xi_t2 = dp * f2 ;
144 xi_p2 = dp * g2 ;
145 xi_l2 = dp * h2 ;
146
147 // Third
148 f3 = - (xi_l0+0.5*xi_l2) * rah
149 * confo2.val_point(rah, M_PI/2., phi0+0.5*dp)
150 + 2. * (xi_p0+0.5*xi_p2)
151 * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
152 g3 = -2. * (xi_t0+0.5*xi_t2)
153 * dlnconfo.val_point(rah, M_PI/2., phi0+0.5*dp) ;
154 h3 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0+0.5*dp))
155 * (xi_t0+0.5*xi_t2) / rah
156 / confo2.val_point(rah, M_PI/2., phi0+0.5*dp) ;
157
158 xi_t3 = dp * f3 ;
159 xi_p3 = dp * g3 ;
160 xi_l3 = dp * h3 ;
161
162 // Forth
163 f4 = - (xi_l0+xi_l3) * rah * confo2.val_point(rah, M_PI/2., phi0+dp)
164 + 2. * (xi_p0+xi_p3) * dlnconfo.val_point(rah, M_PI/2., phi0+dp) ;
165 g4 = -2. * (xi_t0+xi_t3) * dlnconfo.val_point(rah, M_PI/2., phi0+dp) ;
166 h4 = (1. - 2.*laplnconfo.val_point(rah, M_PI/2., phi0+dp))
167 * (xi_t0+xi_t3) / rah / confo2.val_point(rah, M_PI/2., phi0+dp) ;
168
169 xi_t4 = dp * f4 ;
170 xi_p4 = dp * g4 ;
171 xi_l4 = dp * h4 ;
172
173 // Final results
174 // -------------
175 xi_tf = xi_t0 + (xi_t1 + 2.*xi_t2 + 2.*xi_t3 + xi_t4) / 6. ;
176 xi_pf = xi_p0 + (xi_p1 + 2.*xi_p2 + 2.*xi_p3 + xi_p4) / 6. ;
177 xi_lf = xi_l0 + (xi_l1 + 2.*xi_l2 + 2.*xi_l3 + xi_l4) / 6. ;
178
179 // Final results are put into the initial data
180 // in order for the next step
181 // -------------------------------------------
182 xi_t0 = xi_tf ;
183 xi_p0 = xi_pf ;
184 xi_l0 = xi_lf ;
185
186 } // End of the loop
187
188 xi_f.set(0) = xi_tf ;
189 xi_f.set(1) = xi_pf ;
190 xi_f.set(2) = xi_lf ;
191
192 }
193
194 return xi_f ;
195
196}
197}
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
virtual double rad_ah() const
Radius of the apparent horizon.
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Tbl runge_kutta_phi(const Tbl &xi_i, const double &phi_i, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the phi direction for the solution of the Killing ve...
Scalar confo_tot
Total conformal factor.
Definition hole_bhns.h:169
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
const Scalar & dsdt() const
Returns of *this .
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_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition scalar.C:890
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.