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