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
28char 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
94namespace Lorene {
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:
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 }
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 ;
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) /
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 :
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) {
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
262
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
275
277
278 //------------------------------------
279 // Determination of the fluid velocity
280 //------------------------------------
281
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
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) ;
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 ;
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
375
376 cout << "alpha_r = " << alpha_r << endl ;
377
378 // Rescaling of the grid (no adaptation!)
379 //---------------------------------------
381
382 // Readjustment of logn :
383 // -----------------------
384
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 ;
419 }
420
421 } // End of while ( kepler )
422
423 if ( omega_trop_grand ) { // fact_omega is decreased for the
424 // next step
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
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. :
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
539 for (int l=1; l<nzet; 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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
virtual void homothetie(double lambda)=0
Sets a new radial scale.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
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
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.