LORENE
star_bin_equilibrium_xcts.C
1/*
2 * Method of class Star_bin to compute an equilibrium configuration
3 * (see file star.h for documentation).
4 */
5
6/*
7 * Copyright (c) 2010 Michal Bejger
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License version 2
13 * as published by the Free Software Foundation.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26char star_bin_equilibrium_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $" ;
27
28/*
29 * $Id: star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $
30 * $Log: star_bin_equilibrium_xcts.C,v $
31 * Revision 1.13 2014/10/13 08:53:38 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.12 2014/10/06 15:13:16 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.11 2011/03/25 16:28:12 e_gourgoulhon
38 * Still in progress
39 *
40 * Revision 1.10 2010/12/20 15:42:10 m_bejger
41 * Various rearrangements of fields in Poissson equations
42 *
43 * Revision 1.9 2010/12/14 17:34:42 m_bejger
44 * Improved iteration for beta_auto poisson_vect()
45 *
46 * Revision 1.8 2010/12/09 10:48:06 m_bejger
47 * Testing the main equations
48 *
49 * Revision 1.7 2010/10/26 18:46:28 m_bejger
50 * Added table fact_resize for domain resizing
51 *
52 * Revision 1.6 2010/10/18 19:08:14 m_bejger
53 * Changed to allow for calculations with more than one domain in the star
54 *
55 * Revision 1.5 2010/06/23 20:40:56 m_bejger
56 * Corrections in equations for Psi_auto, chi_auto and beta_auto
57 *
58 * Revision 1.4 2010/06/18 13:28:59 m_bejger
59 * Adjusted the computation of the first integral, radial scale
60 *
61 * Revision 1.3 2010/06/17 17:05:06 m_bejger
62 * Testing version
63 *
64 * Revision 1.2 2010/06/15 08:21:21 m_bejger
65 * Minor changes; still not working properly
66 *
67 * Revision 1.1 2010/05/04 07:51:05 m_bejger
68 * Initial version
69 *
70 * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium_xcts.C,v 1.13 2014/10/13 08:53:38 j_novak Exp $
71 *
72 */
73
74// C headers
75#include <cmath>
76
77// Lorene headers
78#include "cmp.h"
79#include "tenseur.h"
80#include "metrique.h"
81#include "star.h"
82#include "param.h"
83#include "graphique.h"
84#include "utilitaires.h"
85#include "tensor.h"
86#include "nbr_spx.h"
87#include "unites.h"
88
89
90namespace Lorene {
92 int mermax,
93 int mermax_potvit,
95 double relax_poisson,
96 double relax_potvit,
97 double thres_adapt,
98 const Tbl& fact_resize,
99 const Tbl* pent_limit,
100 Tbl& diff) {
101
102 // Fundamental constants and units
103 // -------------------------------
104
105 using namespace Unites ;
106
107 // Initializations
108 // ---------------
109
110 const Mg3d* mg = mp.get_mg() ;
111 int nz = mg->get_nzone() ; // total number of domains
112
113 // The following is required to initialize mp_prev as a Map_et:
114 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
115
116 // Domain and radial indices of points at the surface of the star:
117 int l_b = nzet - 1 ;
118 int i_b = mg->get_nr(l_b) - 1 ;
119 int k_b ;
120 int j_b ;
121
122 // Value of the enthalpy defining the surface of the star
123 double ent_b = 0 ;
124
125 // Error indicators
126 // ----------------
127
128 double& diff_ent = diff.set(0) ;
129 double& diff_vel_pot = diff.set(1) ;
130 double& diff_psi = diff.set(2) ;
131 double& diff_chi = diff.set(3) ;
132 double& diff_beta_x = diff.set(4) ;
133 double& diff_beta_y = diff.set(5) ;
134 double& diff_beta_z = diff.set(6) ;
135
136 // Parameters for the function Map_et::adapt
137 // -----------------------------------------
138
140 int nitermax = 100 ;
141 int niter ;
142 int adapt_flag = 1 ; // 1 = performs the full computation,
143 // 0 = performs only the rescaling by
144 // the factor alpha_r
145 //## int nz_search = nzet + 1 ; // Number of domains for searching the
146 // enthalpy
147
148 int nz_search = nzet ; // Number of domains
149 // for searching
150 // the enthalpy isosurfaces
151
152 double precis_secant = 1.e-14 ;
153 double alpha_r ;
154 double reg_map = 1. ; // 1 = regular mapping,
155 // 0 = contracting mapping
156
157 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
158 // locate zeros by the secant method
159
160 par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
161 // to the isosurfaces of ent is to be performed
162
163 par_adapt.add_int(nz_search, 2) ; // number of domains to search
164 // the enthalpy isosurface
165
166 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
167 // 0 = performs only the rescaling by
168 // the factor alpha_r
169
170 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
171 // (theta_*, phi_*)
172
173 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
174 // (theta_*, phi_*)
175
176 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
177 // the secant method
178
179 par_adapt.add_double(precis_secant, 0) ;// required absolute precision in
180 // the determination of zeros by
181 // the secant method
182 par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
183
184 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
185 // distances will be multiplied
186
187 // Enthalpy values for the adaptation
189 if (pent_limit != 0x0) ent_limit = *pent_limit ;
190
191 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
192 // to define the isosurfaces.
193
194 double precis_poisson = 1.e-16 ;
195
198
199 // Parameters for the function Scalar::poisson for Psi_auto
200 // ---------------------------------------------------------------
201
202 Param par_psi ;
203
204 par_psi.add_int(mermax_poisson, 0) ; // maximum number of iterations
205 par_psi.add_double(relax_poisson, 0) ; // relaxation parameter
206 par_psi.add_double(precis_poisson, 1) ; // required precision
207 par_psi.add_int_mod(niter, 0) ; // number of iterations actually used
208 par_psi.add_cmp_mod( ssjm1psi ) ;
209
210 // Parameters for the function Scalar::poisson for chi_auto
211 // ---------------------------------------------------------------
212
213 Param par_chi ;
214
215 par_chi.add_int(mermax_poisson, 0) ; // maximum number of iterations
216 par_chi.add_double(relax_poisson, 0) ; // relaxation parameter
217 par_chi.add_double(precis_poisson, 1) ; // required precision
218 par_chi.add_int_mod(niter, 0) ; // number of iterations actually used
219 par_chi.add_cmp_mod( ssjm1chi ) ;
220
221 // Parameters for the function Vector::poisson for beta
222 // ----------------------------------------------------
223
225
226 par_beta.add_int(mermax_poisson, 0) ; // maximum number of
227 // iterations
228 par_beta.add_double(relax_poisson, 0) ; // relaxation parameter
229 par_beta.add_double(precis_poisson, 1) ; // required precision
230 par_beta.add_int_mod(niter, 0) ; // number of iterations
231 // actually used
232
233 // Sources at the previous step, for a poisson_vect() solver
235
237 ssjm1wbeta.set_etat_qcq() ;
238 for (int i=0; i<3; i++) ssjm1wbeta.set(i) = Cmp(ssjm1_wbeta(i+1)) ;
239
240 par_beta.add_cmp_mod(ssjm1khi) ;
241 par_beta.add_tenseur_mod(ssjm1wbeta) ;
242
243 // Redefinition of external potential
244 // See Eq (99) from Gourgoulhon et al. (2001)
245 // logN = logN_auto + logn_ac_rest = log(chi_auto + 1.)
246 // - log(Psi_auto + 1.) + logn_ac_rest
247 //------------------------------------
248
251
252 Scalar logn_auto = log(chi_auto_p1) - log(Psi_auto_p1) ;
253 logn_auto.std_spectral_base() ;
254
256 - log(1. + Psi_comp/Psi_auto_p1) ;
257
258 logn_ac_rest.std_spectral_base() ;
259
260 //cout << "logn_auto" << norme(logn_auto) << endl ;
261 //cout << "logn_ac_rest" << norme(logn_ac_rest) << endl ;
262 //cout << "pot_centri" << norme(pot_centri) << endl ;
263 //cout << "loggam" << norme(loggam) << endl ;
264
266 //cout << "pot_ext" << norme(pot_ext) << endl ;
267
268 Scalar ent_jm1 = ent ; // Enthalpy at previous step
269
270 Scalar source_tot(mp) ; // source term in equations Psi_auto
271 // and chi_auto
272
273 Vector source_beta(mp, CON, mp.get_bvect_cart()) ; // source term
274 // in the equation
275 // for beta_auto
276
277 //==================================================================
278 // Start of iteration
279 //==================================================================
280
281 for(int mer=0 ; mer<mermax ; mer++ ) {
282
283 cout << "-----------------------------------------------" << endl ;
284 cout << "step: " << mer << endl ;
285 cout << "diff_ent = " << diff_ent << endl ;
286
287 //-----------------------------------------------------
288 // Resolution of the elliptic equation for the velocity
289 // scalar potential
290 //-----------------------------------------------------
291
292 if (irrotational) {
293
295 relax_potvit) ;
296 }
297
298 diff_vel_pot = 0. ; // to avoid the warning
299
300 //-----------------------------------------------------
301 // Computation of the new radial scale
302 //--------------------------------------------------
303
304 // alpha_r (r = alpha_r r') is determined so that the enthalpy
305 // takes the requested value ent_b at the stellar surface
306
307 // Values at the center of the star:
308 double logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
309 double pot_ext_c = pot_ext.val_grid_point(0, 0, 0, 0) ;
310
311 // Search for the reference point (theta_*, phi_*) [notation of
312 // Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
313 // at the surface of the star
314 // ------------------------------------------------------------
315 double alpha_r2 = 0 ;
316 for (int k=0; k<mg->get_np(l_b); k++) {
317 for (int j=0; j<mg->get_nt(l_b); j++) {
318
319 double pot_ext_b = pot_ext.val_grid_point(l_b, k, j, i_b) ;
320 double logn_auto_b = logn_auto.val_grid_point(l_b, k, j, i_b) ;
321 // See Eq (100) from Gourgoulhon et al. (2001)
322 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
324
325 // cout << "k, j, alpha_r2_jk : " << k << " " << j << " " << alpha_r2_jk << endl ;
326
327 if (alpha_r2_jk > alpha_r2) {
329 k_b = k ;
330 j_b = j ;
331 }
332
333 }
334 }
335
337
338 cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
339 << alpha_r << endl ;
340
341 // New value of logn_auto
342 // ----------------------
343
344 Psi_auto = pow(Psi_auto +1.,alpha_r2) - 1. ;
345 chi_auto = pow(chi_auto +1.,alpha_r2) - 1. ;
348
349 //------------------------------------------------------------
350 // Change the values of the inner points of the domain adjascent
351 // to the surface of the star by those of the outer points of
352 // the domain under the surface
353 //------------------------------------------------------------
354
357
358 logn_auto = log(chi_auto + 1.) - log(Psi_auto + 1.) ;
359 logn_auto.std_spectral_base() ;
360
361 logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
362
363 //------------------------------------------
364 // First integral --> enthalpy in all space
365 // See Eq (98) from Gourgoulhon et al. (2001)
366 //-------------------------------------------
367
368 ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
369 cout.precision(8) ;
370 cout << "pot" << norme(pot_ext) << endl ;
371
373
374 //----------------------------------------------------
375 // Adaptation of the mapping to the new enthalpy field
376 //----------------------------------------------------
377
378 // Shall the adaptation be performed (cusp) ?
379 // ------------------------------------------
380
381 double dent_eq = ent.dsdr().val_point(ray_eq(),M_PI/2.,0.) ;
382 double dent_pole = ent.dsdr().val_point(ray_pole(),0.,0.) ;
383 double rap_dent = fabs( dent_eq / dent_pole ) ;
384 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
385
386 if ( rap_dent < thres_adapt ) {
387 adapt_flag = 0 ; // No adaptation of the mapping
388 cout << "******* FROZEN MAPPING *********" << endl ;
389 }
390 else{
391 adapt_flag = 1 ; // The adaptation of the mapping is to be
392 // performed
393 }
394
395 ent_limit.set_etat_qcq() ;
396 for (int l=0; l<nzet; l++) { // loop on domains inside the star
397 ent_limit.set(l) = ent.val_grid_point(l, k_b, j_b, i_b) ;
398
399 } ent_limit.set(nzet-1) = ent_b ;
400
402
403 Cmp ent_cmp(ent) ;
405 ent = ent_cmp ;
406
407 // Readjustment of the external boundary of domain l=nzet
408 // to keep a fixed ratio with respect to star's surface
409
410 double rr_in_1 = mp.val_r(nzet, -1., M_PI/2., 0.) ;
411
412 // Resizes the outer boundary of the shell including the comp. NS
413 double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
414
415 mp.resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
416
417 // Resizes the inner boundary of the shell including the comp. NS
418 double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
419
420 mp.resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
421
422 if (nz > nzet+3) {
423
424 // Resize of the domain from 1(nzet) to N-4
425 double rr_out_nm3_new = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
426
427 for (int i=nzet-1; i<nz-4; i++) {
428
429 double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
430
431 double rr_mid = rr_out_i
432 + (rr_out_nm3_new - rr_out_i) / double(nz - 3 - i) ;
433
434 double rr_2timesi = 2. * rr_out_i ;
435
436 if (rr_2timesi < rr_mid) {
437
438 double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
439
441
442 } else {
443
444 double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
445
446 mp.resize(i+1, rr_mid / rr_out_ip1) ;
447
448 } // End of else
449
450 } // End of i loop
451
452 } // End of (nz > nzet+3) loop
453
454// if (nz>= 5) {
455//
456// double separation = 2. * fabs(mp.get_ori_x()) ;
457// double ray_eqq = ray_eq() ;
458// double ray_eqq_pi = ray_eq_pi() ;
459// double new_rr_out_2 = (separation - ray_eqq) * 0.95 ;
460// double new_rr_out_3 = (separation + ray_eqq_pi) * 1.05 ;
461//
462// double rr_in_1 = mp.val_r(1,-1., M_PI/2, 0.) ;
463// double rr_out_1 = mp.val_r(1, 1., M_PI/2, 0.) ;
464// double rr_out_2 = mp.val_r(2, 1., M_PI/2, 0.) ;
465// double rr_out_3 = mp.val_r(3, 1., M_PI/2, 0.) ;
466//
467// mp.resize(1, 0.5*(new_rr_out_2 + rr_in_1) / rr_out_1) ;
468// mp.resize(2, new_rr_out_2 / rr_out_2) ;
469// mp.resize(3, new_rr_out_3 / rr_out_3) ;
470//
471// for (int dd=4; dd<=nz-2; dd++) {
472// mp.resize(dd, new_rr_out_3 * pow(4., dd-3) /
473// mp.val_r(dd, 1., M_PI/2, 0.)) ;
474// }
475//
476// }
477
478 //else cout << "too small number of domains" << endl ;
479
480
481 //------------------------------------------------------------------
482 // Computation of the enthalpy at the new grid points
483 //------------------------------------------------------------------
484
485 mp_prev.homothetie(alpha_r) ;
486
487 //Cmp ent_cmp2 (ent) ;
488 //mp.reevaluate_symy(&mp_prev, nzet+1, ent_cmp2) ;
489 //ent = ent_cmp2 ;
490
492
493 double ent_s_max = -1 ;
494 int k_s_max = -1 ;
495 int j_s_max = -1 ;
496 for (int k=0; k<mg->get_np(l_b); k++) {
497 for (int j=0; j<mg->get_nt(l_b); j++) {
498 double xx = fabs( ent.val_grid_point(l_b, k, j, i_b) ) ;
499 if (xx > ent_s_max) {
500 ent_s_max = xx ;
501 k_s_max = k ;
502 j_s_max = j ;
503 }
504 }
505 }
506 cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
507 << " and nzet : " << endl ;
508 cout << " " << ent_s_max << " reached for k = " << k_s_max <<
509 " and j = " << j_s_max << endl ;
510
511 //------------------------------------------------------------------
512 // Equation of state
513 //------------------------------------------------------------------
514
515 equation_of_state() ; // computes new values for nbar (n), ener (e)
516 // and press (p) from the new ent (H)
517
518 //------------------------------------------------------------------
519 // Matter source terms in the gravitational field equations
520 //------------------------------------------------------------------
521
522 hydro_euler() ; // computes new values for ener_euler (E),
523 // s_euler (S) and u_euler (U^i)
524
525 // -------------------------------
526 // AUXILIARY QUANTITIES
527 // -------------------------------
528
529 int nr = mp.get_mg()->get_nr(0) ;
530 int nt = mp.get_mg()->get_nt(0) ;
531 int np = mp.get_mg()->get_np(0) ;
532
533 Scalar Psi3 = psi4 / Psi ;
534 Psi3.std_spectral_base() ;
535
536 Scalar sPsi7 = Psi * pow(psi4, -2.) ;
537 sPsi7.std_spectral_base() ;
538
539 //------------------------------------------------------------------
540 // Poisson equation for Psi_auto (Eq. 8.127 of arXiv:gr-qc/0703035)
541 //------------------------------------------------------------------
542
543 // Source
544 //--------
545
546 source_tot = - 0.5 * qpig * psi4 % Psi % ener_euler
547 - 0.125 * sPsi7 * (hacar_auto + hacar_comp) ;
548 source_tot.std_spectral_base() ;
549
550 // Resolution of the Poisson equation
551 // ----------------------------------
552
553 cout << "Resolution of the Poisson equation for Psi_auto:" << endl ;
554 source_tot.poisson(par_psi, Psi_auto) ;
556
557 cout << "Psi_auto: " << endl << norme(Psi_auto/(nr*nt*np)) << endl ;
558
559 // Check: has the Poisson equation been correctly solved ?
560 // -----------------------------------------------------
561
563 cout <<
564 "Relative error in the resolution of the equation for Psi_auto : "
565 << endl ;
566 for (int l=0; l<nz; l++) {
567 cout << tdiff_psi(l) << " " ;
568 }
569 cout << endl
570 << "==========================================================="
571 << endl ;
572
574
575 //------------------------------------------------------------------
576 // Poisson equation for chi_auto (Eq. 8.129 of arXiv:gr-qc/0703035)
577 //------------------------------------------------------------------
578
579 // Source
580 //--------
581
582 source_tot = chi * 0.5 * qpig * psi4 % (ener_euler + 2.*s_euler)
583 + 0.875 * nn % sPsi7 * (hacar_auto + hacar_comp) ;
584 source_tot.std_spectral_base() ;
585
586 // Resolution of the Poisson equation
587 // ----------------------------------
588
589 cout << "Resolution of the Poisson equation for chi_auto:" << endl ;
590 source_tot.poisson(par_chi, chi_auto) ;
592
593 cout << "chi_auto: " << endl << norme(chi_auto/(nr*nt*np)) << endl ;
594
595 // Check: has the Poisson equation been correctly solved ?
596 // -----------------------------------------------------
597
599 cout <<
600 "Relative error in the resolution of the equation for chi_auto : "
601 << endl ;
602
603 for (int l=0; l<nz; l++) {
604 cout << tdiff_chi(l) << " " ;
605 }
606 cout << endl
607 << "==========================================================="
608 << endl ;
609
611
612 //------------------------------------------------------------------
613 // Vector Poisson equation for beta_auto
614 // (Eq. 8.128 of arXiv:gr-qc/0703035)
615 //------------------------------------------------------------------
616
617 // Source
618 //--------
619
620 source_beta = 4.* qpig * chi % Psi3 * (ener_euler + press) * u_euler
621 + 2. * sPsi7 * contract(haij_auto, 1, dcov_chi, 0)
622 -14. * nn % sPsi7 * contract(haij_auto, 1, dcov_Psi, 0) ;
623 source_beta.std_spectral_base() ;
624
625 // Resolution of the Poisson equation
626 // ----------------------------------
627
628 Tenseur source_p(mp, 1, CON, mp.get_bvect_cart() ) ;
629 source_p.set_etat_qcq() ;
630 for (int i=0; i<3; i++) {
631
632 source_p.set(i) = Cmp(source_beta(i+1)) ;
633 //source_p.set(i).filtre(4) ;
634
635 }
636
637 Tenseur vect_auxi (mp, 1, CON, mp.get_bvect_cart()) ;
638 vect_auxi.set_etat_qcq() ;
639 for (int i=0; i<3 ;i++) vect_auxi.set(i) = Cmp(w_beta(i+1)) ;
640 vect_auxi.set_std_base() ;
641
642
644 scal_auxi.set_etat_qcq() ;
645 scal_auxi.set() = Cmp(khi) ;
646 scal_auxi.set_std_base() ;
647
648 Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart()) ;
649 resu_p.set_etat_qcq() ;
650 for (int i=0; i<3 ;i++) resu_p.set(i) = Cmp(beta_auto.set(i+1));
651 resu_p.set_std_base() ;
652
653 // Resolution of the vector Poisson equation
654 // -----------------------------------------
655
656 double lambda = double(1) / double(3) ;
658 //source_p.poisson_vect_oohara(lambda, par_beta, resu_p, scal_auxi) ;
659
660 for (int i=1; i<=3; i++) {
661 beta_auto.set(i) = resu_p(i-1) ;
662 w_beta.set(i) = vect_auxi(i-1) ;
663 }
664 khi = scal_auxi() ;
665
666 // Values of sources from the previous step
668 for (int i=0; i<3; i++) ssjm1_wbeta.set(i+1) = ssjm1wbeta(i) ;
669
670 cout << "beta_auto_x" << endl << norme(beta_auto(1)/(nr*nt*np))
671 << endl ;
672 cout << "beta_auto_y" << endl << norme(beta_auto(2)/(nr*nt*np))
673 << endl ;
674 cout << "beta_auto_z" << endl << norme(beta_auto(3)/(nr*nt*np))
675 << endl ;
676
677 // Check: has the equation for beta_auto been correctly solved ?
678 // --------------------------------------------------------------
679
681 + lambda* beta_auto.divergence(flat).derive_con(flat) ;
682
683 source_beta.dec_dzpuis() ;
684 Tbl tdiff_beta_x = diffrel(lap_beta(1), source_beta(1)) ;
685 Tbl tdiff_beta_y = diffrel(lap_beta(2), source_beta(2)) ;
686 Tbl tdiff_beta_z = diffrel(lap_beta(3), source_beta(3)) ;
687
688 cout <<
689 "Relative error in the resolution of the equation for beta_auto : "
690 << endl ;
691 cout << "x component : " ;
692 for (int l=0; l<nz; l++) {
693 cout << tdiff_beta_x(l) << " " ;
694 }
695 cout << endl ;
696 cout << "y component : " ;
697 for (int l=0; l<nz; l++) {
698 cout << tdiff_beta_y(l) << " " ;
699 }
700 cout << endl ;
701 cout << "z component : " ;
702 for (int l=0; l<nz; l++) {
703 cout << tdiff_beta_z(l) << " " ;
704 }
705 cout << endl ;
706
710
711 // End of relativistic equations
712
713 //-------------------------------------------------
714 // Relative change in enthalpy
715 //-------------------------------------------------
716
719 for (int l=1; l<nzet; l++) diff_ent += diff_ent_tbl(l) ;
720
721 diff_ent /= nzet ;
722
723
724 ent_jm1 = ent ;
725
726
727 } // End of main loop
728
729 //==================================================================
730 // End of iteration
731 //==================================================================
732
733}
734
735
736
737
738
739
740
741}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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 resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
virtual void reevaluate_symy(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.
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
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 Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian 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
Scalar Psi
Total conformal factor .
Definition star.h:1152
Scalar pot_centri
Centrifugal potential.
Definition star.h:1129
Sym_tensor haij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition star.h:1193
Scalar psi4
Conformal factor .
Definition star.h:1158
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition star.h:1177
void equilibrium(double ent_c, int mermax, int mermax_potvit, int mermax_poisson, double relax_poisson, double relax_potvit, double thres_adapt, const Tbl &fact, const Tbl *pent_limit, Tbl &diff)
Computes an equilibrium configuration.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi.
Definition star.h:1227
Scalar khi
Solution for the scalar part of the vector Poisson equation for .
Definition star.h:1168
Vector ssjm1_wbeta
Effective source at the previous step for wbeta of the vector Poisson equation for wbeta (needed for ...
Definition star.h:1234
Scalar Psi_comp
Scalar field generated principally by the companion star.
Definition star.h:1149
Scalar ssjm1_chi
Effective source at the previous step for the resolution of the Poisson equation for \chi_auto .
Definition star.h:1216
Vector dcov_Psi
Covariant derivative of the conformal factor .
Definition star.h:1171
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition star.h:1182
Scalar Psi_auto
Scalar field generated principally by the star.
Definition star.h:1144
Scalar hacar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Definition star.h:1211
Scalar chi
Total function .
Definition star.h:1155
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star.h:1120
Vector dcov_chi
Covariant derivative of the function .
Definition star.h:1174
Scalar chi_auto
Scalar field generated principally by the star.
Definition star.h:1134
Vector w_beta
Solution for the vector part of the vector Poisson equation for .
Definition star.h:1163
Scalar chi_comp
Scalar field generated principally by the companion star.
Definition star.h:1139
Scalar ssjm1_psi
Effective source at the previous step for the resolution of the Poisson equation for \Psi_auto .
Definition star.h:1221
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
bool irrotational
true for an irrotational star, false for a corotating one
Definition star.h:1099
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Scalar hacar_auto
Part of the scalar generated by beta_auto, i.e.
Definition star.h:1205
Scalar nn
Lapse function N .
Definition star.h:225
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
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
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
void smooth(int nzet, Valeur &uuva) const
Changes the function *this as a smooth one when there exists a discontinuity between the nucleus and ...
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
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 dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
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.