LORENE
bhole_pseudo_kerr.C
1/*
2 * Copyright (c) 2000-2001 Philippe Grandclement
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char bhole_pseudo_kerr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole/bhole_pseudo_kerr.C,v 1.7 2014/10/13 08:52:40 j_novak Exp $" ;
24
25/*
26 * $Id: bhole_pseudo_kerr.C,v 1.7 2014/10/13 08:52:40 j_novak Exp $
27 * $Log: bhole_pseudo_kerr.C,v $
28 * Revision 1.7 2014/10/13 08:52:40 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.6 2014/10/06 15:12:58 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.5 2008/08/19 06:41:59 j_novak
35 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
36 * cast-type operations, and constant strings that must be defined as const char*
37 *
38 * Revision 1.4 2004/03/25 10:28:57 j_novak
39 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
40 *
41 * Revision 1.3 2003/10/03 15:58:43 j_novak
42 * Cleaning of some headers
43 *
44 * Revision 1.2 2002/10/16 14:36:32 j_novak
45 * Reorganization of #include instructions of standard C++, in order to
46 * use experimental version 3 of gcc.
47 *
48 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
49 * LORENE
50 *
51 * Revision 2.9 2001/05/16 14:51:36 phil
52 * correction calcul de kij
53 *
54 * Revision 2.8 2001/05/07 12:24:30 phil
55 * correction calcul de J
56 *
57 * Revision 2.7 2001/05/07 09:28:33 phil
58 * *** empty log message ***
59 *
60 * Revision 2.6 2001/02/12 15:36:58 phil
61 * ajout calcul de J a l infini
62 *
63 * Revision 2.5 2000/12/14 14:11:54 phil
64 * correction cl sur psi
65 *
66 * Revision 2.4 2000/12/14 12:41:54 phil
67 * on met les bases dans les sources
68 *
69 * Revision 2.3 2000/12/14 10:57:47 phil
70 * corections diverses et sans importances
71 *
72 * Revision 2.2 2000/12/14 10:45:00 phil
73 * ATTENTION : PASSAGE DE PHI A PSI
74 *
75 * Revision 2.1 2000/11/03 12:56:33 phil
76 * ajout de const
77 *
78 * Revision 2.0 2000/10/20 09:19:09 phil
79 * *** empty log message ***
80 *
81 *
82 * $Header: /cvsroot/Lorene/C++/Source/Bhole/bhole_pseudo_kerr.C,v 1.7 2014/10/13 08:52:40 j_novak Exp $
83 *
84 */
85
86//standard
87#include <cstdlib>
88#include <cmath>
89
90// Lorene
91#include "tenseur.h"
92#include "bhole.h"
93#include "proto.h"
94
95//Resolution pour le lapse pour 1 seul trou
96namespace Lorene {
97void Bhole::solve_lapse_seul (double relax) {
98
99 assert ((relax>0) && (relax<=1)) ;
100
101 cout << "Resolution LAPSE" << endl ;
102
103 // Pour la relaxation ...
104 Cmp lapse_old (n_auto()) ;
106 Tenseur kk (mp) ;
107 kk = 0 ;
108 Tenseur work(mp) ;
109 work.set_etat_qcq() ;
110 for (int i=0 ; i<3 ; i++) {
111 work.set() = auxi(i, i) ;
112 kk = kk + work ;
113 }
114
115 // La source
116 Cmp source
118 +pow(psi_auto(), 4.)*n_auto()*kk()) ;
119 source.std_base_scal() ;
120
121 // On resout pour N-1 :
122 Valeur limite (mp.get_mg()->get_angu()) ;
123 limite = -1 ;
124 limite.std_base_scal() ;
125
126 Cmp soluce (source.poisson_dirichlet(limite, 0)) ;
127 soluce = soluce + 1 ; // Permet de trouver N
128 soluce.raccord(1) ;
129
130 n_auto.set() = relax*soluce + (1-relax)*lapse_old ;
131}
132
133
134// Resolution sur Psi :
135void Bhole::solve_psi_seul (double relax) {
136
137 assert ((relax>0) && (relax<=1)) ;
138
139 cout << "Resolution PSI" << endl ;
140
141 Cmp psi_old (psi_auto()) ;
143 Tenseur kk (mp) ;
144 kk = 0 ;
145 Tenseur work(mp) ;
146 work.set_etat_qcq() ;
147 for (int i=0 ; i<3 ; i++) {
148 work.set() = auxi(i, i) ;
149 kk = kk + work ;
150 }
151
152 // La source :
153 Cmp source (-pow(psi_auto(), 5.)*kk()/8.) ;
154 source.std_base_scal() ;
155
156 // La condition limite de type neumann :
157 int np = mp.get_mg()->get_np(1) ;
158 int nt = mp.get_mg()->get_nt(1) ;
159 Valeur limite (mp.get_mg()->get_angu()) ;
160 limite = 1 ;
161 for (int k=0 ; k<np ; k++)
162 for (int j=0 ; j<nt ; j++)
163 limite.set(0, k, j, 0) = -0.5/rayon*psi_auto()(1, k, j, 0) ;
164 limite.std_base_scal() ;
165
166 Cmp soluce (source.poisson_neumann(limite, 0)) ;
167 soluce = soluce + 1 ;
168 soluce.raccord(1) ;
169
170 psi_auto.set() = relax*soluce + (1-relax)*psi_old ;
171
172}
173
174
175// Le shift. Processus iteratif pour cause de CL.
176void Bhole::solve_shift_seul (double precision, double relax) {
177
178 assert (precision > 0) ;
179 assert ((relax>0) && (relax<=1)) ;
180
181 cout << "resolution SHIFT" << endl ;
182
183 Tenseur shift_old (shift_auto) ;
184
187 source.set_std_base() ;
188
189 // On verifie si les 3 composantes ne sont pas nulles :
190 if (source.get_etat() == ETATQCQ) {
191 int indic = 0 ;
192 for (int i=0 ; i<3 ; i++)
193 if (source(i).get_etat() == ETATQCQ)
194 indic = 1 ;
195 if (indic ==0)
196 source.set_etat_zero() ;
197 }
198
199 // On filtre les hautes frequences pour raison de stabilite :
200 if (source.get_etat() == ETATQCQ)
201 for (int i=0 ; i<3 ; i++)
202 source.set(i).filtre(4) ;
203
204
205 // On determine les conditions limites en fonction du omega et du boost :
206 int np = mp.get_mg()->get_np(1) ;
207 int nt = mp.get_mg()->get_nt(1) ;
208
209 Mtbl x_mtbl (mp.get_mg()) ;
210 x_mtbl.set_etat_qcq() ;
211 Mtbl y_mtbl (mp.get_mg()) ;
212 y_mtbl.set_etat_qcq() ;
213 x_mtbl = mp.x ;
214 y_mtbl = mp.y ;
215
216 // Les bases pour les conditions limites :
217 Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
218
219 Valeur lim_x (mp.get_mg()->get_angu()) ;
220 lim_x = 1 ;
221 for (int k=0 ; k<np ; k++)
222 for (int j=0 ; j<nt ; j++)
223 lim_x.set(0, k, j, 0) = omega*y_mtbl(1, k, j, 0)-boost[0] ;
224 lim_x.base = *bases[0] ;
225
226 Valeur lim_y (mp.get_mg()->get_angu()) ;
227 lim_y = 1 ;
228 for (int k=0 ; k<np ; k++)
229 for (int j=0 ; j<nt ; j++)
230 lim_y.set(0, k, j, 0) = - omega*x_mtbl(1, k, j, 0)-boost[1] ;
231 lim_y.base = *bases[1] ;
232
233 Valeur lim_z (mp.get_mg()->get_angu()) ;
234 lim_z = 1 ;
235 for (int k=0 ; k<np ; k++)
236 for (int j=0 ; j<nt ; j++)
237 lim_z.set(0, k, j, 0) = -boost[2] ;
238 lim_z.base = *bases[2] ;
239
240 // On n'en a plus besoin
241 for (int i=0 ; i<3 ; i++)
242 delete bases[i] ;
243 delete [] bases ;
244
245 // On resout :
246 poisson_vect_frontiere(1./3., source, shift_auto, lim_x, lim_y,
247 lim_z, 0, precision, 20) ;
248
249 shift_auto = relax*shift_auto + (1-relax)*shift_old ;
250}
251
252
253// La regularisation si un seul trou noir :
255
256 // Le vecteur B (non tournant et non boostant)
257 Tenseur tbi (shift_auto) ;
258 for (int i=0 ; i<3 ; i++) {
259 tbi.set(i).va.coef_i() ;
260 tbi.set(i).va.set_etat_c_qcq() ;
261 }
262
263 for (int i=0 ; i<3 ; i++)
264 shift_auto(i).va.coef_i() ;
265
266 tbi.set(0) = *shift_auto(0).va.c - omega*shift_auto.get_mp()->y + boost[0];
267 tbi.set(1) = *shift_auto(1).va.c + omega*shift_auto.get_mp()->x + boost[1];
268 tbi.set(2) = *shift_auto(2).va.c + boost[2];
269 tbi.set_std_base() ;
270
271 // On evite soucis a l'infini (on a besoin que de la valeur sur horizon
272 tbi.set(0).annule(mp.get_mg()->get_nzone()-1) ;
273 tbi.set(1).annule(mp.get_mg()->get_nzone()-1) ;
274
275 Tenseur derive_r (mp, 1, CON, mp.get_bvect_cart()) ;
276 derive_r.set_etat_qcq() ;
277 for (int i=0 ; i<3 ; i++)
278 derive_r.set(i) = tbi(i).dsdr() ;
279
280 Valeur val_hor (mp.get_mg()) ;
281 Valeur fonction_radiale (mp.get_mg()) ;
282 Cmp enleve (mp) ;
283
284 double erreur = 0 ;
285 int nz = mp.get_mg()->get_nzone() ;
286 int np = mp.get_mg()->get_np(1) ;
287 int nt = mp.get_mg()->get_nt(1) ;
288 int nr = mp.get_mg()->get_nr(1) ;
289
290 double r_0 = mp.val_r(1, -1, 0, 0) ;
291 double r_1 = mp.val_r(1, 1, 0, 0) ;
292
293 for (int comp=0 ; comp<3 ; comp++) {
294 val_hor.annule_hard() ;
295 for (int k=0 ; k<np ; k++)
296 for (int j=0 ; j<nt ; j++)
297 for (int i=0 ; i<nr ; i++)
298 val_hor.set(1, k, j, i) = derive_r(comp)(1, k, j, 0) ;
299
300 fonction_radiale = pow(r_1-mp.r, 3.)* (mp.r-r_0)/pow(r_1-r_0, 3.) ;
301 fonction_radiale.annule(0) ;
302 fonction_radiale.annule(2, nz-1) ;
303
304 enleve = fonction_radiale*val_hor ;
305 enleve.va.base = shift_auto(comp).va.base ;
306
307 // Ca devrai annuler la derivee de B sur H et donc rendre K regulier
308 Cmp copie (shift_auto(comp)) ;
309 shift_auto.set(comp) = shift_auto(comp)-enleve ;
310
311 assert (shift_auto(comp).check_dzpuis(0)) ;
312
313 // On regarde l'intensite de la correction si non nul !
314 Tbl norm (norme(shift_auto(comp))) ;
315 if (norm(1) > 1e-5) {
316 Tbl diff (diffrelmax (copie, shift_auto(comp))) ;
317 if (erreur<diff(1))
318 erreur = diff(1) ;
319 }
320 }
321 regul = erreur ;
322}
323
324// On calcul Kij sachant que la regulatisation sur le shift doit
325// etre faite avant.
327
329
330 Cmp lapse_non_sing (division_xpun(n_auto(), 0)) ;
331 Cmp auxi (mp) ;
332
334 for (int i=0 ; i<3 ; i++)
335 for (int j=i ; j<3 ; j++) {
336 auxi = taij_auto(i, j) ;
337 auxi = division_xpun (auxi, 0) ;
338 tkij_auto.set(i, j) = auxi/lapse_non_sing/2. ;
339 tkij_auto.set(i, j).raccord(1) ;
340 }
341}
342
343// Calcul masse ADM via valeur asymptotiques
344double Bhole::masse_adm_seul () const {
345
346 Cmp integrant (psi_auto().dsdr()) ;
347 double masse = mp.integrale_surface_infini (integrant) ;
348 masse /= -2*M_PI ;
349 return masse ;
350}
351
353
354 Cmp integrant (n_auto().dsdr()) ;
355 double masse = mp.integrale_surface_infini (integrant) ;
356 masse /= 4*M_PI ;
357 return masse ;
358}
359
360// Calcul du moment angulaire via integrale a l infini
361// Non valable si le boost != 0 ;
362
364
365 // On verifie si le boost est bien nul :
366 double indic = 0 ;
367 for (int i=0 ; i<3 ; i++)
368 if (boost[i] != 0)
369 indic = 1 ;
370 if (indic == 1) {
371 cout << "Calcul du moment non valable pour un boost != 0" << endl ;
372 abort() ;
373 }
374
375 if (omega == 0)
376 return 0 ;
377 else {
378 Tenseur vecteur_un (mp, 1, CON, mp.get_bvect_cart()) ;
379 vecteur_un.set_etat_qcq() ;
380 for (int i=0 ; i<3 ; i++)
381 vecteur_un.set(i) = tkij_auto(0, i) ;
382 vecteur_un.change_triad (mp.get_bvect_spher()) ;
383 Cmp integrant_un (vecteur_un(0)) ;
384 integrant_un.mult_r_zec() ;
385
386 Tenseur vecteur_deux (mp, 1, CON, mp.get_bvect_cart()) ;
387 vecteur_deux.set_etat_qcq() ;
388 for (int i=0 ; i<3 ; i++)
389 vecteur_deux.set(i) = tkij_auto(1, i) ;
390 vecteur_deux.change_triad (mp.get_bvect_spher()) ;
391 Cmp integrant_deux (vecteur_deux(0)) ;
392 integrant_deux.mult_r_zec() ;
393
394 // On multiplie par y et x :
395 integrant_un.va = integrant_un.va.mult_st() ;
396 integrant_un.va = integrant_un.va.mult_sp() ;
397
398 integrant_deux.va = integrant_deux.va.mult_st() ;
399 integrant_deux.va = integrant_deux.va.mult_cp() ;
400
401 double moment = mp.integrale_surface_infini (-integrant_un+integrant_deux) ;
402 moment /= 8*M_PI ;
403 return moment ;
404 }
405}
406
407// Calcul du moment angulaire via integrale sur l horizon
408// Non valable si le boost != 0 ;
409
411
412 // On verifie si le boost est bien nul :
413 double indic = 0 ;
414 for (int i=0 ; i<3 ; i++)
415 if (boost[i] != 0)
416 indic = 1 ;
417 if (indic == 1) {
418 cout << "Calcul du moment non valable pour un boost != 0" << endl ;
419 abort() ;
420 }
421
422 if (omega == 0)
423 return 0 ;
424 else {
425 // Integrale sur l horizon :
426 Cmp xa (mp) ;
427 xa = mp.xa ;
428 xa.std_base_scal() ;
429
430 Cmp ya (mp) ;
431 ya = mp.ya ;
432 ya.std_base_scal() ;
433
434 Tenseur vecteur (mp, 1, CON, mp.get_bvect_cart()) ;
435 vecteur.set_etat_qcq() ;
436 for (int i=0 ; i<3 ; i++)
437 vecteur.set(i) = (-ya*tkij_auto(0, i)+xa * tkij_auto(1, i)) ;
438 vecteur.set_std_base() ;
439 vecteur.annule(mp.get_mg()->get_nzone()-1) ;
440 vecteur.change_triad (mp.get_bvect_spher()) ;
441
442 Cmp integrant (pow(psi_auto(), 6)*vecteur(0)) ;
443 integrant.std_base_scal() ;
444 double moment = mp.integrale_surface (integrant, rayon)/8/M_PI ;
445 return moment ;
446 }
447}
448
449}
Bases of the spectral expansions.
Definition base_val.h:322
Tenseur psi_auto
Part of generated by the hole.
Definition bhole.h:290
Tenseur shift_auto
Part of generated by the hole.
Definition bhole.h:297
Tenseur_sym tkij_auto
Auto .
Definition bhole.h:307
double omega
Angular velocity in LORENE's units.
Definition bhole.h:275
double regul
Intensity of the correction on the shift vector.
Definition bhole.h:284
Tenseur_sym taij_auto
Part of generated by the hole.
Definition bhole.h:299
void solve_shift_seul(double precis, double relax)
Solves the equation for ~:
void regularise_seul()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon.
double masse_adm_seul() const
Calculates the ADM mass of the black hole using : .
Tenseur n_auto
Part of N generated by the hole.
Definition bhole.h:286
double * boost
Lapse on the horizon.
Definition bhole.h:283
void solve_psi_seul(double relax)
Solves the equation for ~:
double moment_seul_hor() const
Calculates the angular momentum of the black hole using the formula on the horizon : where and H de...
double rayon
Radius of the horizon in LORENE's units.
Definition bhole.h:274
Map_af & mp
Affine mapping.
Definition bhole.h:273
void fait_tkij()
Calculates the total .
double masse_komar_seul() const
Calculates the Komar mass of the black hole using : .
void solve_lapse_seul(double relax)
Solves the equation for N ~:
double moment_seul_inf() const
Calculates the angular momentum of the black hole using the formula at infinity : where .
void fait_taij_auto()
Calculates the part of generated by shift_auto .
Definition bhole.C:379
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void mult_r_zec()
Multiplication by r in the external compactified domain (ZEC)
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Cmp poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Cmp::poisson() .
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition cmp_manip.C:74
Cmp poisson_neumann(const Valeur &, int) const
Idem as Cmp::poisson_dirichlet , the boundary condition being on the radial derivative of the solutio...
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord ya
Absolute y coordinate.
Definition map.h:731
Coord r
r coordinate centered on the grid
Definition map.h:718
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Coord x
x coordinate centered on the grid
Definition map.h:726
Coord xa
Absolute x coordinate.
Definition map.h:730
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
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
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Multi-domain array.
Definition mtbl.h:118
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Basic array class.
Definition tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Cmp ** c
The components.
Definition tenseur.h:322
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition tenseur.C:900
const Map * get_mp() const
Returns pointer on the mapping.
Definition tenseur.h:699
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tenseur.C:668
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Values and coefficients of a (real-value) function.
Definition valeur.h:287
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition valeur.C:744
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
void coef_i() const
Computes the physical value of *this.
const Valeur & mult_st() const
Returns applied to *this.
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
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
const Valeur & mult_cp() const
Returns applied to *this.
const Valeur & mult_sp() const
Returns applied to *this.
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Lorene prototypes.
Definition app_hor.h:64