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
28char strot_dirac_diff_equil_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_diff_equil.C,v 1.3 2014/10/13 08:53:39 j_novak Exp $" ;
29
30/*
31 * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_diff_equil.C,v 1.3 2014/10/13 08:53:39 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
46namespace Lorene {
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:
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 }
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) /
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 :
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) {
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
205
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
218
220
221 //------------------------------------
222 // Determination of the fluid velocity
223 //------------------------------------
224
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
257 }
258
259 // New fluid velocity :
260 //
261
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) ;
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' ;
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) ;
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
339
340 cout << "alpha_r = " << alpha_r << '\n' ;
341
342 // Rescaling of the grid (no adaptation!)
343 //---------------------------------------
345
346 // Readjustment of logn :
347 // -----------------------
348
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
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' ;
383 }
384
385 } // End of while ( kepler )
386
387 if ( omega_trop_grand ) { // fact_omega is decreased for the
388 // next step
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
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
465 for (int l=1; l<nzet; 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}
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.
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
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.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tbl par_frot
Parameters of the function .
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.
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
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
double ray_eq() const
Coordinate radius at , [r_unit].
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
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 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.