LORENE
et_magnetisation.C
1 /*
2  * Methods for axisymmetric rotating neutron stars with magnetisation in the stress-energy tensor.
3  *
4  * See the file et_rot_mag.h for documentation
5  *
6  */
7 
8 /*
9  * Copyright (c) 2013-2014 Jerome Novak, Debarti Chatterjee, Micaela Oertel
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 char et_magnetisation_C[] = "$Header $" ;
30 
31 /*
32  * $Id: et_magnetisation.C,v 1.10 2014/10/21 09:23:53 j_novak Exp $
33  * $Log: et_magnetisation.C,v $
34  * Revision 1.10 2014/10/21 09:23:53 j_novak
35  * Addition of global functions mass_g(), angu_mom(), grv2/3() and mom_quad().
36  *
37  * Revision 1.9 2014/10/13 08:52:56 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.8 2014/05/28 07:46:06 j_novak
41  * Minor modifications.
42  *
43  * Revision 1.7 2014/05/27 12:32:28 j_novak
44  * Added possibility to converge to a given magnetic moment.
45  *
46  * Revision 1.6 2014/05/14 15:19:05 j_novak
47  * The magnetisation field is now filtered.
48  *
49  * Revision 1.5 2014/04/29 13:46:07 j_novak
50  * Addition of switches 'use_B_in_eos' and 'include_magnetisation' to control the model.
51  *
52  * Revision 1.4 2014/04/28 12:48:13 j_novak
53  * Minor modifications.
54  *
55  * Revision 1.3 2013/12/19 17:05:40 j_novak
56  * Corrected a dzpuis problem.
57  *
58  * Revision 1.2 2013/12/13 16:36:51 j_novak
59  * Addition and computation of magnetisation terms in the Einstein equations.
60  *
61  * Revision 1.1 2013/11/25 13:52:11 j_novak
62  * New class Et_magnetisation to include magnetization terms in the stress energy tensor.
63  *
64  *
65  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_magnetisation.C,v 1.10 2014/10/21 09:23:53 j_novak Exp $
66  *
67  */
68 
69 // Headers C
70 #include <cmath>
71 
72 // Headers Lorene
73 #include "et_rot_mag.h"
74 #include "utilitaires.h"
75 #include "unites.h"
76 #include "param.h"
77 #include "eos.h"
78 
79  //--------------//
80  // Constructors //
81  //--------------//
82 // Standard constructor
83 // --------------------
84 
85 
86 namespace Lorene {
87 Et_magnetisation::Et_magnetisation(Map& mp_i, int nzet_i, bool relat,
88  const Eos& eos_i, bool magnetisation,
89  bool B_eos)
90  : Et_rot_mag(mp_i, nzet_i, relat, eos_i, 1),
91  use_B_in_eos(B_eos), include_magnetisation(magnetisation),
92  xmag(mp_i),
93  E_I(mp_i),
94  J_I(mp_i, COV, mp_i.get_bvect_spher()),
95  Sij_I(mp_i, COV, mp_i.get_bvect_spher())
96 {
97  const Eos_mag* tmp_p_eos = dynamic_cast<const Eos_mag*>(&eos_i) ;
98  if (tmp_p_eos == 0x0) {
99  cerr << "Et_magnetisation::Et_magnetisation : " << endl ;
100  cerr << "Only magnetised EoS is admitted for this class!" << endl ;
101  cerr << "Aborting ... " << endl ;
102  abort() ;
103  }
104  if ( include_magnetisation && (!use_B_in_eos) ) {
105  cerr << "Et_magnetisation::Et_magnetisation : " << endl ;
106  cerr << "Magnetisation terms can be included only if "
107  << "the magnetic field is used in the EoS!" << endl;
108  cerr << "Aborting ... " << endl ;
109  abort() ;
110  }
111  xmag = 0 ;
112  E_I = 0 ;
113  J_I.set_etat_zero() ;
114  Sij_I.set_etat_zero() ;
115 
116 }
117 
118 
119 
120 Et_magnetisation::Et_magnetisation(Map& mp_i, const Eos& eos_i, FILE* fich)
121  : Et_rot_mag(mp_i, eos_i, fich, 0),
122  use_B_in_eos(true), include_magnetisation(true),
123  xmag(mp_i),
124  E_I(mp_i),
125  J_I(mp_i, COV, mp_i.get_bvect_spher()),
126  Sij_I(mp_i, COV, mp_i.get_bvect_spher())
127 {
128 
129  // Read of the saved fields:
130  // ------------------------
131 
132  fread(&use_B_in_eos, sizeof(bool), 1, fich) ;
133 
134  fread(&include_magnetisation, sizeof(bool), 1, fich) ;
135 
136  Scalar xmag_file(mp, *mp.get_mg(), fich) ;
137  xmag = xmag_file ;
138 
139  Scalar E_I_file(mp, *mp.get_mg(), fich) ;
140  E_I = E_I_file ;
141 
142  Vector J_I_file(mp, mp.get_bvect_spher(), fich) ;
143  J_I = J_I_file ;
144 
145  Sym_tensor Sij_I_file(mp, mp.get_bvect_spher(), fich) ;
146  Sij_I = Sij_I_file ;
147 
148 }
149 
150 
151 // Copy constructor
152 // ----------------
153 
155  : Et_rot_mag(et),
157  xmag(et.xmag),
158  E_I(et.E_I),
159  J_I(et.J_I),
160  Sij_I(et.Sij_I)
161 { }
162 
163 
164  //------------//
165  // Destructor //
166  //------------//
167 
169 
170 
171 // Assignment to another Et_magnetisation
172 // --------------------------------
173 
175 
176  // Assignement of quantities common to all the derived classes of Et_rot_mag
180  xmag = et.xmag ;
181  E_I = et.E_I ;
182  J_I = et.J_I ;
183  Sij_I = et.Sij_I ;
184 
185 }
186 
187 
189 
190  Cmp ent_eos = ent() ;
191  const Eos_mag* mageos = dynamic_cast<const Eos_mag*>(&eos) ;
192  assert (mageos != 0x0) ;
193 
194  // Slight rescale of the enthalpy field in case of 2 domains inside the
195  // star
196 
197  double epsilon = 1.e-12 ;
198 
199  const Mg3d* mg = mp.get_mg() ;
200  int nz = mg->get_nzone() ;
201 
202  Mtbl xi(mg) ;
203  xi.set_etat_qcq() ;
204  for (int l=0; l<nz; l++) {
205  xi.t[l]->set_etat_qcq() ;
206  for (int k=0; k<mg->get_np(l); k++) {
207  for (int j=0; j<mg->get_nt(l); j++) {
208  for (int i=0; i<mg->get_nr(l); i++) {
209  xi.set(l,k,j,i) = mg->get_grille3d(l)->x[i] ;
210  }
211  }
212  }
213 
214  }
215 
216  Cmp fact_ent(mp) ;
217  fact_ent.allocate_all() ;
218 
219  fact_ent.set(0) = 1 + epsilon * xi(0) * xi(0) ;
220  fact_ent.set(1) = 1 - 0.25 * epsilon * (xi(1) - 1) * (xi(1) - 1) ;
221 
222  for (int l=nzet; l<nz; l++) {
223  fact_ent.set(l) = 1 ;
224  }
225 
226  if (nzet > 1) {
227 
228  if (nzet == 3) {
229  fact_ent.set(1) = 1 - 0.5 * epsilon * (xi(1) - 0.5) * (xi(1) - 0.5) ;
230  fact_ent.set(2) = 1 - 0.25 * epsilon * (xi(2) - 1) * (xi(2) - 1) ;
231  }
232 
233  if (nzet > 3) {
234  cerr << "Et_magnetisation::equation_of_state: "
235  << "not ready yet for nzet > 3 !" << endl ;
236  }
237 
238  ent_eos = fact_ent * ent_eos ;
239  ent_eos.std_base_scal() ;
240  }
241 
242 
243 
244  // Call to a magnetized EOS
245  // with the norm of the magnetic field passed as an argument to the EoS.
246 
247  double magb0 = 0 ;
248  Cmp norm_b(mp) ;
249  norm_b.set_etat_zero() ;
250  if (use_B_in_eos) {
251  Tenseur Bmag = Magn() ;
252  norm_b = sqrt(a_car() * ( Bmag(0)*Bmag(0) + Bmag(1)*Bmag(1) ))
253  / gam_euler() ; // Use of b^\mu in the fluid frame
254  norm_b.std_base_scal() ;
255  }
256  Param par ;
257  par.add_double_mod(magb0) ;
258 
265  xmag.allocate_all() ;
266 
267  for (int l=0; l< nz; l++) {
268  Tbl* tent = ent_eos.va.c->t[l] ;
269  if ( (tent->get_etat() == ETATZERO) || (l >= nzet) ) {
270  nbar.set().set(l).set_etat_zero() ;
271  ener.set().set(l).set_etat_zero() ;
272  press.set().set(l).set_etat_zero() ;
273  xmag.annule_domain(l) ;
274  }
275  else {
276  nbar.set().set(l).annule_hard() ;
277  ener.set().set(l).annule_hard() ;
278  press.set().set(l).annule_hard() ;
280  for (int k=0; k < mg->get_np(l); k++) {
281  for (int j=0; j < mg->get_nt(l); j++) {
282  for (int i=0; i < mg->get_nr(l); i++) {
283  magb0 = norm_b(l, k, j, i) ;
284  double ent0 = ent_eos(l, k, j, i) ;
285  nbar.set().set(l, k, j, i) = mageos->nbar_ent_p(ent0, &par) ;
286  ener.set().set(l, k, j, i) = mageos->ener_ent_p(ent0, &par) ;
287  press.set().set(l, k, j, i) = mageos->press_ent_p(ent0, &par) ;
288  xmag.set_grid_point(l, k, j, i) = mageos->mag_ent_p(ent0, &par) ;
289  }
290  }
291  }
292  }
293  }
294 
296  xmag.set_etat_zero() ;
297 
298  // Set the bases for spectral expansion
299  nbar.set_std_base() ;
300  ener.set_std_base() ;
301  press.set_std_base() ;
303 
304  Scalar tmp_scal(xmag) ;
305  tmp_scal.exponential_filter_r(0, 0, 1) ;
306  xmag = tmp_scal ;
307 
308  // The derived quantities are obsolete
309  del_deriv() ;
310 
311 }
312 
313 
314 
315  //--------------//
316  // Outputs //
317  //--------------//
318 
319 // Save in a file
320 // --------------
321 void Et_magnetisation::sauve(FILE* fich) const {
322 
323  Et_rot_mag::sauve(fich) ;
324 
325  fwrite(&use_B_in_eos, sizeof(bool), 1, fich) ;
326 
327  fwrite(&include_magnetisation, sizeof(bool), 1, fich) ;
328 
329  xmag.sauve(fich) ;
330  E_I.sauve(fich) ;
331  J_I.sauve(fich) ;
332  Sij_I.sauve(fich) ;
333 
334 }
335 
336 
337 // Printing
338 // --------
339 
340 
341 ostream& Et_magnetisation::operator>>(ostream& ost) const {
342 
343  using namespace Unites_mag ;
344 
346  // int theta_eq = mp.get_mg()->get_nt(nzet-1)-1 ;
347  ost << endl ;
348  ost << "Rotating magnetized neutron star"
349  << endl ;
350  if (use_B_in_eos)
351  ost << "Using magnetic field in the EoS" << endl ;
352  if (include_magnetisation) {
353  ost << "Including magnetisation terms in the equations" << endl ;
354  ost << "Maximal value of the magnetization scalar x : "
355  << max(maxabs(xmag)) << endl ;
356  }
357  return ost ;
358 }
359 
360 //=================================================================================
361 //
362 // Computation of equilibrium configuration
363 //
364 //=================================================================================
365 
366 
367 void Et_magnetisation::equilibrium_mag(double ent_c, double omega0,
368  double fact_omega, int nzadapt, const Tbl& ent_limit,
369  const Itbl& icontrol, const Tbl& control, double mbar_wanted,
370  double magmom_wanted,
371  double aexp_mass, Tbl& diff, double Q0, double a_j0,
372  Cmp (*f_j)(const Cmp&, const double),
373  Cmp (*M_j)(const Cmp& x, const double)) {
374 
375  // Fundamental constants and units
376  // -------------------------------
377  using namespace Unites_mag ;
378 
379  // For the display
380  // ---------------
381  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
382  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
383 
384  // Grid parameters
385  // ---------------
386 
387  const Mg3d* mg = mp.get_mg() ;
388  int nz = mg->get_nzone() ; // total number of domains
389  int nzm1 = nz - 1 ;
390 
391  // The following is required to initialize mp_prev as a Map_et:
392  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
393 
394  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
395  assert(mg->get_type_t() == SYM) ;
396  int l_b = nzet - 1 ;
397  int i_b = mg->get_nr(l_b) - 1 ;
398  int j_b = mg->get_nt(l_b) - 1 ;
399  int k_b = 0 ;
400 
401  // Value of the enthalpy defining the surface of the star
402  double ent_b = ent_limit(nzet-1) ;
403 
404  // Parameters to control the iteration
405  // -----------------------------------
406 
407  int mer_max = icontrol(0) ;
408  int mer_rot = icontrol(1) ;
409  int mer_change_omega = icontrol(2) ;
410  int mer_fix_omega = icontrol(3) ;
411  int mer_mass = icontrol(4) ;
412  int mermax_poisson = icontrol(5) ;
413  int delta_mer_kep = icontrol(6) ;
414  int mer_mag = icontrol(7) ;
415  int mer_change_mag = icontrol(8) ;
416  int mer_fix_mag = icontrol(9) ;
417  int mer_magmom = icontrol(10) ;
418 
419  // Protections:
420  if (mer_change_omega < mer_rot) {
421  cerr << "Et_magnetisation::equilibrium_mag: mer_change_omega < mer_rot !"
422  << endl ;
423  cerr << " mer_change_omega = " << mer_change_omega << endl ;
424  cerr << " mer_rot = " << mer_rot << endl ;
425  abort() ;
426  }
427  if (mer_fix_omega < mer_change_omega) {
428  cerr << "Et_magnetisation::equilibrium_mag: mer_fix_omega < mer_change_omega !"
429  << endl ;
430  cerr << " mer_fix_omega = " << mer_fix_omega << endl ;
431  cerr << " mer_change_omega = " << mer_change_omega << endl ;
432  abort() ;
433  }
434 
435  // In order to converge to a given baryon mass, shall the central
436  // enthalpy be varied or Omega ?
437  bool change_ent = true ;
438  if (mer_mass < 0) {
439  change_ent = false ;
440  mer_mass = abs(mer_mass) ;
441  }
442 
443  double precis = control(0) ;
444  double omega_ini = control(1) ;
445  double relax = control(2) ;
446  double relax_prev = double(1) - relax ;
447  double relax_poisson = control(3) ;
448  double thres_adapt = control(4) ;
449  double precis_adapt = control(5) ;
450  double Q_ini = control(6) ;
451  double a_j_ini = control (7) ;
452 
453  // Error indicators
454  // ----------------
455 
456  diff.set_etat_qcq() ;
457  double& diff_ent = diff.set(0) ;
458 
459  // Parameters for the function Map_et::adapt
460  // -----------------------------------------
461 
462  Param par_adapt ;
463  int nitermax = 100 ;
464  int niter ;
465  int adapt_flag = 1 ; // 1 = performs the full computation,
466  // 0 = performs only the rescaling by
467  // the factor alpha_r
468  int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
469  // isosurfaces
470  double alpha_r ;
471  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
472 
473  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
474  // locate zeros by the secant method
475  par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
476  // to the isosurfaces of ent is to be
477  // performed
478  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
479  // the enthalpy isosurface
480  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
481  // 0 = performs only the rescaling by
482  // the factor alpha_r
483  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
484  // (theta_*, phi_*)
485  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
486  // (theta_*, phi_*)
487 
488  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
489  // the secant method
490 
491  par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
492  // the determination of zeros by
493  // the secant method
494  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping,
495  // 0 = contracting mapping
496 
497  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
498  // distances will be multiplied
499 
500  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
501  // to define the isosurfaces.
502 
503 
504  // Parameters for the function Map_et::poisson for nuf
505  // ----------------------------------------------------
506 
507  double precis_poisson = 1.e-16 ;
508 
509  Param par_poisson_nuf ;
510  par_poisson_nuf.add_int(mermax_poisson, 0) ; // maximum number of iterations
511  par_poisson_nuf.add_double(relax_poisson, 0) ; // relaxation parameter
512  par_poisson_nuf.add_double(precis_poisson, 1) ; // required precision
513  par_poisson_nuf.add_int_mod(niter, 0) ; // number of iterations actually used
514  par_poisson_nuf.add_cmp_mod( ssjm1_nuf ) ;
515 
516  Param par_poisson_nuq ;
517  par_poisson_nuq.add_int(mermax_poisson, 0) ; // maximum number of iterations
518  par_poisson_nuq.add_double(relax_poisson, 0) ; // relaxation parameter
519  par_poisson_nuq.add_double(precis_poisson, 1) ; // required precision
520  par_poisson_nuq.add_int_mod(niter, 0) ; // number of iterations actually used
521  par_poisson_nuq.add_cmp_mod( ssjm1_nuq ) ;
522 
523  Param par_poisson_tggg ;
524  par_poisson_tggg.add_int(mermax_poisson, 0) ; // maximum number of iterations
525  par_poisson_tggg.add_double(relax_poisson, 0) ; // relaxation parameter
526  par_poisson_tggg.add_double(precis_poisson, 1) ; // required precision
527  par_poisson_tggg.add_int_mod(niter, 0) ; // number of iterations actually used
528  par_poisson_tggg.add_cmp_mod( ssjm1_tggg ) ;
529  double lambda_tggg ;
530  par_poisson_tggg.add_double_mod( lambda_tggg ) ;
531 
532  Param par_poisson_dzeta ;
533  double lbda_grv2 ;
534  par_poisson_dzeta.add_double_mod( lbda_grv2 ) ;
535 
536  // Parameters for the function Tenseur::poisson_vect
537  // -------------------------------------------------
538 
539  Param par_poisson_vect ;
540 
541  par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of iterations
542  par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
543  par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
544  par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
545  par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
546  par_poisson_vect.add_int_mod(niter, 0) ;
547 
548  // Parameters for the Maxwell equations
549  // -------------------------------------
550 
551  Param par_poisson_At ; // For scalar At Poisson equation
552  Cmp ssjm1_At(mp) ;
553  ssjm1_At.set_etat_zero() ;
554  par_poisson_At.add_int(mermax_poisson, 0) ; // maximum number of iterations
555  par_poisson_At.add_double(relax_poisson, 0) ; // relaxation parameter
556  par_poisson_At.add_double(precis_poisson, 1) ; // required precision
557  par_poisson_At.add_int_mod(niter, 0) ; // number of iterations actually used
558  par_poisson_At.add_cmp_mod( ssjm1_At ) ;
559 
560  Param par_poisson_Avect ; // For vector Aphi Poisson equation
561 
562  Cmp ssjm1_khi_mag(ssjm1_khi) ;
563  Tenseur ssjm1_w_mag(ssjm1_wshift) ;
564 
565  par_poisson_Avect.add_int(mermax_poisson, 0) ; // maximum number of iterations
566  par_poisson_Avect.add_double(relax_poisson, 0) ; // relaxation parameter
567  par_poisson_Avect.add_double(precis_poisson, 1) ; // required precision
568  par_poisson_Avect.add_cmp_mod( ssjm1_khi_mag ) ;
569  par_poisson_Avect.add_tenseur_mod( ssjm1_w_mag ) ;
570  par_poisson_Avect.add_int_mod(niter, 0) ;
571 
572  // Initializations
573  // ---------------
574 
575  // Initial angular velocity / magnetic quantities
576  omega = 0 ;
577  Q = 0 ;
578  a_j = 0 ;
579 
580  double accrois_omega = (omega0 - omega_ini)
581  / double(mer_fix_omega - mer_change_omega) ;
582  double accrois_Q = (Q0 - Q_ini)
583  / double(mer_fix_mag - mer_change_mag);
584  double accrois_a_j = (a_j0 - a_j_ini)
585  / double(mer_fix_mag - mer_change_mag);
586 
587  update_metric() ; // update of the metric coefficients
588 
589  equation_of_state() ; // update of the density, pressure, etc...
590 
591  hydro_euler() ; // update of the hydro quantities relative to the
592  // Eulerian observer
593 
594  MHD_comput() ; // update of EM contributions to stress-energy tensor
595 
596 
597  // Quantities at the previous step :
598  Map_et mp_prev = mp_et ;
599  Tenseur ent_prev = ent ;
600  Tenseur logn_prev = logn ;
601  Tenseur dzeta_prev = dzeta ;
602 
603  // Creation of uninitialized tensors:
604  Tenseur source_nuf(mp) ; // source term in the equation for nuf
605  Tenseur source_nuq(mp) ; // source term in the equation for nuq
606  Tenseur source_dzf(mp) ; // matter source term in the eq. for dzeta
607  Tenseur source_dzq(mp) ; // quadratic source term in the eq. for dzeta
608  Tenseur source_tggg(mp) ; // source term in the eq. for tggg
609  Tenseur source_shift(mp, 1, CON, mp.get_bvect_cart()) ; // source term for shift
610  Tenseur mlngamma(mp) ; // centrifugal potential
611 
612  // Preparations for the Poisson equations:
613  // --------------------------------------
614  if (nuf.get_etat() == ETATZERO) {
615  nuf.set_etat_qcq() ;
616  nuf.set() = 0 ;
617  }
618 
619  if (relativistic) {
620  if (nuq.get_etat() == ETATZERO) {
621  nuq.set_etat_qcq() ;
622  nuq.set() = 0 ;
623  }
624 
625  if (tggg.get_etat() == ETATZERO) {
626  tggg.set_etat_qcq() ;
627  tggg.set() = 0 ;
628  }
629 
630  if (dzeta.get_etat() == ETATZERO) {
631  dzeta.set_etat_qcq() ;
632  dzeta.set() = 0 ;
633  }
634  }
635 
636  ofstream fichconv("convergence.d") ; // Output file for diff_ent
637  fichconv << "# diff_ent GRV2 " << endl ;
638 
639  ofstream fichfreq("frequency.d") ; // Output file for omega
640  fichfreq << "# f [Hz]" << endl ;
641 
642  ofstream fichevol("evolution.d") ; // Output file for various quantities
643  fichevol <<
644  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
645  << endl ;
646 
647  diff_ent = 1 ;
648  double err_grv2 = 1 ;
649 
650  //=========================================================================
651  // Start of iteration
652  //=========================================================================
653 
654  for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
655 
656  cout << "-----------------------------------------------" << endl ;
657  cout << "step: " << mer << endl ;
658  cout << "diff_ent = " << display_bold << diff_ent << display_normal
659  << endl ;
660  cout << "err_grv2 = " << err_grv2 << endl ;
661  fichconv << mer ;
662  fichfreq << mer ;
663  fichevol << mer ;
664 
665  if (mer >= mer_rot) {
666 
667  if (mer < mer_change_omega) {
668  omega = omega_ini ;
669  }
670  else {
671  if (mer <= mer_fix_omega) {
672  omega = omega_ini + accrois_omega *
673  (mer - mer_change_omega) ;
674  }
675  }
676 
677  }
678 
679  if (mer >= mer_mag) {
680  if (mer < mer_change_mag) {
681  Q = Q_ini ;
682  a_j = a_j_ini ;
683  }
684  else {
685  if (mer <= mer_fix_mag) {
686  Q = Q_ini + accrois_Q * (mer - mer_change_mag) ;
687  a_j = a_j_ini + accrois_a_j * (mer - mer_change_mag) ;
688  }
689  }
690  }
691 
692  //-----------------------------------------------
693  // Computation of electromagnetic potentials :
694  // -------------------------------------------
695  magnet_comput(adapt_flag, f_j, par_poisson_At, par_poisson_Avect) ;
696 
697  MHD_comput() ; // computes EM contributions to T_{mu,nu}
698 
699  // S_{rr} + S_{\theta\theta}
700  Tenseur SrrplusStt( Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
701  SrrplusStt = SrrplusStt / a_car ; // S^r_r + S^\theta_\theta
702 
703  Tenseur Spp (Cmp(Sij_I(3, 3))) ; //S_{\phi\phi}
704  Spp = Spp / b_car ; // S^\phi_\phi
705 
706  Cmp temp(E_I) ;
707  Tenseur E_int(temp) ;
708 
709  //-----------------------------------------------
710  // Sources of the Poisson equations
711  //-----------------------------------------------
712 
713  // Source for nu
714  // -------------
715  Tenseur beta = log(bbb) ;
716  beta.set_std_base() ;
717 
718  if (relativistic) {
719  source_nuf =
720  qpig * a_car *( ener_euler + s_euler + E_int + SrrplusStt + Spp ) ;
721 
722  source_nuq = ak_car
724  + beta.gradient_spher())
725  + qpig * a_car * 2*E_em ;
726  }
727  else {
728  source_nuf = qpig * nbar ;
729 
730  source_nuq = 0 ;
731  }
732  source_nuf.set_std_base() ;
733  source_nuq.set_std_base() ;
734 
735  // Source for dzeta
736  // ----------------
737  source_dzf = 2 * qpig * a_car
738  * (press + (ener_euler+press) * uuu*uuu + Spp ) ;
739  source_dzf.set_std_base() ;
740 
741  source_dzq = 2 * qpig * a_car * E_em + 1.5 * ak_car -
743  source_dzq.set_std_base() ;
744 
745 
746  // Source for tggg
747  // ---------------
748 
749  source_tggg = 2 * qpig * nnn * a_car * bbb
750  * ( 2*press + SrrplusStt ) ;
751  source_tggg.set_std_base() ;
752 
753  (source_tggg.set()).mult_rsint() ;
754 
755 
756  // Source for shift
757  // ----------------
758 
759  // Matter term:
760 
761  Cmp tjpem(J_I(3)) ;
762  tjpem += Jp_em() ;
763  tjpem.div_rsint() ;
764 
765  source_shift = (-4*qpig) * nnn * a_car * (ener_euler + press) * u_euler ;
766 
767  // Quadratic terms:
768  Tenseur vtmp = 6 * beta.gradient_spher() - 2 * logn.gradient_spher() ;
769  Tenseur mtmp(mp, 1, COV, mp.get_bvect_spher()) ;
770  if (tjpem.get_etat() == ETATZERO) mtmp.set_etat_zero() ;
771  else {
772  mtmp.set_etat_qcq() ;
773  mtmp.set(0) = 0 ;
774  mtmp.set(1) = 0 ;
775  mtmp.set(2) = (-4*qpig)*tjpem*nnn()*a_car()/b_car() ;
776  }
777  mtmp.change_triad(mp.get_bvect_cart()) ;
778 
779  vtmp.change_triad(mp.get_bvect_cart()) ;
780 
781  Tenseur squad = nnn * flat_scalar_prod(tkij, vtmp) ;
782 
783  // The addition of matter terms and quadratic terms is performed
784  // component by component because u_euler is contravariant,
785  // while squad is covariant.
786 
787  if (squad.get_etat() == ETATQCQ) {
788  for (int i=0; i<3; i++) {
789  source_shift.set(i) += squad(i) ;
790  }
791  }
792  if (mtmp.get_etat() == ETATQCQ) {
793  if (source_shift.get_etat() == ETATZERO) {
794  source_shift.set_etat_qcq() ;
795  for (int i=0; i<3; i++) {
796  source_shift.set(i) = mtmp(i) ;
797  source_shift.set(i).va.coef_i() ;
798  }
799  }
800  else
801  for (int i=0; i<3; i++)
802  source_shift.set(i) += mtmp(i) ;
803  }
804 
805  source_shift.set_std_base() ;
806 
807  //----------------------------------------------
808  // Resolution of the Poisson equation for nuf
809  //----------------------------------------------
810 
811  source_nuf().poisson(par_poisson_nuf, nuf.set()) ;
812 
813  if (relativistic) {
814 
815  //----------------------------------------------
816  // Resolution of the Poisson equation for nuq
817  //----------------------------------------------
818 
819  source_nuq().poisson(par_poisson_nuq, nuq.set()) ;
820 
821  //---------------------------------------------------------
822  // Resolution of the vector Poisson equation for the shift
823  //---------------------------------------------------------
824 
825 
826  if (source_shift.get_etat() != ETATZERO) {
827 
828  for (int i=0; i<3; i++) {
829  if(source_shift(i).dz_nonzero()) {
830  assert( source_shift(i).get_dzpuis() == 4 ) ;
831  }
832  else{
833  (source_shift.set(i)).set_dzpuis(4) ;
834  }
835  }
836 
837  }
838  //##
839  // source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
840 
841  double lambda_shift = double(1) / double(3) ;
842 
843  if ( mg->get_np(0) == 1 ) {
844  lambda_shift = 0 ;
845  }
846 
847  source_shift.poisson_vect(lambda_shift, par_poisson_vect,
848  shift, w_shift, khi_shift) ;
849 
850  // Computation of tnphi and nphi from the Cartesian components
851  // of the shift
852  // -----------------------------------------------------------
853 
854  fait_nphi() ;
855 
856  }
857 
858  //-----------------------------------------
859  // Determination of the fluid velociy U
860  //-----------------------------------------
861 
862  if (mer > mer_fix_omega + delta_mer_kep) {
863 
864  omega *= fact_omega ; // Increase of the angular velocity if
865  } // fact_omega != 1
866 
867  bool omega_trop_grand = false ;
868  bool kepler = true ;
869 
870  while ( kepler ) {
871 
872  // Possible decrease of Omega to ensure a velocity < c
873 
874  bool superlum = true ;
875 
876  while ( superlum ) {
877 
878  // New fluid velocity U :
879 
880  Cmp tmp = omega - nphi() ;
881  tmp.annule(nzm1) ;
882  tmp.std_base_scal() ;
883 
884  tmp.mult_rsint() ; // Multiplication by r sin(theta)
885 
886  uuu = bbb() / nnn() * tmp ;
887 
888  if (uuu.get_etat() == ETATQCQ) {
889  // Same basis as (Omega -N^phi) r sin(theta) :
890  ((uuu.set()).va).set_base( (tmp.va).base ) ;
891  }
892 
893  // Is the new velocity larger than c in the equatorial plane ?
894 
895  superlum = false ;
896 
897  for (int l=0; l<nzet; l++) {
898  for (int i=0; i<mg->get_nr(l); i++) {
899 
900  double u1 = uuu()(l, 0, j_b, i) ;
901  if (u1 >= 1.) { // superluminal velocity
902  superlum = true ;
903  cout << "U > c for l, i : " << l << " " << i
904  << " U = " << u1 << endl ;
905  }
906  }
907  }
908  if ( superlum ) {
909  cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
910  omega /= fact_omega ; // Decrease of Omega
911  cout << "New rotation frequency : "
912  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
913  omega_trop_grand = true ;
914  }
915  } // end of while ( superlum )
916 
917 
918  // New computation of U (which this time is not superluminal)
919  // as well as of gam_euler, ener_euler, etc...
920  // -----------------------------------
921 
922  hydro_euler() ;
923 
924  //------------------------------------------------------
925  // First integral of motion
926  //------------------------------------------------------
927 
928  // Centrifugal potential :
929  if (relativistic) {
930  mlngamma = - log( gam_euler ) ;
931  }
932  else {
933  mlngamma = - 0.5 * uuu*uuu ;
934  }
935 
936  Tenseur mag(mp) ;
937  if (is_conduct()) {
938  mag = mu0*M_j(A_phi, a_j) ;}
939  else{
940  mag = mu0*M_j(omega*A_phi-A_t, a_j) ;}
941 
942  // Equatorial values of various potentials :
943  double nuf_b = nuf()(l_b, k_b, j_b, i_b) ;
944  double nuq_b = nuq()(l_b, k_b, j_b, i_b) ;
945  double mlngamma_b = mlngamma()(l_b, k_b, j_b, i_b) ;
946  double mag_b = mag()(l_b, k_b, j_b, i_b) ;
947 
948  // Central values of various potentials :
949  double nuf_c = nuf()(0,0,0,0) ;
950  double nuq_c = nuq()(0,0,0,0) ;
951  double mlngamma_c = 0 ;
952  double mag_c = mag()(0,0,0,0) ;
953 
954  // Scale factor to ensure that the enthalpy is equal to ent_b at
955  // the equator
956  double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
957  + nuq_c - nuq_b + mag_c - mag_b)
958  / ( nuf_b - nuf_c ) ;
959  alpha_r = sqrt(alpha_r2) ;
960  cout << "alpha_r = " << alpha_r << endl ;
961 
962  // Readjustment of nu :
963  // -------------------
964 
965  logn = alpha_r2 * nuf + nuq ;
966  double nu_c = logn()(0,0,0,0) ;
967 
968  // First integral --> enthalpy in all space
969  //-----------------
970  ent = (ent_c + nu_c + mlngamma_c + mag_c) - logn - mlngamma - mag ;
971 
972  // Test: is the enthalpy negative somewhere in the equatorial plane
973  // inside the star ? If yes, this means that the Keplerian velocity
974  // has been overstep.
975 
976  kepler = false ;
977  for (int l=0; l<nzet; l++) {
978  int imax = mg->get_nr(l) - 1 ;
979  if (l == l_b) imax-- ; // The surface point is skipped
980  for (int i=0; i<imax; i++) {
981  if ( ent()(l, 0, j_b, i) < 0. ) {
982  kepler = true ;
983  cout << "ent < 0 for l, i : " << l << " " << i
984  << " ent = " << ent()(l, 0, j_b, i) << endl ;
985  }
986  }
987  }
988 
989  if ( kepler ) {
990  cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ;
991  omega /= fact_omega ; // Omega is decreased
992  cout << "New rotation frequency : "
993  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
994  omega_trop_grand = true ;
995  }
996 
997  } // End of while ( kepler )
998 
999  if ( omega_trop_grand ) { // fact_omega is decreased for the
1000  // next step
1001  fact_omega = sqrt( fact_omega ) ;
1002  cout << "**** New fact_omega : " << fact_omega << endl ;
1003  }
1004 
1005  //----------------------------------------------------
1006  // Adaptation of the mapping to the new enthalpy field
1007  //----------------------------------------------------
1008 
1009  // Shall the adaptation be performed (cusp) ?
1010  // ------------------------------------------
1011 
1012  double dent_eq = ent().dsdr()(l_b, k_b, j_b, i_b) ;
1013  double dent_pole = ent().dsdr()(l_b, k_b, 0, i_b) ;
1014  double rap_dent = fabs( dent_eq / dent_pole ) ;
1015  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
1016 
1017  if ( rap_dent < thres_adapt ) {
1018  adapt_flag = 0 ; // No adaptation of the mapping
1019  cout << "******* FROZEN MAPPING *********" << endl ;
1020  }
1021  else{
1022  adapt_flag = 1 ; // The adaptation of the mapping is to be
1023  // performed
1024  }
1025 
1026  mp_prev = mp_et ;
1027 
1028  mp.adapt(ent(), par_adapt) ;
1029 
1030  //----------------------------------------------------
1031  // Computation of the enthalpy at the new grid points
1032  //----------------------------------------------------
1033 
1034  mp_prev.homothetie(alpha_r) ;
1035 
1036  mp.reevaluate(&mp_prev, nzet+1, ent.set()) ;
1037 
1038  //----------------------------------------------------
1039  // Equation of state
1040  //----------------------------------------------------
1041 
1042  equation_of_state() ; // computes new values for nbar (n), ener (e)
1043  // and press (p) from the new ent (H)
1044 
1045  //---------------------------------------------------------
1046  // Matter source terms in the gravitational field equations
1047  //---------------------------------------------------------
1048 
1049  //## Computation of tnphi and nphi from the Cartesian components
1050  // of the shift for the test in hydro_euler():
1051 
1052  fait_nphi() ;
1053 
1054  hydro_euler() ; // computes new values for ener_euler (E),
1055  // s_euler (S) and u_euler (U^i)
1056 
1057  if (relativistic) {
1058 
1059  //-------------------------------------------------------
1060  // 2-D Poisson equation for tggg
1061  //-------------------------------------------------------
1062 
1063  mp.poisson2d(source_tggg(), mp.cmp_zero(), par_poisson_tggg, tggg.set()) ;
1064 
1065  //-------------------------------------------------------
1066  // 2-D Poisson equation for dzeta
1067  //-------------------------------------------------------
1068 
1069  mp.poisson2d(source_dzf(), source_dzq(), par_poisson_dzeta, dzeta.set()) ;
1070 
1071  err_grv2 = lbda_grv2 - 1;
1072  cout << "GRV2: " << err_grv2 << endl ;
1073 
1074  }
1075  else {
1076  err_grv2 = grv2() ;
1077  }
1078 
1079  //---------------------------------------
1080  // Computation of the metric coefficients (except for N^phi)
1081  //---------------------------------------
1082 
1083  // Relaxations on nu and dzeta :
1084 
1085  if (mer >= 10) {
1086  logn = relax * logn + relax_prev * logn_prev ;
1087 
1088  dzeta = relax * dzeta + relax_prev * dzeta_prev ;
1089  }
1090 
1091  // Update of the metric coefficients N, A, B and computation of K_ij :
1092 
1093  update_metric() ;
1094 
1095  //-----------------------
1096  // Informations display
1097  //-----------------------
1098 
1099  // partial_display(cout) ;
1100  fichfreq << " " << omega / (2*M_PI) * f_unit ;
1101  fichevol << " " << rap_dent ;
1102  fichevol << " " << ray_pole() / ray_eq() ;
1103  fichevol << " " << ent_c ;
1104 
1105  //-----------------------------------------
1106  // Convergence towards a given baryon mass
1107  //-----------------------------------------
1108 
1109  if (mer > mer_mass) {
1110 
1111  double xx ;
1112  if (mbar_wanted > 0.) {
1113  xx = mass_b() / mbar_wanted - 1. ;
1114  cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx
1115  << endl ;
1116  }
1117  else{
1118  xx = mass_g() / fabs(mbar_wanted) - 1. ;
1119  cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx
1120  << endl ;
1121  }
1122  double xprog = ( mer > 2*mer_mass) ? 1. :
1123  double(mer-mer_mass)/double(mer_mass) ;
1124  xx *= xprog ;
1125  double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
1126  double fact = pow(ax, aexp_mass) ;
1127  cout << " xprog, xx, ax, fact : " << xprog << " " <<
1128  xx << " " << ax << " " << fact << endl ;
1129 
1130  if ( change_ent ) {
1131  ent_c *= fact ;
1132  }
1133  else {
1134  if (mer%4 == 0) omega *= fact ;
1135  }
1136  }
1137 
1138  //---------------------------------------------
1139  // Convergence towards a given magnetic moment
1140  //---------------------------------------------
1141 
1142  if (mer > mer_magmom) {
1143 
1144  // if(mer > mer_mass) {
1145  // cerr << "et_magnetisation::equilibrium()" << endl ;
1146  // cerr << "Impossible to converge to a given baryon mass" << endl ;
1147  // cerr << "and a given magnetic moment!" << endl ;
1148  // cerr << "Aborting..." << endl ;
1149  // abort() ;
1150  // }
1151  double xx = MagMom() / magmom_wanted - 1. ;
1152  cout << "Discrep. mag. moment <-> wanted mag. moment : " << xx
1153  << endl ;
1154 
1155  double xprog = ( mer > 2*mer_magmom) ? 1. :
1156  double(mer-mer_magmom)/double(mer_magmom) ;
1157  xx *= xprog ;
1158  double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
1159  double fact = pow(ax, aexp_mass) ;
1160  cout << " xprog, xx, ax, fact : " << xprog << " " <<
1161  xx << " " << ax << " " << fact << endl ;
1162 
1163  a_j *= fact ;
1164  }
1165 
1166  //------------------------------------------------------------
1167  // Relative change in enthalpy with respect to previous step
1168  //------------------------------------------------------------
1169 
1170  Tbl diff_ent_tbl = diffrel( ent(), ent_prev() ) ;
1171  diff_ent = diff_ent_tbl(0) ;
1172  for (int l=1; l<nzet; l++) {
1173  diff_ent += diff_ent_tbl(l) ;
1174  }
1175  diff_ent /= nzet ;
1176 
1177  fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
1178  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
1179 
1180  //------------------------------
1181  // Recycling for the next step
1182  //------------------------------
1183 
1184  ent_prev = ent ;
1185  logn_prev = logn ;
1186  dzeta_prev = dzeta ;
1187 
1188  fichconv << endl ;
1189  fichfreq << endl ;
1190  fichevol << endl ;
1191  fichconv.flush() ;
1192  fichfreq.flush() ;
1193  fichevol.flush() ;
1194 
1195  } // End of main loop
1196 
1197  //=========================================================================
1198  // End of iteration
1199  //=========================================================================
1200 
1201  fichconv.close() ;
1202  fichfreq.close() ;
1203  fichevol.close() ;
1204 }
1205 
1206 
1207 
1208 
1209 
1210 
1211 
1212 
1213 
1214 
1215 
1216 
1217 
1218 
1219 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:500
Tenseur khi_shift
Scalar used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: etoile.h:1560
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:372
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1142
virtual double nbar_ent_p(double ent, const Param *par=0x0) const
Computes the baryon density from the log-enthalpy.
Definition: eos_mag.C:330
Tenseur E_em
electromagnetic energy density in the Eulerian frame
Definition: et_rot_mag.h:161
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: etoile.h:1616
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Radial mapping of rather general form.
Definition: map.h:2752
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition: tenseur.C:1548
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
Definition: map.h:807
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Et_magnetisation(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i, bool include_mag=true, bool use_B=true)
Standard constructor.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Multi-domain array.
Definition: mtbl.h:118
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift ...
Definition: etoile_rot.C:781
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
double Q
In the case of a perfect conductor, the requated baryonic charge.
Definition: et_rot_mag.h:179
Lorene prototypes.
Definition: app_hor.h:64
Tenseur w_shift
Vector used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: etoile.h:1550
Scalar E_I
Interaction (magnetisation) energy density.
Definition: et_rot_mag.h:624
Equation of state base class.
Definition: eos.h:190
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for ...
Definition: etoile.h:1625
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
int get_etat() const
Returns the logical state.
Definition: cmp.h:896
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:468
double ray_eq() const
Coordinate radius at , [r_unit].
Base class for coordinate mappings.
Definition: map.h:670
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1507
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Tenseur Magn() const
Computes the magnetic field spherical components in Lorene&#39;s units.
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
Basic integer array class.
Definition: itbl.h:122
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the radial direction.
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
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:906
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:365
Tenseur press
Fluid pressure.
Definition: etoile.h:461
double * x
Array of values of at the nr collocation points.
Definition: grilles.h:209
Tensor field of valence 1.
Definition: vector.h:188
Tenseur shift
Total shift vector.
Definition: etoile.h:512
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition: param.C:453
Cmp A_phi
-component of the electromagnetic potential 1-form divided by .
Definition: et_rot_mag.h:155
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
void update_metric()
Computes metric coefficients from known potentials.
Definition: et_rot_upmetr.C:69
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:686
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:474
Cmp ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition: etoile.h:1598
void operator=(const Et_rot_mag &)
Assignment to another Et_rot_mag.
Definition: et_rot_mag.C:283
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
virtual double press_ent_p(double ent, const Param *par=0x0) const
Computes the pressure from the log-enthalpy.
Definition: eos_mag.C:419
double a_j
Amplitude of the curent/charge function.
Definition: et_rot_mag.h:180
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.
bool is_conduct() const
Tells if the star is made of conducting or isolating material.
Definition: et_rot_mag.h:241
Vector J_I
Interaction momentum density 3-vector.
Definition: et_rot_mag.h:627
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
Scalar xmag
The magnetisation scalar.
Definition: et_rot_mag.h:622
bool use_B_in_eos
Flag : true if the value of the magnetic field is used in the Eos.
Definition: et_rot_mag.h:617
Class for magnetized (isolator or perfect conductor), rigidly rotating stars.
Definition: et_rot_mag.h:143
double mag_ent_p(double ent, const Param *par=0x0) const
Computes the magnetisation.
Definition: eos_mag.C:460
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: et_rot_mag.C:336
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:471
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
bool include_magnetisation
Flag : true if magnetisation terms are included in the equations.
Definition: et_rot_mag.h:620
virtual double mass_b() const
Baryon mass.
Parameter storage.
Definition: param.h:125
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:522
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:116
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
const Eos & eos
Equation of state of the stellar matter.
Definition: etoile.h:451
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:905
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Cmp ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition: etoile.h:1592
Tenseur ak_car
Scalar .
Definition: etoile.h:1586
virtual void sauve(FILE *) const
Save in a file.
int get_etat() const
Returns the logical state.
Definition: tenseur.h:704
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: et_rot_mag.C:258
void operator=(const Et_magnetisation &)
Assignment to another Et_rot_mag.
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition: et_rot_hydro.C:83
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Class for a magnetized (tabulated) equation of state.
Definition: eos_mag.h:78
Tenseur uuu
Norm of u_euler.
Definition: etoile.h:1518
Tbl & set(int l)
Read/write of the Tbl in a given domain.
Definition: mtbl.h:221
Tenseur tggg
Metric potential .
Definition: etoile.h:1537
double omega
Rotation angular velocity ([f_unit] )
Definition: etoile.h:1501
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: tbl.C:347
Tenseur nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: etoile.h:1526
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
virtual double grv2() const
Error on the virial identity GRV2.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: cmp.C:323
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:437
Multi-domain grid.
Definition: grilles.h:273
double ray_pole() const
Coordinate radius at [r_unit].
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:721
Tenseur ener
Total energy density in the fluid frame.
Definition: etoile.h:460
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:684
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
void equilibrium_mag(double ent_c, double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, double mbar_wanted, double magmom_wanted, double aexp_mass, Tbl &diff, double Q0, double a_j0, Cmp(*f_j)(const Cmp &x, const double), Cmp(*M_j)(const Cmp &x, const double))
Computes an equilibrium configuration.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:457
virtual void sauve(FILE *) const
Save in a file.
Definition: et_rot_mag.C:311
virtual double ener_ent_p(double ent, const Param *par=0x0) const
Computes the total energy density from the log-enthalpy.
Definition: eos_mag.C:373
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
virtual void magnet_comput(const int adapt_flag, Cmp(*f_j)(const Cmp &x, const double), Param &par_poisson_At, Param &par_poisson_Avect)
Computes the electromagnetic quantities solving the Maxwell equations (6) and (7) of [Bocquet...
Sym_tensor Sij_I
Interaction stress 3-tensor.
Definition: et_rot_mag.h:630
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.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
virtual double mass_g() const
Gravitational mass.
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
Definition: et_rot_mag.h:150
Tenseur & logn
Metric potential = logn_auto.
Definition: etoile.h:1521
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:701
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
virtual ~Et_magnetisation()
Destructor.
double MagMom() const
Magnetic Momentum in SI units.
Tenseur nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: etoile.h:1531
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:465
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:372
Tenseur & dzeta
Metric potential = beta_auto.
Definition: etoile.h:1534
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:461
Cmp ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg ...
Definition: etoile.h:1608
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition: etoile.h:1567
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Standard electro-magnetic units.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
Tenseur Jp_em
component of the electromagnetic momentum density 3-vector, as measured in the Eulerian frame...
Definition: et_rot_mag.h:167
virtual void MHD_comput()
Computes the electromagnetic part of the stress-energy tensor.