LORENE
strot_dirac_equilibrium.C
1 /*
2  * Function Star_rot_Dirac::equilibrium
3  *
4  * (see file star_rot_dirac.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
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 char strot_dirac_equilibrium_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_equilibrium.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $" ;
29 
30 /*
31  * $Id: strot_dirac_equilibrium.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $
32  * $Log: strot_dirac_equilibrium.C,v $
33  * Revision 1.13 2014/10/13 08:53:40 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.12 2014/10/06 15:13:18 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.11 2009/10/26 10:54:33 j_novak
40  * Added the case of a NONSYM base in theta.
41  *
42  * Revision 1.10 2008/02/25 10:40:52 j_novak
43  * Added the flag mer_hij to control the step from which the equation for h^{ij}
44  * is being solved.
45  *
46  * Revision 1.9 2006/03/14 15:18:21 lm_lin
47  *
48  * Add convergence to a given baryon mass.
49  *
50  * Revision 1.8 2005/11/28 14:45:16 j_novak
51  * Improved solution of the Poisson tensor equation in the case of a transverse
52  * tensor.
53  *
54  * Revision 1.7 2005/09/16 14:04:49 j_novak
55  * The equation for hij is now solved only for mer > mer_fix_omega. It uses the
56  * Poisson solver of the class Sym_tensor_trans.
57  *
58  * Revision 1.6 2005/04/20 14:26:29 j_novak
59  * Removed some outputs.
60  *
61  * Revision 1.5 2005/04/05 09:24:05 j_novak
62  * minor modifs
63  *
64  * Revision 1.4 2005/03/10 09:39:19 j_novak
65  * The order of resolution has been changed in the iteration step.
66  *
67  * Revision 1.3 2005/02/17 17:30:09 f_limousin
68  * Change the name of some quantities to be consistent with other classes
69  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
70  *
71  * Revision 1.2 2005/02/09 13:36:42 lm_lin
72  *
73  * Calculate GRV2 during iterations.
74  *
75  * Revision 1.1 2005/01/31 08:51:48 j_novak
76  * New files for rotating stars in Dirac gauge (still under developement).
77  *
78  *
79  * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_equilibrium.C,v 1.13 2014/10/13 08:53:40 j_novak Exp $
80  *
81  */
82 
83 
84 // C headers
85 #include <cmath>
86 #include <cassert>
87 
88 // Lorene headers
89 #include "star_rot_dirac.h"
90 
91 #include "utilitaires.h"
92 #include "unites.h"
93 
94 namespace Lorene {
95 void Star_rot_Dirac::equilibrium(double ent_c, double omega0,
96  double fact_omega, int , const Tbl& ent_limit,
97  const Itbl& icontrol, const Tbl& control,
98  double mbar_wanted, double aexp_mass, Tbl& diff){
99 
100 
101  // Fundamental constants and units
102  // --------------------------------
103  using namespace Unites ;
104 
105  // For the display
106  // ---------------
107  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
108  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
109 
110 
111  // Grid parameters
112  // ----------------
113 
114  const Mg3d* mg = mp.get_mg() ;
115  int nz = mg->get_nzone() ; // total number of domains
116  int nzm1 = nz - 1 ;
117 
118  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
119  int type_t = mg->get_type_t() ;
120  assert( ( type_t == SYM) || (type_t == NONSYM) ) ;
121  int l_b = nzet - 1 ;
122  int i_b = mg->get_nr(l_b) - 1 ;
123  int j_b = (type_t == SYM ? mg->get_nt(l_b) - 1 : mg->get_nt(l_b)/2 ) ;
124  int k_b = 0 ;
125 
126  // Value of the enthalpy defining the surface of the star
127  double ent_b = ent_limit(nzet-1) ;
128 
129  // Parameters to control the iteration
130  // -----------------------------------
131 
132  int mer_max = icontrol(0) ;
133  int mer_rot = icontrol(1) ;
134  int mer_change_omega = icontrol(2) ;
135  int mer_fix_omega = icontrol(3) ;
136  int mer_mass = icontrol(4) ;
137  int delta_mer_kep = icontrol(5) ;
138  int mer_hij = icontrol(6) ;
139 
140  // Protections:
141  if (mer_change_omega < mer_rot) {
142  cout << "Star_rot_Dirac::equilibrium: mer_change_omega < mer_rot !"
143  << endl ;
144  cout << " mer_change_omega = " << mer_change_omega << endl ;
145  cout << " mer_rot = " << mer_rot << endl ;
146  abort() ;
147  }
148  if (mer_fix_omega < mer_change_omega) {
149  cout << "Star_rot_Dirac::equilibrium: mer_fix_omega < mer_change_omega !"
150  << endl ;
151  cout << " mer_fix_omega = " << mer_fix_omega << endl ;
152  cout << " mer_change_omega = " << mer_change_omega << endl ;
153  abort() ;
154  }
155 
156  // In order to converge to a given baryon mass, shall the central
157  // enthalpy be varied or Omega ?
158  bool change_ent = true ;
159  if (mer_mass < 0) {
160  change_ent = false ;
161  mer_mass = abs(mer_mass) ;
162  }
163 
164 
165  double precis = control(0) ;
166  double omega_ini = control(1) ;
167  double relax = control(2) ;
168  double relax_prev = double(1) - relax ;
169 
170  // Error indicators
171  // ----------------
172 
173  diff.annule_hard() ;
174  double& diff_ent = diff.set(0) ;
175 
176  double alpha_r = 1 ;
177 
178  // Initializations
179  // ---------------
180 
181  // Initial angular velocities
182  omega = 0 ;
183 
184  double accrois_omega = (omega0 - omega_ini) /
185  double(mer_fix_omega - mer_change_omega) ;
186 
187 
188  update_metric() ; //update of the metric quantities
189 
190  equation_of_state() ; // update of the density, pressure,...etc
191 
192  hydro_euler() ; //update of the hydro quantities relative to the
193  // Eulerian observer
194 
195  // Quantities at the previous step :
196  Scalar ent_prev = ent ;
197  Scalar logn_prev = logn ;
198  Scalar qqq_prev = qqq ;
199  // Vector beta_prev = beta ;
200  // Sym_tensor_trans hh_prev = hh ;
201 
202  // Output files
203  // -------------
204 
205  ofstream fichconv("convergence.d") ; // Output file for diff_ent
206  fichconv << "# diff_ent GRV2 max_triax vit_triax" << endl ;
207 
208  ofstream fichfreq("frequency.d") ; // Output file for omega
209  fichfreq << "# f [Hz]" << endl ;
210 
211  ofstream fichevol("evolution.d") ; // Output file for various quantities
212  fichevol <<
213  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
214  << endl ;
215 
216  diff_ent = 1 ;
217  double err_grv2 = 1 ;
218 
219 
220 
221 //=========================================================================
222 // Start of iteration
223 //=========================================================================
224 
225  for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++) {
226 
227  cout << "-----------------------------------------------" << endl ;
228  cout << "step: " << mer << endl ;
229  cout << "diff_ent = " << display_bold << diff_ent << display_normal
230  << endl ;
231  cout << "err_grv2 = " << err_grv2 << endl ;
232  fichconv << mer ;
233  fichfreq << mer ;
234  fichevol << mer ;
235 
236 
237  // switch on rotation
238  if (mer >= mer_rot) {
239 
240  if (mer < mer_change_omega) {
241  omega = omega_ini ;
242  }
243  else {
244  if (mer <= mer_fix_omega) {
245  omega = omega_ini + accrois_omega *
246  (mer - mer_change_omega) ;
247  }
248  }
249 
250 
251  }
252 
253 
254  //---------------------------------------------------//
255  // Resolution of the Poisson equation for logn //
256  // Note: ln_f is due to the fluid part //
257  // ln_q is due to the quadratic metric part //
258  //---------------------------------------------------//
259 
260  Scalar ln_f_new(mp) ;
261  Scalar ln_q_new(mp) ;
262 
263  solve_logn_f( ln_f_new ) ;
264  solve_logn_q( ln_q_new ) ;
265 
266  ln_f_new.std_spectral_base() ;
267  ln_q_new.std_spectral_base() ;
268 
269 
270  //--------------------------------------------------//
271  // Resolution of the Poisson equation for shift //
272  //--------------------------------------------------//
273 
274  Vector beta_new(mp, CON, mp.get_bvect_spher()) ;
275 
276  solve_shift( beta_new ) ;
277 
278  //------------------------------------
279  // Determination of the fluid velocity
280  //------------------------------------
281 
282  if (mer > mer_fix_omega + delta_mer_kep) {
283 
284  omega *= fact_omega ; // Increase of the angular velocity if
285  } // fact_omega != 1
286 
287  bool omega_trop_grand = false ;
288  bool kepler = true ;
289 
290  while ( kepler ) {
291 
292  // Possible decrease of Omega to ensure a velocity < c
293 
294  bool superlum = true ;
295 
296  while ( superlum ){
297 
298  // New fluid velocity :
299  //
300 
301  u_euler.set(1).set_etat_zero() ;
302  u_euler.set(2).set_etat_zero() ;
303 
304  u_euler.set(3) = omega ;
305  u_euler.set(3).annule(nzet,nzm1) ; // nzet is defined in class Star
307  u_euler.set(3).mult_rsint() ;
308  u_euler.set(3) += beta(3) ;
309  u_euler.set(3).annule(nzet,nzm1) ;
310 
311  u_euler = u_euler / nn ;
312 
313 
314  // v2 (square of norm of u_euler)
315  // -------------------------------
316 
317  v2 = contract(contract(gamma.cov(), 0, u_euler, 0), 0, u_euler, 0) ;
318 
319  // Is the new velocity larger than c in the equatorial plane ?
320 
321  superlum = false ;
322 
323  for (int l=0; l<nzet; l++) {
324  for (int i=0; i<mg->get_nr(l); i++) {
325 
326  double u2 = v2.val_grid_point(l, 0, j_b, i) ;
327  if (u2 >= 1.) { // superluminal velocity
328  superlum = true ;
329  cout << "U > c for l, i : " << l << " " << i
330  << " U = " << sqrt( u2 ) << endl ;
331  }
332  }
333  }
334  if ( superlum ) {
335  cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
336  omega /= fact_omega ; // Decrease of Omega
337  cout << "New rotation frequency : "
338  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
339  omega_trop_grand = true ;
340  }
341  } // end of while ( superlum )
342 
343 
344  // New computation of U (this time is not superluminal)
345  // as well as of gam_euler, ener_euler,...etc
346 
347  hydro_euler() ;
348 
349 
350 
351  //--------------------------------//
352  // First integral of motion //
353  //--------------------------------//
354 
355  Scalar mlngamma(mp) ; // -log( gam_euler )
356 
357  mlngamma = - log( gam_euler ) ;
358 
359  // Equatorial values of various potentials :
360  double ln_f_b = ln_f_new.val_grid_point(l_b, k_b, j_b, i_b) ;
361  double ln_q_b = ln_q_new.val_grid_point(l_b, k_b, j_b, i_b) ;
362  double mlngamma_b = mlngamma.val_grid_point(l_b, k_b, j_b, i_b) ;
363 
364  // Central values of various potentials :
365  double ln_f_c = ln_f_new.val_grid_point(0,0,0,0) ;
366  double ln_q_c = ln_q_new.val_grid_point(0,0,0,0) ;
367  double mlngamma_c = 0 ;
368 
369  // Scale factor to ensure that the (log of) enthalpy is equal to
370  // ent_b at the equator
371  double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
372  + ln_q_c - ln_q_b) / ( ln_f_b - ln_f_c ) ;
373 
374  alpha_r = sqrt(alpha_r2) ;
375 
376  cout << "alpha_r = " << alpha_r << endl ;
377 
378  // Rescaling of the grid (no adaptation!)
379  //---------------------------------------
380  mp.homothetie(alpha_r) ;
381 
382  // Readjustment of logn :
383  // -----------------------
384 
385  logn = alpha_r2 * ln_f_new + ln_q_new ;
386 
387  double logn_c = logn.val_grid_point(0,0,0,0) ;
388 
389  // First integral of motion -> (log of) enthalpy in all space
390  // ----------------------------------------------------------
391 
392  ent = (ent_c + logn_c + mlngamma_c) - logn - mlngamma ;
393 
394 
395  // --------------------------------------------------------------
396  // Test: is the enthalpy negative somewhere in the equatorial plane
397  // inside the star?
398  // --------------------------------------------------------
399 
400  kepler = false ;
401  for (int l=0; l<nzet; l++) {
402  int imax = mg->get_nr(l) - 1 ;
403  if (l == l_b) imax-- ; // The surface point is skipped
404  for (int i=0; i<imax; i++) {
405  if ( ent.val_grid_point(l, 0, j_b, i) < 0. ) {
406  kepler = true ;
407  cout << "ent < 0 for l, i : " << l << " " << i
408  << " ent = " << ent.val_grid_point(l, 0, j_b, i) << endl ;
409  }
410  }
411  }
412 
413  if ( kepler ) {
414  cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ;
415  omega /= fact_omega ; // Omega is decreased
416  cout << "New rotation frequency : "
417  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
418  omega_trop_grand = true ;
419  }
420 
421  } // End of while ( kepler )
422 
423  if ( omega_trop_grand ) { // fact_omega is decreased for the
424  // next step
425  fact_omega = sqrt( fact_omega ) ;
426  cout << "**** New fact_omega : " << fact_omega << endl ;
427  }
428 
429 
430 //---------------------------------
431  // Equation of state
432  //---------------------------------
433 
434  equation_of_state() ; // computes new values for nbar (n), ener (e),
435  // and press (p) from the new ent (H)
436 
437  hydro_euler() ;
438 
439  //---------------------------------------------//
440  // Resolution of the Poisson equation for qqq //
441  //---------------------------------------------//
442 
443  Scalar q_new(mp) ;
444 
445  solve_qqq( q_new ) ;
446 
447  q_new.std_spectral_base() ;
448 
449  //----------------------------------------------//
450  // Resolution of the Poisson equation for hh //
451  //----------------------------------------------//
452 
453  Sym_tensor_trans hij_new(mp, mp.get_bvect_spher(), flat) ;
454 
455  if (mer > mer_hij )
456  solve_hij( hij_new ) ;
457  else
458  hij_new.set_etat_zero() ;
459 
460  hh = hij_new ;
461  qqq = q_new ;
462  beta = beta_new ;
463 
464  //---------------------------------------
465  // Calculate error of the GRV2 identity
466  //---------------------------------------
467 
468  err_grv2 = grv2() ;
469 
470 
471  //--------------------------------------
472  // Relaxations on some quantities....?
473  //
474  // ** On logn and qqq?
475  //--------------------------------------
476 
477  if (mer >= 10) {
478  logn = relax * logn + relax_prev * logn_prev ;
479 
480  qqq = relax * qqq + relax_prev * qqq_prev ;
481 
482  }
483 
484  // Update of the metric quantities :
485 
486  update_metric() ;
487 
488  //---------------------------
489  // Informations display
490  // More to come later......
491  //---------------------------
492 
493  // partial_display(cout) ; // What is partial_display(cout) ?
494  fichfreq << " " << omega / (2*M_PI) * f_unit ;
495  fichevol << " " << ent_c ;
496 
497 
498  //-----------------------------------------
499  // Convergence towards a given baryon mass
500  //-----------------------------------------
501 
502  if (mer > mer_mass) {
503 
504  double xx ;
505  if (mbar_wanted > 0.) {
506  xx = mass_b() / mbar_wanted - 1. ;
507  cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx
508  << endl ;
509  }
510  else{
511  xx = mass_g() / fabs(mbar_wanted) - 1. ;
512  cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx
513  << endl ;
514  }
515  double xprog = ( mer > 2*mer_mass) ? 1. :
516  double(mer-mer_mass)/double(mer_mass) ;
517  xx *= xprog ;
518  double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
519  double fact = pow(ax, aexp_mass) ;
520  cout << " xprog, xx, ax, fact : " << xprog << " " <<
521  xx << " " << ax << " " << fact << endl ;
522 
523  if ( change_ent ) {
524  ent_c *= fact ;
525  }
526  else {
527  if (mer%4 == 0) omega *= fact ;
528  }
529  }
530 
531 
532  //-----------------------------------------------------------
533  // Relative change in enthalpy with respect to previous step
534  // ** Check: Is diffrel(ent, ent_prev) ok?
535  //-----------------------------------------------------------
536 
537  Tbl diff_ent_tbl = diffrel( ent, ent_prev ) ;
538  diff_ent = diff_ent_tbl(0) ;
539  for (int l=1; l<nzet; l++) {
540  diff_ent += diff_ent_tbl(l) ;
541  }
542  diff_ent /= nzet ;
543 
544  fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
545  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
546 
547  //------------------------------
548  // Recycling for the next step
549  //------------------------------
550 
551  ent_prev = ent ;
552  logn_prev = logn ;
553  qqq_prev = qqq ;
554 
555  fichconv << endl ;
556  fichfreq << endl ;
557  fichevol << endl ;
558  fichconv.flush() ;
559  fichfreq.flush() ;
560  fichevol.flush() ;
561 
562 
563  } // End of main loop
564 
565  //=================================================
566  // End of iteration
567  //=================================================
568 
569  fichconv.close() ;
570  fichfreq.close() ;
571  fichevol.close() ;
572 
573 
574 }
575 }
Basic integer array class.
Definition: itbl.h:122
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
virtual void homothetie(double lambda)=0
Sets a new radial scale.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
Multi-domain grid.
Definition: grilles.h:273
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:485
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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
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 set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
virtual double mass_g() const
Gravitational mass.
Sym_tensor_trans hh
is defined by .
double omega
Rotation angular velocity ([f_unit] )
const Metric_flat & flat
flat metric (spherical components)
void solve_logn_f(Scalar &ln_f_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.
void update_metric()
Computes metric quantities from known potentials.
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
virtual double mass_b() const
Baryonic mass.
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.
void solve_shift(Vector &shift_new) const
Solution of the shift equation for rotating stars in Dirac gauge.
virtual double grv2() const
Error on the virial identity GRV2.
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Scalar nn
Lapse function N .
Definition: star.h:225
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
Metric gamma
3-metric
Definition: star.h:235
Scalar ent
Log-enthalpy.
Definition: star.h:190
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Map & mp
Mapping associated with the star.
Definition: star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Vector beta
Shift vector.
Definition: star.h:228
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:608
Basic array class.
Definition: tbl.h:161
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 1.
Definition: vector.h:188
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.