LORENE
bhole_equations_bin.C
1/*
2 * Copyright (c) 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_equations_bin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $" ;
24
25/*
26 * $Id: bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $
27 * $Log: bhole_equations_bin.C,v $
28 * Revision 1.5 2014/10/13 08:52:40 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.4 2014/10/06 15:12:58 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.3 2006/04/27 09:12:32 p_grandclement
35 * First try at irrotational black holes
36 *
37 * Revision 1.2 2002/10/16 14:36:33 j_novak
38 * Reorganization of #include instructions of standard C++, in order to
39 * use experimental version 3 of gcc.
40 *
41 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
42 * LORENE
43 *
44 * Revision 2.6 2001/05/07 09:11:56 phil
45 * *** empty log message ***
46 *
47 * Revision 2.5 2001/04/30 09:30:50 phil
48 * filtre pour poisson vectoriel
49 *
50 * Revision 2.4 2001/04/26 12:04:34 phil
51 * *** empty log message ***
52 *
53 * Revision 2.3 2001/04/25 15:08:48 phil
54 * vire fait_tkij
55 *
56 * Revision 2.2 2001/04/02 12:16:11 phil
57 * *** empty log message ***
58 *
59 * Revision 2.1 2001/03/22 10:41:36 phil
60 * modification prototypage
61 *
62 * Revision 2.0 2001/03/01 08:18:04 phil
63 * *** empty log message ***
64 *
65 *
66 * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_equations_bin.C,v 1.5 2014/10/13 08:52:40 j_novak Exp $
67 *
68 */
69
70//standard
71#include <cstdlib>
72#include <cmath>
73
74// Lorene
75#include "nbr_spx.h"
76#include "tenseur.h"
77#include "bhole.h"
78#include "proto.h"
79#include "utilitaires.h"
80#include "graphique.h"
81
82// Resolution pour le lapse
83namespace Lorene {
84void Bhole_binaire::solve_lapse (double precision, double relax) {
85
86 assert ((relax >0) && (relax<=1)) ;
87
88 cout << "-----------------------------------------------" << endl ;
89 cout << "Resolution LAPSE" << endl ;
90
91 Tenseur lapse_un_old (hole1.n_auto) ;
92 Tenseur lapse_deux_old (hole2.n_auto) ;
93
95 Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
96
98 Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
99
100 // Les sources
101
102 Cmp source_un
104 +hole1.n_tot()*pow(hole1.psi_tot(), 4.)*kk_un) ;
105 source_un.std_base_scal() ;
106
107 Cmp source_deux
109 +hole2.n_tot()*pow(hole2.psi_tot(), 4.)*kk_deux) ;
110 source_deux.std_base_scal() ;
111
112 //On resout
113 hole1.n_auto.set() = hole1.n_auto() - 1./2. ;
114 hole2.n_auto.set() = hole2.n_auto() - 1./2. ;
115
116 dirichlet_binaire (source_un, source_deux, -1., -1., hole1.n_auto.set(),
117 hole2.n_auto.set(), 0, precision) ;
118
119 hole1.n_auto.set() = hole1.n_auto() + 1./2. ;
120 hole2.n_auto.set() = hole2.n_auto() + 1./2. ;
121
122 hole1.n_auto.set().raccord(1) ;
123 hole2.n_auto.set().raccord(1) ;
124
125 // La relaxation :
126 hole1.n_auto.set() = relax*hole1.n_auto() + (1-relax)*lapse_un_old() ;
127 hole2.n_auto.set() = relax*hole2.n_auto() + (1-relax)*lapse_deux_old() ;
128
131}
132
133//Resolution sur Psi
134void Bhole_binaire::solve_psi (double precision, double relax) {
135
136 assert ((relax>0) && (relax<=1)) ;
137
138 cout << "-----------------------------------------------" << endl ;
139 cout << "Resolution PSI" << endl ;
140
141 Tenseur psi_un_old (hole1.psi_auto) ;
142 Tenseur psi_deux_old (hole2.psi_auto) ;
143
145 Cmp kk_un (auxi_un(0, 0)+auxi_un(1, 1)+auxi_un(2, 2)) ;
146
148 Cmp kk_deux (auxi_deux(0, 0)+auxi_deux(1, 1)+auxi_deux(2, 2)) ;
149
150 // Les sources
151 Cmp source_un (-pow(hole1.psi_tot(), 5.)*kk_un/8.) ;
152 source_un.std_base_scal() ;
153
154 Cmp source_deux (-pow(hole2.psi_tot(), 5.)*kk_deux/8.) ;
155 source_deux.std_base_scal() ;
156
157 // Les valeurs limites :
158 int np_un = hole1.mp.get_mg()->get_np(1) ;
159 int nt_un = hole1.mp.get_mg()->get_nt(1) ;
160 Valeur lim_un (hole1.mp.get_mg()->get_angu()) ;
161 lim_un = 1 ;
162 for (int k=0 ; k<np_un ; k++)
163 for (int j=0 ; j<nt_un ; j++)
164 lim_un.set(0, k, j, 0) = -0.5/hole1.rayon*hole1.psi_tot()(1, k, j, 0) ;
165 lim_un.std_base_scal() ;
166
167 int np_deux = hole2.mp.get_mg()->get_np(1) ;
168 int nt_deux = hole2.mp.get_mg()->get_nt(1) ;
169 Valeur lim_deux (hole2.mp.get_mg()->get_angu()) ;
170 lim_deux = 1 ;
171 for (int k=0 ; k<np_deux ; k++)
172 for (int j=0 ; j<nt_deux ; j++)
173 lim_deux.set(0, k, j, 0) = -0.5/hole2.rayon*hole2.psi_tot()(1, k, j, 0) ;
174 lim_deux.std_base_scal() ;
175
176 //On resout
177 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
178 hole1.psi_auto.set(), hole2.psi_auto.set(), 0, precision) ;
179
180 hole1.psi_auto.set() = hole1.psi_auto() + 1./2. ;
181 hole2.psi_auto.set() = hole2.psi_auto() + 1./2. ;
182
185
186 // La relaxation :
187 hole1.psi_auto.set() = relax*hole1.psi_auto() + (1-relax)*psi_un_old() ;
188 hole2.psi_auto.set() = relax*hole2.psi_auto() + (1-relax)*psi_deux_old() ;
189
192}
193
194
195// Resolution sur shift a omega bloque.
196void Bhole_binaire::solve_shift (double precision, double relax) {
197
198 cout << "------------------------------------------------" << endl ;
199 cout << "Resolution shift : Omega = " << omega << endl ;
200
201 // On determine les sources
202 Tenseur source_un (
205 if (source_un.get_etat() == ETATZERO) {
206 source_un.set_etat_qcq() ;
207 for (int i=0 ; i<3 ; i++)
208 source_un.set(i).set_etat_zero() ;
209 }
210 source_un.set_std_base() ;
211
212 Tenseur source_deux (
215 if (source_deux.get_etat() == ETATZERO) {
216 source_deux.set_etat_qcq() ;
217 for (int i=0 ; i<3 ; i++)
218 source_deux.set(i).set_etat_zero() ;
219 }
220 source_deux.set_std_base() ;
221
222 // On filtre les hautes frequences.
223 for (int i=0 ; i<3 ; i++) {
224 if (source_un(i).get_etat() != ETATZERO)
225 source_un.set(i).filtre(4) ;
226 if (source_deux(i).get_etat() != ETATZERO)
227 source_deux.set(i).filtre(4) ;
228 }
229
230 // Les alignemenents pour le signe des CL.
231 double orientation_un = hole1.mp.get_rot_phi() ;
232 assert ((orientation_un==0) || (orientation_un == M_PI)) ;
233
234 double orientation_deux = hole2.mp.get_rot_phi() ;
235 assert ((orientation_deux==0) || (orientation_deux == M_PI)) ;
236
237 int aligne_un = (orientation_un == 0) ? 1 : -1 ;
238 int aligne_deux = (orientation_deux == 0) ? 1 : -1 ;
239
240 // On regarde si toutes les composantes sont nulles :
241 int ind1 = 0 ;
242 int ind2 = 0 ;
243 for (int i=0 ; i<3 ; i++) {
244 if (source_un(i).get_etat() == ETATQCQ)
245 ind1 = 1 ;
246 if (source_deux(i).get_etat() == ETATQCQ)
247 ind2 = 1 ;
248 }
249
250 if (ind1==0)
251 source_un.set_etat_zero() ;
252 if (ind2==0)
253 source_deux.set_etat_zero() ;
254
255 // On determine les Cl en fonction de omega :
256 int np_un = hole1.mp.get_mg()->get_np (1) ;
257 int nt_un = hole1.mp.get_mg()->get_nt (1) ;
258
259 Mtbl xa_mtbl_un (source_un.get_mp()->get_mg()) ;
260 xa_mtbl_un.set_etat_qcq() ;
261 Mtbl ya_mtbl_un (source_un.get_mp()->get_mg()) ;
262 ya_mtbl_un.set_etat_qcq() ;
263
264 xa_mtbl_un = source_un.get_mp()->xa ;
265 ya_mtbl_un = source_un.get_mp()->ya ;
266
267 Mtbl x_mtbl_un (source_un.get_mp()->get_mg()) ;
268 x_mtbl_un.set_etat_qcq() ;
269 Mtbl y_mtbl_un (source_un.get_mp()->get_mg()) ;
270 y_mtbl_un.set_etat_qcq() ;
271
272 x_mtbl_un = source_un.get_mp()->x ;
273 y_mtbl_un = source_un.get_mp()->y ;
274
275 // Les bases
276 Base_val** bases_un = hole1.mp.get_mg()->std_base_vect_cart() ;
277 Base_val** bases_deux = hole2.mp.get_mg()->std_base_vect_cart() ;
278
279 Valeur lim_x_un (*hole1.mp.get_mg()->get_angu()) ;
280 lim_x_un = 1 ; // Juste pour affecter dans espace des configs ;
281 lim_x_un.set_etat_c_qcq() ;
282 for (int k=0 ; k<np_un ; k++)
283 for (int j=0 ; j<nt_un ; j++)
284 lim_x_un.set(0, k, j, 0) = aligne_un*omega*ya_mtbl_un(0, 0, 0, 0) + aligne_un*hole1.omega_local*y_mtbl_un(1,k,j,0) ;
285 lim_x_un.base = *bases_un[0] ;
286
287 Valeur lim_y_un (*hole1.mp.get_mg()->get_angu()) ;
288 lim_y_un = 1 ; // Juste pour affecter dans espace des configs ;
289 lim_y_un.set_etat_c_qcq() ;
290 for (int k=0 ; k<np_un ; k++)
291 for (int j=0 ; j<nt_un ; j++)
292 lim_y_un.set(0, k, j, 0) = -aligne_un*omega*xa_mtbl_un(0, 0, 0, 0) - aligne_un*hole1.omega_local*x_mtbl_un(1,k,j,0) ;;
293 lim_y_un.base = *bases_un[1] ;
294
295 Valeur lim_z_un (*hole1.mp.get_mg()->get_angu()) ;
296 lim_z_un = 1 ;
297 for (int k=0 ; k<np_un ; k++)
298 for (int j=0 ; j<nt_un ; j++)
299 lim_z_un.set(0, k, j, 0) = 0 ;
300 lim_z_un.base = *bases_un[2] ;
301
302 // On determine les Cl en fonction de omega :
303 int np_deux = hole2.mp.get_mg()->get_np (1) ;
304 int nt_deux = hole2.mp.get_mg()->get_nt (1) ;
305
306 Mtbl xa_mtbl_deux (source_deux.get_mp()->get_mg()) ;
307 xa_mtbl_deux.set_etat_qcq() ;
308 Mtbl ya_mtbl_deux (source_deux.get_mp()->get_mg()) ;
309 ya_mtbl_deux.set_etat_qcq() ;
310
311 xa_mtbl_deux = source_deux.get_mp()->xa ;
312 ya_mtbl_deux = source_deux.get_mp()->ya ;
313
314 Mtbl x_mtbl_deux (source_deux.get_mp()->get_mg()) ;
315 x_mtbl_deux.set_etat_qcq() ;
316 Mtbl y_mtbl_deux (source_deux.get_mp()->get_mg()) ;
317 y_mtbl_deux.set_etat_qcq() ;
318
319 x_mtbl_deux = source_deux.get_mp()->x ;
320 y_mtbl_deux = source_deux.get_mp()->y ;
321
322 Valeur lim_x_deux (*hole2.mp.get_mg()->get_angu()) ;
323 lim_x_deux = 1 ; // Juste pour affecter dans espace des configs ;
324 lim_x_deux.set_etat_c_qcq() ;
325 for (int k=0 ; k<np_deux ; k++)
326 for (int j=0 ; j<nt_deux ; j++)
327 lim_x_deux.set(0, k, j, 0) = aligne_deux*omega*ya_mtbl_deux(0, 0, 0, 0) + aligne_deux*hole2.omega_local*y_mtbl_deux(1,k,j,0) ;
328 lim_x_deux.base = *bases_deux[0] ;
329
330 Valeur lim_y_deux (*hole2.mp.get_mg()->get_angu()) ;
331 lim_y_deux = 1 ; // Juste pour affecter dans espace des configs ;
332 lim_y_deux.set_etat_c_qcq() ;
333 for (int k=0 ; k<np_deux ; k++)
334 for (int j=0 ; j<nt_deux ; j++)
335 lim_y_deux.set(0, k, j, 0) = -aligne_deux*omega*xa_mtbl_deux(0, 0, 0, 0) - aligne_deux*hole2.omega_local*x_mtbl_deux(1,k,j,0) ;
336 lim_y_deux.base = *bases_deux[1] ;
337
338 Valeur lim_z_deux (*hole2.mp.get_mg()->get_angu()) ;
339 lim_z_deux = 1 ;
340 for (int k=0 ; k<np_deux ; k++)
341 for (int j=0 ; j<nt_deux ; j++)
342 lim_z_deux.set(0, k, j, 0) = 0 ;
343 lim_z_deux.base = *bases_deux[2] ;
344
345 for (int i=0 ; i<3 ; i++) {
346 delete bases_un[i] ;
347 delete bases_deux[i] ;
348 }
349 delete [] bases_un ;
350 delete [] bases_deux ;
351
352 // On resout le truc :
353 Tenseur shift_un_old (hole1.shift_auto) ;
354 Tenseur shift_deux_old (hole2.shift_auto) ;
355
356 poisson_vect_binaire (1./3., source_un, source_deux,
357 lim_x_un, lim_y_un, lim_z_un,
358 lim_x_deux, lim_y_deux, lim_z_deux,
359 hole1.shift_auto, hole2.shift_auto, 0, precision) ;
360
361 for (int i=0 ; i<3 ; i++) {
364 }
365
366 // On regularise les shift.
368 (1-relax)*shift_un_old ;
370 (1-relax)*shift_deux_old ;
371
372 double diff_un = regle (hole1.shift_auto, hole2.shift_auto, omega, hole1.omega_local) ;
373 double diff_deux = regle (hole2.shift_auto, hole1.shift_auto, omega, hole2.omega_local) ;
374 hole1.regul = diff_un ;
375 hole2.regul = diff_deux ;
376
377 cout << "Difference relatives : " << diff_un << " " << diff_deux << endl ;
378}
379}
Bases of the spectral expansions.
Definition base_val.h:322
double omega
Position of the axis of rotation.
Definition bhole.h:769
Bhole hole1
Black hole one.
Definition bhole.h:762
void solve_shift(double precis, double relax)
Solves the equation for the shift, using the Oohara-Nakarmure scheme~:
Bhole hole2
Black hole two.
Definition bhole.h:763
void solve_psi(double precis, double relax)
Solves the equation for the conformal factor~:
void solve_lapse(double precis, double relax)
Solves the equation for the lapse~:
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
double omega_local
local angular velocity
Definition bhole.h:276
Tenseur psi_tot
Total .
Definition bhole.h:292
Tenseur_sym tkij_auto
Auto .
Definition bhole.h:307
double regul
Intensity of the correction on the shift vector.
Definition bhole.h:284
Tenseur grad_n_tot
Total gradient of N .
Definition bhole.h:294
Tenseur n_auto
Part of N generated by the hole.
Definition bhole.h:286
double rayon
Radius of the horizon in LORENE's units.
Definition bhole.h:274
Map_af & mp
Affine mapping.
Definition bhole.h:273
Tenseur_sym taij_tot
Total , which must be zero on the horizon of the regularisation on the shift has been done.
Definition bhole.h:305
Tenseur n_tot
Total N .
Definition bhole.h:288
void fait_psi_comp(const Bhole &comp)
Imports the part of due to the companion hole comp .
Definition bhole.C:280
Tenseur_sym tkij_tot
Total .
Definition bhole.h:308
void fait_n_comp(const Bhole &comp)
Imports the part of N due to the companion hole comp .
Definition bhole.C:254
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition cmp_manip.C:74
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord ya
Absolute y coordinate.
Definition map.h:731
Coord x
x coordinate centered on the grid
Definition map.h:726
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
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
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
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
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
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
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
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