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