LORENE
star_rot_equil.C
1/*
2 * Method Star_rot::equilibrium
3 *
4 * (see file star_h.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2010 Eric Gourgoulhon
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 as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char star_rot_equil_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_rot_equil.C,v 1.6 2014/10/13 08:53:39 j_novak Exp $" ;
31
32/*
33 * $Id: star_rot_equil.C,v 1.6 2014/10/13 08:53:39 j_novak Exp $
34 * $Log: star_rot_equil.C,v $
35 * Revision 1.6 2014/10/13 08:53:39 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.5 2014/10/06 15:13:17 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.4 2010/01/26 16:49:53 e_gourgoulhon
42 * Reformated some outputs to the screen.
43 *
44 * Revision 1.3 2010/01/26 10:49:43 e_gourgoulhon
45 * First working version.
46 *
47 * Revision 1.2 2010/01/25 22:32:38 e_gourgoulhon
48 * Start of real implementation
49 *
50 * Revision 1.1 2010/01/25 18:15:52 e_gourgoulhon
51 * First version.
52 *
53 *
54 * $Header: /cvsroot/Lorene/C++/Source/Star/star_rot_equil.C,v 1.6 2014/10/13 08:53:39 j_novak Exp $
55 *
56 */
57
58// Headers C
59#include <cmath>
60
61// Headers Lorene
62#include "star_rot.h"
63#include "param.h"
64#include "tenseur.h"
65
66#include "graphique.h"
67#include "utilitaires.h"
68#include "unites.h"
69
70namespace Lorene {
71void Star_rot::equilibrium(double ent_c, double omega0, double fact_omega,
72 int nzadapt, const Tbl& ent_limit, const Itbl& icontrol,
73 const Tbl& control, double mbar_wanted,
74 double aexp_mass, Tbl& diff, Param*) {
75
76 // Fundamental constants and units
77 // -------------------------------
78
79 using namespace Unites ;
80
81 // For the display
82 // ---------------
83 char display_bold[]="x[1m" ; display_bold[0] = 27 ;
84 char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
85
86 // Grid parameters
87 // ---------------
88
89 const Mg3d* mg = mp.get_mg() ;
90 int nz = mg->get_nzone() ; // total number of domains
91 int nzm1 = nz - 1 ;
92
93 // The following is required to initialize mp_prev as a Map_et:
94 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
95
96 // Index of the point at phi=0, theta=pi/2 at the surface of the star:
97 assert(mg->get_type_t() == SYM) ;
98 int l_b = nzet - 1 ;
99 int i_b = mg->get_nr(l_b) - 1 ;
100 int j_b = mg->get_nt(l_b) - 1 ;
101 int k_b = 0 ;
102
103 // Value of the enthalpy defining the surface of the star
104 double ent_b = ent_limit(nzet-1) ;
105
106 // Parameters to control the iteration
107 // -----------------------------------
108
109 int mer_max = icontrol(0) ;
110 int mer_rot = icontrol(1) ;
111 int mer_change_omega = icontrol(2) ;
112 int mer_fix_omega = icontrol(3) ;
113 int mer_mass = icontrol(4) ;
114 int mermax_poisson = icontrol(5) ;
115 int mer_triax = icontrol(6) ;
116 int delta_mer_kep = icontrol(7) ;
117
118 // Protections:
120 cout << "Star_rot::equilibrium: mer_change_omega < mer_rot !" << endl ;
121 cout << " mer_change_omega = " << mer_change_omega << endl ;
122 cout << " mer_rot = " << mer_rot << endl ;
123 abort() ;
124 }
126 cout << "Star_rot::equilibrium: mer_fix_omega < mer_change_omega !"
127 << endl ;
128 cout << " mer_fix_omega = " << mer_fix_omega << endl ;
129 cout << " mer_change_omega = " << mer_change_omega << endl ;
130 abort() ;
131 }
132
133 // In order to converge to a given baryon mass, shall the central
134 // enthalpy be varied or Omega ?
135 bool change_ent = true ;
136 if (mer_mass < 0) {
137 change_ent = false ;
139 }
140
141 double precis = control(0) ;
142 double omega_ini = control(1) ;
143 double relax = control(2) ;
144 double relax_prev = double(1) - relax ;
145 double relax_poisson = control(3) ;
146 double thres_adapt = control(4) ;
147 double ampli_triax = control(5) ;
148 double precis_adapt = control(6) ;
149
150
151 // Error indicators
152 // ----------------
153
154 diff.set_etat_qcq() ;
155 double& diff_ent = diff.set(0) ;
156 double& diff_nuf = diff.set(1) ;
157 double& diff_nuq = diff.set(2) ;
158// double& diff_dzeta = diff.set(3) ;
159// double& diff_ggg = diff.set(4) ;
160 double& diff_shift_x = diff.set(5) ;
161 double& diff_shift_y = diff.set(6) ;
162 double& vit_triax = diff.set(7) ;
163
164 // Parameters for the function Map_et::adapt
165 // -----------------------------------------
166
168 int nitermax = 100 ;
169 int niter ;
170 int adapt_flag = 1 ; // 1 = performs the full computation,
171 // 0 = performs only the rescaling by
172 // the factor alpha_r
173 int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
174 // isosurfaces
175 double alpha_r ;
176 double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
177
178 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
179 // locate zeros by the secant method
180 par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
181 // to the isosurfaces of ent is to be
182 // performed
183 par_adapt.add_int(nz_search, 2) ; // number of domains to search for
184 // the enthalpy isosurface
185 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
186 // 0 = performs only the rescaling by
187 // the factor alpha_r
188 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
189 // (theta_*, phi_*)
190 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
191 // (theta_*, phi_*)
192
193 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
194 // the secant method
195
196 par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
197 // the determination of zeros by
198 // the secant method
199 par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
200
201 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
202 // distances will be multiplied
203
204 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
205 // to define the isosurfaces.
206
207 // Parameters for the function Map_et::poisson for nuf
208 // ----------------------------------------------------
209
210 double precis_poisson = 1.e-16 ;
211
217 cssjm1_wshift.set_etat_qcq() ;
218 for (int i=1; i<=3; i++) {
219 cssjm1_wshift.set(i-1) = ssjm1_wshift(i) ;
220 }
221
222 Tenseur cshift(mp, 1, CON, mp.get_bvect_cart() ) ;
223 cshift.set_etat_qcq() ;
224 for (int i=1; i<=3; i++) {
225 cshift.set(i-1) = -beta(i) ;
226 }
227
228 Tenseur cw_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
229 cw_shift.set_etat_qcq() ;
230 for (int i=1; i<=3; i++) {
231 cw_shift.set(i-1) = w_shift(i) ;
232 }
233
235 ckhi_shift.set_etat_qcq() ;
236 ckhi_shift.set() = khi_shift ;
237
239 par_poisson_nuf.add_int(mermax_poisson, 0) ; // maximum number of iterations
240 par_poisson_nuf.add_double(relax_poisson, 0) ; // relaxation parameter
241 par_poisson_nuf.add_double(precis_poisson, 1) ; // required precision
242 par_poisson_nuf.add_int_mod(niter, 0) ; // number of iterations actually used
243 par_poisson_nuf.add_cmp_mod( cssjm1_nuf ) ;
244
246 par_poisson_nuq.add_int(mermax_poisson, 0) ; // maximum number of iterations
247 par_poisson_nuq.add_double(relax_poisson, 0) ; // relaxation parameter
248 par_poisson_nuq.add_double(precis_poisson, 1) ; // required precision
249 par_poisson_nuq.add_int_mod(niter, 0) ; // number of iterations actually used
250 par_poisson_nuq.add_cmp_mod( cssjm1_nuq ) ;
251
253 par_poisson_tggg.add_int(mermax_poisson, 0) ; // maximum number of iterations
254 par_poisson_tggg.add_double(relax_poisson, 0) ; // relaxation parameter
255 par_poisson_tggg.add_double(precis_poisson, 1) ; // required precision
256 par_poisson_tggg.add_int_mod(niter, 0) ; // number of iterations actually used
257 par_poisson_tggg.add_cmp_mod( cssjm1_tggg ) ;
258 double lambda_tggg ;
259 par_poisson_tggg.add_double_mod( lambda_tggg ) ;
260
262 double lbda_grv2 ;
263 par_poisson_dzeta.add_double_mod( lbda_grv2 ) ;
264
265 // Parameters for the function Scalar::poisson_vect
266 // -------------------------------------------------
267
269
270 par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of iterations
271 par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
272 par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
273 par_poisson_vect.add_cmp_mod( cssjm1_khi ) ;
274 par_poisson_vect.add_tenseur_mod( cssjm1_wshift ) ;
275 par_poisson_vect.add_int_mod(niter, 0) ;
276
277
278 // Initializations
279 // ---------------
280
281 // Initial angular velocity
282 omega = 0 ;
283
284 double accrois_omega = (omega0 - omega_ini) /
286
287
288 update_metric() ; // update of the metric coefficients
289
290 equation_of_state() ; // update of the density, pressure, etc...
291
292 hydro_euler() ; // update of the hydro quantities relative to the
293 // Eulerian observer
294
295 // Quantities at the previous step :
296 Map_et mp_prev = mp_et ;
297 Scalar ent_prev = ent ;
300
301 // Creation of uninitialized tensors:
302 Scalar source_nuf(mp) ; // source term in the equation for nuf
303 Scalar source_nuq(mp) ; // source term in the equation for nuq
304 Scalar source_dzf(mp) ; // matter source term in the eq. for dzeta
305 Scalar source_dzq(mp) ; // quadratic source term in the eq. for dzeta
306 Scalar source_tggg(mp) ; // source term in the eq. for tggg
308 // source term for shift
309 Scalar mlngamma(mp) ; // centrifugal potential
310
311
312 ofstream fichconv("convergence.d") ; // Output file for diff_ent
313 fichconv << "# diff_ent GRV2 max_triax vit_triax" << endl ;
314
315 ofstream fichfreq("frequency.d") ; // Output file for omega
316 fichfreq << "# f [Hz]" << endl ;
317
318 ofstream fichevol("evolution.d") ; // Output file for various quantities
319 fichevol <<
320 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
321 << endl ;
322
323 diff_ent = 1 ;
324 double err_grv2 = 1 ;
325 double max_triax_prev = 0 ; // Triaxial amplitude at previous step
326
327 //=========================================================================
328 // Start of iteration
329 //=========================================================================
330
331 for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
332
333 cout << "-----------------------------------------------" << endl ;
334 cout << "step: " << mer << endl ;
335 cout << "diff_ent = " << display_bold << diff_ent << display_normal
336 << endl ;
337 cout << "err_grv2 = " << err_grv2 << endl ;
338 fichconv << mer ;
339 fichfreq << mer ;
340 fichevol << mer ;
341
342 if (mer >= mer_rot) {
343
344 if (mer < mer_change_omega) {
345 omega = omega_ini ;
346 }
347 else {
348 if (mer <= mer_fix_omega) {
351 }
352 }
353
354 }
355
356 //-----------------------------------------------
357 // Sources of the Poisson equations
358 //-----------------------------------------------
359
360 // Source for nu
361 // -------------
362 Scalar bet = log(bbb) ;
363 bet.std_spectral_base() ;
364
366 Vector d_bet = bet.derive_cov( mp.flat_met_spher() ) ;
367
368 if (relativistic) {
369 source_nuf = qpig * a_car *( ener_euler + s_euler ) ;
370
371 source_nuq = ak_car - d_logn(1)*(d_logn(1)+d_bet(1))
372 - d_logn(2)*(d_logn(2)+d_bet(2))
373 - d_logn(3)*(d_logn(3)+d_bet(3)) ;
374 }
375 else {
376 source_nuf = qpig * nbar ;
377
378 source_nuq = 0 ;
379 }
380 source_nuf.std_spectral_base() ;
381 source_nuq.std_spectral_base() ;
382
383 // Source for dzeta
384 // ----------------
385 source_dzf = 2 * qpig * a_car * (press + (ener_euler+press) * uuu*uuu ) ;
386 source_dzf.std_spectral_base() ;
387
388 source_dzq = 1.5 * ak_car
389 - d_logn(1)*d_logn(1) - d_logn(2)*d_logn(2) - d_logn(3)*d_logn(3) ;
390 source_dzq.std_spectral_base() ;
391
392 // Source for tggg
393 // ---------------
394
395 source_tggg = 4 * qpig * nn * a_car * bbb * press ;
396 source_tggg.std_spectral_base() ;
397
398 source_tggg.mult_rsint() ;
399
400
401 // Source for shift
402 // ----------------
403
404 // Matter term:
405 source_shift = (-4*qpig) * nn * a_car * (ener_euler + press)
406 * u_euler ;
407
408 // Quadratic terms:
409 Vector vtmp = 6 * bet.derive_con( mp.flat_met_spher() )
410 - 2 * logn.derive_con( mp.flat_met_spher() ) ;
411 vtmp.change_triad(mp.get_bvect_cart()) ;
412
413 Vector squad = nn * contract(tkij, 1, vtmp, 0) ;
414
416
417 //----------------------------------------------
418 // Resolution of the Poisson equation for nuf
419 //----------------------------------------------
420
421 source_nuf.poisson(par_poisson_nuf, nuf) ;
422
423 cout << "Test of the Poisson equation for nuf :" << endl ;
424 Tbl err = source_nuf.test_poisson(nuf, cout, true) ;
425 diff_nuf = err(0, 0) ;
426
427 //---------------------------------------
428 // Triaxial perturbation of nuf
429 //---------------------------------------
430
431 if (mer == mer_triax) {
432
433 if ( mg->get_np(0) == 1 ) {
434 cout <<
435 "Star_rot::equilibrium: np must be stricly greater than 1"
436 << endl << " to set a triaxial perturbation !" << endl ;
437 abort() ;
438 }
439
440 const Coord& phi = mp.phi ;
441 const Coord& sint = mp.sint ;
442 Scalar perturb(mp) ;
443 perturb = 1 + ampli_triax * sint*sint * cos(2*phi) ;
444 nuf = nuf * perturb ;
445
446 nuf.std_spectral_base() ; // set the bases for spectral expansions
447 // to be the standard ones for a
448 // scalar field
449
450 }
451
452 // Monitoring of the triaxial perturbation
453 // ---------------------------------------
454
455 const Valeur& va_nuf = nuf.get_spectral_va() ;
456 va_nuf.coef() ; // Computes the spectral coefficients
457 double max_triax = 0 ;
458
459 if ( mg->get_np(0) > 1 ) {
460
461 for (int l=0; l<nz; l++) { // loop on the domains
462 for (int j=0; j<mg->get_nt(l); j++) {
463 for (int i=0; i<mg->get_nr(l); i++) {
464
465 // Coefficient of cos(2 phi) :
466 double xcos2p = (*(va_nuf.c_cf))(l, 2, j, i) ;
467
468 // Coefficient of sin(2 phi) :
469 double xsin2p = (*(va_nuf.c_cf))(l, 3, j, i) ;
470
471 double xx = sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ;
472
473 max_triax = ( xx > max_triax ) ? xx : max_triax ;
474 }
475 }
476 }
477
478 }
479
480 cout << "Triaxial part of nuf : " << max_triax << endl ;
481
482 if (relativistic) {
483
484 //----------------------------------------------
485 // Resolution of the Poisson equation for nuq
486 //----------------------------------------------
487
488 source_nuq.poisson(par_poisson_nuq, nuq) ;
489
490 cout << "Test of the Poisson equation for nuq :" << endl ;
491 err = source_nuq.test_poisson(nuq, cout, true) ;
492 diff_nuq = err(0, 0) ;
493
494 //---------------------------------------------------------
495 // Resolution of the vector Poisson equation for the shift
496 //---------------------------------------------------------
497
498
499 for (int i=1; i<=3; i++) {
500 if(source_shift(i).get_etat() != ETATZERO) {
501 if(source_shift(i).dz_nonzero()) {
502 assert( source_shift(i).get_dzpuis() == 4 ) ;
503 }
504 else{
505 (source_shift.set(i)).set_dzpuis(4) ;
506 }
507 }
508 }
509
510 double lambda_shift = double(1) / double(3) ;
511
512 if ( mg->get_np(0) == 1 ) {
513 lambda_shift = 0 ;
514 }
515
517 csource_shift.set_etat_qcq() ;
518 for (int i=1; i<=3; i++) {
519 csource_shift.set(i-1) = source_shift(i) ;
520 }
521 csource_shift.set(2).set_etat_zero() ; //## bizarre...
522
525
526 for (int i=1; i<=3; i++) {
527 beta.set(i) = - cshift(i-1) ;
528 beta.set(i).set_dzpuis(0) ; //## bizarre...
529 w_shift.set(i) = cw_shift(i-1) ;
530 }
531 khi_shift = ckhi_shift() ;
532
533 cout << "Test of the Poisson equation for shift_x :" << endl ;
534 err = source_shift(1).test_poisson(-beta(1), cout, true) ;
535 diff_shift_x = err(0, 0) ;
536
537 cout << "Test of the Poisson equation for shift_y :" << endl ;
538 err = source_shift(2).test_poisson(-beta(2), cout, true) ;
539 diff_shift_y = err(0, 0) ;
540
541 // Computation of tnphi and nphi from the Cartesian components
542 // of the shift
543 // -----------------------------------------------------------
544
545 fait_nphi() ;
546
547 }
548
549 //-----------------------------------------
550 // Determination of the fluid velociy U
551 //-----------------------------------------
552
554
555 omega *= fact_omega ; // Increase of the angular velocity if
556 } // fact_omega != 1
557
558 bool omega_trop_grand = false ;
559 bool kepler = true ;
560
561 while ( kepler ) {
562
563 // Possible decrease of Omega to ensure a velocity < c
564
565 bool superlum = true ;
566
567 while ( superlum ) {
568
569 // New fluid velocity U :
570
571 Scalar tmp = omega - nphi ;
572 tmp.annule_domain(nzm1) ;
573 tmp.std_spectral_base() ;
574
575 tmp.mult_rsint() ; // Multiplication by r sin(theta)
576
577 uuu = bbb / nn * tmp ;
578
579 if (uuu.get_etat() == ETATQCQ) {
580 // Same basis as (Omega -N^phi) r sin(theta) :
581 (uuu.set_spectral_va()).set_base( tmp.get_spectral_va().get_base() ) ;
582 }
583
584 // Is the new velocity larger than c in the equatorial plane ?
585
586 superlum = false ;
587
588 for (int l=0; l<nzet; l++) {
589 for (int i=0; i<mg->get_nr(l); i++) {
590
591 double u1 = uuu.val_grid_point(l, 0, j_b, i) ;
592 if (u1 >= 1.) { // superluminal velocity
593 superlum = true ;
594 cout << "U > c for l, i : " << l << " " << i
595 << " U = " << u1 << endl ;
596 }
597 }
598 }
599 if ( superlum ) {
600 cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
601 omega /= fact_omega ; // Decrease of Omega
602 cout << "New rotation frequency : "
603 << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
605 }
606 } // end of while ( superlum )
607
608
609 // New computation of U (which this time is not superluminal)
610 // as well as of gam_euler, ener_euler, etc...
611 // -----------------------------------
612
613 hydro_euler() ;
614
615
616 //------------------------------------------------------
617 // First integral of motion
618 //------------------------------------------------------
619
620 // Centrifugal potential :
621 if (relativistic) {
622 mlngamma = - log( gam_euler ) ;
623 }
624 else {
625 mlngamma = - 0.5 * uuu*uuu ;
626 }
627
628 // Equatorial values of various potentials :
629 double nuf_b = nuf.val_grid_point(l_b, k_b, j_b, i_b) ;
630 double nuq_b = nuq.val_grid_point(l_b, k_b, j_b, i_b) ;
631 double mlngamma_b = mlngamma.val_grid_point(l_b, k_b, j_b, i_b) ;
632
633 // Central values of various potentials :
634 double nuf_c = nuf.val_grid_point(0,0,0,0) ;
635 double nuq_c = nuq.val_grid_point(0,0,0,0) ;
636 double mlngamma_c = 0 ;
637
638 // Scale factor to ensure that the enthalpy is equal to ent_b at
639 // the equator
640 double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
641 + nuq_c - nuq_b) / ( nuf_b - nuf_c ) ;
643 cout << "alpha_r = " << alpha_r << endl ;
644
645 // Readjustment of nu :
646 // -------------------
647
648 logn = alpha_r2 * nuf + nuq ;
649 double nu_c = logn.val_grid_point(0,0,0,0) ;
650
651 // First integral --> enthalpy in all space
652 //-----------------
653
654 ent = (ent_c + nu_c + mlngamma_c) - logn - mlngamma ;
655
656 // Test: is the enthalpy negative somewhere in the equatorial plane
657 // inside the star ? If yes, this means that the Keplerian velocity
658 // has been overstep.
659
660 kepler = false ;
661 for (int l=0; l<nzet; l++) {
662 int imax = mg->get_nr(l) - 1 ;
663 if (l == l_b) imax-- ; // The surface point is skipped
664 for (int i=0; i<imax; i++) {
665 if ( ent.val_grid_point(l, 0, j_b, i) < 0. ) {
666 kepler = true ;
667 cout << "ent < 0 for l, i : " << l << " " << i
668 << " ent = " << ent.val_grid_point(l, 0, j_b, i) << endl ;
669 }
670 }
671 }
672
673 if ( kepler ) {
674 cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ;
675 omega /= fact_omega ; // Omega is decreased
676 cout << "New rotation frequency : "
677 << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
679 }
680
681 } // End of while ( kepler )
682
683 if ( omega_trop_grand ) { // fact_omega is decreased for the
684 // next step
686 cout << "**** New fact_omega : " << fact_omega << endl ;
687 }
688
689 //----------------------------------------------------
690 // Adaptation of the mapping to the new enthalpy field
691 //----------------------------------------------------
692
693 // Shall the adaptation be performed (cusp) ?
694 // ------------------------------------------
695
696 double dent_eq = ent.dsdr().val_grid_point(l_b, k_b, j_b, i_b) ;
697 double dent_pole = ent.dsdr().val_grid_point(l_b, k_b, 0, i_b) ;
698 double rap_dent = fabs( dent_eq / dent_pole ) ;
699 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
700
701 if ( rap_dent < thres_adapt ) {
702 adapt_flag = 0 ; // No adaptation of the mapping
703 cout << "******* FROZEN MAPPING *********" << endl ;
704 }
705 else{
706 adapt_flag = 1 ; // The adaptation of the mapping is to be
707 // performed
708 }
709
710 mp_prev = mp_et ;
711
712 Cmp cent(ent) ;
714
715 //----------------------------------------------------
716 // Computation of the enthalpy at the new grid points
717 //----------------------------------------------------
718
719 mp_prev.homothetie(alpha_r) ;
720
722 ent = cent ;
723
724 //----------------------------------------------------
725 // Equation of state
726 //----------------------------------------------------
727
728 equation_of_state() ; // computes new values for nbar (n), ener (e)
729 // and press (p) from the new ent (H)
730
731 //---------------------------------------------------------
732 // Matter source terms in the gravitational field equations
733 //---------------------------------------------------------
734
735 //## Computation of tnphi and nphi from the Cartesian components
736 // of the shift for the test in hydro_euler():
737
738 fait_nphi() ;
739
740 hydro_euler() ; // computes new values for ener_euler (E),
741 // s_euler (S) and u_euler (U^i)
742
743 if (relativistic) {
744
745 //-------------------------------------------------------
746 // 2-D Poisson equation for tggg
747 //-------------------------------------------------------
748
750 Cmp ctggg(tggg) ;
752 ctggg) ;
753 tggg = ctggg ;
754
755
756 //-------------------------------------------------------
757 // 2-D Poisson equation for dzeta
758 //-------------------------------------------------------
759
762 Cmp cdzeta(dzeta) ;
764 cdzeta) ;
765 dzeta = cdzeta ;
766
767 err_grv2 = lbda_grv2 - 1;
768 cout << "GRV2: " << err_grv2 << endl ;
769
770 }
771 else {
772 err_grv2 = grv2() ;
773 }
774
775
776 //---------------------------------------
777 // Computation of the metric coefficients (except for N^phi)
778 //---------------------------------------
779
780 // Relaxations on nu and dzeta :
781
782 if (mer >= 10) {
783 logn = relax * logn + relax_prev * logn_prev ;
784
785 dzeta = relax * dzeta + relax_prev * dzeta_prev ;
786 }
787
788 // Update of the metric coefficients N, A, B and computation of K_ij :
789
790 update_metric() ;
791
792 //-----------------------
793 // Informations display
794 //-----------------------
795
797 fichfreq << " " << omega / (2*M_PI) * f_unit ;
798 fichevol << " " << rap_dent ;
799 fichevol << " " << ray_pole() / ray_eq() ;
800 fichevol << " " << ent_c ;
801
802 //-----------------------------------------
803 // Convergence towards a given baryon mass
804 //-----------------------------------------
805
806 if (mer > mer_mass) {
807
808 double xx ;
809 if (mbar_wanted > 0.) {
810 xx = mass_b() / mbar_wanted - 1. ;
811 cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx
812 << endl ;
813 }
814 else{
815 xx = mass_g() / fabs(mbar_wanted) - 1. ;
816 cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx
817 << endl ;
818 }
819 double xprog = ( mer > 2*mer_mass) ? 1. :
821 xx *= xprog ;
822 double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
823 double fact = pow(ax, aexp_mass) ;
824 cout << " xprog, xx, ax, fact : " << xprog << " " <<
825 xx << " " << ax << " " << fact << endl ;
826
827 if ( change_ent ) {
828 ent_c *= fact ;
829 }
830 else {
831 if (mer%4 == 0) omega *= fact ;
832 }
833 }
834
835
836 //------------------------------------------------------------
837 // Relative change in enthalpy with respect to previous step
838 //------------------------------------------------------------
839
841 diff_ent = diff_ent_tbl(0) ;
842 for (int l=1; l<nzet; l++) {
843 diff_ent += diff_ent_tbl(l) ;
844 }
845 diff_ent /= nzet ;
846
847 fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
848 fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
849 fichconv << " " << log10( fabs(max_triax) + 1.e-16 ) ;
850
851 vit_triax = 0 ;
852 if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
854 }
855
856 fichconv << " " << vit_triax ;
857
858 //------------------------------
859 // Recycling for the next step
860 //------------------------------
861
862 ent_prev = ent ;
863 logn_prev = logn ;
864 dzeta_prev = dzeta ;
866
867 fichconv << endl ;
868 fichfreq << endl ;
869 fichevol << endl ;
870 fichconv.flush() ;
871 fichfreq.flush() ;
872 fichevol.flush() ;
873
874 } // End of main loop
875
876 //=========================================================================
877 // End of iteration
878 //=========================================================================
879
884 for (int i=1; i<=3; i++) {
886 }
887
888 fichconv.close() ;
889 fichfreq.close() ;
890 fichevol.close() ;
891
892}
893}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Radial mapping of rather general form.
Definition map.h:2752
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
virtual void reevaluate(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.
Coord sint
Definition map.h:721
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
Definition map.h:807
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition map.C:331
virtual void poisson2d(const Cmp &source_mat, const Cmp &source_quad, Param &par, Cmp &uu) const =0
Computes the solution of a 2-D Poisson equation.
Coord phi
coordinate centered on the grid
Definition map.h:720
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321
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_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
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
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
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
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition star_rot.h:208
Sym_tensor tkij
Tensor related to the extrinsic curvature tensor by .
Definition star_rot.h:167
Scalar tggg
Metric potential .
Definition star_rot.h:137
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition star_rot.h:150
Scalar bbb
Metric factor B.
Definition star_rot.h:107
virtual double mass_b() const
Baryon mass.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
double omega
Rotation angular velocity ([f_unit] )
Definition star_rot.h:101
Scalar ak_car
Scalar .
Definition star_rot.h:186
Scalar uuu
Norm of u_euler.
Definition star_rot.h:121
Scalar nphi
Metric coefficient .
Definition star_rot.h:113
void update_metric()
Computes metric coefficients from known potentials.
virtual double mass_g() const
Gravitational mass.
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition star_rot.h:198
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Definition star_rot.h:126
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition star_rot.h:94
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition star_rot.h:216
Scalar dzeta
Metric potential .
Definition star_rot.h:134
Scalar a_car
Square of the metric factor A.
Definition star_rot.h:104
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition star_rot.h:192
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition star_rot.h:131
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition star_rot.h:225
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, Param *=0x0)
Computes an equilibrium configuration.
virtual double grv2() const
Error on the virial identity GRV2.
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Definition star_rot.C:590
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition star_rot.h:160
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition star_rot.C:748
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Scalar nn
Lapse function N .
Definition star.h:225
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
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].
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Scalar press
Fluid pressure.
Definition star.h:194
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
double ray_pole() const
Coordinate radius at [r_unit].
Basic array class.
Definition tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Values and coefficients of a (real-value) function.
Definition valeur.h:287
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 cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
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.