LORENE
et_bin_bhns_extr_equil.C
1 /*
2  * Method of class Et_bin_bhns_extr to compute an equilibrium configuration
3  * of a BH-NS binary system with an extreme mass ratio
4  *
5  * (see file et_bin_bhns_extr.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2004-2005 Keisuke Taniguchi
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: et_bin_bhns_extr_equil.C,v 1.12 2016/12/05 16:17:52 j_novak Exp $
33  * $Log: et_bin_bhns_extr_equil.C,v $
34  * Revision 1.12 2016/12/05 16:17:52 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.11 2014/10/13 08:52:54 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.10 2014/10/06 15:13:07 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.9 2005/02/28 23:11:05 k_taniguchi
44  * Modification to include the case of the conformally flat background metric
45  *
46  * Revision 1.8 2005/01/25 17:33:19 k_taniguchi
47  * Suppression of the filter for the source term of the shift vector.
48  *
49  * Revision 1.7 2005/01/03 18:01:12 k_taniguchi
50  * Addition of the method to fix the position of the neutron star
51  * in the coordinate system.
52  *
53  * Revision 1.6 2004/12/29 16:29:55 k_taniguchi
54  * Suppression of "dzpius" for the shift vector.
55  *
56  * Revision 1.5 2004/12/22 18:26:53 k_taniguchi
57  * Change an argument of poisson_vect_falloff.
58  *
59  * Revision 1.4 2004/12/06 17:59:50 k_taniguchi
60  * Change the position of resize.
61  *
62  * Revision 1.3 2004/12/02 21:31:56 k_taniguchi
63  * Set a filter for the shift vector.
64  *
65  * Revision 1.2 2004/12/02 15:05:36 k_taniguchi
66  * Modification of the procedure for resize.
67  *
68  * Revision 1.1 2004/11/30 20:48:45 k_taniguchi
69  * *** empty log message ***
70  *
71  *
72  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_equil.C,v 1.12 2016/12/05 16:17:52 j_novak Exp $
73  *
74  */
75 
76 // C headers
77 #include <cmath>
78 
79 // Lorene headers
80 #include "et_bin_bhns_extr.h"
81 #include "etoile.h"
82 #include "map.h"
83 #include "coord.h"
84 #include "param.h"
85 #include "eos.h"
86 #include "graphique.h"
87 #include "utilitaires.h"
88 #include "unites.h"
89 
90 namespace Lorene {
91 void Et_bin_bhns_extr::equil_bhns_extr_ks(double ent_c, const double& mass,
92  const double& sepa, int mermax,
93  int mermax_poisson,
94  double relax_poisson,
95  int mermax_potvit,
96  double relax_potvit, int np_filter,
97  double thres_adapt,
98  Tbl& diff) {
99 
100  // Fundamental constants and units
101  // -------------------------------
102  using namespace Unites ;
103 
104  assert( kerrschild == true ) ;
105 
106  // Initializations
107  // ---------------
108 
109  const Mg3d* mg = mp.get_mg() ;
110  int nz = mg->get_nzone() ; // total number of domains
111 
112  // The following is required to initialize mp_prev as a Map_et:
113  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
114 
115  // Domain and radial indices of points at the surface of the star:
116  int l_b = nzet - 1 ;
117  int i_b = mg->get_nr(l_b) - 1 ;
118  int k_b ;
119  int j_b ;
120 
121  // Value of the enthalpy defining the surface of the star
122  double ent_b = 0 ;
123 
124  // Error indicators
125  // ----------------
126 
127  double& diff_ent = diff.set(0) ;
128  double& diff_vel_pot = diff.set(1) ;
129  double& diff_logn = diff.set(2) ;
130  double& diff_beta = diff.set(3) ;
131  double& diff_shift_x = diff.set(4) ;
132  double& diff_shift_y = diff.set(5) ;
133  double& diff_shift_z = diff.set(6) ;
134 
135  // Parameters for the function Map_et::adapt
136  // -----------------------------------------
137 
138  Param par_adapt ;
139  int nitermax = 100 ;
140  int niter ;
141  int adapt_flag = 1 ; // 1 = performs the full computation,
142  // 0 = performs only the rescaling by
143  // the factor alpha_r
144  int nz_search = nzet ; // Number of domains for searching
145  // the enthalpy isosurfaces
146  double precis_secant = 1.e-14 ;
147  double alpha_r ;
148  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
149 
150  Tbl ent_limit(nz) ;
151 
152  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
153  // locate zeros by the secant method
154  par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
155  // to the isosurfaces of ent is to be
156  // performed
157  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
158  // the enthalpy isosurface
159  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
160  // 0 = performs only the rescaling by
161  // the factor alpha_r
162  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
163  // (theta_*, phi_*)
164  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
165  // (theta_*, phi_*)
166  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
167  // used in the secant method
168  par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
169  // the determination of zeros by
170  // the secant method
171  par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
172  // 0 = contracting mapping
173  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
174  // distances will be multiplied
175  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
176  // to define the isosurfaces
177 
178  // Parameters for the function Map_et::poisson for logn_auto
179  // ---------------------------------------------------------
180 
181  double precis_poisson = 1.e-16 ;
182 
183  Param par_poisson1 ;
184 
185  par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
186  par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
187  par_poisson1.add_double(precis_poisson, 1) ; // required precision
188  par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
189  // used
190  par_poisson1.add_cmp_mod( ssjm1_logn ) ;
191 
192  // Parameters for the function Map_et::poisson for beta_auto
193  // ---------------------------------------------------------
194 
195  Param par_poisson2 ;
196 
197  par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
198  par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
199  par_poisson2.add_double(precis_poisson, 1) ; // required precision
200  par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
201  // used
202  par_poisson2.add_cmp_mod( ssjm1_beta ) ;
203 
204  // Parameters for the function Tenseur::poisson_vect
205  // -------------------------------------------------
206 
207  Param par_poisson_vect ;
208 
209  par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
210  // iterations
211  par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
212  par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
213  par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
214  par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
215  par_poisson_vect.add_int_mod(niter, 0) ;
216 
217  // External potential
218  // See Eq (99) from Gourgoulhon et al. (2001)
219  // -----------------------------------------
220 
221  Tenseur pot_ext = logn_comp + pot_centri + loggam ;
222 
223  Tenseur ent_jm1 = ent ; // Enthalpy at previous step
224 
225  Tenseur source(mp) ; // source term in the equation for logn_auto
226  // and beta_auto
227 
228  Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
229  // equation for shift_auto
230 
231  //==========================================================//
232  // Start of iteration //
233  //==========================================================//
234 
235  for(int mer=0 ; mer<mermax ; mer++ ) {
236 
237  cout << "-----------------------------------------------" << endl ;
238  cout << "step: " << mer << endl ;
239  cout << "diff_ent = " << diff_ent << endl ;
240 
241  //------------------------------------------------------
242  // Resolution of the elliptic equation for the velocity
243  // scalar potential
244  //------------------------------------------------------
245 
246  if (irrotational) {
247  diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
248  precis_poisson, relax_potvit) ;
249  }
250 
251  //-------------------------------------
252  // Computation of the new radial scale
253  //-------------------------------------
254 
255  // alpha_r (r = alpha_r r') is determined so that the enthalpy
256  // takes the requested value ent_b at the stellar surface
257 
258  // Values at the center of the star:
259  double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
260  double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
261 
262  // Search for the reference point (theta_*, phi_*)
263  // [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
264  // at the surface of the star
265  // ------------------------------------------------------------------
266  double alpha_r2 = 0 ;
267  for (int k=0; k<mg->get_np(l_b); k++) {
268  for (int j=0; j<mg->get_nt(l_b); j++) {
269 
270  double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
271  double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
272 
273  // See Eq (100) from Gourgoulhon et al. (2001)
274  double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
275  / ( logn_auto_b - logn_auto_c ) ;
276 
277  if (alpha_r2_jk > alpha_r2) {
278  alpha_r2 = alpha_r2_jk ;
279  k_b = k ;
280  j_b = j ;
281  }
282 
283  }
284  }
285 
286  alpha_r = sqrt(alpha_r2) ;
287 
288  cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
289  << alpha_r << endl ;
290 
291  // New value of logn_auto
292  // ----------------------
293 
294  logn_auto = alpha_r2 * logn_auto ;
295  logn_auto_regu = alpha_r2 * logn_auto_regu ;
296  logn_auto_c = logn_auto()(0, 0, 0, 0) ;
297 
298  //------------------------------------------------------------
299  // Change the values of the inner points of the second domain
300  // by those of the outer points of the first domain
301  //------------------------------------------------------------
302 
303  (logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
304 
305  //--------------------------------------------
306  // First integral --> enthalpy in all space
307  // See Eq (98) from Gourgoulhon et al. (2001)
308  //--------------------------------------------
309 
310  ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
311 
312  //----------------------------------------------------------
313  // Change the enthalpy field to be set its maximum position
314  // at the coordinate center
315  //----------------------------------------------------------
316 
317  double dentdx = ent().dsdx()(0, 0, 0, 0) ;
318  double dentdy = ent().dsdy()(0, 0, 0, 0) ;
319 
320  cout << "dH/dx|_center = " << dentdx << endl ;
321  cout << "dH/dy|_center = " << dentdy << endl ;
322 
323  double dec_fact = 1. ;
324 
325  Tenseur func_in(mp) ;
326  func_in.set_etat_qcq() ;
327  func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x
328  - dec_fact * (dentdy/ent_c) * mp.y ;
329  func_in.set().annule(nzet, nz-1) ;
330  func_in.set_std_base() ;
331 
332  Tenseur func_ex(mp) ;
333  func_ex.set_etat_qcq() ;
334  func_ex.set() = 1. ;
335  func_ex.set().annule(0, nzet-1) ;
336  func_ex.set_std_base() ;
337 
338  // New enthalpy field
339  // ------------------
340  ent.set() = ent() * (func_in() + func_ex()) ;
341 
342  (ent().va).smooth(nzet, (ent.set()).va) ;
343 
344  double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
345  double dentdy_new = ent().dsdy()(0, 0, 0, 0) ;
346  cout << "dH/dx|_new = " << dentdx_new << endl ;
347  cout << "dH/dy|_new = " << dentdy_new << endl ;
348 
349  //-----------------------------------------------------
350  // Adaptation of the mapping to the new enthalpy field
351  //----------------------------------------------------
352 
353  // Shall the adaptation be performed (cusp) ?
354  // ------------------------------------------
355 
356  double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
357  double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
358  double rap_dent = fabs( dent_eq / dent_pole ) ;
359  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
360 
361  if ( rap_dent < thres_adapt ) {
362  adapt_flag = 0 ; // No adaptation of the mapping
363  cout << "******* FROZEN MAPPING *********" << endl ;
364  }
365  else{
366  adapt_flag = 1 ; // The adaptation of the mapping is to be
367  // performed
368  }
369 
370  ent_limit.set_etat_qcq() ;
371  for (int l=0; l<nzet; l++) { // loop on domains inside the star
372  ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
373  }
374 
375  ent_limit.set(nzet-1) = ent_b ;
376  Map_et mp_prev = mp_et ;
377 
378  mp.adapt(ent(), par_adapt) ;
379 
380  //----------------------------------------------------
381  // Computation of the enthalpy at the new grid points
382  //----------------------------------------------------
383 
384  mp_prev.homothetie(alpha_r) ;
385 
386  mp.reevaluate(&mp_prev, nzet+1, ent.set()) ;
387 
388  double fact_resize = 1. / alpha_r ;
389  for (int l=nzet; l<nz-1; l++) {
390  mp_et.resize(l, fact_resize) ;
391  }
392  mp_et.resize_extr(fact_resize) ;
393 
394  double ent_s_max = -1 ;
395  int k_s_max = -1 ;
396  int j_s_max = -1 ;
397  for (int k=0; k<mg->get_np(l_b); k++) {
398  for (int j=0; j<mg->get_nt(l_b); j++) {
399  double xx = fabs( ent()(l_b, k, j, i_b) ) ;
400  if (xx > ent_s_max) {
401  ent_s_max = xx ;
402  k_s_max = k ;
403  j_s_max = j ;
404  }
405  }
406  }
407  cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
408  << " and nzet : " << endl ;
409  cout << " " << ent_s_max << " reached for k = " << k_s_max
410  << " and j = " << j_s_max << endl ;
411 
412  //-------------------
413  // Equation of state
414  //-------------------
415 
416  equation_of_state() ; // computes new values for nbar (n), ener (e)
417  // and press (p) from the new ent (H)
418 
419  //----------------------------------------------------------
420  // Matter source terms in the gravitational field equations
421  //---------------------------------------------------------
422 
423  hydro_euler_extr(mass, sepa) ; // computes new values for
424  // ener_euler (E), s_euler (S)
425  // and u_euler (U^i)
426 
427 
428  //----------------------------------------------------
429  // Auxiliary terms for the source of Poisson equation
430  //----------------------------------------------------
431 
432  const Coord& xx = mp.x ;
433  const Coord& yy = mp.y ;
434  const Coord& zz = mp.z ;
435 
436  Tenseur r_bh(mp) ;
437  r_bh.set_etat_qcq() ;
438  r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
439  r_bh.set_std_base() ;
440 
441  Tenseur xx_cov(mp, 1, COV, ref_triad) ;
442  xx_cov.set_etat_qcq() ;
443  xx_cov.set(0) = xx + sepa ;
444  xx_cov.set(1) = yy ;
445  xx_cov.set(2) = zz ;
446  xx_cov.set_std_base() ;
447 
448  Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
449  xsr_cov = xx_cov / r_bh ;
450  xsr_cov.set_std_base() ;
451 
452  Tenseur msr(mp) ;
453  msr = ggrav * mass / r_bh ;
454  msr.set_std_base() ;
455 
456  Tenseur lapse_bh(mp) ;
457  lapse_bh = 1. / sqrt( 1.+2.*msr ) ;
458  lapse_bh.set_std_base() ;
459 
460  Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
461  lapse_bh2 = 1. / (1.+2.*msr) ;
462  lapse_bh2.set_std_base() ;
463 
464  Tenseur ldnu(mp) ;
465  ldnu = flat_scalar_prod_desal(xsr_cov, d_logn_auto) ;
466  ldnu.set_std_base() ;
467 
468  Tenseur ldbeta(mp) ;
469  ldbeta = flat_scalar_prod_desal(xsr_cov, d_beta_auto) ;
470  ldbeta.set_std_base() ;
471 
472  Tenseur lltkij(mp) ;
473  lltkij.set_etat_qcq() ;
474  lltkij.set() = 0 ;
475  lltkij.set_std_base() ;
476 
477  for (int i=0; i<3; i++)
478  for (int j=0; j<3; j++)
479  lltkij.set() += xsr_cov(i) * xsr_cov(j) * tkij_auto(i, j) ;
480 
481  Tenseur lshift(mp) ;
482  lshift = flat_scalar_prod_desal(xsr_cov, shift_auto) ;
483  lshift.set_std_base() ;
484 
485  Tenseur d_ldnu(mp, 1, COV, ref_triad) ;
486  d_ldnu = ldnu.gradient() ; // (d/dx, d/dy, d/dz)
487  d_ldnu.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
488 
489  Tenseur ldldnu(mp) ;
490  ldldnu = flat_scalar_prod_desal(xsr_cov, d_ldnu) ;
491  ldldnu.set_std_base() ;
492 
493  Tenseur d_ldbeta(mp, 1, COV, ref_triad) ;
494  d_ldbeta = ldbeta.gradient() ; // (d/dx, d/dy, d/dz)
495  d_ldbeta.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
496 
497  Tenseur ldldbeta(mp) ;
498  ldldbeta = flat_scalar_prod_desal(xsr_cov, d_ldbeta) ;
499  ldldbeta.set_std_base() ;
500 
501  //------------------------------------------
502  // Poisson equation for logn_auto (nu_auto)
503  //------------------------------------------
504 
505  // Source
506  // ------
507 
508  if (relativistic) {
509 
510  source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
512  + 2.*lapse_bh2 % msr % (ldnu % ldbeta + ldldnu)
513  + lapse_bh2 % lapse_bh2 % msr % (2.*(ldnu + 4.*msr % ldnu)
514  - ldbeta) / r_bh
515  - (4.*a_car % lapse_bh2 % lapse_bh2 % msr / 3. / nnn / r_bh)
516  * (2.+3.*msr) * (3.+4.*msr) * lltkij
517  + (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
518  / nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
519  + (4.*pow(lapse_bh2, 3.) % msr % msr / 3. / r_bh / r_bh)
520  % (2.*(a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
521  + (a_car - 1.) % pow(1.+3.*msr, 2.)
522  - 3.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
523 
524  }
525  else {
526  cout << "The computation of BH-NS binary systems"
527  << " should be relativistic !!!" << endl ;
528  abort() ;
529  }
530 
531  source.set_std_base() ;
532 
533  // Resolution of the Poisson equation
534  // ----------------------------------
535 
536  int k_falloff = 1 ;
537 
538  source().poisson_falloff(par_poisson1, logn_auto.set(), k_falloff) ;
539 
540  // Construct logn_auto_regu for et_bin_upmetr_extr.C
541  // -------------------------------------------------
542 
544 
545  // Check: has the Poisson equation been correctly solved ?
546  // -----------------------------------------------------
547 
548  Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
549 
550  cout <<
551  "Relative error in the resolution of the equation for logn_auto : "
552  << endl ;
553 
554  for (int l=0; l<nz; l++) {
555  cout << tdiff_logn(l) << " " ;
556  }
557  cout << endl ;
558  diff_logn = max(abs(tdiff_logn)) ;
559 
560  if (relativistic) {
561 
562  //--------------------------------
563  // Poisson equation for beta_auto
564  //--------------------------------
565 
566  // Source
567  // ------
568 
569  source = qpig * a_car % s_euler + 0.75 * akcar_auto
572  + lapse_bh2 % msr % (ldnu%ldnu + ldbeta%ldbeta + 2.*ldldbeta)
573  + lapse_bh2 % lapse_bh2 % msr % (2.*(1.+4.*msr) * ldbeta
574  - ldnu) / r_bh
575  - (a_car % lapse_bh2 % lapse_bh2 % msr / nnn / r_bh)
576  * (2.+3.*msr) * (3.+4.*msr) * lltkij
577  + (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
578  / nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
579  + (2.*pow(lapse_bh2, 3.) % msr % msr / r_bh / r_bh)
580  % ((a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
581  + (a_car - 1.) * pow(1.+3.*msr, 2.)
582  - 2.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
583 
584  source.set_std_base() ;
585 
586  // Resolution of the Poisson equation
587  // ----------------------------------
588 
589  source().poisson_falloff(par_poisson2, beta_auto.set(),
590  k_falloff) ;
591 
592  // Check: has the Poisson equation been correctly solved ?
593  // -----------------------------------------------------
594 
595  Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
596 
597  cout << "Relative error in the resolution of the equation for "
598  << "beta_auto : " << endl ;
599  for (int l=0; l<nz; l++) {
600  cout << tdiff_beta(l) << " " ;
601  }
602  cout << endl ;
603  diff_beta = max(abs(tdiff_beta)) ;
604 
605  //----------------------------------------
606  // Vector Poisson equation for shift_auto
607  //----------------------------------------
608 
609  // Some auxiliary terms for the source
610  // -----------------------------------
611 
612  Tenseur xx_con(mp, 1, CON, ref_triad) ;
613  xx_con.set_etat_qcq() ;
614  xx_con.set(0) = xx + sepa ;
615  xx_con.set(1) = yy ;
616  xx_con.set(2) = zz ;
617  xx_con.set_std_base() ;
618 
619  Tenseur xsr_con(mp, 1, CON, ref_triad) ;
620  xsr_con = xx_con / r_bh ;
621  xsr_con.set_std_base() ;
622 
623  // Components of shift_auto with respect to the Cartesian triad
624  // (d/dx, d/dy, d/dz) of the mapping :
625 
626  Tenseur shift_auto_local = shift_auto ;
627  shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
628 
629  // Gradient (partial derivatives with respect to the Cartesian
630  // coordinates of the mapping)
631  // dn(i, j) = D_i N^j
632 
633  Tenseur dn = shift_auto_local.gradient() ;
634 
635  // Return to the absolute reference frame
636  dn.change_triad(ref_triad) ;
637 
638  // Trace of D_i N^j = divergence of N^j :
639  Tenseur divn = contract(dn, 0, 1) ;
640 
641  // l^j D_j N^i
642  Tenseur ldn_con = contract(xsr_con, 0, dn, 0) ;
643 
644  // D_j (l^k D_k N^i): dldn(j, i)
645  Tenseur ldn_local = ldn_con ;
646  ldn_local.change_triad( mp.get_bvect_cart() ) ;
647  Tenseur dldn = ldn_local.gradient() ;
648  dldn.change_triad(ref_triad) ;
649 
650  // l^j D_j (l^k D_k N^i)
651  Tenseur ldldn = contract(xsr_con, 0, dldn, 0) ;
652 
653  // l_k D_j N^k
654  Tenseur ldn_cov = contract(xsr_cov, 0, dn, 1) ;
655 
656  // l^j l_k D_j N^k
657  Tenseur lldn_cov = contract(xsr_con, 0, ldn_cov, 0) ;
658 
659  // eta^{ij} l_k D_j N^k
660  Tenseur eldn(mp, 1, CON, ref_triad) ;
661  eldn.set_etat_qcq() ;
662  eldn.set(0) = ldn_cov(0) ;
663  eldn.set(1) = ldn_cov(1) ;
664  eldn.set(2) = ldn_cov(2) ;
665  eldn.set_std_base() ;
666 
667  // l^i D_j N^j
668  Tenseur ldivn = xsr_con % divn ;
669 
670  // D_j (l^i D_k N^k): dldivn(j, i)
671  Tenseur ldivn_local = ldivn ;
672  ldivn_local.change_triad( mp.get_bvect_cart() ) ;
673  Tenseur dldivn = ldivn_local.gradient() ;
674  dldivn.change_triad(ref_triad) ;
675 
676  // l^j D_j (l^i D_k N^k)
677  Tenseur ldldivn = contract(xsr_con, 0, dldivn, 0) ;
678 
679  // l_j N^j
680  Tenseur ln = contract(xsr_cov, 0, shift_auto, 0) ;
681 
682  Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
683 
684  Tenseur lvtmp = contract(xsr_con, 0, vtmp, 0) ;
685 
686  // eta^{ij} vtmp_j
687  Tenseur evtmp(mp, 1, CON, ref_triad) ;
688  evtmp.set_etat_qcq() ;
689  evtmp.set(0) = vtmp(0) ;
690  evtmp.set(1) = vtmp(1) ;
691  evtmp.set(2) = vtmp(2) ;
692  evtmp.set_std_base() ;
693 
694  // lapse_ns
695  Tenseur lapse_ns(mp) ;
696  lapse_ns = exp(logn_auto) ;
697  lapse_ns.set_std_base() ;
698 
699  // Source
700  // ------
701 
702  source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
703  % u_euler
705  - 2.*nnn % lapse_bh2 * msr / r_bh
707  + 2.*lapse_bh2 * msr * (3.*ldldn + ldldivn) / 3.
708  - lapse_bh2 * msr / r_bh
709  * (4.*ldivn - lapse_bh2 % (3.*ldn_con + 8.*msr * ldn_con)
710  - (eldn + 2.*lapse_bh2*(9.+11.*msr)*lldn_cov%xsr_con) / 3.)
711  - 2.*lapse_bh2 % lapse_bh2 * msr / r_bh / r_bh
712  * ( (4.+11.*msr) * shift_auto
713  - lapse_bh2 * (12.+51.*msr+46.*msr*msr) * ln % xsr_con )
714  / 3.
715  + 8.*pow(lapse_bh2, 4.) * msr / r_bh / r_bh
716  % (lapse_ns - 1.) * (2.+10.*msr+9.*msr*msr) * xsr_con / 3.
717  + 2.*pow(lapse_bh2, 3.) * msr / r_bh * (2.+3.*msr)
718  * ( (1.+2.*msr) * evtmp - (3.+2.*msr) * lvtmp * xsr_con) / 3. ;
719 
720  source_shift.set_std_base() ;
721 
722  // Resolution of the Poisson equation
723  // ----------------------------------
724 
725  // Filter for the source of shift vector :
726 
727  for (int i=0; i<3; i++) {
728  for (int l=0; l<nz; l++) {
729  if (source_shift(i).get_etat() != ETATZERO)
730  source_shift.set(i).filtre_phi(np_filter, l) ;
731  }
732  }
733 
734  // For Tenseur::poisson_vect, the triad must be the mapping
735  // triad, not the reference one:
736 
737  source_shift.change_triad( mp.get_bvect_cart() ) ;
738  /*
739  for (int i=0; i<3; i++) {
740  if(source_shift(i).dz_nonzero()) {
741  assert( source_shift(i).get_dzpuis() == 4 ) ;
742  }
743  else {
744  (source_shift.set(i)).set_dzpuis(4) ;
745  }
746  }
747 
748  source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
749  */
750  double lambda_shift = double(1) / double(3) ;
751 
752  int* shift_falloff ;
753  shift_falloff = new int[4] ;
754  shift_falloff[0] = 1 ;
755  shift_falloff[1] = 1 ;
756  shift_falloff[2] = 2 ;
757  shift_falloff[3] = 1 ;
758 
759  source_shift.poisson_vect_falloff(lambda_shift, par_poisson_vect,
761  khi_shift, shift_falloff) ;
762 
763  delete[] shift_falloff ;
764 
765  // Check: has the equation for shift_auto been correctly solved ?
766  // --------------------------------------------------------------
767 
768  // Divergence of shift_auto :
769  Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
770  // divna.dec2_dzpuis() ; // dzpuis 2 -> 0
771 
772  // Grad(div) :
773  Tenseur graddivn = divna.gradient() ;
774  // graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
775 
776  // Full operator :
777  Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
778  lap_shift.set_etat_qcq() ;
779  for (int i=0; i<3; i++) {
780  lap_shift.set(i) = shift_auto(i).laplacien()
781  + lambda_shift * graddivn(i) ;
782  }
783 
784  Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
785  Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
786  Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
787 
788  cout << "Relative error in the resolution of the equation for "
789  << "shift_auto : " << endl ;
790  cout << "x component : " ;
791  for (int l=0; l<nz; l++) {
792  cout << tdiff_shift_x(l) << " " ;
793  }
794  cout << endl ;
795  cout << "y component : " ;
796  for (int l=0; l<nz; l++) {
797  cout << tdiff_shift_y(l) << " " ;
798  }
799  cout << endl ;
800  cout << "z component : " ;
801  for (int l=0; l<nz; l++) {
802  cout << tdiff_shift_z(l) << " " ;
803  }
804  cout << endl ;
805 
806  diff_shift_x = max(abs(tdiff_shift_x)) ;
807  diff_shift_y = max(abs(tdiff_shift_y)) ;
808  diff_shift_z = max(abs(tdiff_shift_z)) ;
809 
810  // Final result
811  // ------------
812  // The output of Tenseur::poisson_vect_falloff is on the mapping
813  // triad, it should therefore be transformed to components on the
814  // reference triad :
815 
817 
818  } // End of relativistic equations
819 
820  //------------------------------
821  // Relative change in enthalpy
822  //------------------------------
823 
824  Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
825  diff_ent = diff_ent_tbl(0) ;
826  for (int l=1; l<nzet; l++) {
827  diff_ent += diff_ent_tbl(l) ;
828  }
829  diff_ent /= nzet ;
830 
831  ent_jm1 = ent ;
832 
833  } // End of main loop
834 
835  //========================================================//
836  // End of iteration //
837  //========================================================//
838 
839 }
840 
841 
842 void Et_bin_bhns_extr::equil_bhns_extr_cf(double ent_c, const double& mass,
843  const double& sepa, int mermax,
844  int mermax_poisson,
845  double relax_poisson,
846  int mermax_potvit,
847  double relax_potvit, int np_filter,
848  double thres_adapt,
849  Tbl& diff) {
850 
851  // Fundamental constants and units
852  // -------------------------------
853  using namespace Unites ;
854 
855  assert( kerrschild == false ) ;
856 
857  // Initializations
858  // ---------------
859 
860  const Mg3d* mg = mp.get_mg() ;
861  int nz = mg->get_nzone() ; // total number of domains
862 
863  // The following is required to initialize mp_prev as a Map_et:
864  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
865 
866  // Domain and radial indices of points at the surface of the star:
867  int l_b = nzet - 1 ;
868  int i_b = mg->get_nr(l_b) - 1 ;
869  int k_b ;
870  int j_b ;
871 
872  // Value of the enthalpy defining the surface of the star
873  double ent_b = 0 ;
874 
875  // Error indicators
876  // ----------------
877 
878  double& diff_ent = diff.set(0) ;
879  double& diff_vel_pot = diff.set(1) ;
880  double& diff_logn = diff.set(2) ;
881  double& diff_beta = diff.set(3) ;
882  double& diff_shift_x = diff.set(4) ;
883  double& diff_shift_y = diff.set(5) ;
884  double& diff_shift_z = diff.set(6) ;
885 
886  // Parameters for the function Map_et::adapt
887  // -----------------------------------------
888 
889  Param par_adapt ;
890  int nitermax = 100 ;
891  int niter ;
892  int adapt_flag = 1 ; // 1 = performs the full computation,
893  // 0 = performs only the rescaling by
894  // the factor alpha_r
895  int nz_search = nzet ; // Number of domains for searching
896  // the enthalpy isosurfaces
897  double precis_secant = 1.e-14 ;
898  double alpha_r ;
899  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
900 
901  Tbl ent_limit(nz) ;
902 
903  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
904  // locate zeros by the secant method
905  par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
906  // to the isosurfaces of ent is to be
907  // performed
908  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
909  // the enthalpy isosurface
910  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
911  // 0 = performs only the rescaling by
912  // the factor alpha_r
913  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
914  // (theta_*, phi_*)
915  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
916  // (theta_*, phi_*)
917  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
918  // used in the secant method
919  par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
920  // the determination of zeros by
921  // the secant method
922  par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
923  // 0 = contracting mapping
924  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
925  // distances will be multiplied
926  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
927  // to define the isosurfaces
928 
929  // Parameters for the function Map_et::poisson for logn_auto
930  // ---------------------------------------------------------
931 
932  double precis_poisson = 1.e-16 ;
933 
934  Param par_poisson1 ;
935 
936  par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
937  par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
938  par_poisson1.add_double(precis_poisson, 1) ; // required precision
939  par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
940  // used
941  par_poisson1.add_cmp_mod( ssjm1_logn ) ;
942 
943  // Parameters for the function Map_et::poisson for beta_auto
944  // ---------------------------------------------------------
945 
946  Param par_poisson2 ;
947 
948  par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
949  par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
950  par_poisson2.add_double(precis_poisson, 1) ; // required precision
951  par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
952  // used
953  par_poisson2.add_cmp_mod( ssjm1_beta ) ;
954 
955  // Parameters for the function Tenseur::poisson_vect
956  // -------------------------------------------------
957 
958  Param par_poisson_vect ;
959 
960  par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
961  // iterations
962  par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
963  par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
964  par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
965  par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
966  par_poisson_vect.add_int_mod(niter, 0) ;
967 
968  // External potential
969  // See Eq (99) from Gourgoulhon et al. (2001)
970  // -----------------------------------------
971 
972  Tenseur pot_ext = logn_comp + pot_centri + loggam ;
973 
974  Tenseur ent_jm1 = ent ; // Enthalpy at previous step
975 
976  Tenseur source(mp) ; // source term in the equation for logn_auto
977  // and beta_auto
978 
979  Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
980  // equation for shift_auto
981 
982  //==========================================================//
983  // Start of iteration //
984  //==========================================================//
985 
986  for(int mer=0 ; mer<mermax ; mer++ ) {
987 
988  cout << "-----------------------------------------------" << endl ;
989  cout << "step: " << mer << endl ;
990  cout << "diff_ent = " << diff_ent << endl ;
991 
992  //------------------------------------------------------
993  // Resolution of the elliptic equation for the velocity
994  // scalar potential
995  //------------------------------------------------------
996 
997  if (irrotational) {
998  diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
999  precis_poisson, relax_potvit) ;
1000  }
1001 
1002  //-------------------------------------
1003  // Computation of the new radial scale
1004  //-------------------------------------
1005 
1006  // alpha_r (r = alpha_r r') is determined so that the enthalpy
1007  // takes the requested value ent_b at the stellar surface
1008 
1009  // Values at the center of the star:
1010  double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
1011  double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
1012 
1013  // Search for the reference point (theta_*, phi_*)
1014  // [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
1015  // at the surface of the star
1016  // ------------------------------------------------------------------
1017  double alpha_r2 = 0 ;
1018  for (int k=0; k<mg->get_np(l_b); k++) {
1019  for (int j=0; j<mg->get_nt(l_b); j++) {
1020 
1021  double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
1022  double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
1023 
1024  // See Eq (100) from Gourgoulhon et al. (2001)
1025  double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
1026  / ( logn_auto_b - logn_auto_c ) ;
1027 
1028  if (alpha_r2_jk > alpha_r2) {
1029  alpha_r2 = alpha_r2_jk ;
1030  k_b = k ;
1031  j_b = j ;
1032  }
1033 
1034  }
1035  }
1036 
1037  alpha_r = sqrt(alpha_r2) ;
1038 
1039  cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
1040  << alpha_r << endl ;
1041 
1042  // New value of logn_auto
1043  // ----------------------
1044 
1045  logn_auto = alpha_r2 * logn_auto ;
1046  logn_auto_regu = alpha_r2 * logn_auto_regu ;
1047  logn_auto_c = logn_auto()(0, 0, 0, 0) ;
1048 
1049  //------------------------------------------------------------
1050  // Change the values of the inner points of the second domain
1051  // by those of the outer points of the first domain
1052  //------------------------------------------------------------
1053 
1054  (logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
1055 
1056  //--------------------------------------------
1057  // First integral --> enthalpy in all space
1058  // See Eq (98) from Gourgoulhon et al. (2001)
1059  //--------------------------------------------
1060 
1061  ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
1062 
1063  //---------------------------------------------------------
1064  // Change the enthalpy field to accelerate the convergence
1065  //---------------------------------------------------------
1066 
1067  double dentdx = ent().dsdx()(0, 0, 0, 0) ;
1068  double dentdy = ent().dsdy()(0, 0, 0, 0) ;
1069 
1070  cout << "dH/dx|_center = " << dentdx << endl ;
1071  cout << "dH/dy|_center = " << dentdy << endl ;
1072 
1073  double dec_fact = 1. ;
1074 
1075  Tenseur func_in(mp) ;
1076  func_in.set_etat_qcq() ;
1077  func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x ;
1078  func_in.set().annule(nzet, nz-1) ;
1079  func_in.set_std_base() ;
1080 
1081  Tenseur func_ex(mp) ;
1082  func_ex.set_etat_qcq() ;
1083  func_ex.set() = 1. ;
1084  func_ex.set().annule(0, nzet-1) ;
1085  func_ex.set_std_base() ;
1086 
1087  // New enthalpy field
1088  // ------------------
1089  ent.set() = ent() * (func_in() + func_ex()) ;
1090 
1091  (ent().va).smooth(nzet, (ent.set()).va) ;
1092 
1093  double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
1094 
1095  cout << "dH/dx|_new = " << dentdx_new << endl ;
1096 
1097  //-----------------------------------------------------
1098  // Adaptation of the mapping to the new enthalpy field
1099  //----------------------------------------------------
1100 
1101  // Shall the adaptation be performed (cusp) ?
1102  // ------------------------------------------
1103 
1104  double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
1105  double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
1106  double rap_dent = fabs( dent_eq / dent_pole ) ;
1107  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
1108 
1109  if ( rap_dent < thres_adapt ) {
1110  adapt_flag = 0 ; // No adaptation of the mapping
1111  cout << "******* FROZEN MAPPING *********" << endl ;
1112  }
1113  else{
1114  adapt_flag = 1 ; // The adaptation of the mapping is to be
1115  // performed
1116  }
1117 
1118  ent_limit.set_etat_qcq() ;
1119  for (int l=0; l<nzet; l++) { // loop on domains inside the star
1120  ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
1121  }
1122 
1123  ent_limit.set(nzet-1) = ent_b ;
1124  Map_et mp_prev = mp_et ;
1125 
1126  mp.adapt(ent(), par_adapt) ;
1127 
1128  //----------------------------------------------------
1129  // Computation of the enthalpy at the new grid points
1130  //----------------------------------------------------
1131 
1132  mp_prev.homothetie(alpha_r) ;
1133 
1134  mp.reevaluate_symy(&mp_prev, nzet+1, ent.set()) ;
1135 
1136  double fact_resize = 1. / alpha_r ;
1137  for (int l=nzet; l<nz-1; l++) {
1138  mp_et.resize(l, fact_resize) ;
1139  }
1140  mp_et.resize_extr(fact_resize) ;
1141 
1142  double ent_s_max = -1 ;
1143  int k_s_max = -1 ;
1144  int j_s_max = -1 ;
1145  for (int k=0; k<mg->get_np(l_b); k++) {
1146  for (int j=0; j<mg->get_nt(l_b); j++) {
1147  double xx = fabs( ent()(l_b, k, j, i_b) ) ;
1148  if (xx > ent_s_max) {
1149  ent_s_max = xx ;
1150  k_s_max = k ;
1151  j_s_max = j ;
1152  }
1153  }
1154  }
1155  cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
1156  << " and nzet : " << endl ;
1157  cout << " " << ent_s_max << " reached for k = " << k_s_max
1158  << " and j = " << j_s_max << endl ;
1159 
1160  //-------------------
1161  // Equation of state
1162  //-------------------
1163 
1164  equation_of_state() ; // computes new values for nbar (n), ener (e)
1165  // and press (p) from the new ent (H)
1166 
1167  //----------------------------------------------------------
1168  // Matter source terms in the gravitational field equations
1169  //---------------------------------------------------------
1170 
1171  hydro_euler_extr(mass, sepa) ; // computes new values for
1172  // ener_euler (E), s_euler (S)
1173  // and u_euler (U^i)
1174 
1175 
1176  //----------------------------------------------------
1177  // Auxiliary terms for the source of Poisson equation
1178  //----------------------------------------------------
1179 
1180  const Coord& xx = mp.x ;
1181  const Coord& yy = mp.y ;
1182  const Coord& zz = mp.z ;
1183 
1184  Tenseur r_bh(mp) ;
1185  r_bh.set_etat_qcq() ;
1186  r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
1187  r_bh.set_std_base() ;
1188 
1189  Tenseur xx_cov(mp, 1, COV, ref_triad) ;
1190  xx_cov.set_etat_qcq() ;
1191  xx_cov.set(0) = xx + sepa ;
1192  xx_cov.set(1) = yy ;
1193  xx_cov.set(2) = zz ;
1194  xx_cov.set_std_base() ;
1195 
1196  Tenseur msr(mp) ;
1197  msr = ggrav * mass / r_bh ;
1198  msr.set_std_base() ;
1199 
1200  Tenseur tmp(mp) ;
1201  tmp = 1. / ( 1. - 0.25*msr*msr ) ;
1202  tmp.set_std_base() ;
1203 
1204  Tenseur xdnu(mp) ;
1205  xdnu = flat_scalar_prod_desal(xx_cov, d_logn_auto) ;
1206  xdnu.set_std_base() ;
1207 
1208  Tenseur xdbeta(mp) ;
1209  xdbeta = flat_scalar_prod_desal(xx_cov, d_beta_auto) ;
1210  xdbeta.set_std_base() ;
1211 
1212  //------------------------------------------
1213  // Poisson equation for logn_auto (nu_auto)
1214  //------------------------------------------
1215 
1216  // Source
1217  // ------
1218 
1219  if (relativistic) {
1220 
1221  source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
1223  - 0.5 * tmp % msr % msr % xdnu / r_bh / r_bh
1224  - tmp % msr % xdbeta / r_bh / r_bh ;
1225 
1226  }
1227  else {
1228  cout << "The computation of BH-NS binary systems"
1229  << " should be relativistic !!!" << endl ;
1230  abort() ;
1231  }
1232 
1233  source.set_std_base() ;
1234 
1235  // Resolution of the Poisson equation
1236  // ----------------------------------
1237 
1238  int k_falloff = 1 ;
1239 
1240  source().poisson_falloff(par_poisson1, logn_auto.set(), k_falloff) ;
1241 
1242  // Construct logn_auto_regu for et_bin_upmetr_extr.C
1243  // -------------------------------------------------
1244 
1246 
1247  // Check: has the Poisson equation been correctly solved ?
1248  // -----------------------------------------------------
1249 
1250  Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
1251 
1252  cout <<
1253  "Relative error in the resolution of the equation for logn_auto : "
1254  << endl ;
1255 
1256  for (int l=0; l<nz; l++) {
1257  cout << tdiff_logn(l) << " " ;
1258  }
1259  cout << endl ;
1260  diff_logn = max(abs(tdiff_logn)) ;
1261 
1262  if (relativistic) {
1263 
1264  //--------------------------------
1265  // Poisson equation for beta_auto
1266  //--------------------------------
1267 
1268  // Source
1269  // ------
1270 
1271  source = qpig * a_car % s_euler + 0.75 * akcar_auto
1274  - tmp % msr % xdnu / r_bh / r_bh
1275  - 0.5 * tmp % msr %msr % xdbeta / r_bh / r_bh ;
1276 
1277  source.set_std_base() ;
1278 
1279  // Resolution of the Poisson equation
1280  // ----------------------------------
1281 
1282  source().poisson_falloff(par_poisson2, beta_auto.set(),
1283  k_falloff) ;
1284 
1285  // Check: has the Poisson equation been correctly solved ?
1286  // -----------------------------------------------------
1287 
1288  Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
1289 
1290  cout << "Relative error in the resolution of the equation for "
1291  << "beta_auto : " << endl ;
1292  for (int l=0; l<nz; l++) {
1293  cout << tdiff_beta(l) << " " ;
1294  }
1295  cout << endl ;
1296  diff_beta = max(abs(tdiff_beta)) ;
1297 
1298  //----------------------------------------
1299  // Vector Poisson equation for shift_auto
1300  //----------------------------------------
1301 
1302  // Some auxiliary terms for the source
1303  // -----------------------------------
1304 
1305  Tenseur bhtmp(mp, 1, COV, ref_triad) ;
1306  bhtmp.set_etat_qcq() ;
1307  bhtmp = tmp % msr % (3.*msr-8.) % xx_cov / r_bh / r_bh ;
1308  bhtmp.set_std_base() ;
1309 
1310  Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
1311 
1312  // Source
1313  // ------
1314 
1315  source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
1316  % u_euler
1317  + nnn % flat_scalar_prod_desal(tkij_auto, vtmp+bhtmp) ;
1318 
1319  source_shift.set_std_base() ;
1320 
1321  // Resolution of the Poisson equation
1322  // ----------------------------------
1323 
1324  // Filter for the source of shift vector :
1325 
1326  for (int i=0; i<3; i++) {
1327  for (int l=0; l<nz; l++) {
1328  if (source_shift(i).get_etat() != ETATZERO)
1329  source_shift.set(i).filtre_phi(np_filter, l) ;
1330  }
1331  }
1332 
1333  // For Tenseur::poisson_vect, the triad must be the mapping
1334  // triad, not the reference one:
1335 
1336  source_shift.change_triad( mp.get_bvect_cart() ) ;
1337  /*
1338  for (int i=0; i<3; i++) {
1339  if(source_shift(i).dz_nonzero()) {
1340  assert( source_shift(i).get_dzpuis() == 4 ) ;
1341  }
1342  else {
1343  (source_shift.set(i)).set_dzpuis(4) ;
1344  }
1345  }
1346 
1347  source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
1348  */
1349  double lambda_shift = double(1) / double(3) ;
1350 
1351  int* shift_falloff ;
1352  shift_falloff = new int[4] ;
1353  shift_falloff[0] = 1 ;
1354  shift_falloff[1] = 1 ;
1355  shift_falloff[2] = 2 ;
1356  shift_falloff[3] = 1 ;
1357 
1358  source_shift.poisson_vect_falloff(lambda_shift, par_poisson_vect,
1360  khi_shift, shift_falloff) ;
1361 
1362  delete[] shift_falloff ;
1363 
1364  // Check: has the equation for shift_auto been correctly solved ?
1365  // --------------------------------------------------------------
1366 
1367  // Divergence of shift_auto :
1368  Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
1369  // divna.dec2_dzpuis() ; // dzpuis 2 -> 0
1370 
1371  // Grad(div) :
1372  Tenseur graddivn = divna.gradient() ;
1373  // graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
1374 
1375  // Full operator :
1376  Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
1377  lap_shift.set_etat_qcq() ;
1378  for (int i=0; i<3; i++) {
1379  lap_shift.set(i) = shift_auto(i).laplacien()
1380  + lambda_shift * graddivn(i) ;
1381  }
1382 
1383  Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
1384  Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
1385  Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
1386 
1387  cout << "Relative error in the resolution of the equation for "
1388  << "shift_auto : " << endl ;
1389  cout << "x component : " ;
1390  for (int l=0; l<nz; l++) {
1391  cout << tdiff_shift_x(l) << " " ;
1392  }
1393  cout << endl ;
1394  cout << "y component : " ;
1395  for (int l=0; l<nz; l++) {
1396  cout << tdiff_shift_y(l) << " " ;
1397  }
1398  cout << endl ;
1399  cout << "z component : " ;
1400  for (int l=0; l<nz; l++) {
1401  cout << tdiff_shift_z(l) << " " ;
1402  }
1403  cout << endl ;
1404 
1405  diff_shift_x = max(abs(tdiff_shift_x)) ;
1406  diff_shift_y = max(abs(tdiff_shift_y)) ;
1407  diff_shift_z = max(abs(tdiff_shift_z)) ;
1408 
1409  // Final result
1410  // ------------
1411  // The output of Tenseur::poisson_vect_falloff is on the mapping
1412  // triad, it should therefore be transformed to components on the
1413  // reference triad :
1414 
1416 
1417  } // End of relativistic equations
1418 
1419  //------------------------------
1420  // Relative change in enthalpy
1421  //------------------------------
1422 
1423  Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
1424  diff_ent = diff_ent_tbl(0) ;
1425  for (int l=1; l<nzet; l++) {
1426  diff_ent += diff_ent_tbl(l) ;
1427  }
1428  diff_ent /= nzet ;
1429 
1430  ent_jm1 = ent ;
1431 
1432  } // End of main loop
1433 
1434  //========================================================//
1435  // End of iteration //
1436  //========================================================//
1437 
1438 }
1439 
1440 }
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1145
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur &#39;s...
Definition: etoile.h:831
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
Radial mapping of rather general form.
Definition: map.h:2783
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
void equil_bhns_extr_cf(double ent_c, const double &mass, const double &sepa, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the conformally flat background met...
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1186
Tenseur pot_centri
Centrifugal potential.
Definition: etoile.h:956
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Definition: etoile.h:962
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star...
Definition: etoile.h:494
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tenseur nnn
Total lapse function.
Definition: etoile.h:512
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:471
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Tenseur press
Fluid pressure.
Definition: etoile.h:464
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Definition: etoile.h:882
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
Definition: cmp_manip.C:104
void hydro_euler_extr(const double &mass, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition: etoile.h:892
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
bool irrotational
true for an irrotational star, false for a corotating one
Definition: etoile.h:825
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata&#39;s prescription [Prog...
Definition: etoile.h:911
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:477
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:684
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Definition: etoile.h:857
double velocity_pot_extr(const double &mass, const double &sepa, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Definition: etoile.h:862
void equil_bhns_extr_ks(double ent_c, const double &mass, const double &sepa, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the Kerr-Schild background metric u...
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
Definition: map_et.C:951
Parameter storage.
Definition: param.h:125
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:525
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:569
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:928
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata&#39;s prescription [Prog...
Definition: etoile.h:921
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition: etoile.h:928
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:435
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Tenseur a_car
Total conformal factor .
Definition: etoile.h:518
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:440
Multi-domain grid.
Definition: grilles.h:279
double ray_pole() const
Coordinate radius at [r_unit].
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Definition: etoile.h:968
virtual void reevaluate_symy(const Map *mp_prev, int nzet, Cmp &uu) const =0
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition: etoile.h:487
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: etoile.h:852
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:460
Coord y
y coordinate centered on the grid
Definition: map.h:745
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
virtual void reevaluate(const Map *mp_prev, int nzet, Cmp &uu) const =0
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
void resize_extr(double lambda)
Rescales the outer boundary of the outermost domain in the case of non-compactified external domain...
Coord x
x coordinate centered on the grid
Definition: map.h:744
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition: etoile.h:509
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1007
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition: etoile.h:941
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:468
Coord z
z coordinate centered on the grid
Definition: map.h:746
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1558
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: etoile.h:976
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition: etoile.h:986