LORENE
star_bhns_equilibrium.C
1 /*
2  * Method of class Star_bhns to compute an equilibrium configuration
3  * of a neutron star in a black hole-neutron star binary
4  *
5  * (see file star_bhns.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2005-2007 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 star_bhns_equilibrium_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $" ;
30 
31 /*
32  * $Id: star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $
33  * $Log: star_bhns_equilibrium.C,v $
34  * Revision 1.4 2014/10/13 08:53:40 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.3 2014/10/06 15:13:16 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.2 2008/05/15 19:13:45 k_taniguchi
41  * Change of some parameters.
42  *
43  * Revision 1.1 2007/06/22 01:30:45 k_taniguchi
44  * *** empty log message ***
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $
48  *
49  */
50 
51 // C++ headers
52 //#include <>
53 
54 // C headers
55 #include <cmath>
56 
57 // Lorene headers
58 #include "star_bhns.h"
59 #include "param.h"
60 #include "cmp.h"
61 #include "tenseur.h"
62 #include "utilitaires.h"
63 #include "unites.h"
64 //#include "graphique.h"
65 
66 namespace Lorene {
67 void Star_bhns::equilibrium_bhns(double ent_c, const double& mass_bh,
68  const double& sepa, bool kerrschild,
69  int mer, int mermax_ns, int mermax_potvit,
70  int mermax_poisson, int filter_r,
71  int filter_r_s, int filter_p_s,
72  double relax_poisson, double relax_potvit,
73  double thres_adapt, double resize_ns,
74  const Tbl& fact_resize, Tbl& diff) {
75 
76  // Fundamental constants and units
77  // -------------------------------
78  using namespace Unites ;
79 
80  // Initializations
81  // ---------------
82 
83  const Mg3d* mg = mp.get_mg() ;
84  int nz = mg->get_nzone() ; // Total number of domain
85  // int nt = mg->get_nt(0) ;
86 
87  // The following is required to initialize mp_prev as a Map_et:
88  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
89 
90  // Domain and radial indices of points at the surface of the star:
91  int l_b = nzet - 1 ;
92  int i_b = mg->get_nr(l_b) - 1 ;
93  int k_b ;
94  int j_b ;
95 
96  // Value of the enthalpy defining the surface of the star
97  double ent_b = 0 ;
98 
99  // Error indicators
100  // ----------------
101 
102  double& diff_ent = diff.set(0) ;
103  double& diff_vel_pot = diff.set(1) ;
104  double& diff_lapconf = diff.set(2) ;
105  double& diff_confo = diff.set(3) ;
106  double& diff_shift_x = diff.set(4) ;
107  double& diff_shift_y = diff.set(5) ;
108  double& diff_shift_z = diff.set(6) ;
109  double& diff_dHdr = diff.set(7) ;
110  double& diff_dHdr_min = diff.set(8) ;
111  double& diff_phi_min = diff.set(9) ;
112  double& diff_radius = diff.set(10) ;
113 
114  // Parameters for the function Map_et::adapt
115  // -----------------------------------------
116 
117  Param par_adapt ;
118  int nitermax = 100 ;
119  int niter ;
120  int adapt_flag = 1 ; // 1 = performs the full computation,
121  // 0 = performs only the rescaling by
122  // the factor alpha_r
123 
124  // Number of domains for searching the enthalpy isosurfaces
125  int nz_search = nzet + 1 ;
126  // int nz_search = nzet ;
127 
128  double precis_secant = 1.e-14 ;
129  double alpha_r ;
130  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
131 
132  Tbl ent_limit(nz) ;
133 
134  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
135  // locate zeros by the secant method
136  par_adapt.add_int(nzet, 1) ; // number of domains where the
137  // adjustment to the isosurfaces of
138  // ent is to be performed
139  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
140  // the enthalpy isosurface
141  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
142  // 0 = performs only the rescaling by
143  // the factor alpha_r
144  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
145  // (theta_*, phi_*)
146  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
147  // (theta_*, phi_*)
148  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used
149  // in the secant method
150  par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
151  // the determination of zeros by
152  // the secant method
153  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping,
154  // 0 = contracting mapping
155  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
156  // distances will be multiplied
157  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
158  // to define the isosurfaces
159 
160  Cmp ssjm1lapconf(ssjm1_lapconf) ;
161  Cmp ssjm1confo(ssjm1_confo) ;
162 
163  ssjm1lapconf.set_etat_qcq() ;
164  ssjm1confo.set_etat_qcq() ;
165 
166  double precis_poisson = 1.e-14 ;
167 
168  // Parameters for the function Scalar::poisson for lapconf_auto
169  // ------------------------------------------------------------
170 
171  Param par_lapconf ;
172 
173  par_lapconf.add_int(mermax_poisson, 0) ; // maximum number of iterations
174  par_lapconf.add_double(relax_poisson, 0) ; // relaxation parameter
175  par_lapconf.add_double(precis_poisson, 1) ; // required precision
176  par_lapconf.add_int_mod(niter, 0) ; // number of iterations actually used
177  par_lapconf.add_cmp_mod( ssjm1lapconf ) ;
178 
179  // Parameters for the function Scalar::poisson for confo_auto
180  // ----------------------------------------------------------
181 
182  Param par_confo ;
183 
184  par_confo.add_int(mermax_poisson, 0) ; // maximum number of iterations
185  par_confo.add_double(relax_poisson, 0) ; // relaxation parameter
186  par_confo.add_double(precis_poisson, 1) ; // required precision
187  par_confo.add_int_mod(niter, 0) ; // number of iterations actually used
188  par_confo.add_cmp_mod( ssjm1confo ) ;
189 
190  // Parameters for the function Vector::poisson for shift method 2
191  // --------------------------------------------------------------
192 
193  Param par_shift2 ;
194 
195  par_shift2.add_int(mermax_poisson, 0) ; // maximum number of iterations
196  par_shift2.add_double(relax_poisson, 0) ; // relaxation parameter
197  par_shift2.add_double(precis_poisson, 1) ; // required precision
198  par_shift2.add_int_mod(niter, 0) ; // number of iterations actually used
199 
200  Cmp ssjm1khi(ssjm1_khi) ;
201  ssjm1khi.set_etat_qcq() ;
202 
203  Tenseur ssjm1wshift(mp, 1, CON, mp.get_bvect_cart()) ;
204  ssjm1wshift.set_etat_qcq() ;
205  for (int i=0; i<3; i++) {
206  ssjm1wshift.set(i) = Cmp(ssjm1_wshift(i+1)) ;
207  }
208 
209  par_shift2.add_cmp_mod(ssjm1khi) ;
210  par_shift2.add_tenseur_mod(ssjm1wshift) ;
211 
212  Scalar ent_jm1 = ent ; // Enthalpy at previous step
213 
214  Scalar lapconf_m1(mp) ; // = lapconf_auto - 0.5
215  Scalar confo_m1(mp) ; // = confo_auto - 0.5
216 
217  Scalar source_lapconf(mp) ; // Source term in the equation for lapconf_auto
218  source_lapconf.set_etat_qcq() ;
219  Scalar source_confo(mp) ; // Source term in the equation for confo_auto
220  source_confo.set_etat_qcq() ;
221  Vector source_shift(mp, CON, mp.get_bvect_cart()) ; // Source term
222  // in the equation for shift_auto
223  source_shift.set_etat_qcq() ;
224 
225  //===========================================================
226  // Start of iteration
227  //===========================================================
228 
229  for (int mer_ns=0; mer_ns<mermax_ns; mer_ns++) {
230 
231  cout << "-----------------------------------------------" << endl ;
232  cout << "step: " << mer_ns << endl ;
233  cout << "diff_ent = " << diff_ent << endl ;
234 
235  //------------------------------------------------------
236  // Resolution of the elliptic equation for the velocity
237  // scalar potential
238  //------------------------------------------------------
239 
240  if (irrotational) {
241  diff_vel_pot = velo_pot_bhns(mass_bh, sepa, kerrschild,
242  mermax_potvit, precis_poisson,
243  relax_potvit) ;
244  }
245  else {
246  diff_vel_pot = 0. ; // to avoid the warning
247  }
248 
249  //-------------------------------------
250  // Computation of the new radial scale
251  //-------------------------------------
252 
253  // alpha_r (r = alpha_r r') is determined so that the enthalpy
254  // takes the requested value ent_b at the stellar surface
255 
256  // Values at the center of the star:
257  // lapconf_auto = lapconf_auto=0(r->infty) + 0.5
258  // The term lapconf_auto=0(r->infty) should be rescaled.
259  double lapconf_auto_c = lapconf_auto.val_grid_point(0,0,0,0) - 0.5 ;
260  double lapconf_comp_c = lapconf_comp.val_grid_point(0,0,0,0) ;
261 
262  double confo_c = confo_tot.val_grid_point(0,0,0,0) ;
263 
264  double gam_c = gam.val_grid_point(0,0,0,0) ;
265  double gam0_c = gam0.val_grid_point(0,0,0,0) ;
266 
267  double hhh_c = exp(ent_c) ;
268  double hhh_b = exp(ent_b) ;
269 
270  // Search for the reference point (theta_*, phi_*) [notation of
271  // Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
272  // at the surface of the star
273  // ------------------------------------------------------------
274  double alpha_r2 = 0. ;
275 
276  for (int k=0; k<mg->get_np(l_b); k++) {
277  for (int j=0; j<mg->get_nt(l_b); j++) {
278 
279  double lapconf_auto_b =
280  lapconf_auto.val_grid_point(l_b,k,j,i_b) - 0.5 ;
281  double lapconf_comp_b =
282  lapconf_comp.val_grid_point(l_b,k,j,i_b) ;
283 
284  double confo_b = confo_tot.val_grid_point(l_b,k,j,i_b) ;
285 
286  double gam_b = gam.val_grid_point(l_b,k,j,i_b) ;
287  double gam0_b = gam0.val_grid_point(l_b,k,j,i_b) ;
288 
289  double aaa = (gam0_c*gam_b*hhh_b*confo_c)
290  / (gam0_b*gam_c*hhh_c*confo_b) ;
291 
292  // See Eq (100) from Gourgoulhon et al. (2001)
293  double alpha_r2_jk = (aaa * lapconf_comp_b - lapconf_comp_c
294  + 0.5 * (aaa - 1.))
295  / (lapconf_auto_c - aaa * lapconf_auto_b ) ;
296 
297  if (alpha_r2_jk > alpha_r2) {
298  alpha_r2 = alpha_r2_jk ;
299  k_b = k ;
300  j_b = j ;
301  }
302  }
303  }
304 
305  alpha_r = sqrt(alpha_r2) ;
306 
307  cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
308  << alpha_r << endl ;
309 
310  // New value of lapconf_auto
311  // -------------------------
312 
313  lapconf_auto = alpha_r2 * (lapconf_auto - 0.5) + 0.5 ;
314  Scalar lapconf_tot_tmp = lapconf_auto + lapconf_comp ;
315  lapconf_tot_tmp.std_spectral_base() ;
316 
317  /*
318  confo_auto = alpha_r2 * (confo_auto - 0.5) + 0.5 ;
319  Scalar confo_tot_tmp = confo_auto + confo_comp ;
320  confo_tot_tmp.std_spectral_base() ;
321  */
322  //------------------------------------------------------------
323  // Change the values of the inner points of the second domain
324  // by those of the outer points of the first domain
325  //------------------------------------------------------------
326 
328 
329  //--------------------------------------------
330  // First integral --> enthalpy in all space
331  // See Eq (98) from Gourgoulhon et al. (2001)
332  //--------------------------------------------
333 
334  double log_lapconf_c = log(lapconf_tot_tmp.val_grid_point(0,0,0,0)) ;
335  double log_confo_c = log(confo_tot.val_grid_point(0,0,0,0)) ;
336  double loggam_c = loggam.val_grid_point(0,0,0,0) ;
337  double pot_centri_c = pot_centri.val_grid_point(0,0,0,0) ;
338 
339  ent = (ent_c + log_lapconf_c - log_confo_c + loggam_c + pot_centri_c)
340  - log(lapconf_tot_tmp) + log(confo_tot) - loggam - pot_centri ;
342 
343 
344  //----------------------------------------------------------
345  // Change the enthalpy field to be set its maximum position
346  // at the coordinate center
347  //----------------------------------------------------------
348 
349  double dentdx = ent.dsdx().val_grid_point(0,0,0,0) ;
350  double dentdy = ent.dsdy().val_grid_point(0,0,0,0) ;
351 
352  cout << "dH/dx|_center = " << dentdx << endl ;
353  cout << "dH/dy|_center = " << dentdy << endl ;
354 
355  double dec_fact_x = 1. ;
356  double dec_fact_y = 1. ;
357 
358  Scalar func_in(mp) ;
359  func_in = 1. - dec_fact_x * (dentdx/ent_c) * mp.x
360  - dec_fact_y * (dentdy/ent_c) * mp.y ;
361 
362  func_in.annule(nzet, nz-1) ;
363  func_in.std_spectral_base() ;
364 
365  Scalar func_ex(mp) ;
366  func_ex = 1. ;
367  func_ex.annule(0, nzet-1) ;
368  func_ex.std_spectral_base() ;
369 
370  // New enthalpy field
371  // ------------------
372  ent = ent * (func_in + func_ex) ;
373 
374  (ent.set_spectral_va()).smooth(nzet, ent.set_spectral_va()) ;
375 
376  double dentdx_new = ent.dsdx().val_grid_point(0,0,0,0) ;
377  double dentdy_new = ent.dsdy().val_grid_point(0,0,0,0) ;
378  cout << "dH/dx|_new = " << dentdx_new << endl ;
379  cout << "dH/dy|_new = " << dentdy_new << endl ;
380 
381  //-----------------------------------------------------
382  // Adaptation of the mapping to the new enthalpy field
383  //-----------------------------------------------------
384 
385  double dent_eq = ent.dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
386  double dent_pole = ent.dsdr().val_point(ray_pole(),0.,0.) ;
387  double rap_dent = fabs( dent_eq / dent_pole ) ;
388  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
389  diff_dHdr = rap_dent ;
390 
391  if ( rap_dent < thres_adapt ) {
392  adapt_flag = 0 ; // No adaptation of the mapping
393  cout << "******* FROZEN MAPPING *********" << endl ;
394  }
395  else {
396  adapt_flag = 1 ; // The adaptation of the mapping is to be
397  // performed
398  }
399 
400  ent_limit.set_etat_qcq() ;
401  for (int l=0; l<nzet; l++) { // loop on domains inside the star
402  ent_limit.set(l) = ent.val_grid_point(l, k_b, j_b, i_b) ;
403  }
404  ent_limit.set(nzet-1) = ent_b ;
405 
406  Map_et mp_prev = mp_et ;
407 
408  Cmp ent_cmp(ent) ;
409  mp.adapt(ent_cmp, par_adapt) ;
410  ent = ent_cmp ;
411 
412  // Readjustment of the external boundary of domain l=nzet
413  // to keep a fixed ratio with respect to star's surface
414 
415  double rr_in_1 = mp.val_r(1, -1., M_PI/2., 0.) ;
416 
417  // Resize of the outer boundary of the shell including the BH
418  double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
419  mp.resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
420 
421  // Resize of the inner boundary of the shell including the BH
422  double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
423  mp.resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
424 
425  // if (mer % 2 == 0) {
426 
427  if (nz > 4) {
428 
429  // Resize of the domain 1
430  double rr_out_1 = mp.val_r(1, 1., M_PI/2., 0.) ;
431  mp.resize(1, rr_in_1/rr_out_1 * resize_ns) ;
432 
433  if (nz > 5) {
434 
435  // Resize of the domain from 2 to N-4
436  double rr_out_nm3_new = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
437 
438  for (int i=1; i<nz-4; i++) {
439 
440  double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
441 
442  double rr_mid = rr_out_i
443  + (rr_out_nm3_new - rr_out_i) / double(nz - 3 - i) ;
444 
445  double rr_2timesi = 2. * rr_out_i ;
446 
447  if (rr_2timesi < rr_mid) {
448 
449  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
450  mp.resize(i+1, rr_2timesi / rr_out_ip1) ;
451 
452  }
453  else {
454 
455  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
456  mp.resize(i+1, rr_mid / rr_out_ip1) ;
457 
458  } // End of else
459 
460  } // End of i loop
461 
462  } // End of (nz > 5) loop
463 
464  } // End of (nz > 4) loop
465 
466  // } // End of (mer % 2) loop
467 
468  //----------------------------------------------------
469  // Computation of the enthalpy at the new grid points
470  //----------------------------------------------------
471 
472  mp_prev.homothetie(alpha_r) ;
473 
474  Cmp ent_cmp2 (ent) ;
475  mp.reevaluate(&mp_prev, nzet+1, ent_cmp2) ;
476  ent = ent_cmp2 ;
477 
478  double ent_s_max = -1 ;
479  int k_s_max = -1 ;
480  int j_s_max = -1 ;
481 
482  for (int k=0; k<mg->get_np(l_b); k++) {
483  for (int j=0; j<mg->get_nt(l_b); j++) {
484  double xx = fabs( ent.val_grid_point(l_b, k, j, i_b) ) ;
485  if (xx > ent_s_max) {
486  ent_s_max = xx ;
487  k_s_max = k ;
488  j_s_max = j ;
489  }
490  }
491  }
492  cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
493  << " and nzet : " << endl ;
494  cout << " " << ent_s_max << " reached for k = " << k_s_max
495  << " and j = " << j_s_max << endl ;
496 
497  //-------------------
498  // Equation of state
499  //-------------------
500 
501  equation_of_state() ; // computes new values for nbar (n), ener (e)
502  // and press (p) from the new ent (H)
503 
504  //----------------------------------------------------------
505  // Matter source terms in the gravitational field equations
506  //----------------------------------------------------------
507 
508  hydro_euler_bhns(kerrschild, mass_bh, sepa) ;
509  // computes new values for ener_euler (E),
510  // s_euler (S) and u_euler (U^i)
511 
512 
513  //-------------------------------------------------
514  // Computation of the minimum of the indicator chi
515  //-------------------------------------------------
516  double azimu_min = phi_min() ;
517  double rad_chi_min = radius_p(azimu_min) ;
518  double chi_min = chi_rp(rad_chi_min, azimu_min) ;
519 
520  cout << "| dH/dr_eq / dH/dr_pole | (minimum) = " << chi_min << endl ;
521  cout << " phi = " << azimu_min / M_PI << " [M_PI]" << endl ;
522  cout << " radius = " << rad_chi_min / km << " [km]" << endl ;
523 
524  diff_dHdr_min = chi_min ;
525  diff_phi_min = azimu_min ;
526  diff_radius = rad_chi_min ;
527 
528  //-----------------------------------
529  // Poisson equation for lapconf_auto
530  //-----------------------------------
531 
532  // Source
533  //--------
534 
535  Scalar sou_lap1(mp) ; // dzpuis = 0
536  sou_lap1 = qpig * lapconf_tot_tmp * pow(confo_tot,4.)
537  * (0.5*ener_euler + s_euler) ;
538 
539  sou_lap1.std_spectral_base() ;
540  sou_lap1.annule(nzet,nz-1) ;
541  sou_lap1.inc_dzpuis(4) ; // dzpuis : 0 -> 4
542 
543  Scalar sou_lap2(mp) ; // dzpuis = 4
544  sou_lap2 = 0.875 * (lapconf_auto+0.5) * taij_quad_auto
545  / pow(confo_auto+0.5,8.) ;
546  sou_lap2.std_spectral_base() ;
547 
548  source_lapconf = sou_lap1 + sou_lap2 ;
549 
550  source_lapconf.std_spectral_base() ;
551  // source_lapse.annule(nzet,nz-1) ;
552 
553  if (filter_r != 0) {
554  if (source_lapconf.get_etat() != ETATZERO) {
555  source_lapconf.filtre(filter_r) ;
556  // source_lapse.filtre_r(filter_r,0) ;
557  }
558  }
559 
560  assert(source_lapconf.get_dzpuis() == 4) ;
561 
562  // Resolution of the Poisson equation (Outer BC : lapconf_m1 -> 0)
563  // ----------------------------------
564 
565  lapconf_m1.set_etat_qcq() ;
566  lapconf_m1 = lapconf_auto - 0.5 ;
567  source_lapconf.poisson(par_lapconf, lapconf_m1) ;
568  ssjm1_lapconf = ssjm1lapconf ;
569 
570  // Check: has the Poisson equation been correctly solved ?
571  // -------------------------------------------------------
572 
573  Tbl tdiff_lapconf = diffrel(lapconf_m1.laplacian(), source_lapconf) ;
574  cout <<
575  "Relative error in the resolution of the equation for lapconf_auto : "
576  << endl ;
577  for (int l=0; l<nz; l++) {
578  cout << tdiff_lapconf(l) << " " ;
579  }
580  cout << endl ;
581  diff_lapconf = max(abs(tdiff_lapconf)) ;
582 
583  // Re-construction of the lapconf function
584  // ---------------------------------------
585  lapconf_auto = lapconf_m1 + 0.5 ;
586  // lapconf_tot = lapconf_auto + lapconf_comp
587  // lapconf_auto, _comp -> 0.5 (r -> inf)
588  // lapconf_tot -> 1 (r -> inf)
589 
590  //---------------------------------
591  // Poisson equation for confo_auto
592  //---------------------------------
593 
594  // Source
595  //--------
596 
597  Scalar sou_con1(mp) ; // dzpuis = 0
598  sou_con1 = - 0.5 * qpig * pow(confo_tot,5.) * ener_euler ;
599  sou_con1.std_spectral_base() ;
600  sou_con1.annule(nzet,nz-1) ;
601  sou_con1.inc_dzpuis(4) ; // dzpuis : 0 -> 4
602 
603  Scalar sou_con2(mp) ; // dzpuis = 4
604  sou_con2 = - 0.125 * taij_quad_auto / pow(confo_auto+0.5,7.) ;
605  sou_con2.std_spectral_base() ;
606 
607  source_confo = sou_con1 + sou_con2 ;
608 
609  source_confo.std_spectral_base() ;
610  // source_confo.annule(nzet,nz-1) ;
611 
612  if (filter_r != 0) {
613  if (source_confo.get_etat() != ETATZERO) {
614  source_confo.filtre(filter_r) ;
615  // source_confo.filtre_r(filter_r,0) ;
616  }
617  }
618 
619  assert(source_confo.get_dzpuis() == 4) ;
620 
621  // Resolution of the Poisson equation (Outer BC : confo_m1 -> 0)
622  // ----------------------------------
623 
624  confo_m1.set_etat_qcq() ;
625  confo_m1 = confo_auto - 0.5 ;
626  source_confo.poisson(par_confo, confo_m1) ;
627  ssjm1_confo = ssjm1confo ;
628 
629  // Check: has the Poisson equation been correctly solved ?
630  // -------------------------------------------------------
631 
632  Tbl tdiff_confo = diffrel(confo_m1.laplacian(), source_confo) ;
633  cout <<
634  "Relative error in the resolution of the equation for confo_auto : "
635  << endl ;
636  for (int l=0; l<nz; l++) {
637  cout << tdiff_confo(l) << " " ;
638  }
639  cout << endl ;
640  diff_confo = max(abs(tdiff_confo)) ;
641 
642  // Re-construction of the conformal factor
643  // ---------------------------------------
644  confo_auto = confo_m1 + 0.5 ; // confo_tot = confo_auto + confo_comp
645  // confo_auto, _comp -> 0.5 (r -> inf)
646  // confo_tot -> 1 (r -> inf)
647 
648  //----------------------------------------
649  // Vector Poisson equation for shift_auto
650  //----------------------------------------
651 
652  // Source
653  // ------
654 
655  Vector sou_shif1(mp, CON, mp.get_bvect_cart()) ; // dzpuis = 0
656  sou_shif1.set_etat_qcq() ;
657 
658  for (int i=1; i<=3; i++) {
659  sou_shif1.set(i) = 4.*qpig * lapconf_tot_tmp
660  * pow(confo_tot, 3.)
661  * (ener_euler + press) * u_euler(i) ;
662  }
663 
664  sou_shif1.std_spectral_base() ;
665  sou_shif1.annule(nzet, nz-1) ;
666 
667  for (int i=1; i<=3; i++) {
668  (sou_shif1.set(i)).inc_dzpuis(4) ; // dzpuis: 0 -> 4
669  }
670 
671  Vector sou_shif2(mp, CON, mp.get_bvect_cart()) ; // dzpuis = 4
672  sou_shif2.set_etat_qcq() ;
673  for (int i=1; i<=3; i++) {
674  sou_shif2.set(i) = 2. *
675  (taij_auto(i,1)*(d_lapconf_auto(1)
676  -7.*(lapconf_auto+0.5)*d_confo_auto(1)
677  /(confo_auto+0.5))
678  +taij_auto(i,2)*(d_lapconf_auto(2)
679  -7.*(lapconf_auto+0.5)*d_confo_auto(2)
680  /(confo_auto+0.5))
681  +taij_auto(i,3)*(d_lapconf_auto(3)
682  -7.*(lapconf_auto+0.5)*d_confo_auto(3)
683  /(confo_auto+0.5))
684  ) / pow(confo_auto+0.5,7.) ;
685  }
686  sou_shif2.std_spectral_base() ;
687 
688  source_shift = sou_shif1 + sou_shif2 ;
689 
690  source_shift.std_spectral_base() ;
691  // source_shift.annule(nzet, nz-1) ;
692 
693  // Resolution of the Poisson equation
694  // ----------------------------------
695 
696  // Filter for the source of shift vector
697 
698  if (filter_r_s != 0) {
699  for (int i=1; i<=3; i++) {
700  if (source_shift(i).get_etat() != ETATZERO) {
701  source_shift.set(i).filtre(filter_r_s) ;
702  // source_shift.set(i).filtre_r(filter_r_s, 0) ;
703  }
704  }
705  }
706 
707  if (filter_p_s != 0) {
708  for (int i=1; i<=3; i++) {
709  if (source_shift(i).get_etat() != ETATZERO) {
710  (source_shift.set(i)).filtre_phi(filter_p_s, nz-1) ;
711  // (source_shift.set(i)).filtre_phi(filter_p_s, 0) ;
712  }
713  }
714  }
715 
716  for (int i=1; i<=3; i++) {
717  if(source_shift(i).dz_nonzero()) {
718  assert( source_shift(i).get_dzpuis() == 4 ) ;
719  }
720  else {
721  (source_shift.set(i)).set_dzpuis(4) ;
722  }
723  }
724 
725  double lambda = double(1) / double(3) ;
726 
727  Tenseur source_p(mp, 1, CON, mp.get_bvect_cart() ) ;
728  source_p.set_etat_qcq() ;
729  for (int i=0; i<3; i++) {
730  source_p.set(i) = Cmp(source_shift(i+1)) ;
731  }
732 
733  Tenseur vect_auxi(mp, 1, CON, mp.get_bvect_cart()) ;
734  vect_auxi.set_etat_qcq() ;
735  for (int i=0; i<3 ;i++) {
736  vect_auxi.set(i) = 0. ;
737  }
738  Tenseur scal_auxi(mp) ;
739  scal_auxi.set_etat_qcq() ;
740  scal_auxi.set().annule_hard() ;
741  scal_auxi.set_std_base() ;
742 
743  Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart() ) ;
744  resu_p.set_etat_qcq() ;
745 
746  source_p.poisson_vect(lambda, par_shift2, resu_p, vect_auxi,
747  scal_auxi) ;
748 
749  for (int i=1; i<=3; i++)
750  shift_auto.set(i) = resu_p(i-1) ;
751 
752  ssjm1_khi = ssjm1khi ;
753 
754  for (int i=0; i<3; i++) {
755  ssjm1_wshift.set(i+1) = ssjm1wshift(i) ;
756  }
757 
758  // Check: has the equation for shift_auto been correctly solved ?
759  // --------------------------------------------------------------
760 
762  + lambda * shift_auto.divergence(flat).derive_con(flat) ;
763 
764  source_shift.dec_dzpuis() ;
765 
766  Tbl tdiff_shift_x = diffrel(lap_shift(1), source_shift(1)) ;
767  Tbl tdiff_shift_y = diffrel(lap_shift(2), source_shift(2)) ;
768  Tbl tdiff_shift_z = diffrel(lap_shift(3), source_shift(3)) ;
769 
770  cout <<
771  "Relative error in the resolution of the equation for shift_auto : "
772  << endl ;
773  cout << "x component : " ;
774  for (int l=0; l<nz; l++) {
775  cout << tdiff_shift_x(l) << " " ;
776  }
777  cout << endl ;
778  cout << "y component : " ;
779  for (int l=0; l<nz; l++) {
780  cout << tdiff_shift_y(l) << " " ;
781  }
782  cout << endl ;
783  cout << "z component : " ;
784  for (int l=0; l<nz; l++) {
785  cout << tdiff_shift_z(l) << " " ;
786  }
787  cout << endl ;
788 
789  diff_shift_x = max(abs(tdiff_shift_x)) ;
790  diff_shift_y = max(abs(tdiff_shift_y)) ;
791  diff_shift_z = max(abs(tdiff_shift_z)) ;
792 
793 
794  //-----------------------------
795  // Relative change in enthalpy
796  //-----------------------------
797 
798  Tbl diff_ent_tbl = diffrel( ent, ent_jm1 ) ;
799  diff_ent = diff_ent_tbl(0) ;
800  for (int l=0; l<nzet; l++) {
801  diff_ent += diff_ent_tbl(l) ;
802  }
803  diff_ent /= nzet ;
804 
805  ent_jm1 = ent ;
806 
807  /*
808  des_profile( lapconf_auto, 0., 10.,
809  M_PI/2., M_PI, "Self lapconf function of NS",
810  "Lapconf (theta=pi/2, phi=0)" ) ;
811 
812  des_profile( lapconf_tot, 0., 10.,
813  M_PI/2., M_PI, "Total lapconf function seen by NS",
814  "Lapconf (theta=pi/2, phi=0)" ) ;
815 
816  des_profile( confo_auto, 0., 10.,
817  M_PI/2., M_PI, "Self conformal factor of NS",
818  "Confo (theta=pi/2, phi=0)" ) ;
819 
820  des_profile( confo_tot, 0., 10.,
821  M_PI/2., M_PI, "Total conformal factor seen by NS",
822  "Confo (theta=pi/2, phi=0)" ) ;
823 
824  des_coupe_vect_z( shift_auto, 0., -3., 0.5, 3,
825  "Self shift vector of NS") ;
826 
827  des_coupe_vect_z( shift_tot, 0., -3., 0.5, 3,
828  "Total shift vector seen by NS") ;
829  */
830  } // End of main loop
831 
832  //=========================================================
833  // End of iteration
834  //=========================================================
835 
836 
837 }
838 }
Lorene::Map_et
Radial mapping of rather general form.
Definition: map.h:2752
Lorene::Star_bhns::loggam
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: star_bhns.h:93
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::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Lorene::Scalar::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition: scalar_r_manip.C:470
Lorene::Scalar::val_point
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition: scalar.C:890
Lorene::Scalar::dsdy
const Scalar & dsdy() const
Returns of *this , where .
Definition: scalar_deriv.C:297
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::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
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::Star::press
Scalar press
Fluid pressure.
Definition: star.h:194
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Scalar::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition: scalar_manip.C:151
Lorene::Tensor::annule
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition: tensor.C:671
Lorene::Tenseur::poisson_vect
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Definition: tenseur_pde.C:118
Lorene::Star_bhns::lapconf_comp
Scalar lapconf_comp
Lapconf function generated by the companion black hole.
Definition: star_bhns.h:116
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:808
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::Map::val_r
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Lorene::Star_bhns::confo_tot
Scalar confo_tot
Total conformal factor.
Definition: star_bhns.h:163
Lorene::Star_bhns::ssjm1_khi
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: star_bhns.h:210
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Lorene::Star_bhns::pot_centri
Scalar pot_centri
Centrifugal potential.
Definition: star_bhns.h:110
Lorene::Star_bhns::d_lapconf_auto
Vector d_lapconf_auto
Derivative of the lapconf function generated by the star .
Definition: star_bhns.h:130
Lorene::Map::resize
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Star_bhns::radius_p
double radius_p(double phi)
Radius of the star to the direction of and .
Definition: star_bhns_chi.C:74
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Lorene::Star::u_euler
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Unites
Standard units of space, time and mass.
Lorene::Star::ray_eq_pi
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition: star_global.C:186
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Lorene::Star_bhns::gam0
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
Definition: star_bhns.h:107
Lorene::Star_bhns::lapconf_auto
Scalar lapconf_auto
Lapconf function generated by the star.
Definition: star_bhns.h:113
Lorene::Star::equation_of_state
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition: star.h:180
Lorene::Star_bhns::shift_auto
Vector shift_auto
Shift vector generated by the star.
Definition: star_bhns.h:138
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
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::Star_bhns::ssjm1_lapconf
Scalar ssjm1_lapconf
Effective source at the previous step for the resolution of the Poisson equation for lapconf_auto .
Definition: star_bhns.h:197
Lorene::Cmp::annule_hard
void annule_hard()
Sets the Cmp to zero in a hard way.
Definition: cmp.C:338
Lorene::Star::ener_euler
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition: star.h:190
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
Lorene::Star_bhns::equilibrium_bhns
void equilibrium_bhns(double ent_c, const double &mass_bh, const double &sepa, bool kerrschild, int mer, int mermax_ns, int mermax_potvit, int mermax_poisson, int filter_r, int filter_r_s, int filter_p_s, double relax_poisson, double relax_potvit, double thres_adapt, double resize_ns, const Tbl &fact_resize, Tbl &diff)
Computes an equilibrium configuration.
Definition: star_bhns_equilibrium.C:67
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Scalar::dsdx
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
Lorene::Star_bhns::irrotational
bool irrotational
true for an irrotational star, false for a corotating one
Definition: star_bhns.h:72
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::Star_bhns::gam
Scalar gam
Lorentz factor between the fluid and the co-orbiting observer.
Definition: star_bhns.h:102
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
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::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Star_bhns::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star ).
Definition: star_bhns.h:192
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
Lorene::Tensor::divergence
const Tensor & divergence(const Metric &gam) const
Computes the divergence of this with respect to some metric .
Definition: tensor.C:1055
Lorene::Map::x
Coord x
x coordinate centered on the grid
Definition: map.h:726
Lorene::Star_bhns::ssjm1_wshift
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition: star_bhns.h:220
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
Lorene::Star_bhns::ssjm1_confo
Scalar ssjm1_confo
Effective source at the previous step for the resolution of the Poisson equation for confo_auto .
Definition: star_bhns.h:202
Lorene::Scalar::derive_con
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Lorene::Star_bhns::phi_min
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Definition: star_bhns_chi.C:83
Lorene::Star_bhns::hydro_euler_bhns
void hydro_euler_bhns(bool kerrschild, const double &mass_bh, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition: star_bhns_hydro.C:59
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::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
Lorene::Star_bhns::velo_pot_bhns
double velo_pot_bhns(const double &mass_bh, const double &sepa, bool kerrschild, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Definition: star_bhns_vel_pot.C:67
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::Star_bhns::taij_quad_auto
Scalar taij_quad_auto
Part of the scalar generated by .
Definition: star_bhns.h:187
Lorene::Star_bhns::confo_auto
Scalar confo_auto
Conformal factor generated by the star.
Definition: star_bhns.h:157
Lorene::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
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::Star::s_euler
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:637
Lorene::Star_bhns::chi_rp
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
Definition: star_bhns_chi.C:61
Lorene::Cmp::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
Lorene::Star::ray_pole
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
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::Valeur::smooth
void smooth(int nzet, Valeur &uuva) const
Changes the function *this as a smooth one when there exists a discontinuity between the nucleus and ...
Definition: valeur_smooth.C:77
Lorene::Star_bhns::d_confo_auto
Vector d_confo_auto
Derivative of the conformal factor generated by the star .
Definition: star_bhns.h:168
Lorene::Tenseur::set_std_base
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Lorene::Star_bhns::taij_auto
Sym_tensor taij_auto
Part of the extrinsic curvature tensor generated by shift_auto , lapse_auto , and confo_auto .
Definition: star_bhns.h:182
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
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316