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