LORENE
bhole_coal.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 
23 
24 
25 /*
26  * $Id: bhole_coal.C,v 1.9 2016/12/05 16:17:45 j_novak Exp $
27  * $Log: bhole_coal.C,v $
28  * Revision 1.9 2016/12/05 16:17:45 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.8 2014/10/13 08:52:40 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.7 2014/10/06 15:12:58 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.6 2005/08/31 09:48:00 m_saijo
38  * Delete one <math.h>
39  *
40  * Revision 1.5 2005/08/31 09:06:18 p_grandclement
41  * add math.h in bhole_coal.C
42  *
43  * Revision 1.4 2005/08/29 15:10:14 p_grandclement
44  * Addition of things needed :
45  * 1) For BBH with different masses
46  * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
47  * WORKING YET !!!
48  *
49  * Revision 1.3 2003/11/13 13:43:54 p_grandclement
50  * Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
51  *
52  * Revision 1.2 2002/10/16 14:36:32 j_novak
53  * Reorganization of #include instructions of standard C++, in order to
54  * use experimental version 3 of gcc.
55  *
56  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
57  * LORENE
58  *
59  * Revision 2.15 2001/05/07 09:12:17 phil
60  * *** empty log message ***
61  *
62  * Revision 2.14 2001/04/26 12:23:06 phil
63  * *** empty log message ***
64  *
65  * Revision 2.13 2001/04/26 12:04:17 phil
66  * *** empty log message ***
67  *
68  * Revision 2.12 2001/03/22 10:49:42 phil
69  * *** empty log message ***
70  *
71  * Revision 2.11 2001/02/28 13:23:54 phil
72  * vire kk_auto
73  *
74  * Revision 2.10 2001/01/29 14:31:04 phil
75  * ajout tuype rotation
76  *
77  * Revision 2.9 2001/01/22 09:29:34 phil
78  * vire convergence vers bare masse
79  *
80  * Revision 2.8 2001/01/10 09:31:52 phil
81  * ajoute fait_kk_auto
82  *
83  * Revision 2.7 2000/12/20 15:02:57 phil
84  * *** empty log message ***
85  *
86  * Revision 2.6 2000/12/20 09:09:48 phil
87  * ajout set_statiques
88  *
89  * Revision 2.5 2000/12/18 17:43:06 phil
90  * ajout sortie pour le rayon
91  *
92  * Revision 2.4 2000/12/18 16:38:39 phil
93  * ajout convergence vers une masse donneee
94  *
95  * Revision 2.3 2000/12/14 10:45:38 phil
96  * ATTENTION : PASSAGE DE PHI A PSI
97  *
98  * Revision 2.2 2000/12/04 14:29:17 phil
99  * changement nom omega pour eviter interference avec membre prive
100  *
101  * Revision 2.1 2000/11/17 10:07:14 phil
102  * corrections diverses
103  *
104  * Revision 2.0 2000/11/17 10:04:08 phil
105  * *** empty log message ***
106  *
107  *
108  * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/bhole_coal.C,v 1.9 2016/12/05 16:17:45 j_novak Exp $
109  *
110  */
111 
112 //standard
113 #include <cmath>
114 #include <cstdlib>
115 
116 // Lorene
117 #include "tenseur.h"
118 #include "bhole.h"
119 
120 namespace Lorene {
121 void Bhole_binaire::set_statiques (double precis, double relax) {
122 
123  int nz_un = hole1.mp.get_mg()->get_nzone() ;
124  int nz_deux = hole2.mp.get_mg()->get_nzone() ;
125 
126  set_omega(0) ;
128 
129  int indic = 1 ;
130  int conte = 0 ;
131 
132  cout << "TROUS STATIQUES : " << endl ;
133  while (indic == 1) {
134  Cmp lapse_un_old (hole1.get_n_auto()()) ;
135  Cmp lapse_deux_old (hole2.get_n_auto()()) ;
136 
137  solve_psi (precis, relax) ;
138  solve_lapse (precis, relax) ;
139 
140  double erreur = 0 ;
141  Tbl diff_un (diffrelmax (lapse_un_old, hole1.get_n_auto()())) ;
142  for (int i=1 ; i<nz_un ; i++)
143  if (diff_un(i) > erreur)
144  erreur = diff_un(i) ;
145 
146  Tbl diff_deux (diffrelmax (lapse_deux_old, hole2.get_n_auto()())) ;
147  for (int i=1 ; i<nz_deux ; i++)
148  if (diff_deux(i) > erreur)
149  erreur = diff_deux(i) ;
150 
151 
152  cout << "PAS TOTAL : " << conte << " DIFFERENCE : " << erreur << endl ;
153 
154  if (erreur < precis)
155  indic = -1 ;
156  conte ++ ;
157  }
158 }
159 
160 void Bhole_binaire::coal (double precis, double relax, int nbre_ome, double seuil_search, double m1, double m2, const int sortie) {
161 
162  assert (omega == 0) ;
163  int nz1 = hole1.mp.get_mg()->get_nzone() ;
164  int nz2 = hole1.mp.get_mg()->get_nzone() ;
165 
166  // Distance initiale
167  double distance = hole1.mp.get_ori_x()-hole2.mp.get_ori_x() ;
168  set_pos_axe (distance*(hole2.rayon/(hole1.rayon+hole2.rayon))) ;
169  double scale_linear = (hole1.rayon + hole2.rayon)/2.*distance ;
170 
171  // Omega initial :
172  double angulaire = sqrt((hole1.rayon+hole2.rayon)/distance/distance/distance) ;
173 
174  int indic = 1 ;
175  int conte = 0 ;
176 
177  char name_iteration[40] ;
178  char name_correction[40] ;
179  char name_viriel[40] ;
180  char name_ome [40] ;
181  char name_linear[40] ;
182  char name_axe[40] ;
183  char name_error_m1[40] ;
184  char name_error_m2[40] ;
185  char name_r1[40] ;
186  char name_r2[40] ;
187 
188  sprintf(name_iteration, "ite.dat") ;
189  sprintf(name_correction, "cor.dat") ;
190  sprintf(name_viriel, "vir.dat") ;
191  sprintf(name_ome, "ome.dat") ;
192  sprintf(name_linear, "linear.dat") ;
193  sprintf(name_axe, "axe.dat") ;
194  sprintf(name_error_m1, "error_m1.dat") ;
195  sprintf(name_error_m2, "error_m2.dat") ;
196  sprintf(name_r1, "r1.dat") ;
197  sprintf(name_r2, "r2.dat") ;
198 
199  ofstream fiche_iteration(name_iteration) ;
200  fiche_iteration.precision(8) ;
201 
202  ofstream fiche_correction(name_correction) ;
203  fiche_correction.precision(8) ;
204 
205  ofstream fiche_viriel(name_viriel) ;
206  fiche_viriel.precision(8) ;
207 
208  ofstream fiche_ome(name_ome) ;
209  fiche_ome.precision(8) ;
210 
211  ofstream fiche_linear(name_linear) ;
212  fiche_linear.precision(8) ;
213 
214  ofstream fiche_axe(name_axe) ;
215  fiche_axe.precision(8) ;
216 
217  ofstream fiche_error_m1 (name_error_m1) ;
218  fiche_error_m1.precision(8) ;
219 
220  ofstream fiche_error_m2 (name_error_m2) ;
221  fiche_error_m2.precision(8) ;
222 
223  ofstream fiche_r1 (name_r1) ;
224  fiche_r1.precision(8) ;
225 
226  ofstream fiche_r2 (name_r2) ;
227  fiche_r2.precision(8) ;
228 
229  // LA BOUCLE EN AUGMENTANT OMEGA A LA MAIN PROGRESSIVEMENT :
230  cout << "OMEGA AUGMENTE A LA MAIN." << endl ;
231  double homme = 0 ;
232  for (int pas = 0 ; pas <nbre_ome ; pas ++) {
233 
234  homme += angulaire/nbre_ome ;
235  set_omega (homme) ;
236 
237  Cmp shift_un_old (hole1.get_shift_auto()(0)) ;
238  Cmp shift_deux_old (hole2.get_shift_auto()(0)) ;
239 
240  solve_shift (precis, relax) ;
241  fait_tkij() ;
242 
243  solve_psi (precis, relax) ;
244  solve_lapse (precis, relax) ;
245 
246  double erreur = 0 ;
247  Tbl diff_un (diffrelmax (shift_un_old, hole1.get_shift_auto()(0))) ;
248  for (int i=1 ; i<nz1 ; i++)
249  if (diff_un(i) > erreur)
250  erreur = diff_un(i) ;
251 
252  Tbl diff_deux (diffrelmax (shift_deux_old, hole2.get_shift_auto()(0))) ;
253  for (int i=1 ; i<nz2 ; i++)
254  if (diff_deux(i) > erreur)
255  erreur = diff_deux(i) ;
256 
257  double error_viriel = viriel() ;
258  double error_linear = linear_momentum_systeme_inf()(1)/scale_linear ;
259  double error_m1 = 1.-sqrt(hole1.area()/16./M_PI)/m1 ;
260  double error_m2 = 1.-sqrt(hole2.area()/16./M_PI)/m2 ;
261  double r1 = hole1.mp.val_r(0, 1, 0, 0) ;
262  double r2 = hole2.mp.val_r(0, 1, 0, 0) ;
263 
264  if (sortie != 0) {
265  fiche_iteration << conte << " " << erreur << endl ;
266  fiche_correction << conte << " " << hole1.get_regul() << " " << hole2.get_regul() << endl ;
267  fiche_viriel << conte << " " << error_viriel << endl ;
268  fiche_ome << conte << " " << homme << endl ;
269  fiche_linear << conte << " " << error_linear << endl ;
270  fiche_axe << conte << " " << pos_axe << endl ;
271  fiche_error_m1 << conte << " " << error_m1 << endl ;
272  fiche_error_m2 << conte << " " << error_m2 << endl ;
273  fiche_r1 << conte << " " << r1 << endl ;
274  fiche_r2 << conte << " " << r2 << endl ;
275  }
276 
277  cout << "PAS TOTAL : " << conte << " DIFFERENCE : " << erreur << endl ;
278  conte ++ ;
279  }
280 
281  // BOUCLE AVEC BLOQUE :
282  cout << "OMEGA VARIABLE" << endl ;
283  indic = 1 ;
284  bool scale = false ;
285 
286  while (indic == 1) {
287 
288  Cmp shift_un_old (hole1.get_shift_auto()(0)) ;
289  Cmp shift_deux_old (hole2.get_shift_auto()(0)) ;
290 
291  solve_shift (precis, relax) ;
292  fait_tkij() ;
293 
294  solve_psi (precis, relax) ;
295  solve_lapse (precis, relax) ;
296 
297  double erreur = 0 ;
298  Tbl diff_un (diffrelmax (shift_un_old, hole1.get_shift_auto()(0))) ;
299  for (int i=1 ; i<nz1 ; i++)
300  if (diff_un(i) > erreur)
301  erreur = diff_un(i) ;
302 
303  Tbl diff_deux (diffrelmax (shift_deux_old, hole2.get_shift_auto()(0))) ;
304  for (int i=1 ; i<nz2 ; i++)
305  if (diff_deux(i) > erreur)
306  erreur = diff_deux(i) ;
307 
308  double error_viriel = viriel() ;
309  double error_linear = linear_momentum_systeme_inf()(1)/scale_linear ;
310  double error_m1 = 1.-sqrt(hole1.area()/16./M_PI)/m1 ;
311  double error_m2 = 1.-sqrt(hole2.area()/16./M_PI)/m2 ;
312  double r1 = hole1.mp.val_r(0, 1, 0, 0) ;
313  double r2 = hole2.mp.val_r(0, 1, 0, 0) ;
314 
315  if (sortie != 0) {
316  fiche_iteration << conte << " " << erreur << endl ;
317  fiche_correction << conte << " " << hole1.regul << " " << hole2.regul << endl ;
318  fiche_viriel << conte << " " << error_viriel << endl ;
319  fiche_ome << conte << " " << omega << endl ;
320  fiche_linear << conte << " " << error_linear << endl ;
321  fiche_axe << conte << " " << pos_axe << endl ;
322  fiche_error_m1 << conte << " " << error_m1 << endl ;
323  fiche_error_m2 << conte << " " << error_m2 << endl ;
324  fiche_r1 << conte << " " << r1 << endl ;
325  fiche_r2 << conte << " " << r2 << endl ;
326  }
327 
328  // On modifie omega, position de l'axe et les masses !
329  if (erreur <= seuil_search)
330  scale = true ;
331  if (scale) {
332  double scaling_ome = pow((2-error_viriel)/(2-2*error_viriel), 1.) ;
333  set_omega (omega*scaling_ome) ;
334 
335  double scaling_axe = pow((2-error_linear)/(2-2*error_linear), 0.1) ;
336  set_pos_axe (pos_axe*scaling_axe) ;
337 
338  double scaling_r1 = pow((2-error_m1)/(2-2*error_m1), 0.1) ;
339  hole1.mp.homothetie_interne(scaling_r1) ;
340 
341  double scaling_r2 = pow((2-error_m2)/(2-2*error_m2), 0.1) ;
342  hole2.mp.homothetie_interne(scaling_r2) ;
343  }
344 
345  cout << "PAS TOTAL : " << conte << " DIFFERENCE : " << erreur << endl ;
346  if (erreur < precis)
347  indic = -1 ;
348  conte ++ ;
349  }
350 
351  fiche_iteration.close() ;
352  fiche_correction.close() ;
353  fiche_viriel.close() ;
354  fiche_ome.close() ;
355  fiche_linear.close() ;
356  fiche_axe.close() ;
357  fiche_error_m1.close() ;
358  fiche_error_m2.close() ;
359  fiche_r1.close() ;
360  fiche_r2.close() ;
361 }
362 }
void fait_tkij()
Calculation af the extrinsic curvature tensor.
Definition: bhole_kij.C:91
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
void solve_shift(double precis, double relax)
Solves the equation for the shift, using the Oohara-Nakarmure scheme~: using the current for the ...
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:786
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.
Definition: map_af_radius.C:99
const Tenseur & get_n_auto() const
Returns the part of N generated by the hole.
Definition: bhole.h:395
void init_bhole_binaire()
Initialisation of the system.
void coal(double precis, double relax, int nbre_ome, double search_ome, double m1, double m2, const int sortie=0)
Solves the equation for a particular angular velocity, the systeme being initialized to Misner-Lindqu...
Definition: bhole_coal.C:160
double get_regul() const
Returns the intensity of the regularisation on .
Definition: bhole.h:390
double omega
Position of the axis of rotation.
Definition: bhole.h:769
double area() const
Computes the area of the throat.
Definition: bhole.C:548
Map_af & mp
Affine mapping.
Definition: bhole.h:273
Bhole hole1
Black hole one.
Definition: bhole.h:762
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void solve_psi(double precis, double relax)
Solves the equation for the conformal factor~: The fields are the total values excpet those with s...
double viriel() const
Computes the viriel error, that is the difference between the ADM and the Komar masses, calculated by the asymptotic behaviours of respectively and N .
double rayon
Radius of the horizon in LORENE&#39;s units.
Definition: bhole.h:274
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void set_statiques(double precis, double relax)
Initialize the systeme to Misner Lindquist solution, that is solving for N and in the case ...
Definition: bhole_coal.C:121
void set_omega(double ome)
Sets the orbital velocity to ome.
Definition: bhole.h:791
void solve_lapse(double precis, double relax)
Solves the equation for the lapse~: The fields are the total values excpet those with subscript ...
const Tenseur & get_shift_auto() const
Returns the part of generated by the hole.
Definition: bhole.h:440
Basic array class.
Definition: tbl.h:164
void homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
Definition: map_af.C:744
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542
double regul
Intensity of the correction on the shift vector.
Definition: bhole.h:284
Bhole hole2
Black hole two.
Definition: bhole.h:763