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