LORENE
strot_dirac_diff_equil.C
1 /*
2  * Function Star_rot_Dirac_diff::equilibrium
3  *
4  * (see file star_rot_dirac_diff.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Motoyuki Saijo
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_diff_equil.C,v 1.4 2016/12/05 16:18:15 j_novak Exp $
32  *
33  */
34 
35 
36 // C headers
37 #include <cmath>
38 #include <cassert>
39 
40 // Lorene headers
41 #include "star_rot_dirac_diff.h"
42 #include "param.h"
43 #include "utilitaires.h"
44 #include "unites.h"
45 
46 namespace Lorene {
47 void Star_rot_Dirac_diff::equilibrium(double ent_c, double omega_c0,
48  double fact_omega, int , const Tbl& ent_limit,
49  const Itbl& icontrol, const Tbl& control,
50  double, double, Tbl& diff){
51 
52 
53  // Fundamental constants and units
54  // --------------------------------
55  using namespace Unites ;
56 
57  // For the display
58  // ---------------
59  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
60  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
61 
62 
63  // Grid parameters
64  // ----------------
65 
66  const Mg3d* mg = mp.get_mg() ;
67  int nz = mg->get_nzone() ; // total number of domains
68  int nzm1 = nz - 1 ;
69 
70  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
71  assert(mg->get_type_t() == SYM) ;
72  int l_b = nzet - 1 ;
73  int i_b = mg->get_nr(l_b) - 1 ;
74  int j_b = mg->get_nt(l_b) - 1 ;
75  int k_b = 0 ;
76 
77  // Value of the enthalpy defining the surface of the star
78  double ent_b = ent_limit(nzet-1) ;
79 
80  // Parameters to control the iteration
81  // -----------------------------------
82 
83  int mer_max = icontrol(0) ;
84  int mer_rot = icontrol(1) ;
85  int mer_change_omega = icontrol(2) ;
86  int mer_fix_omega = icontrol(3) ;
87  //int mer_mass = icontrol(4) ;
88  int delta_mer_kep = icontrol(5) ;
89 
90  // Protections:
91  if (mer_change_omega < mer_rot) {
92  cout << "Star_rot_Dirac_diff::equilibrium: mer_change_omega < mer_rot !"
93  << '\n' ;
94  cout << " mer_change_omega = " << mer_change_omega << '\n' ;
95  cout << " mer_rot = " << mer_rot << '\n' ;
96  abort() ;
97  }
98  if (mer_fix_omega < mer_change_omega) {
99  cout << "Star_rot_Dirac_diff::equilibrium: mer_fix_omega < mer_change_omega !"
100  << '\n' ;
101  cout << " mer_fix_omega = " << mer_fix_omega << '\n' ;
102  cout << " mer_change_omega = " << mer_change_omega << '\n' ;
103  abort() ;
104  }
105 
106  double precis = control(0) ;
107  double omega_ini = control(1) ;
108  double relax = control(2) ;
109  double relax_prev = double(1) - relax ;
110 
111  // Error indicators
112  // ----------------
113 
114  diff.set_etat_qcq() ;
115  double& diff_ent = diff.set(0) ;
116 
117  double alpha_r = 1 ;
118 
119  // Initializations
120  // ---------------
121 
122  // Initial angular velocities
123  double omega_c = 0 ;
124 
125  double accrois_omega = (omega_c0 - omega_ini) /
126  double(mer_fix_omega - mer_change_omega) ;
127 
128 
129  update_metric() ; //update of the metric quantities
130 
131  equation_of_state() ; // update of the density, pressure,...etc
132 
133  hydro_euler() ; //update of the hydro quantities relative to the
134  // Eulerian observer
135 
136  // Quantities at the previous step :
137  Scalar ent_prev = ent ;
138  Scalar logn_prev = logn ;
139  Scalar qqq_prev = qqq ;
140  // Vector beta_prev = beta ;
141  // Sym_tensor_trans hh_prev = hh ;
142 
143  // Output files
144  // -------------
145 
146  ofstream fichconv("convergence.d") ; // Output file for diff_ent
147  fichconv << "# diff_ent GRV2 max_triax vit_triax" << '\n' ;
148 
149  ofstream fichfreq("frequency.d") ; // Output file for omega
150  fichfreq << "# f [Hz]" << '\n' ;
151 
152  ofstream fichevol("evolution.d") ; // Output file for various quantities
153  fichevol <<
154  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
155  << '\n' ;
156 
157  diff_ent = 1 ;
158  double err_grv2 = 1 ;
159 
160 
161 
162 //=========================================================================
163 // Start of iteration
164 //=========================================================================
165 
166  for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++) {
167 
168  cout << "-----------------------------------------------" << '\n' ;
169  cout << "step: " << mer << '\n' ;
170  cout << "ent_c = " << display_bold << ent_c << display_normal
171  << '\n' ;
172  cout << "diff_ent = " << display_bold << diff_ent << display_normal
173  << '\n' ;
174  cout << "err_grv2 = " << err_grv2 << '\n' ;
175  fichconv << mer ;
176  fichfreq << mer ;
177  fichevol << mer ;
178 
179 
180  // switch on rotation
181  if (mer >= mer_rot) {
182 
183  if (mer < mer_change_omega) {
184  omega_c = omega_ini ;
185  }
186  else {
187  if (mer <= mer_fix_omega) {
188  omega_c = omega_ini + accrois_omega *
189  (mer - mer_change_omega) ;
190  }
191  }
192 
193 
194  }
195 
196 
197  //---------------------------------------------------//
198  // Resolution of the Poisson equation for logn //
199  // Note: ln_f is due to the fluid part //
200  // ln_q is due to the quadratic metric part //
201  //---------------------------------------------------//
202 
203  Scalar ln_f_new(mp) ;
204  Scalar ln_q_new(mp) ;
205 
206  solve_logn_f( ln_f_new ) ;
207  solve_logn_q( ln_q_new ) ;
208 
209  ln_f_new.std_spectral_base() ;
210  ln_q_new.std_spectral_base() ;
211 
212 
213  //--------------------------------------------------//
214  // Resolution of the Poisson equation for shift //
215  //--------------------------------------------------//
216 
217  Vector beta_new(mp, CON, mp.get_bvect_spher()) ;
218 
219  solve_shift( beta_new ) ;
220 
221  //------------------------------------
222  // Determination of the fluid velocity
223  //------------------------------------
224 
225  if (mer > mer_fix_omega + delta_mer_kep) {
226 
227  omega_c *= fact_omega ; // Increase of the angular velocity if
228  } // fact_omega != 1
229 
230  bool omega_trop_grand = false ;
231  bool kepler = true ;
232 
233  while ( kepler ) {
234 
235  // Possible decrease of Omega to ensure a velocity < c
236 
237  bool superlum = true ;
238 
239  while ( superlum ){
240 
241  // Computation of Omega(r,theta)
242 
243  if (omega_c == 0.) {
244  omega_field = 0 ;
245  }
246  else {
247  par_frot.set(0) = omega_c ;
248  if (par_frot(2) != double(0)) { // fixed a = R_eq / R_0
249  par_frot.set(1) = ray_eq() / par_frot(2) ;
250  }
251  double omeg_min = 0 ;
252  double omeg_max = omega_c ;
253  double precis1 = 1.e-14 ;
254  int nitermax1 = 200 ;
255 
256  fait_omega_field(omeg_min, omeg_max, precis1, nitermax1) ;
257  }
258 
259  // New fluid velocity :
260  //
261 
262  u_euler.set(1).set_etat_zero() ;
263  u_euler.set(2).set_etat_zero() ;
264 
265  u_euler.set(3) = omega_field ;
266  u_euler.set(3).annule(nzet,nzm1) ; // nzet is defined in class Star
268  u_euler.set(3).mult_rsint() ;
269  u_euler.set(3) += beta(3) ;
270  u_euler.set(3).annule(nzet,nzm1) ;
271 
272  u_euler = u_euler / nn ;
273 
274 
275  // v2 (square of norm of u_euler)
276  // -------------------------------
277 
278  v2 = contract(contract(gamma.cov(), 0, u_euler, 0), 0, u_euler, 0) ;
279 
280  // Is the new velocity larger than c in the equatorial plane ?
281 
282  superlum = false ;
283 
284  for (int l=0; l<nzet; l++) {
285  for (int i=0; i<mg->get_nr(l); i++) {
286 
287  double u2 = v2.val_grid_point(l, 0, j_b, i) ;
288  if (u2 >= 1.) { // superluminal velocity
289  superlum = true ;
290  cout << "U > c for l, i : " << l << " " << i
291  << " U = " << sqrt( u2 ) << '\n' ;
292  }
293  }
294  }
295  if ( superlum ) {
296  cout << "**** VELOCITY OF LIGHT REACHED ****" << '\n' ;
297  omega /= fact_omega ; // Decrease of Omega
298  cout << "New rotation frequency : "
299  << omega/(2.*M_PI) * f_unit << " Hz" << '\n' ;
300  omega_trop_grand = true ;
301  }
302  } // end of while ( superlum )
303 
304 
305  // New computation of U (this time is not superluminal)
306  // as well as of gam_euler, ener_euler,...etc
307 
308  hydro_euler() ;
309 
310 
311 
312  //--------------------------------//
313  // First integral of motion //
314  //--------------------------------//
315 
316  Scalar mlngamma(mp) ; // -log( gam_euler )
317 
318  mlngamma = - log( gam_euler ) ;
319 
320  // Equatorial values of various potentials :
321  double ln_f_b = ln_f_new.val_grid_point(l_b, k_b, j_b, i_b) ;
322  double ln_q_b = ln_q_new.val_grid_point(l_b, k_b, j_b, i_b) ;
323  double mlngamma_b = mlngamma.val_grid_point(l_b, k_b, j_b, i_b) ;
324  double primf_b = prim_field.val_grid_point(l_b, k_b, j_b, i_b) ;
325 
326  // Central values of various potentials :
327  double ln_f_c = ln_f_new.val_grid_point(0,0,0,0) ;
328  double ln_q_c = ln_q_new.val_grid_point(0,0,0,0) ;
329  double mlngamma_c = 0 ;
330  double primf_c = prim_field.val_grid_point(0,0,0,0) ;
331 
332  // Scale factor to ensure that the (log of) enthalpy is equal to
333  // ent_b at the equator
334  double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
335  + ln_q_c - ln_q_b + primf_c - primf_b)
336  / ( ln_f_b - ln_f_c ) ;
337 
338  alpha_r = sqrt(alpha_r2) ;
339 
340  cout << "alpha_r = " << alpha_r << '\n' ;
341 
342  // Rescaling of the grid (no adaptation!)
343  //---------------------------------------
344  mp.homothetie(alpha_r) ;
345 
346  // Readjustment of logn :
347  // -----------------------
348 
349  logn = alpha_r2 * ln_f_new + ln_q_new ;
350 
351  double logn_c = logn.val_grid_point(0,0,0,0) ;
352 
353  // First integral of motion -> (log of) enthalpy in all space
354  // ----------------------------------------------------------
355 
356  ent = (ent_c + logn_c + mlngamma_c) - logn - mlngamma - prim_field;
357 
358 
359  // --------------------------------------------------------------
360  // Test: is the enthalpy negative somewhere in the equatorial plane
361  // inside the star?
362  // --------------------------------------------------------
363 
364  kepler = false ;
365  for (int l=0; l<nzet; l++) {
366  int imax = mg->get_nr(l) - 1 ;
367  if (l == l_b) imax-- ; // The surface point is skipped
368  for (int i=0; i<imax; i++) {
369  if ( ent.val_grid_point(l, 0, j_b, i) < 0. ) {
370  kepler = true ;
371  cout << "ent < 0 for l, i : " << l << " " << i
372  << " ent = " << ent.val_grid_point(l, 0, j_b, i) << '\n' ;
373  }
374  }
375  }
376 
377  if ( kepler ) {
378  cout << "**** KEPLERIAN VELOCITY REACHED ****" << '\n' ;
379  omega /= fact_omega ; // Omega is decreased
380  cout << "New central rotation frequency : "
381  << omega_c/(2.*M_PI) * f_unit << " Hz" << '\n' ;
382  omega_trop_grand = true ;
383  }
384 
385  } // End of while ( kepler )
386 
387  if ( omega_trop_grand ) { // fact_omega is decreased for the
388  // next step
389  fact_omega = sqrt( fact_omega ) ;
390  cout << "**** New fact_omega : " << fact_omega << '\n' ;
391  }
392 
393 
394 //---------------------------------
395  // Equation of state
396  //---------------------------------
397 
398  equation_of_state() ; // computes new values for nbar (n), ener (e),
399  // and press (p) from the new ent (H)
400 
401  hydro_euler() ;
402 
403  //---------------------------------------------//
404  // Resolution of the Poisson equation for qqq //
405  //---------------------------------------------//
406 
407  Scalar q_new(mp) ;
408 
409  solve_qqq( q_new ) ;
410 
411  q_new.std_spectral_base() ;
412 
413  //----------------------------------------------//
414  // Resolution of the Poisson equation for hh //
415  //----------------------------------------------//
416 
417  Sym_tensor_trans hij_new(mp, mp.get_bvect_spher(), flat) ;
418 
419  solve_hij( hij_new ) ;
420 
421  hh = hij_new ;
422  qqq = q_new ;
423  beta = beta_new ;
424 
425  //---------------------------------------
426  // Calculate error of the GRV2 identity
427  //---------------------------------------
428 
429  err_grv2 = grv2() ;
430 
431 
432  //--------------------------------------
433  // Relaxations on some quantities....?
434  //
435  // ** On logn and qqq?
436  //--------------------------------------
437 
438  if (mer >= 10) {
439  logn = relax * logn + relax_prev * logn_prev ;
440 
441  qqq = relax * qqq + relax_prev * qqq_prev ;
442 
443  }
444 
445  // Update of the metric quantities :
446 
447  update_metric() ;
448 
449  //---------------------------
450  // Informations display
451  // More to come later......
452  //---------------------------
453 
454  // partial_display(cout) ; // What is partial_display(cout) ?
455  fichfreq << " " << omega_c / (2*M_PI) * f_unit ;
456  fichevol << " " << ent_c ;
457 
458  //-----------------------------------------------------------
459  // Relative change in enthalpy with respect to previous step
460  // ** Check: Is diffrel(ent, ent_prev) ok?
461  //-----------------------------------------------------------
462 
463  Tbl diff_ent_tbl = diffrel( ent, ent_prev ) ;
464  diff_ent = diff_ent_tbl(0) ;
465  for (int l=1; l<nzet; l++) {
466  diff_ent += diff_ent_tbl(l) ;
467  }
468  diff_ent /= nzet ;
469 
470  fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
471  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
472 
473  //------------------------------
474  // Recycling for the next step
475  //------------------------------
476 
477  ent_prev = ent ;
478  logn_prev = logn ;
479  qqq_prev = qqq ;
480 
481  fichconv << '\n' ;
482  fichfreq << '\n' ;
483  fichevol << '\n' ;
484  fichconv.flush() ;
485  fichfreq.flush() ;
486  fichevol.flush() ;
487 
488 
489  } // End of main loop
490 
491  //=================================================
492  // End of iteration
493  //=================================================
494 
495  fichconv.close() ;
496  fichfreq.close() ;
497  fichevol.close() ;
498 
499 
500 }
501 }
void solve_logn_f(Scalar &ln_f_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
double omega
Rotation angular velocity ([f_unit] )
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Map & mp
Mapping associated with the star.
Definition: star.h:180
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Metric gamma
3-metric
Definition: star.h:235
virtual void equilibrium(double ent_c, double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, double mbar_wanted, double aexp_mass, Tbl &diff)
Computes an equilibrium configuration.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Sym_tensor_trans hh
is defined by .
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
Basic integer array class.
Definition: itbl.h:122
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Scalar ent
Log-enthalpy.
Definition: star.h:190
Vector beta
Shift vector.
Definition: star.h:228
Tensor field of valence 1.
Definition: vector.h:188
void update_metric()
Computes metric quantities from known potentials.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Tbl par_frot
Parameters of the function .
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:643
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
virtual void homothetie(double lambda)=0
Sets a new radial scale.
virtual double grv2() const
Error on the virial identity GRV2.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:611
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Multi-domain grid.
Definition: grilles.h:279
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
Scalar nn
Lapse function N .
Definition: star.h:225
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:111
void solve_shift(Vector &shift_new) const
Solution of the shift equation for rotating stars in Dirac gauge.
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
void solve_qqq(Scalar &q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:465
const Metric_flat & flat
flat metric (spherical components)