LORENE
star_bin_equilibrium_xcts.C
1 /*
2  * Method of class Star_bin to compute an equilibrium configuration
3  * (see file star.h for documentation).
4  */
5 
6 /*
7  * Copyright (c) 2010 Michal Bejger
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License version 2
13  * as published by the Free Software Foundation.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 char star_bin_equilibrium_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $" ;
27 
28 /*
29  * $Id: star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $
30  * $Log: star_bin_equilibrium_xcts.C,v $
31  * Revision 1.13 2014/10/13 08:53:38 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.12 2014/10/06 15:13:16 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.11 2011/03/25 16:28:12 e_gourgoulhon
38  * Still in progress
39  *
40  * Revision 1.10 2010/12/20 15:42:10 m_bejger
41  * Various rearrangements of fields in Poissson equations
42  *
43  * Revision 1.9 2010/12/14 17:34:42 m_bejger
44  * Improved iteration for beta_auto poisson_vect()
45  *
46  * Revision 1.8 2010/12/09 10:48:06 m_bejger
47  * Testing the main equations
48  *
49  * Revision 1.7 2010/10/26 18:46:28 m_bejger
50  * Added table fact_resize for domain resizing
51  *
52  * Revision 1.6 2010/10/18 19:08:14 m_bejger
53  * Changed to allow for calculations with more than one domain in the star
54  *
55  * Revision 1.5 2010/06/23 20:40:56 m_bejger
56  * Corrections in equations for Psi_auto, chi_auto and beta_auto
57  *
58  * Revision 1.4 2010/06/18 13:28:59 m_bejger
59  * Adjusted the computation of the first integral, radial scale
60  *
61  * Revision 1.3 2010/06/17 17:05:06 m_bejger
62  * Testing version
63  *
64  * Revision 1.2 2010/06/15 08:21:21 m_bejger
65  * Minor changes; still not working properly
66  *
67  * Revision 1.1 2010/05/04 07:51:05 m_bejger
68  * Initial version
69  *
70  * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $
71  *
72  */
73 
74 // C headers
75 #include <cmath>
76 
77 // Lorene headers
78 #include "cmp.h"
79 #include "tenseur.h"
80 #include "metrique.h"
81 #include "star.h"
82 #include "param.h"
83 #include "graphique.h"
84 #include "utilitaires.h"
85 #include "tensor.h"
86 #include "nbr_spx.h"
87 #include "unites.h"
88 
89 
90 namespace Lorene {
91 void Star_bin_xcts::equilibrium(double ent_c,
92  int mermax,
93  int mermax_potvit,
94  int mermax_poisson,
95  double relax_poisson,
96  double relax_potvit,
97  double thres_adapt,
98  const Tbl& fact_resize,
99  const Tbl* pent_limit,
100  Tbl& diff) {
101 
102  // Fundamental constants and units
103  // -------------------------------
104 
105  using namespace Unites ;
106 
107  // Initializations
108  // ---------------
109 
110  const Mg3d* mg = mp.get_mg() ;
111  int nz = mg->get_nzone() ; // total number of domains
112 
113  // The following is required to initialize mp_prev as a Map_et:
114  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
115 
116  // Domain and radial indices of points at the surface of the star:
117  int l_b = nzet - 1 ;
118  int i_b = mg->get_nr(l_b) - 1 ;
119  int k_b ;
120  int j_b ;
121 
122  // Value of the enthalpy defining the surface of the star
123  double ent_b = 0 ;
124 
125  // Error indicators
126  // ----------------
127 
128  double& diff_ent = diff.set(0) ;
129  double& diff_vel_pot = diff.set(1) ;
130  double& diff_psi = diff.set(2) ;
131  double& diff_chi = diff.set(3) ;
132  double& diff_beta_x = diff.set(4) ;
133  double& diff_beta_y = diff.set(5) ;
134  double& diff_beta_z = diff.set(6) ;
135 
136  // Parameters for the function Map_et::adapt
137  // -----------------------------------------
138 
139  Param par_adapt ;
140  int nitermax = 100 ;
141  int niter ;
142  int adapt_flag = 1 ; // 1 = performs the full computation,
143  // 0 = performs only the rescaling by
144  // the factor alpha_r
145  //## int nz_search = nzet + 1 ; // Number of domains for searching the
146  // enthalpy
147 
148  int nz_search = nzet ; // Number of domains
149  // for searching
150  // the enthalpy isosurfaces
151 
152  double precis_secant = 1.e-14 ;
153  double alpha_r ;
154  double reg_map = 1. ; // 1 = regular mapping,
155  // 0 = contracting mapping
156 
157  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
158  // locate zeros by the secant method
159 
160  par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
161  // to the isosurfaces of ent is to be performed
162 
163  par_adapt.add_int(nz_search, 2) ; // number of domains to search
164  // the enthalpy isosurface
165 
166  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
167  // 0 = performs only the rescaling by
168  // the factor alpha_r
169 
170  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
171  // (theta_*, phi_*)
172 
173  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
174  // (theta_*, phi_*)
175 
176  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
177  // the secant method
178 
179  par_adapt.add_double(precis_secant, 0) ;// required absolute precision in
180  // the determination of zeros by
181  // the secant method
182  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
183 
184  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
185  // distances will be multiplied
186 
187  // Enthalpy values for the adaptation
188  Tbl ent_limit(nzet) ;
189  if (pent_limit != 0x0) ent_limit = *pent_limit ;
190 
191  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
192  // to define the isosurfaces.
193 
194  double precis_poisson = 1.e-16 ;
195 
196  Cmp ssjm1psi (ssjm1_psi) ;
197  Cmp ssjm1chi (ssjm1_chi) ;
198 
199  // Parameters for the function Scalar::poisson for Psi_auto
200  // ---------------------------------------------------------------
201 
202  Param par_psi ;
203 
204  par_psi.add_int(mermax_poisson, 0) ; // maximum number of iterations
205  par_psi.add_double(relax_poisson, 0) ; // relaxation parameter
206  par_psi.add_double(precis_poisson, 1) ; // required precision
207  par_psi.add_int_mod(niter, 0) ; // number of iterations actually used
208  par_psi.add_cmp_mod( ssjm1psi ) ;
209 
210  // Parameters for the function Scalar::poisson for chi_auto
211  // ---------------------------------------------------------------
212 
213  Param par_chi ;
214 
215  par_chi.add_int(mermax_poisson, 0) ; // maximum number of iterations
216  par_chi.add_double(relax_poisson, 0) ; // relaxation parameter
217  par_chi.add_double(precis_poisson, 1) ; // required precision
218  par_chi.add_int_mod(niter, 0) ; // number of iterations actually used
219  par_chi.add_cmp_mod( ssjm1chi ) ;
220 
221  // Parameters for the function Vector::poisson for beta
222  // ----------------------------------------------------
223 
224  Param par_beta ;
225 
226  par_beta.add_int(mermax_poisson, 0) ; // maximum number of
227  // iterations
228  par_beta.add_double(relax_poisson, 0) ; // relaxation parameter
229  par_beta.add_double(precis_poisson, 1) ; // required precision
230  par_beta.add_int_mod(niter, 0) ; // number of iterations
231  // actually used
232 
233  // Sources at the previous step, for a poisson_vect() solver
234  Cmp ssjm1khi (ssjm1_khi) ;
235 
236  Tenseur ssjm1wbeta(mp, 1, CON, mp.get_bvect_cart()) ;
237  ssjm1wbeta.set_etat_qcq() ;
238  for (int i=0; i<3; i++) ssjm1wbeta.set(i) = Cmp(ssjm1_wbeta(i+1)) ;
239 
240  par_beta.add_cmp_mod(ssjm1khi) ;
241  par_beta.add_tenseur_mod(ssjm1wbeta) ;
242 
243  // Redefinition of external potential
244  // See Eq (99) from Gourgoulhon et al. (2001)
245  // logN = logN_auto + logn_ac_rest = log(chi_auto + 1.)
246  // - log(Psi_auto + 1.) + logn_ac_rest
247  //------------------------------------
248 
249  Scalar Psi_auto_p1 = Psi_auto + 1. ;
250  Scalar chi_auto_p1 = chi_auto + 1. ;
251 
252  Scalar logn_auto = log(chi_auto_p1) - log(Psi_auto_p1) ;
253  logn_auto.std_spectral_base() ;
254 
255  Scalar logn_ac_rest = log(1. + chi_comp/chi_auto_p1)
256  - log(1. + Psi_comp/Psi_auto_p1) ;
257 
258  logn_ac_rest.std_spectral_base() ;
259 
260  //cout << "logn_auto" << norme(logn_auto) << endl ;
261  //cout << "logn_ac_rest" << norme(logn_ac_rest) << endl ;
262  //cout << "pot_centri" << norme(pot_centri) << endl ;
263  //cout << "loggam" << norme(loggam) << endl ;
264 
265  Scalar pot_ext = logn_ac_rest + pot_centri + loggam ;
266  //cout << "pot_ext" << norme(pot_ext) << endl ;
267 
268  Scalar ent_jm1 = ent ; // Enthalpy at previous step
269 
270  Scalar source_tot(mp) ; // source term in equations Psi_auto
271  // and chi_auto
272 
273  Vector source_beta(mp, CON, mp.get_bvect_cart()) ; // source term
274  // in the equation
275  // for beta_auto
276 
277  //==================================================================
278  // Start of iteration
279  //==================================================================
280 
281  for(int mer=0 ; mer<mermax ; mer++ ) {
282 
283  cout << "-----------------------------------------------" << endl ;
284  cout << "step: " << mer << endl ;
285  cout << "diff_ent = " << diff_ent << endl ;
286 
287  //-----------------------------------------------------
288  // Resolution of the elliptic equation for the velocity
289  // scalar potential
290  //-----------------------------------------------------
291 
292  if (irrotational) {
293 
294  diff_vel_pot = velocity_potential(mermax_potvit, precis_poisson,
295  relax_potvit) ;
296  }
297 
298  diff_vel_pot = 0. ; // to avoid the warning
299 
300  //-----------------------------------------------------
301  // Computation of the new radial scale
302  //--------------------------------------------------
303 
304  // alpha_r (r = alpha_r r') is determined so that the enthalpy
305  // takes the requested value ent_b at the stellar surface
306 
307  // Values at the center of the star:
308  double logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
309  double pot_ext_c = pot_ext.val_grid_point(0, 0, 0, 0) ;
310 
311  // Search for the reference point (theta_*, phi_*) [notation of
312  // Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
313  // at the surface of the star
314  // ------------------------------------------------------------
315  double alpha_r2 = 0 ;
316  for (int k=0; k<mg->get_np(l_b); k++) {
317  for (int j=0; j<mg->get_nt(l_b); j++) {
318 
319  double pot_ext_b = pot_ext.val_grid_point(l_b, k, j, i_b) ;
320  double logn_auto_b = logn_auto.val_grid_point(l_b, k, j, i_b) ;
321  // See Eq (100) from Gourgoulhon et al. (2001)
322  double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
323  ( logn_auto_b - logn_auto_c ) ;
324 
325  // cout << "k, j, alpha_r2_jk : " << k << " " << j << " " << alpha_r2_jk << endl ;
326 
327  if (alpha_r2_jk > alpha_r2) {
328  alpha_r2 = alpha_r2_jk ;
329  k_b = k ;
330  j_b = j ;
331  }
332 
333  }
334  }
335 
336  alpha_r = sqrt(alpha_r2) ;
337 
338  cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
339  << alpha_r << endl ;
340 
341  // New value of logn_auto
342  // ----------------------
343 
344  Psi_auto = pow(Psi_auto +1.,alpha_r2) - 1. ;
345  chi_auto = pow(chi_auto +1.,alpha_r2) - 1. ;
348 
349  //------------------------------------------------------------
350  // Change the values of the inner points of the domain adjascent
351  // to the surface of the star by those of the outer points of
352  // the domain under the surface
353  //------------------------------------------------------------
354 
357 
358  logn_auto = log(chi_auto + 1.) - log(Psi_auto + 1.) ;
359  logn_auto.std_spectral_base() ;
360 
361  logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
362 
363  //------------------------------------------
364  // First integral --> enthalpy in all space
365  // See Eq (98) from Gourgoulhon et al. (2001)
366  //-------------------------------------------
367 
368  ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
369  cout.precision(8) ;
370  cout << "pot" << norme(pot_ext) << endl ;
371 
372  (ent.set_spectral_va()).smooth(nzet, ent.set_spectral_va()) ;
373 
374  //----------------------------------------------------
375  // Adaptation of the mapping to the new enthalpy field
376  //----------------------------------------------------
377 
378  // Shall the adaptation be performed (cusp) ?
379  // ------------------------------------------
380 
381  double dent_eq = ent.dsdr().val_point(ray_eq(),M_PI/2.,0.) ;
382  double dent_pole = ent.dsdr().val_point(ray_pole(),0.,0.) ;
383  double rap_dent = fabs( dent_eq / dent_pole ) ;
384  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
385 
386  if ( rap_dent < thres_adapt ) {
387  adapt_flag = 0 ; // No adaptation of the mapping
388  cout << "******* FROZEN MAPPING *********" << endl ;
389  }
390  else{
391  adapt_flag = 1 ; // The adaptation of the mapping is to be
392  // performed
393  }
394 
395  ent_limit.set_etat_qcq() ;
396  for (int l=0; l<nzet; l++) { // loop on domains inside the star
397  ent_limit.set(l) = ent.val_grid_point(l, k_b, j_b, i_b) ;
398 
399  } ent_limit.set(nzet-1) = ent_b ;
400 
401  Map_et mp_prev = mp_et ;
402 
403  Cmp ent_cmp(ent) ;
404  mp.adapt(ent_cmp, par_adapt) ;
405  ent = ent_cmp ;
406 
407  // Readjustment of the external boundary of domain l=nzet
408  // to keep a fixed ratio with respect to star's surface
409 
410  double rr_in_1 = mp.val_r(nzet, -1., M_PI/2., 0.) ;
411 
412  // Resizes the outer boundary of the shell including the comp. NS
413  double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
414 
415  mp.resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
416 
417  // Resizes the inner boundary of the shell including the comp. NS
418  double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
419 
420  mp.resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
421 
422  if (nz > nzet+3) {
423 
424  // Resize of the domain from 1(nzet) to N-4
425  double rr_out_nm3_new = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
426 
427  for (int i=nzet-1; i<nz-4; i++) {
428 
429  double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
430 
431  double rr_mid = rr_out_i
432  + (rr_out_nm3_new - rr_out_i) / double(nz - 3 - i) ;
433 
434  double rr_2timesi = 2. * rr_out_i ;
435 
436  if (rr_2timesi < rr_mid) {
437 
438  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
439 
440  mp.resize(i+1, rr_2timesi / rr_out_ip1) ;
441 
442  } else {
443 
444  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
445 
446  mp.resize(i+1, rr_mid / rr_out_ip1) ;
447 
448  } // End of else
449 
450  } // End of i loop
451 
452  } // End of (nz > nzet+3) loop
453 
454 // if (nz>= 5) {
455 //
456 // double separation = 2. * fabs(mp.get_ori_x()) ;
457 // double ray_eqq = ray_eq() ;
458 // double ray_eqq_pi = ray_eq_pi() ;
459 // double new_rr_out_2 = (separation - ray_eqq) * 0.95 ;
460 // double new_rr_out_3 = (separation + ray_eqq_pi) * 1.05 ;
461 //
462 // double rr_in_1 = mp.val_r(1,-1., M_PI/2, 0.) ;
463 // double rr_out_1 = mp.val_r(1, 1., M_PI/2, 0.) ;
464 // double rr_out_2 = mp.val_r(2, 1., M_PI/2, 0.) ;
465 // double rr_out_3 = mp.val_r(3, 1., M_PI/2, 0.) ;
466 //
467 // mp.resize(1, 0.5*(new_rr_out_2 + rr_in_1) / rr_out_1) ;
468 // mp.resize(2, new_rr_out_2 / rr_out_2) ;
469 // mp.resize(3, new_rr_out_3 / rr_out_3) ;
470 //
471 // for (int dd=4; dd<=nz-2; dd++) {
472 // mp.resize(dd, new_rr_out_3 * pow(4., dd-3) /
473 // mp.val_r(dd, 1., M_PI/2, 0.)) ;
474 // }
475 //
476 // }
477 
478  //else cout << "too small number of domains" << endl ;
479 
480 
481  //------------------------------------------------------------------
482  // Computation of the enthalpy at the new grid points
483  //------------------------------------------------------------------
484 
485  mp_prev.homothetie(alpha_r) ;
486 
487  //Cmp ent_cmp2 (ent) ;
488  //mp.reevaluate_symy(&mp_prev, nzet+1, ent_cmp2) ;
489  //ent = ent_cmp2 ;
490 
491  mp.reevaluate_symy(&mp_prev, nzet+1, ent) ;
492 
493  double ent_s_max = -1 ;
494  int k_s_max = -1 ;
495  int j_s_max = -1 ;
496  for (int k=0; k<mg->get_np(l_b); k++) {
497  for (int j=0; j<mg->get_nt(l_b); j++) {
498  double xx = fabs( ent.val_grid_point(l_b, k, j, i_b) ) ;
499  if (xx > ent_s_max) {
500  ent_s_max = xx ;
501  k_s_max = k ;
502  j_s_max = j ;
503  }
504  }
505  }
506  cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
507  << " and nzet : " << endl ;
508  cout << " " << ent_s_max << " reached for k = " << k_s_max <<
509  " and j = " << j_s_max << endl ;
510 
511  //------------------------------------------------------------------
512  // Equation of state
513  //------------------------------------------------------------------
514 
515  equation_of_state() ; // computes new values for nbar (n), ener (e)
516  // and press (p) from the new ent (H)
517 
518  //------------------------------------------------------------------
519  // Matter source terms in the gravitational field equations
520  //------------------------------------------------------------------
521 
522  hydro_euler() ; // computes new values for ener_euler (E),
523  // s_euler (S) and u_euler (U^i)
524 
525  // -------------------------------
526  // AUXILIARY QUANTITIES
527  // -------------------------------
528 
529  int nr = mp.get_mg()->get_nr(0) ;
530  int nt = mp.get_mg()->get_nt(0) ;
531  int np = mp.get_mg()->get_np(0) ;
532 
533  Scalar Psi3 = psi4 / Psi ;
534  Psi3.std_spectral_base() ;
535 
536  Scalar sPsi7 = Psi * pow(psi4, -2.) ;
537  sPsi7.std_spectral_base() ;
538 
539  //------------------------------------------------------------------
540  // Poisson equation for Psi_auto (Eq. 8.127 of arXiv:gr-qc/0703035)
541  //------------------------------------------------------------------
542 
543  // Source
544  //--------
545 
546  source_tot = - 0.5 * qpig * psi4 % Psi % ener_euler
547  - 0.125 * sPsi7 * (hacar_auto + hacar_comp) ;
548  source_tot.std_spectral_base() ;
549 
550  // Resolution of the Poisson equation
551  // ----------------------------------
552 
553  cout << "Resolution of the Poisson equation for Psi_auto:" << endl ;
554  source_tot.poisson(par_psi, Psi_auto) ;
555  ssjm1_psi = ssjm1psi ;
556 
557  cout << "Psi_auto: " << endl << norme(Psi_auto/(nr*nt*np)) << endl ;
558 
559  // Check: has the Poisson equation been correctly solved ?
560  // -----------------------------------------------------
561 
562  Tbl tdiff_psi = diffrel(Psi_auto.laplacian(), source_tot) ;
563  cout <<
564  "Relative error in the resolution of the equation for Psi_auto : "
565  << endl ;
566  for (int l=0; l<nz; l++) {
567  cout << tdiff_psi(l) << " " ;
568  }
569  cout << endl
570  << "==========================================================="
571  << endl ;
572 
573  diff_psi = max(abs(tdiff_psi)) ;
574 
575  //------------------------------------------------------------------
576  // Poisson equation for chi_auto (Eq. 8.129 of arXiv:gr-qc/0703035)
577  //------------------------------------------------------------------
578 
579  // Source
580  //--------
581 
582  source_tot = chi * 0.5 * qpig * psi4 % (ener_euler + 2.*s_euler)
583  + 0.875 * nn % sPsi7 * (hacar_auto + hacar_comp) ;
584  source_tot.std_spectral_base() ;
585 
586  // Resolution of the Poisson equation
587  // ----------------------------------
588 
589  cout << "Resolution of the Poisson equation for chi_auto:" << endl ;
590  source_tot.poisson(par_chi, chi_auto) ;
591  ssjm1_chi = ssjm1chi ;
592 
593  cout << "chi_auto: " << endl << norme(chi_auto/(nr*nt*np)) << endl ;
594 
595  // Check: has the Poisson equation been correctly solved ?
596  // -----------------------------------------------------
597 
598  Tbl tdiff_chi = diffrel(chi_auto.laplacian(), source_tot) ;
599  cout <<
600  "Relative error in the resolution of the equation for chi_auto : "
601  << endl ;
602 
603  for (int l=0; l<nz; l++) {
604  cout << tdiff_chi(l) << " " ;
605  }
606  cout << endl
607  << "==========================================================="
608  << endl ;
609 
610  diff_chi = max(abs(tdiff_chi)) ;
611 
612  //------------------------------------------------------------------
613  // Vector Poisson equation for beta_auto
614  // (Eq. 8.128 of arXiv:gr-qc/0703035)
615  //------------------------------------------------------------------
616 
617  // Source
618  //--------
619 
620  source_beta = 4.* qpig * chi % Psi3 * (ener_euler + press) * u_euler
621  + 2. * sPsi7 * contract(haij_auto, 1, dcov_chi, 0)
622  -14. * nn % sPsi7 * contract(haij_auto, 1, dcov_Psi, 0) ;
623  source_beta.std_spectral_base() ;
624 
625  // Resolution of the Poisson equation
626  // ----------------------------------
627 
628  Tenseur source_p(mp, 1, CON, mp.get_bvect_cart() ) ;
629  source_p.set_etat_qcq() ;
630  for (int i=0; i<3; i++) {
631 
632  source_p.set(i) = Cmp(source_beta(i+1)) ;
633  //source_p.set(i).filtre(4) ;
634 
635  }
636 
637  Tenseur vect_auxi (mp, 1, CON, mp.get_bvect_cart()) ;
638  vect_auxi.set_etat_qcq() ;
639  for (int i=0; i<3 ;i++) vect_auxi.set(i) = Cmp(w_beta(i+1)) ;
640  vect_auxi.set_std_base() ;
641 
642 
643  Tenseur scal_auxi (mp) ;
644  scal_auxi.set_etat_qcq() ;
645  scal_auxi.set() = Cmp(khi) ;
646  scal_auxi.set_std_base() ;
647 
648  Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart()) ;
649  resu_p.set_etat_qcq() ;
650  for (int i=0; i<3 ;i++) resu_p.set(i) = Cmp(beta_auto.set(i+1));
651  resu_p.set_std_base() ;
652 
653  // Resolution of the vector Poisson equation
654  // -----------------------------------------
655 
656  double lambda = double(1) / double(3) ;
657  source_p.poisson_vect(lambda, par_beta, resu_p, vect_auxi, scal_auxi) ;
658  //source_p.poisson_vect_oohara(lambda, par_beta, resu_p, scal_auxi) ;
659 
660  for (int i=1; i<=3; i++) {
661  beta_auto.set(i) = resu_p(i-1) ;
662  w_beta.set(i) = vect_auxi(i-1) ;
663  }
664  khi = scal_auxi() ;
665 
666  // Values of sources from the previous step
667  ssjm1_khi = ssjm1khi ;
668  for (int i=0; i<3; i++) ssjm1_wbeta.set(i+1) = ssjm1wbeta(i) ;
669 
670  cout << "beta_auto_x" << endl << norme(beta_auto(1)/(nr*nt*np))
671  << endl ;
672  cout << "beta_auto_y" << endl << norme(beta_auto(2)/(nr*nt*np))
673  << endl ;
674  cout << "beta_auto_z" << endl << norme(beta_auto(3)/(nr*nt*np))
675  << endl ;
676 
677  // Check: has the equation for beta_auto been correctly solved ?
678  // --------------------------------------------------------------
679 
680  Vector lap_beta = (beta_auto.derive_con(flat)).divergence(flat)
681  + lambda* beta_auto.divergence(flat).derive_con(flat) ;
682 
683  source_beta.dec_dzpuis() ;
684  Tbl tdiff_beta_x = diffrel(lap_beta(1), source_beta(1)) ;
685  Tbl tdiff_beta_y = diffrel(lap_beta(2), source_beta(2)) ;
686  Tbl tdiff_beta_z = diffrel(lap_beta(3), source_beta(3)) ;
687 
688  cout <<
689  "Relative error in the resolution of the equation for beta_auto : "
690  << endl ;
691  cout << "x component : " ;
692  for (int l=0; l<nz; l++) {
693  cout << tdiff_beta_x(l) << " " ;
694  }
695  cout << endl ;
696  cout << "y component : " ;
697  for (int l=0; l<nz; l++) {
698  cout << tdiff_beta_y(l) << " " ;
699  }
700  cout << endl ;
701  cout << "z component : " ;
702  for (int l=0; l<nz; l++) {
703  cout << tdiff_beta_z(l) << " " ;
704  }
705  cout << endl ;
706 
707  diff_beta_x = max(abs(tdiff_beta_x)) ;
708  diff_beta_y = max(abs(tdiff_beta_y)) ;
709  diff_beta_z = max(abs(tdiff_beta_z)) ;
710 
711  // End of relativistic equations
712 
713  //-------------------------------------------------
714  // Relative change in enthalpy
715  //-------------------------------------------------
716 
717  Tbl diff_ent_tbl = diffrel( ent, ent_jm1 ) ;
718  diff_ent = diff_ent_tbl(0) ;
719  for (int l=1; l<nzet; l++) diff_ent += diff_ent_tbl(l) ;
720 
721  diff_ent /= nzet ;
722 
723 
724  ent_jm1 = ent ;
725 
726 
727  } // End of main loop
728 
729  //==================================================================
730  // End of iteration
731  //==================================================================
732 
733 }
734 
735 
736 
737 
738 
739 
740 
741 }
Scalar chi
Total function .
Definition: star.h:1155
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1142
Vector dcov_chi
Covariant derivative of the function .
Definition: star.h:1174
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Radial mapping of rather general form.
Definition: map.h:2752
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
bool irrotational
true for an irrotational star, false for a corotating one
Definition: star.h:1099
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Vector dcov_Psi
Covariant derivative of the conformal factor .
Definition: star.h:1171
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
Scalar ent
Log-enthalpy.
Definition: star.h:190
void equilibrium(double ent_c, int mermax, int mermax_potvit, int mermax_poisson, double relax_poisson, double relax_potvit, double thres_adapt, const Tbl &fact, const Tbl *pent_limit, Tbl &diff)
Computes an equilibrium configuration.
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
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
Tensor field of valence 1.
Definition: vector.h:188
Sym_tensor haij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition: star.h:1193
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
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Scalar ssjm1_chi
Effective source at the previous step for the resolution of the Poisson equation for ...
Definition: star.h:1216
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
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
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.
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: star.h:1120
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Scalar pot_centri
Centrifugal potential.
Definition: star.h:1129
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Scalar Psi_auto
Scalar field generated principally by the star.
Definition: star.h:1144
Vector ssjm1_wbeta
Effective source at the previous step for wbeta of the vector Poisson equation for wbeta (needed for ...
Definition: star.h:1234
Scalar press
Fluid pressure.
Definition: star.h:194
Vector w_beta
Solution for the vector part of the vector Poisson equation for .
Definition: star.h:1163
Scalar chi_comp
Scalar field generated principally by the companion star.
Definition: star.h:1139
Parameter storage.
Definition: param.h:125
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:522
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:905
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Scalar chi_auto
Scalar field generated principally by the star.
Definition: star.h:1134
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi...
Definition: star.h:1227
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
Scalar psi4
Conformal factor .
Definition: star.h:1158
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) ...
Definition: star.h:1177
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Scalar Psi
Total conformal factor .
Definition: star.h:1152
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
Multi-domain grid.
Definition: grilles.h:273
Scalar hacar_auto
Part of the scalar generated by beta_auto, i.e.
Definition: star.h:1205
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
virtual void reevaluate_symy(const Map *mp_prev, int nzet, Cmp &uu) const =0
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
Scalar nn
Lapse function N .
Definition: star.h:225
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:1182
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Scalar Psi_comp
Scalar field generated principally by the companion star.
Definition: star.h:1149
Scalar hacar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Definition: star.h:1211
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Scalar khi
Solution for the scalar part of the vector Poisson equation for .
Definition: star.h:1168
Scalar ssjm1_psi
Effective source at the previous step for the resolution of the Poisson equation for ...
Definition: star.h:1221
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385