LORENE
bin_ns_bh_coal.C
1/*
2 * Copyright (c) 2004 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 bin_ns_bh_coal_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $" ;
24
25/*
26 * $Id: bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $
27 * $Log: bin_ns_bh_coal.C,v $
28 * Revision 1.12 2014/10/13 08:52:43 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.11 2014/10/06 15:13:01 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.10 2007/04/24 20:13:53 f_limousin
35 * Implementation of Dirichlet and Neumann BC for the lapse
36 *
37 * Revision 1.9 2006/09/05 13:39:43 p_grandclement
38 * update of the bin_ns_bh project
39 *
40 * Revision 1.8 2006/06/23 07:09:24 p_grandclement
41 * Addition of spinning black hole
42 *
43 * Revision 1.7 2006/06/01 12:47:52 p_grandclement
44 * update of the Bin_ns_bh project
45 *
46 * Revision 1.6 2006/04/27 09:12:33 p_grandclement
47 * First try at irrotational black holes
48 *
49 * Revision 1.5 2006/04/25 07:21:57 p_grandclement
50 * Various changes for the NS_BH project
51 *
52 * Revision 1.4 2006/03/30 07:33:45 p_grandclement
53 * *** empty log message ***
54 *
55 * Revision 1.3 2005/12/01 12:59:10 p_grandclement
56 * Files for bin_ns_bh project
57 *
58 * Revision 1.2 2005/10/18 13:12:32 p_grandclement
59 * update of the mixted binary codes
60 *
61 * Revision 1.1 2005/08/29 15:10:15 p_grandclement
62 * Addition of things needed :
63 * 1) For BBH with different masses
64 * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
65 * WORKING YET !!!
66 *
67 *
68 * $Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $
69 *
70 */
71
72//standard
73#include <cstdlib>
74
75// Lorene
76#include "tenseur.h"
77#include "bin_ns_bh.h"
78#include "unites.h"
79#include "graphique.h"
80
81namespace Lorene {
82void Bin_ns_bh::coal (double precis, double relax, int itemax_equil,
83 int itemax_mp_et, double ent_c_init, double seuil_masses,
84 double dist, double m1, double m2, double spin_cible,
85 double scale_ome_local, const int sortie, int bound_nn,
86 double lim_nn) {
87
88 using namespace Unites ;
89
90 int nz_bh = hole.mp.get_mg()->get_nzone() ;
91 int nz_ns = star.mp.get_mg()->get_nzone() ;
92
93 // Distance initiale
94 double distance = fabs(hole.mp.get_ori_x()-star.mp.get_ori_x()) ;
95 double mass_ns = star.mass_g() * ggrav;
96 double mass_bh = hole.masse_adm_seul() ;
97 double axe_pos = star.mp.get_ori_x() ;
98 double scale_linear = (mass_ns+mass_bh)/2.*distance*omega ;
99
100 char name_iteration[40] ;
101 char name_correction[40] ;
102 char name_viriel[40] ;
103 char name_ome [40] ;
104 char name_linear[40] ;
105 char name_axe[40] ;
106 char name_error_m1[40] ;
107 char name_error_m2[40] ;
108 char name_hor[40] ;
109 char name_entc[40] ;
110 char name_dist[40] ;
111 char name_spin[40] ;
112 char name_ome_local[40] ;
113
114 sprintf(name_iteration, "ite.dat") ;
115 sprintf(name_correction, "cor.dat") ;
116 sprintf(name_viriel, "vir.dat") ;
117 sprintf(name_ome, "ome.dat") ;
118 sprintf(name_linear, "linear.dat") ;
119 sprintf(name_axe, "axe.dat") ;
120 sprintf(name_error_m1, "error_m_bh.dat") ;
121 sprintf(name_error_m2, "error_m_ns.dat") ;
122 sprintf(name_hor, "hor.dat") ;
123 sprintf(name_entc, "entc.dat") ;
124 sprintf(name_dist, "distance.dat") ;
125 sprintf(name_spin, "spin.dat") ;
126 sprintf(name_ome_local, "ome_local.dat") ;
127
128 ofstream fiche_iteration(name_iteration) ;
129 fiche_iteration.precision(8) ;
130
131 ofstream fiche_correction(name_correction) ;
132 fiche_correction.precision(8) ;
133
134 ofstream fiche_viriel(name_viriel) ;
135 fiche_viriel.precision(8) ;
136
137 ofstream fiche_ome(name_ome) ;
138 fiche_ome.precision(8) ;
139
140 ofstream fiche_linear(name_linear) ;
141 fiche_linear.precision(8) ;
142
143 ofstream fiche_axe(name_axe) ;
144 fiche_axe.precision(8) ;
145
146 ofstream fiche_error_m1 (name_error_m1) ;
147 fiche_error_m1.precision(8) ;
148
149 ofstream fiche_error_m2 (name_error_m2) ;
150 fiche_error_m2.precision(8) ;
151
152 ofstream fiche_hor (name_hor) ;
153 fiche_hor.precision(8) ;
154
155 ofstream fiche_entc (name_entc) ;
156 fiche_entc.precision(8) ;
157
158 ofstream fiche_dist (name_dist) ;
159 fiche_dist.precision(8) ;
160
161 ofstream fiche_spin (name_spin) ;
162 fiche_spin.precision(8) ;
163
164 ofstream fiche_ome_local (name_ome_local) ;
165 fiche_ome_local.precision(8) ;
166
167 bool loop = true ;
168 bool search_masses = false ;
169 double ent_c = ent_c_init ;
170
171 Cmp shift_bh_old (hole.mp) ;
172 Cmp shift_ns_old (star.mp) ;
173
174 double erreur = 1 ;
175
176 int conte = 0 ;
177
178 double old_mass_ns ;
179 mass_ns = star.mass_b() ;
180 double spin_old ;
181 double spin = 1;
182 bool adapt = true ;
183
184 while (loop) {
185
186 spin_old = spin ;
187 spin = hole.local_momentum() ;
188 if (sortie !=0) {
189 fiche_ome_local << conte << " " << hole.omega_local << endl ;
190 fiche_spin << conte << " " << spin/m1/m1 << endl ;
191 }
192
193 double conv_spin = fabs(spin-spin_old) ;
194 double error_spin = spin - spin_cible ;
195 double rel_diff_spin = (spin_cible==0) ? fabs(error_spin) :
196 fabs(error_spin)/spin_cible ;
197 if ((conv_spin*2<rel_diff_spin) && (search_masses) &&
198 hole.get_rot_state() != COROT) {
199 double func = scale_ome_local*log((2+error_spin)/(2+2*error_spin));
201 }
202
203 old_mass_ns = mass_ns ;
204
205 if (hole.get_shift_auto().get_etat() != ETATZERO)
206 shift_bh_old = hole.get_shift_auto()(0) ;
207 else
208 shift_bh_old = 0 ;
209
210 if (star.get_shift_auto().get_etat() != ETATZERO)
211 shift_ns_old = star.get_shift_auto()(0) ;
212 else
213 shift_ns_old = 0 ;
214
216 star.fait_d_psi() ;
217 star.hydro_euler() ;
218
219 Tbl diff (7) ;
220 diff.set_etat_qcq() ;
221 int ite ;
222
223
224 star.equilibrium_nsbh (adapt, ent_c, ite, itemax_equil, itemax_mp_et,
225 relax, itemax_mp_et, relax, diff) ;
226
227 hole.update_metric(star) ;
228
229 hole.equilibrium (star, precis, relax, bound_nn, lim_nn) ;
230 cout << "Apres equilibrium" << endl ;
231
233 cout << "Apres star::update_metric" << endl ;
234
236 cout << "Apres star::update_metric_der_comp" << endl ;
237 fait_tkij(bound_nn, lim_nn) ;
238 cout << "Apres Bin_ns_bh::fait_tkij" << endl ;
239
240 erreur = 0 ;
241 Tbl diff_bh (diffrelmax (shift_bh_old, hole.get_shift_auto()(0))) ;
242 for (int i=1 ; i<nz_bh ; i++)
243 if (diff_bh(i) > erreur)
244 erreur = diff_bh(i) ;
245 Tbl diff_ns (diffrelmax (shift_ns_old, star.get_shift_auto()(0))) ;
246 for (int i=0 ; i<nz_ns ; i++)
247 if (diff_ns(i) > erreur)
248 erreur = diff_ns(i) ;
249
250 if (erreur<seuil_masses)
251 search_masses = true ;
252
253 mass_ns = star.mass_b() ;
254
255 cout << "Avant viriel" << endl ;
256 double error_viriel = viriel() ;
257 cout << "Apres viriel" << endl ;
258 double error_linear = linear_momentum_systeme_inf()(1)/scale_linear ;
259 cout << "Apres linear" << endl ;
260 double error_m1 = 1.-sqrt(hole.area()/16./M_PI)/m1 ;
261 cout << "Apres Mbh" << endl ;
262 double error_m2 = 1 - mass_ns/m2 ;
263 cout << "Apres Mns" << endl ;
264
265 if (sortie != 0) {
266 fiche_iteration << conte << " " << erreur << endl ;
267 fiche_correction << conte << " " << hole.regul << endl ;
268 fiche_viriel << conte << " " << error_viriel << endl ;
269 fiche_linear << conte << " " << error_linear << endl ;
270 fiche_error_m1 << conte << " " << error_m1 << endl ;
271 fiche_error_m2 << conte << " " << error_m2 << endl ;
272 fiche_hor << conte << " " << hole.mp.val_r(0, 1, 0,0) << endl ;
273 fiche_entc << conte << " " << ent_c << endl ;
274 fiche_dist << conte << " " << distance << endl ;
275 fiche_ome << conte << " " << omega << endl ;
276 fiche_axe << conte << " " << axe_pos << endl ;
277 }
278
279 // The axis position
280 double scaling_axe = pow((2+error_linear)/(2+2*error_linear), 0.1) ;
281 axe_pos *= scaling_axe ;
282 star.set_mp().set_ori (axe_pos, 0, 0) ;
283 hole.set_mp().set_ori (-distance + axe_pos, 0, 0) ;
284
285 // Value of omega
286 double new_ome = star.compute_angul() ;
287 if (new_ome !=0) {
288 set_omega(new_ome) ;
289 if (hole.get_rot_state() == COROT)
290 hole.set_omega_local(new_ome) ;
291 }
292
293 // The right distance
294 double error_dist = (distance-dist)/dist ;
295 double scale_d = pow((2+error_dist)/(2+2*error_dist), 0.2) ;
296 distance *= scale_d ;
297
298
299 // Converge to the right masses :
300 if (search_masses) {
301
302 double scaling_r = pow((2-error_m1)/(2-2*error_m1), 0.25) ;
303 hole.mp.homothetie_interne(scaling_r) ;
304 hole.set_rayon(hole.get_rayon()*scaling_r) ;
305 }
306
307 // The case of the NS :
308 double convergence = fabs(mass_ns - old_mass_ns)/mass_ns ;
309 double rel_diff = fabs(error_m2) ;
310 if ((search_masses) && (convergence*2 < rel_diff)) {
311 double scaling_ent = pow((2-error_m2)/(2-2*error_m2), 1) ;
312 ent_c *= scaling_ent ;
313
314 }
315
316
317
318 cout << "PAS TOTAL : " << conte << " DIFFERENCE : " << erreur << endl ;
319 //if (erreur < 1e-4)
320 //adapt = false ;
321
322 if (erreur < precis)
323 loop = false ;
324 conte ++ ;
325 }
326
327
328 fiche_iteration.close() ;
329 fiche_correction.close() ;
330 fiche_viriel.close() ;
331 fiche_ome.close() ;
332 fiche_linear.close() ;
333 fiche_axe.close() ;
334 fiche_error_m1.close() ;
335 fiche_error_m2.close() ;
336 fiche_hor.close() ;
337 fiche_entc.close() ;
338 fiche_dist.close() ;
339 fiche_spin.close() ;
340 fiche_ome_local.close() ;
341}
342}
double omega_local
local angular velocity
Definition bhole.h:276
double regul
Intensity of the correction on the shift vector.
Definition bhole.h:284
double masse_adm_seul() const
Calculates the ADM mass of the black hole using : .
Map_af & set_mp()
Read/write of the mapping.
Definition bhole.h:342
double area() const
Computes the area of the throat.
Definition bhole.C:545
double get_rot_state() const
Returns the state of rotation.
Definition bhole.h:373
Map_af & mp
Affine mapping.
Definition bhole.h:273
const Tenseur & get_shift_auto() const
Returns the part of generated by the hole.
Definition bhole.h:440
double get_rayon() const
Returns the radius of the horizon.
Definition bhole.h:347
void set_omega_local(double ome)
Sets the local angular velocity to ome .
Definition bhole.h:369
void set_rayon(double ray)
Sets the radius of the horizon to ray .
Definition bhole.h:352
void fait_tkij(int bound_nn=-1, double lim_nn=0)
Computation of the extrinsic curvature tensor for both {\tt star} and {\tt bhole}.
void set_omega(double)
Sets the orbital angular velocity [{\tt f_unit}].
Definition bin_ns_bh.C:218
double x_axe
Absolute X coordinate of the rotation axis.
Definition bin_ns_bh.h:140
Et_bin_nsbh star
The neutron star.
Definition bin_ns_bh.h:128
Bhole hole
The black hole.
Definition bin_ns_bh.h:131
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition bin_ns_bh.h:136
void update_metric(const Bhole &comp)
Computes metric coefficients from known potentials, when the companion is a black hole.
virtual void equilibrium_nsbh(double ent_c, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, double thres_adapt, const Tbl &fact, Tbl &diff)
Computes an equilibrium configuration in a NS-BH binary system.
virtual void kinematics(double omega, double x_axe)
Computes the quantities bsn and pot_centri .
void update_metric_der_comp(const Bhole &comp)
Computes the derivative of metric functions related to the companion black hole.
virtual double mass_b() const
Baryon mass.
const Tenseur & get_shift_auto() const
Returns the part of the shift vector generated principaly by the star.
Definition etoile.h:1152
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
virtual double mass_g() const
Gravitational mass.
void fait_d_psi()
Computes the gradient of the total velocity potential .
Definition etoile_bin.C:764
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Map & set_mp()
Read/write of the mapping.
Definition etoile.h:601
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
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.
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
Definition map.C:253
double get_ori_x() const
Returns the x coordinate of the origin.
Definition map.h:768
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.