LORENE
binhor_coal.C
1/*
2 * Copyright (c) 2004-2005 Francois Limousin
3 * Jose-Luis Jaramillo
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24char binhor_coal_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $" ;
25
26/*
27 * $Id: binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $
28 * $Log: binhor_coal.C,v $
29 * Revision 1.15 2014/10/13 08:52:42 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.14 2014/10/06 15:13:01 j_novak
33 * Modified #include directives to use c++ syntax.
34 *
35 * Revision 1.13 2007/04/13 15:28:55 f_limousin
36 * Lots of improvements, generalisation to an arbitrary state of
37 * rotation, implementation of the spatial metric given by Samaya.
38 *
39 * Revision 1.12 2006/08/01 14:37:19 f_limousin
40 * New version
41 *
42 * Revision 1.11 2006/06/28 13:36:09 f_limousin
43 * Convergence to a given irreductible mass
44 *
45 * Revision 1.10 2006/05/24 16:56:37 f_limousin
46 * Many small modifs.
47 *
48 * Revision 1.9 2005/09/13 18:33:15 f_limousin
49 * New function vv_bound_cart_bin(double) for computing binaries with
50 * berlin condition for the shift vector.
51 * Suppress all the symy and asymy in the importations.
52 *
53 * Revision 1.8 2005/07/11 08:21:57 f_limousin
54 * Implementation of a new boundary condition for the lapse in the binary
55 * case : boundary_nn_Dir_lapl().
56 *
57 * Revision 1.7 2005/03/10 17:21:52 f_limousin
58 * Add the Berlin boundary condition for the shift.
59 * Some changes to avoid warnings.
60 *
61 * Revision 1.6 2005/03/10 17:09:05 f_limousin
62 * Display the logarithm of viriel and convergence.
63 *
64 * Revision 1.5 2005/03/10 16:57:00 f_limousin
65 * Improve the convergence of the code coal_bh.
66 *
67 * Revision 1.4 2005/02/24 17:24:26 f_limousin
68 * The boundary conditions for psi, N and beta are now parameters in
69 * par_init.d and par_coal.d.
70 *
71 * Revision 1.3 2005/02/07 10:43:36 f_limousin
72 * Add the printing of the regularisation of the shift in the case N=0
73 * on the horizon.
74 *
75 * Revision 1.2 2004/12/31 15:40:21 f_limousin
76 * Improve the initialisation of several quantities in set_statiques().
77 *
78 * Revision 1.1 2004/12/29 16:11:19 f_limousin
79 * First version
80 *
81 *
82 * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $
83 *
84 */
85
86//standard
87#include <cstdlib>
88
89// Lorene
90#include "tensor.h"
91#include "isol_hor.h"
92#include "graphique.h"
93
94
95namespace Lorene {
96void Bin_hor::set_statiques (double precis, double relax, int bound_nn,
97 double lim_nn, int bound_psi) {
98
99 int nz = hole1.mp.get_mg()->get_nzone() ;
100
101 set_omega(0) ;
104 init_bin_hor() ;
106
107 int indic = 1 ;
108 int conte = 0 ;
109
110 cout << "Static black holes : " << endl ;
111 while (indic == 1) {
112 Scalar lapse_un_old (hole1.n_auto) ;
113
114 solve_psi (precis, relax, bound_psi) ;
115 solve_lapse (precis, relax, bound_nn, lim_nn) ;
116
117 // des_profile(hole1.nn(), 0, 20, M_PI/2, M_PI) ;
118
119 double erreur = 0 ;
120 Tbl diff (diffrelmax (lapse_un_old, hole1.n_auto)) ;
121 for (int i=1 ; i<nz ; i++)
122 if (diff(i) > erreur)
123 erreur = diff(i) ;
124
125 cout << "Step : " << conte << " Difference : " << erreur << endl ;
126
127 if (erreur < precis)
128 indic = -1 ;
129 conte ++ ;
130 }
131}
132
133double Bin_hor::coal (double angu_vel, double relax, int nb_ome,
134 int nb_it, int bound_nn, double lim_nn,
135 int bound_psi, int bound_beta, double omega_eff,
136 double alpha,
137 ostream& fich_iteration, ostream& fich_correction,
138 ostream& fich_viriel, ostream& fich_kss,
139 int step, int search_mass, double mass_irr,
140 const int sortie) {
141
142 int nz = hole1.mp.get_mg()->get_nzone() ;
143
144 double precis = 1e-7 ;
145
146 // LOOP INCREASING OMEGA :
147 cout << "OMEGA INCREASED STEP BY STEP." << endl ;
148 double homme = get_omega() ;
149 double inc_homme = (angu_vel - homme)/nb_ome ;
150 for (int pas = 0 ; pas <nb_ome ; pas ++) {
151
152 bool verif = false ;
153 if (omega_eff == alpha*homme ) verif = true ;
154
155 homme += inc_homme ;
156 set_omega (homme) ;
157 if (verif)
158 omega_eff = alpha*homme ;
159 Scalar beta_un_old (hole1.beta_auto(1)) ;
160
161 solve_shift (precis, relax, bound_beta, omega_eff) ;
163
164 solve_psi (precis, relax, bound_psi) ;
165 solve_lapse (precis, relax, bound_nn, lim_nn) ;
166
167 // Convergence to the given irreductible mass
168 if (search_mass == 1 && step >= 30) {
169 double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
170 sqrt(hole2.area_hor()/16/M_PI) ;
171 double error_m = (mass_irr - mass_area) / mass_irr ;
172 double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
173 hole1.mp.homothetie_interne(scaling_r) ;
174 hole1.radius = hole1.radius *scaling_r ;
175 hole2.mp.homothetie_interne(scaling_r) ;
176 hole2.radius = hole2.radius *scaling_r ;
177
178 // Update of the different metrics (another possibility would
179 // be to set all derived quantities to 0x0, especially
180 // the connection p_connect
185
186 }
187
188 cout << "Angular momentum computed at the horizon : " << ang_mom_hor()
189 << endl ;
190
191 double erreur = 0 ;
192 Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
193 for (int i=1 ; i<nz ; i++)
194 if (diff(i) > erreur)
195 erreur = diff(i) ;
196
197 // Saving ok K_{ij}s^is^j
198 // -----------------------
199
200 Scalar kkss (contract(hole1.get_k_dd(), 0, 1,
202 hole1.get_gam().radial_vect(), 0, 1)) ;
203 double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
204 double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
205 int nnp = hole1.mp.get_mg()->get_np(1) ;
206 int nnt = hole2.mp.get_mg()->get_nt(1) ;
207 for (int k=0 ; k<nnp ; k++)
208 for (int j=0 ; j<nnt ; j++){
209 if (kkss.val_grid_point(1, k, j, 0) > max_kss)
210 max_kss = kkss.val_grid_point(1, k, j, 0) ;
211 if (kkss.val_grid_point(1, k, j, 0) < min_kss)
212 min_kss = kkss.val_grid_point(1, k, j, 0) ;
213 }
214
215 if (sortie != 0) {
216 fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
217 fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
218 // fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
219 fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
220 fich_kss << step << " " << max_kss << " " << min_kss << endl ;
221 }
222
223 cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
224 step ++ ;
225 }
226
227 // LOOP WITH FIXED OMEGA :
228
229 if (nb_it !=0)
230 cout << "OMEGA FIXED" << endl ;
231 double erreur ;
232
233 for (int pas = 0 ; pas <nb_it ; pas ++) {
234
235 Scalar beta_un_old (hole1.beta_auto(1)) ;
236
237 solve_shift (precis, relax, bound_beta, omega_eff) ;
239
240 solve_psi (precis, relax, bound_psi) ;
241 solve_lapse (precis, relax, bound_nn, lim_nn) ;
242
243 // Convergence to the given irreductible mass
244 if (search_mass == 1 && step >= 30) {
245 double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
246 sqrt(hole2.area_hor()/16/M_PI) ;
247 double error_m = (mass_irr - mass_area) / mass_irr ;
248 double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
249
250 hole1.mp.homothetie_interne(scaling_r) ;
251 hole1.radius = hole1.radius *scaling_r ;
252 hole2.mp.homothetie_interne(scaling_r) ;
253 hole2.radius = hole2.radius *scaling_r ;
254 }
255
256 erreur = 0 ;
257 Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
258 for (int i=1 ; i<nz ; i++)
259 if (diff(i) > erreur)
260 erreur = diff(i) ;
261
262 // Saving ok K_{ij}s^is^j
263 // -----------------------
264
265 Scalar kkss (contract(hole1.get_k_dd(), 0, 1,
267 hole1.get_gam().radial_vect(), 0, 1)) ;
268 double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
269 double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
270 int nnp = hole1.mp.get_mg()->get_np(1) ;
271 int nnt = hole2.mp.get_mg()->get_nt(1) ;
272 for (int k=0 ; k<nnp ; k++)
273 for (int j=0 ; j<nnt ; j++){
274 if (kkss.val_grid_point(1, k, j, 0) > max_kss)
275 max_kss = kkss.val_grid_point(1, k, j, 0) ;
276 if (kkss.val_grid_point(1, k, j, 0) < min_kss)
277 min_kss = kkss.val_grid_point(1, k, j, 0) ;
278 }
279
280
281 if (sortie != 0) {
282 fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
283 fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
284 // fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
285 fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
286 fich_kss << step << " " << max_kss << " " << min_kss << endl ;
287 }
288
289 cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
290 step ++ ;
291 }
292
293 if (nb_it != 0){
294 fich_iteration << "#----------------------------" << endl ;
295 fich_correction << "#-----------------------------" << endl ;
296 fich_viriel << "#------------------------------" << endl ;
297 }
298
299 return viriel() ;
300}
301}
double viriel() const
Computes the viriel error, that is the difference between the ADM and the Komar masses,...
void solve_psi(double precis, double relax, int bound_psi)
Solves the equation for the conformal factor : The fields are the total values excpet those with subs...
double coal(double ang_vel, double relax, int nb_om, int nb_it, int bound_nn, double lim_nn, int bound_psi, int bound_beta, double omega_eff, double alpha, ostream &fich_iteration, ostream &fich_correction, ostream &fich_viriel, ostream &fich_kss, int step, int search_mass, double mass_irr, const int sortie=0)
Solves the equation for a particular angular velocity, the systeme being initialized to Misner-Lindqu...
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
void init_bin_hor()
Initialisation of the system.
Definition bin_hor.C:161
void solve_shift(double precis, double relax, int bound_beta, double omega_eff)
Solves the equation for the shift, using the Oohara-Nakarmure scheme : The fields are the total value...
void set_statiques(double precis, double relax, int bound_nn, double lim_nn, int bound_psi)
Initialize the systeme to Misner Lindquist solution, that is solving for N and in the case .
Definition binhor_coal.C:96
void solve_lapse(double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for the lapse : The fields are the total values except those with subscript ,...
void set_omega(double ome)
Sets the orbital velocity to ome.
Definition isol_hor.h:1409
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition binhor_kij.C:88
double get_omega() const
Returns the angular velocity.
Definition isol_hor.h:1422
void homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
Definition map_af.C:614
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition metric.C:362
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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
double radius
Radius of the horizon in LORENE's units.
Definition isol_hor.h:906
const Metric & get_gam() const
metric
Definition single_hor.C:339
Vector beta_auto
Shift function .
Definition isol_hor.h:944
const Sym_tensor & get_k_dd() const
k_dd
Definition single_hor.C:348
Metric_flat ff
3 metric flat
Definition isol_hor.h:980
Map_af & mp
Affine mapping.
Definition isol_hor.h:900
void init_met_trK()
Sets the 3-metric tilde to the flat metric and gamt_point, trK and trK_point to zero.
Definition single_hor.C:536
double area_hor() const
Area of the horizon.
double regul
Intensity of the correction on the shift vector.
Definition isol_hor.h:912
double omega_hor() const
Orbital velocity
Scalar n_auto
Lapse function .
Definition isol_hor.h:915
Metric tgam
3 metric tilde
Definition isol_hor.h:977
Basic array class.
Definition tbl.h:161
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
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 contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64