LORENE
bhole_with_ns.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
23char bhole_with_ns_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole/bhole_with_ns.C,v 1.12 2014/10/13 08:52:40 j_novak Exp $" ;
24
25/*
26 * $Id: bhole_with_ns.C,v 1.12 2014/10/13 08:52:40 j_novak Exp $
27 * $Log: bhole_with_ns.C,v $
28 * Revision 1.12 2014/10/13 08:52:40 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.11 2014/10/06 15:12:58 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.10 2007/04/24 20:14:04 f_limousin
35 * Implementation of Dirichlet and Neumann BC for the lapse
36 *
37 * Revision 1.9 2007/02/03 07:46:30 p_grandclement
38 * Addition of term kss for psi BC
39 *
40 * Revision 1.8 2006/04/27 09:12:31 p_grandclement
41 * First try at irrotational black holes
42 *
43 * Revision 1.7 2006/04/25 07:21:57 p_grandclement
44 * Various changes for the NS_BH project
45 *
46 * Revision 1.6 2005/08/29 15:10:13 p_grandclement
47 * Addition of things needed :
48 * 1) For BBH with different masses
49 * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
50 * WORKING YET !!!
51 *
52 * Revision 1.5 2004/03/25 10:28:57 j_novak
53 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
54 *
55 * Revision 1.4 2003/11/25 07:11:09 k_taniguchi
56 * Change some arguments from the class Etolie_bin to Et_bin_nsbh.
57 *
58 * Revision 1.3 2003/11/13 13:43:53 p_grandclement
59 * Addition of things needed for Bhole::update_metric (const Etoile_bin&, double, double)
60 *
61 * Revision 1.2 2003/10/24 13:05:49 p_grandclement
62 * correction of the equations for Bin_ns_bh...
63 *
64 * Revision 1.1 2003/02/13 16:40:25 p_grandclement
65 * Addition of various things for the Bin_ns_bh project, non of them being
66 * completely tested
67 *
68 *
69 *
70 * $Header: /cvsroot/Lorene/C++/Source/Bhole/bhole_with_ns.C,v 1.12 2014/10/13 08:52:40 j_novak Exp $
71 *
72 */
73
74//standard
75#include <cstdlib>
76#include <cmath>
77
78// Lorene
79#include "tenseur.h"
80#include "bhole.h"
81#include "proto.h"
82#include "utilitaires.h"
83#include "et_bin_nsbh.h"
84#include "graphique.h"
85#include "scalar.h"
86
87//Resolution pour le lapse pour 1 seul trou
88namespace Lorene {
89void Bhole::solve_lapse_with_ns (double relax, int bound_nn, double lim_nn) {
90
91 assert ((relax>0) && (relax<=1)) ;
92
93 cout << "Resolution LAPSE" << endl ;
94
95 // Pour la relaxation ...
96 Cmp lapse_old (n_auto()) ;
98 Tenseur kk (mp) ;
99 kk = 0 ;
100 Tenseur work(mp) ;
101 work.set_etat_qcq() ;
102 for (int i=0 ; i<3 ; i++) {
103 work.set() = auxi(i, i) ;
104 kk = kk + work ;
105 }
106
107 // La source
108 Cmp psiq (pow(psi_tot(), 4.)) ;
109 psiq.std_base_scal() ;
110 Cmp source
112 +psiq*n_tot()*kk()) ;
113 source.std_base_scal() ;
114
115 Cmp soluce(n_auto()) ;
116
117 if (bound_nn == 0){
118 // Dirichlet
119 Valeur limite (mp.get_mg()->get_angu()) ;
120 limite = -0.5 + lim_nn ;
121 int np = mp.get_mg()->get_np(1) ;
122 int nt = mp.get_mg()->get_nt(1) ;
123 for (int k=0 ; k<np ; k++)
124 for (int j=0 ; j<nt ; j++)
125 limite.set(0,k,j,0) -= n_comp() (1, k, j, 0) ;
126 limite.std_base_scal() ;
127
128 soluce = source.poisson_dirichlet(limite, 0) ;
129 }
130 else {
131 assert(bound_nn == 1);
132 // Neumann
133 Valeur limite (mp.get_mg()->get_angu()) ;
134 limite.annule_hard() ;
135 int np = mp.get_mg()->get_np(1) ;
136 int nt = mp.get_mg()->get_nt(1) ;
137 for (int k=0 ; k<np ; k++)
138 for (int j=0 ; j<nt ; j++)
139 limite.set(0,k,j,0) -= n_tot()(1, k, j, 0)/psi_tot()(1,k,j,0)*
140 psi_tot().dsdr()(1,k,j,0) ;
141 limite.std_base_scal() ;
142
143 soluce = source.poisson_neumann(limite, 0) ;
144 }
145
146 soluce = soluce + 0.5 ;
147
148 n_auto.set() = relax*soluce + (1-relax)*lapse_old ;
149 n_auto.set().raccord(3) ;
150}
151
152// Resolution sur Psi :
153void Bhole::solve_psi_with_ns (double relax) {
154
155 assert ((relax>0) && (relax<=1)) ;
156
157 cout << "Resolution PSI" << endl ;
158
159 Cmp psi_old (psi_auto()) ;
161 Tenseur kk (mp) ;
162 kk = 0 ;
163 Tenseur work(mp) ;
164 work.set_etat_qcq() ;
165 for (int i=0 ; i<3 ; i++) {
166 work.set() = auxi(i, i) ;
167 kk = kk + work ;
168 }
169 Cmp psic (pow(psi_tot(), 5.)) ;
170 psic.std_base_scal() ;
171
172 // La source :
173 Cmp source (-psic*kk()/8.) ;
174 source.std_base_scal() ;
175
176 // Condition limite :
177 Valeur limite (mp.get_mg()->get_angu()) ;
178 limite = 1 ;
179
180
181 int np = mp.get_mg()->get_np(1) ;
182 int nt = mp.get_mg()->get_nt(1) ;
183
184 double* vec_s = new double[3] ;
185 Mtbl tet_mtbl (mp.get_mg()) ;
186 tet_mtbl = mp.tet ;
187 Mtbl phi_mtbl (mp.get_mg()) ;
188 phi_mtbl = mp.phi ;
189
190 for (int k=0 ; k<np ; k++)
191 for (int j=0 ; j<nt ; j++) {
192
193 double tet = tet_mtbl(1,k,j,0) ;
194 double phi = phi_mtbl(1,k,j,0) ;
195 vec_s[0] = cos(phi)*sin(tet) ;
196 vec_s[1] = sin(phi)*sin(tet) ;
197 vec_s[2] = cos(tet) ;
198 double part_ss = 0 ;
199 if (tkij_tot.get_etat()==ETATQCQ)
200 for (int m=0 ; m<3 ; m++)
201 for (int n=0 ; n<3 ; n++)
202 part_ss += vec_s[m]*vec_s[n]*tkij_tot(m,n)(1,k,j,0) ;
203 part_ss *= pow(psi_tot()(1,k,j,0),3.)/4. ;
204
205
206 limite.set(0, k, j, 0) = -0.5/rayon*psi_tot()(1, k, j, 0) -
207 psi_comp().dsdr()(1, k, j, 0) - part_ss ;
208}
209
210
211 limite.std_base_scal() ;
212
213 Cmp soluce (source.poisson_neumann(limite, 0)) ;
214 soluce = soluce + 1./2. ;
215
216 psi_auto.set() = relax*soluce + (1-relax)*psi_old ;
217 psi_auto.set().raccord(3) ;
218
219}
220
221// Le shift. Processus iteratif pour cause de CL.
223 double precision, double relax,
224 int bound_nn, double lim_nn) {
225
226 assert (precision > 0) ;
227 assert ((relax>0) && (relax<=1)) ;
228
229 cout << "resolution SHIFT" << endl ;
230
231 Tenseur shift_old (shift_auto) ;
232
235 source.set_std_base() ;
236
237 // On verifie si les 3 composantes ne sont pas nulles :
238 if (source.get_etat() == ETATQCQ) {
239 int indic = 0 ;
240 for (int i=0 ; i<3 ; i++)
241 if (source(i).get_etat() == ETATQCQ)
242 indic = 1 ;
243 if (indic ==0)
244 for (int i=0 ; i<3 ; i++)
245 source.set_etat_zero() ;
246 }
247
248
249 // On filtre les hautes frequences pour raison de stabilite :
250 if (source.get_etat() == ETATQCQ)
251 for (int i=0 ; i<3 ; i++)
252 source.set(i).filtre(4) ;
253
254
255 // On determine les conditions limites en fonction de omega et de NS :
256 int np = mp.get_mg()->get_np(1) ;
257 int nt = mp.get_mg()->get_nt(1) ;
258
259 Mtbl x_mtbl (mp.get_mg()) ;
260 x_mtbl.set_etat_qcq() ;
261 Mtbl y_mtbl (mp.get_mg()) ;
262 y_mtbl.set_etat_qcq() ;
263 x_mtbl = mp.x ;
264 y_mtbl = mp.y ;
265
266 double air, theta, phi, xabs, yabs, zabs ;
267 Mtbl Xabs (mp.get_mg()) ;
268 Xabs = mp.xa ;
269 Mtbl Yabs (mp.get_mg()) ;
270 Yabs = mp.ya ;
271 Mtbl Zabs (mp.get_mg()) ;
272 Zabs = mp.za ;
273
274 Mtbl tet_mtbl (mp.get_mg()) ;
275 tet_mtbl = mp.tet ;
276 Mtbl phi_mtbl (mp.get_mg()) ;
277 phi_mtbl = mp.phi ;
278
279 // Les bases pour les conditions limites :
280 Base_val** bases = mp.get_mg()->std_base_vect_cart() ;
281
282 Valeur lim_x (mp.get_mg()->get_angu()) ;
283 lim_x = 1 ;
284 Valeur lim_y (mp.get_mg()->get_angu()) ;
285 lim_y = 1 ;
286 Valeur lim_z (mp.get_mg()->get_angu()) ;
287 lim_z = 1 ;
288
289 for (int k=0 ; k<np ; k++)
290 for (int j=0 ; j<nt ; j++) {
291
292 double tet = tet_mtbl(1,k,j,0) ;
293 double phy = phi_mtbl(1,k,j,0) ;
294
295 xabs = Xabs (1, k, j, 0) ;
296 yabs = Yabs (1, k, j, 0) ;
297 zabs = Zabs (1, k, j, 0) ;
298
299 ns.get_mp().convert_absolute (xabs, yabs, zabs, air, theta, phi) ;
300
301 lim_x.set(0, k, j, 0) = omega*Yabs(0, 0, 0, 0) +
302 omega_local*y_mtbl(1,k,j,0) -
303 ns.get_shift_auto()(0).val_point(air, theta, phi) +
304 n_tot()(1,k,j,0)/psi_tot()(1,k,j,0)/psi_tot()(1,k,j,0)*
305 cos(phy)*sin(tet) ;
306 lim_x.base = *bases[0] ;
307
308
309 lim_y.set(0, k, j, 0) = -omega*Xabs(0, 0, 0, 0) -
310 omega_local*x_mtbl(1,k,j,0) -
311 ns.get_shift_auto()(1).val_point(air, theta, phi) +
312 n_tot()(1,k,j,0)/psi_tot()(1,k,j,0)/psi_tot()(1,k,j,0)*
313 sin(phy)*sin(tet) ;
314
315 lim_z.set(0, k, j, 0) = -
316 ns.get_shift_auto()(2).val_point(air, theta, phi) +
317 n_tot()(1,k,j,0)/psi_tot()(1,k,j,0)/psi_tot()(1,k,j,0)*cos(tet) ;
318 }
319
320 lim_x.base = *bases[0] ;
321 lim_y.base = *bases[1] ;
322 lim_z.base = *bases[2] ;
323
324 // On n'en a plus besoin
325 for (int i=0 ; i<3 ; i++)
326 delete bases[i] ;
327 delete [] bases ;
328
329 // On resout :
330 poisson_vect_frontiere(1./3., source, shift_auto, lim_x, lim_y,
331 lim_z, 0, precision, 20) ;
332
333 shift_auto = relax*shift_auto + (1-relax)*shift_old ;
334
335 for (int i=0; i<3; i++)
336 shift_auto.set(i).raccord(3) ;
337
338
339 // Regularisation of the shift if necessary
340 // -----------------------------------------
341 if (bound_nn == 0 && lim_nn == 0)
342 regul = regle (shift_auto, ns.get_shift_auto(), omega, omega_local) ;
343 else
344 regul = 0. ;
345
346}
347
348
349void Bhole::equilibrium (const Et_bin_nsbh& comp,
350 double precision, double relax,
351 int bound_nn, double lim_nn) {
352
353 // Solve for the lapse :
354 solve_lapse_with_ns (relax, bound_nn, lim_nn) ;
355
356 // Solve for the conformal factor :
357 solve_psi_with_ns (relax) ;
358
359 if (omega != 0)
360 // Solve for the shift vector :
361 solve_shift_with_ns (comp, precision, relax, bound_nn, lim_nn) ;
362
363}
364
365
366void Bhole::update_metric (const Et_bin_nsbh& comp) {
367
368 fait_n_comp(comp) ;
369 fait_psi_comp(comp) ;
370 /*
371 Scalar lapse_auto (n_auto()) ;
372 Scalar lapse_tot (n_tot()) ;
373 Scalar lapse_comp (n_comp()) ;
374 des_meridian(lapse_auto, 0, 7, "n_auto", 0) ;
375 des_meridian(lapse_comp, 0, 7, "n_comp", 11) ;
376 des_meridian(lapse_tot, 0, 7, "n_tot", 1) ;
377
378 Scalar psiauto (psi_auto()) ;
379 Scalar psitot (psi_tot()) ;
380 des_meridian(psiauto, 0, 7, "psi_auto", 2) ;
381 des_meridian(psitot, 0, 7, "psi_tot", 3) ;
382 */
383
384}
385
386
387}
Bases of the spectral expansions.
Definition base_val.h:322
Tenseur psi_auto
Part of generated by the hole.
Definition bhole.h:290
Tenseur shift_auto
Part of generated by the hole.
Definition bhole.h:297
double omega_local
local angular velocity
Definition bhole.h:276
Tenseur psi_tot
Total .
Definition bhole.h:292
Tenseur psi_comp
Part of generated by the companion hole.
Definition bhole.h:291
Tenseur_sym tkij_auto
Auto .
Definition bhole.h:307
double omega
Angular velocity in LORENE's units.
Definition bhole.h:275
double regul
Intensity of the correction on the shift vector.
Definition bhole.h:284
Tenseur grad_n_tot
Total gradient of N .
Definition bhole.h:294
Tenseur n_auto
Part of N generated by the hole.
Definition bhole.h:286
double rayon
Radius of the horizon in LORENE's units.
Definition bhole.h:274
void solve_psi_with_ns(double relax)
Solves the equation for ~:
Tenseur n_comp
Part of N generated by the companion hole.
Definition bhole.h:287
Map_af & mp
Affine mapping.
Definition bhole.h:273
Tenseur_sym taij_tot
Total , which must be zero on the horizon of the regularisation on the shift has been done.
Definition bhole.h:305
Tenseur n_tot
Total N .
Definition bhole.h:288
void solve_lapse_with_ns(double relax, int bound_nn, double lim_nn)
Solves the equation for N ~:
void fait_psi_comp(const Bhole &comp)
Imports the part of due to the companion hole comp .
Definition bhole.C:280
Tenseur_sym tkij_tot
Total .
Definition bhole.h:308
void solve_shift_with_ns(const Et_bin_nsbh &ns, double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for ~:
void fait_n_comp(const Bhole &comp)
Imports the part of N due to the companion hole comp .
Definition bhole.C:254
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
Cmp poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Cmp::poisson() .
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition cmp_manip.C:74
Cmp poisson_neumann(const Valeur &, int) const
Idem as Cmp::poisson_dirichlet , the boundary condition being on the radial derivative of the solutio...
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Class for a star in a NS-BH binary system.
Definition et_bin_nsbh.h:79
const Tenseur & get_shift_auto() const
Returns the part of the shift vector generated principaly by the star.
Definition etoile.h:1152
const Map & get_mp() const
Returns the mapping.
Definition etoile.h:659
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord ya
Absolute y coordinate.
Definition map.h:731
void convert_absolute(double xx, double yy, double zz, double &rr, double &theta, double &pphi) const
Determines the coordinates corresponding to given absolute Cartesian coordinates (X,...
Definition map.C:302
Coord za
Absolute z coordinate.
Definition map.h:732
Coord tet
coordinate centered on the grid
Definition map.h:719
Coord x
x coordinate centered on the grid
Definition map.h:726
Coord phi
coordinate centered on the grid
Definition map.h:720
Coord xa
Absolute x coordinate.
Definition map.h:730
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Multi-domain array.
Definition mtbl.h:118
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition mtbl.C:299
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Lorene prototypes.
Definition app_hor.h:64