LORENE
et_rot_mag_equil.C
1 /*
2  * Function et_rot_mag::equilibrium_mag
3  *
4  * Computes rotating equilibirum with a magnetic field
5  * (see file et_rot_mag.h for documentation)
6  *
7  */
8 
9 /*
10  * Copyright (c) 2002 Eric Gourgoulhon
11  * Copyright (c) 2002 Emmanuel Marcq
12  * Copyright (c) 2002 Jerome Novak
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License as published by
18  * the Free Software Foundation; either version 2 of the License, or
19  * (at your option) any later version.
20  *
21  * LORENE is distributed in the hope that it will be useful,
22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24  * GNU General Public License for more details.
25  *
26  * You should have received a copy of the GNU General Public License
27  * along with LORENE; if not, write to the Free Software
28  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29  *
30  */
31 
32 char et_rot_mag_equil_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag_equil.C,v 1.21 2014/10/13 08:52:58 j_novak Exp $" ;
33 
34 /*
35  * $Id: et_rot_mag_equil.C,v 1.21 2014/10/13 08:52:58 j_novak Exp $
36  * $Log: et_rot_mag_equil.C,v $
37  * Revision 1.21 2014/10/13 08:52:58 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.20 2014/10/06 15:13:09 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.19 2014/09/03 15:33:42 j_novak
44  * Filtering of Maxwell sources is now optional.
45  *
46  * Revision 1.18 2012/08/12 17:48:36 p_cerda
47  * Magnetstar: New classes for magnetstar. Allowing for non-equatorial symmetry in Etoile et al. Adding B_phi in Et_rot_mag.
48  *
49  * Revision 1.17 2004/03/25 10:43:04 j_novak
50  * Some units forgotten...
51  *
52  * Revision 1.16 2003/11/19 22:01:57 e_gourgoulhon
53  * -- Relaxation on logn and dzeta performed only if mer >= 10.
54  * -- err_grv2 is now evaluated also in the Newtonian case.
55  *
56  * Revision 1.15 2003/10/27 10:54:43 e_gourgoulhon
57  * Changed local variable name lambda_grv2 to lbda_grv2 in order not
58  * to shadow method name.
59  *
60  * Revision 1.14 2003/10/03 15:58:47 j_novak
61  * Cleaning of some headers
62  *
63  * Revision 1.13 2002/10/16 14:36:36 j_novak
64  * Reorganization of #include instructions of standard C++, in order to
65  * use experimental version 3 of gcc.
66  *
67  * Revision 1.12 2002/10/11 11:47:35 j_novak
68  * Et_rot_mag::MHD_comput is now virtual.
69  * Use of standard constructor for Tenseur mtmp in Et_rot_mag::equilibrium_mag
70  *
71  * Revision 1.11 2002/06/05 15:15:59 j_novak
72  * The case of non-adapted mapping is treated.
73  * parmag.d and parrot.d have been merged.
74  *
75  * Revision 1.10 2002/06/03 13:23:16 j_novak
76  * The case when the mapping is not adapted is now treated
77  *
78  * Revision 1.9 2002/06/03 13:00:45 e_marcq
79  *
80  * conduc parameter read in parmag.d
81  *
82  * Revision 1.6 2002/05/17 15:08:01 e_marcq
83  *
84  * Rotation progressive plug-in, units corrected, Q and a_j new member data
85  *
86  * Revision 1.5 2002/05/16 10:02:09 j_novak
87  * Errors in stress energy tensor corrected
88  *
89  * Revision 1.4 2002/05/15 09:53:59 j_novak
90  * First operational version
91  *
92  * Revision 1.3 2002/05/14 13:38:36 e_marcq
93  *
94  *
95  * Unit update, new outputs
96  *
97  * Revision 1.1 2002/05/10 09:26:52 j_novak
98  * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
99  *
100  *
101  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag_equil.C,v 1.21 2014/10/13 08:52:58 j_novak Exp $
102  *
103  */
104 
105 // Headers C
106 #include <cmath>
107 
108 // Headers Lorene
109 #include "et_rot_mag.h"
110 #include "param.h"
111 #include "unites.h"
112 
113 namespace Lorene {
114 void Et_rot_mag::equilibrium_mag(double ent_c, double omega0,
115  double fact_omega, int nzadapt, const Tbl& ent_limit,
116  const Itbl& icontrol, const Tbl& control, double mbar_wanted,
117  double aexp_mass, Tbl& diff, const double Q0, const double a_j0,
118  Cmp (*f_j)(const Cmp&, const double),
119  Cmp (*M_j)(const Cmp& x, const double)) {
120 
121  // Fundamental constants and units
122  // -------------------------------
123  using namespace Unites_mag ;
124 
125  // For the display
126  // ---------------
127  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
128  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
129 
130  // Grid parameters
131  // ---------------
132 
133  const Mg3d* mg = mp.get_mg() ;
134  int nz = mg->get_nzone() ; // total number of domains
135  int nzm1 = nz - 1 ;
136 
137  // The following is required to initialize mp_prev as a Map_et:
138  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
139 
140  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
141  // assert(mg->get_type_t() == SYM) ;
142  int l_b = nzet - 1 ;
143  int i_b = mg->get_nr(l_b) - 1 ;
144  int j_b = mg->get_nt(l_b) - 1 ;
145  int k_b = 0 ;
146 
147  // Value of the enthalpy defining the surface of the star
148  double ent_b = ent_limit(nzet-1) ;
149 
150  // Parameters to control the iteration
151  // -----------------------------------
152 
153  int mer_max = icontrol(0) ;
154  int mer_rot = icontrol(1) ;
155  int mer_change_omega = icontrol(2) ;
156  int mer_fix_omega = icontrol(3) ;
157  int mer_mass = icontrol(4) ;
158  int mermax_poisson = icontrol(5) ;
159  int delta_mer_kep = icontrol(6) ;
160  int mer_mag = icontrol(7) ;
161  int mer_change_mag = icontrol(8) ;
162  int mer_fix_mag = icontrol(9) ;
163  int mag_filter = icontrol(10) ;
164 
165  // Protections:
166  if (mer_change_omega < mer_rot) {
167  cout << "Etoile_rot::equilibrium: mer_change_omega < mer_rot !" << endl ;
168  cout << " mer_change_omega = " << mer_change_omega << endl ;
169  cout << " mer_rot = " << mer_rot << endl ;
170  abort() ;
171  }
172  if (mer_fix_omega < mer_change_omega) {
173  cout << "Etoile_rot::equilibrium: mer_fix_omega < mer_change_omega !"
174  << endl ;
175  cout << " mer_fix_omega = " << mer_fix_omega << endl ;
176  cout << " mer_change_omega = " << mer_change_omega << endl ;
177  abort() ;
178  }
179 
180  // In order to converge to a given baryon mass, shall the central
181  // enthalpy be varied or Omega ?
182  bool change_ent = true ;
183  if (mer_mass < 0) {
184  change_ent = false ;
185  mer_mass = abs(mer_mass) ;
186  }
187 
188  double precis = control(0) ;
189  double omega_ini = control(1) ;
190  double relax = control(2) ;
191  double relax_prev = double(1) - relax ;
192  double relax_poisson = control(3) ;
193  double thres_adapt = control(4) ;
194  double precis_adapt = control(5) ;
195  double Q_ini = control(6) ;
196  double a_j_ini = control (7) ;
197 
198  // Error indicators
199  // ----------------
200 
201  diff.set_etat_qcq() ;
202  double& diff_ent = diff.set(0) ;
203 
204  // Parameters for the function Map_et::adapt
205  // -----------------------------------------
206 
207  Param par_adapt ;
208  int nitermax = 100 ;
209  int niter ;
210  int adapt_flag = 1 ; // 1 = performs the full computation,
211  // 0 = performs only the rescaling by
212  // the factor alpha_r
213  int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
214  // isosurfaces
215  double alpha_r ;
216  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
217 
218  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
219  // locate zeros by the secant method
220  par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
221  // to the isosurfaces of ent is to be
222  // performed
223  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
224  // the enthalpy isosurface
225  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
226  // 0 = performs only the rescaling by
227  // the factor alpha_r
228  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
229  // (theta_*, phi_*)
230  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
231  // (theta_*, phi_*)
232 
233  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
234  // the secant method
235 
236  par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
237  // the determination of zeros by
238  // the secant method
239  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
240 
241  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
242  // distances will be multiplied
243 
244  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
245  // to define the isosurfaces.
246 
247  // Parameters for the function Map_et::poisson for nuf
248  // ----------------------------------------------------
249 
250  double precis_poisson = 1.e-16 ;
251 
252  Param par_poisson_nuf ;
253  par_poisson_nuf.add_int(mermax_poisson, 0) ; // maximum number of iterations
254  par_poisson_nuf.add_double(relax_poisson, 0) ; // relaxation parameter
255  par_poisson_nuf.add_double(precis_poisson, 1) ; // required precision
256  par_poisson_nuf.add_int_mod(niter, 0) ; // number of iterations actually used
257  par_poisson_nuf.add_cmp_mod( ssjm1_nuf ) ;
258 
259  Param par_poisson_nuq ;
260  par_poisson_nuq.add_int(mermax_poisson, 0) ; // maximum number of iterations
261  par_poisson_nuq.add_double(relax_poisson, 0) ; // relaxation parameter
262  par_poisson_nuq.add_double(precis_poisson, 1) ; // required precision
263  par_poisson_nuq.add_int_mod(niter, 0) ; // number of iterations actually used
264  par_poisson_nuq.add_cmp_mod( ssjm1_nuq ) ;
265 
266  Param par_poisson_tggg ;
267  par_poisson_tggg.add_int(mermax_poisson, 0) ; // maximum number of iterations
268  par_poisson_tggg.add_double(relax_poisson, 0) ; // relaxation parameter
269  par_poisson_tggg.add_double(precis_poisson, 1) ; // required precision
270  par_poisson_tggg.add_int_mod(niter, 0) ; // number of iterations actually used
271  par_poisson_tggg.add_cmp_mod( ssjm1_tggg ) ;
272  double lambda_tggg ;
273  par_poisson_tggg.add_double_mod( lambda_tggg ) ;
274 
275  Param par_poisson_dzeta ;
276  double lbda_grv2 ;
277  par_poisson_dzeta.add_double_mod( lbda_grv2 ) ;
278 
279  // Parameters for the function Tenseur::poisson_vect
280  // -------------------------------------------------
281 
282  Param par_poisson_vect ;
283 
284  par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of iterations
285  par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
286  par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
287  par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
288  par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
289  par_poisson_vect.add_int_mod(niter, 0) ;
290 
291 
292  // Parameters for the Maxwell equations
293  // -------------------------------------
294 
295  Param par_poisson_At ; // For scalar At Poisson equation
296  Cmp ssjm1_At(mp) ;
297  ssjm1_At.set_etat_zero() ;
298  par_poisson_At.add_int(mermax_poisson, 0) ; // maximum number of iterations
299  par_poisson_At.add_double(relax_poisson, 0) ; // relaxation parameter
300  par_poisson_At.add_double(precis_poisson, 1) ; // required precision
301  par_poisson_At.add_int_mod(niter, 0) ; // number of iterations actually used
302  par_poisson_At.add_cmp_mod( ssjm1_At ) ;
303  par_poisson_At.add_int(mag_filter, 1) ; //filtering of Maxwell sources
304 
305  Param par_poisson_Avect ; // For vector Aphi Poisson equation
306 
307  Cmp ssjm1_khi_mag(ssjm1_khi) ;
308  Tenseur ssjm1_w_mag(ssjm1_wshift) ;
309 
310  par_poisson_Avect.add_int(mermax_poisson, 0) ; // maximum number of iterations
311  par_poisson_Avect.add_double(relax_poisson, 0) ; // relaxation parameter
312  par_poisson_Avect.add_double(precis_poisson, 1) ; // required precision
313  par_poisson_Avect.add_cmp_mod( ssjm1_khi_mag ) ;
314  par_poisson_Avect.add_tenseur_mod( ssjm1_w_mag ) ;
315  par_poisson_Avect.add_int_mod(niter, 0) ;
316  par_poisson_Avect.add_int(mag_filter, 1) ; //filtering of Maxwell sources
317 
318 
319  // Initializations
320  // ---------------
321 
322  // Initial angular velocity / magnetic quantities
323  omega = 0 ;
324  Q = 0 ;
325  a_j = 0 ;
326 
327  double accrois_omega = (omega0 - omega_ini) /
328  double(mer_fix_omega - mer_change_omega) ;
329  double accrois_Q = (Q0 - Q_ini) /
330  double(mer_fix_mag - mer_change_mag);
331  double accrois_a_j = (a_j0 - a_j_ini) /
332  double(mer_fix_mag - mer_change_mag);
333 
334  update_metric() ; // update of the metric coefficients
335 
336  equation_of_state() ; // update of the density, pressure, etc...
337 
338  hydro_euler() ; // update of the hydro quantities relative to the
339  // Eulerian observer
340 
341  MHD_comput() ; // update of EM contributions to stress-energy tensor
342 
343 
344  // Quantities at the previous step :
345  Map_et mp_prev = mp_et ;
346  Tenseur ent_prev = ent ;
347  Tenseur logn_prev = logn ;
348  Tenseur dzeta_prev = dzeta ;
349 
350  // Creation of uninitialized tensors:
351  Tenseur source_nuf(mp) ; // source term in the equation for nuf
352  Tenseur source_nuq(mp) ; // source term in the equation for nuq
353  Tenseur source_dzf(mp) ; // matter source term in the eq. for dzeta
354  Tenseur source_dzq(mp) ; // quadratic source term in the eq. for dzeta
355  Tenseur source_tggg(mp) ; // source term in the eq. for tggg
356  Tenseur source_shift(mp, 1, CON, mp.get_bvect_cart()) ;
357  // source term for shift
358  Tenseur mlngamma(mp) ; // centrifugal potential
359 
360  // Preparations for the Poisson equations:
361  // --------------------------------------
362  if (nuf.get_etat() == ETATZERO) {
363  nuf.set_etat_qcq() ;
364  nuf.set() = 0 ;
365  }
366 
367  if (relativistic) {
368  if (nuq.get_etat() == ETATZERO) {
369  nuq.set_etat_qcq() ;
370  nuq.set() = 0 ;
371  }
372 
373  if (tggg.get_etat() == ETATZERO) {
374  tggg.set_etat_qcq() ;
375  tggg.set() = 0 ;
376  }
377 
378  if (dzeta.get_etat() == ETATZERO) {
379  dzeta.set_etat_qcq() ;
380  dzeta.set() = 0 ;
381  }
382  }
383 
384  ofstream fichconv("convergence.d") ; // Output file for diff_ent
385  fichconv << "# diff_ent GRV2 " << endl ;
386 
387  ofstream fichfreq("frequency.d") ; // Output file for omega
388  fichfreq << "# f [Hz]" << endl ;
389 
390  ofstream fichevol("evolution.d") ; // Output file for various quantities
391  fichevol <<
392  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
393  << endl ;
394 
395  diff_ent = 1 ;
396  double err_grv2 = 1 ;
397 
398  //=========================================================================
399  // Start of iteration
400  //=========================================================================
401 
402  for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
403 
404  cout << "-----------------------------------------------" << endl ;
405  cout << "step: " << mer << endl ;
406  cout << "diff_ent = " << display_bold << diff_ent << display_normal
407  << endl ;
408  cout << "err_grv2 = " << err_grv2 << endl ;
409  fichconv << mer ;
410  fichfreq << mer ;
411  fichevol << mer ;
412 
413  if (mer >= mer_rot) {
414 
415  if (mer < mer_change_omega) {
416  omega = omega_ini ;
417  }
418  else {
419  if (mer <= mer_fix_omega) {
420  omega = omega_ini + accrois_omega *
421  (mer - mer_change_omega) ;
422  }
423  }
424 
425  }
426 
427  if (mer >= mer_mag) {
428  if (mer < mer_change_mag) {
429  Q = Q_ini ;
430  a_j = a_j_ini ;
431  }
432  else {
433  if (mer <= mer_fix_mag) {
434  Q = Q_ini + accrois_Q * (mer - mer_change_mag) ;
435  a_j = a_j_ini + accrois_a_j * (mer - mer_change_mag) ;
436  }
437  }
438  }
439 
440 
441  //-----------------------------------------------
442  // Computation of electromagnetic potentials :
443  // -------------------------------------------
444  magnet_comput(adapt_flag,
445  f_j, par_poisson_At, par_poisson_Avect) ;
446 
447  MHD_comput() ; // computes EM contributions to T_{mu,nu}
448 
449  //-----------------------------------------------
450  // Sources of the Poisson equations
451  //-----------------------------------------------
452 
453  // Source for nu
454  // -------------
455  Tenseur beta = log(bbb) ;
456  beta.set_std_base() ;
457 
458  if (relativistic) {
459  source_nuf = qpig * a_car *( ener_euler + s_euler) ;
460 
461  source_nuq = ak_car - flat_scalar_prod(logn.gradient_spher(),
462  logn.gradient_spher() + beta.gradient_spher())
463  + qpig * a_car * 2*E_em ;
464  }
465  else {
466  source_nuf = qpig * nbar ;
467 
468  source_nuq = 0 ;
469  }
470  source_nuf.set_std_base() ;
471  source_nuq.set_std_base() ;
472 
473  // Source for dzeta
474  // ----------------
475  source_dzf = 2 * qpig * a_car * (press + (ener_euler+press) * uuu*uuu ) ;
476  source_dzf.set_std_base() ;
477 
478  source_dzq = 2 * qpig * a_car * E_em + 1.5 * ak_car -
480  source_dzq.set_std_base() ;
481 
482  // Source for tggg
483  // ---------------
484 
485  source_tggg = 4 * qpig * nnn * a_car * bbb * press ;
486  source_tggg.set_std_base() ;
487 
488  (source_tggg.set()).mult_rsint() ;
489 
490 
491  // Source for shift
492  // ----------------
493 
494  // Matter term:
495 
496  Cmp tjpem(Jp_em()) ;
497  tjpem.div_rsint() ;
498 
499  source_shift = (-4*qpig) * nnn * a_car * (ener_euler + press)
500  * u_euler ;
501 
502  // Quadratic terms:
503  Tenseur vtmp = 6 * beta.gradient_spher() - 2 * logn.gradient_spher() ;
504  Tenseur mtmp(mp, 1, COV, mp.get_bvect_spher()) ;
505  if (tjpem.get_etat() == ETATZERO) mtmp.set_etat_zero() ;
506  else {
507  mtmp.set_etat_qcq() ;
508  mtmp.set(0) = 0 ;
509  mtmp.set(1) = 0 ;
510  mtmp.set(2) = (-4*qpig)*tjpem*nnn()*a_car()/b_car() ;
511  }
512  mtmp.change_triad(mp.get_bvect_cart()) ;
513 
514  vtmp.change_triad(mp.get_bvect_cart()) ;
515 
516  Tenseur squad = nnn * flat_scalar_prod(tkij, vtmp) ;
517 
518  // The addition of matter terms and quadratic terms is performed
519  // component by component because u_euler is contravariant,
520  // while squad is covariant.
521 
522  if (squad.get_etat() == ETATQCQ) {
523  for (int i=0; i<3; i++) {
524  source_shift.set(i) += squad(i) ;
525  }
526  }
527  if (mtmp.get_etat() == ETATQCQ) {
528  if (source_shift.get_etat() == ETATZERO) {
529  source_shift.set_etat_qcq() ;
530  for (int i=0; i<3; i++) {
531  source_shift.set(i) = mtmp(i) ;
532  source_shift.set(i).va.coef_i() ;
533  }
534  }
535  else
536  for (int i=0; i<3; i++)
537  source_shift.set(i) += mtmp(i) ;
538  }
539 
540  source_shift.set_std_base() ;
541 
542  //----------------------------------------------
543  // Resolution of the Poisson equation for nuf
544  //----------------------------------------------
545 
546  source_nuf().poisson(par_poisson_nuf, nuf.set()) ;
547 
548  if (relativistic) {
549 
550  //----------------------------------------------
551  // Resolution of the Poisson equation for nuq
552  //----------------------------------------------
553 
554  source_nuq().poisson(par_poisson_nuq, nuq.set()) ;
555 
556  //---------------------------------------------------------
557  // Resolution of the vector Poisson equation for the shift
558  //---------------------------------------------------------
559 
560 
561  if (source_shift.get_etat() != ETATZERO) {
562 
563  for (int i=0; i<3; i++) {
564  if(source_shift(i).dz_nonzero()) {
565  assert( source_shift(i).get_dzpuis() == 4 ) ;
566  }
567  else{
568  (source_shift.set(i)).set_dzpuis(4) ;
569  }
570  }
571 
572  }
573  //##
574  // source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
575 
576  double lambda_shift = double(1) / double(3) ;
577 
578  if ( mg->get_np(0) == 1 ) {
579  lambda_shift = 0 ;
580  }
581 
582  source_shift.poisson_vect(lambda_shift, par_poisson_vect,
583  shift, w_shift, khi_shift) ;
584 
585  // Computation of tnphi and nphi from the Cartesian components
586  // of the shift
587  // -----------------------------------------------------------
588 
589  fait_nphi() ;
590 
591  //## cout.precision(10) ;
592  // cout << "nphi : " << nphi()(0, 0, 0, 0)
593  // << " " << nphi()(l_b, k_b, j_b, i_b/2)
594  // << " " << nphi()(l_b, k_b, j_b, i_b) << endl ;
595 
596  }
597 
598  //-----------------------------------------
599  // Determination of the fluid velociy U
600  //-----------------------------------------
601 
602  if (mer > mer_fix_omega + delta_mer_kep) {
603 
604  omega *= fact_omega ; // Increase of the angular velocity if
605  } // fact_omega != 1
606 
607  bool omega_trop_grand = false ;
608  bool kepler = true ;
609 
610  while ( kepler ) {
611 
612  // Possible decrease of Omega to ensure a velocity < c
613 
614  bool superlum = true ;
615 
616  while ( superlum ) {
617 
618  // New fluid velocity U :
619 
620  Cmp tmp = omega - nphi() ;
621  tmp.annule(nzm1) ;
622  tmp.std_base_scal() ;
623 
624  tmp.mult_rsint() ; // Multiplication by r sin(theta)
625 
626  uuu = bbb() / nnn() * tmp ;
627 
628  if (uuu.get_etat() == ETATQCQ) {
629  // Same basis as (Omega -N^phi) r sin(theta) :
630  ((uuu.set()).va).set_base( (tmp.va).base ) ;
631  }
632 
633 
634  // Is the new velocity larger than c in the equatorial plane ?
635 
636  superlum = false ;
637 
638  for (int l=0; l<nzet; l++) {
639  for (int i=0; i<mg->get_nr(l); i++) {
640 
641  double u1 = uuu()(l, 0, j_b, i) ;
642  if (u1 >= 1.) { // superluminal velocity
643  superlum = true ;
644  cout << "U > c for l, i : " << l << " " << i
645  << " U = " << u1 << endl ;
646  }
647  }
648  }
649  if ( superlum ) {
650  cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
651  omega /= fact_omega ; // Decrease of Omega
652  cout << "New rotation frequency : "
653  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
654  omega_trop_grand = true ;
655  }
656  } // end of while ( superlum )
657 
658 
659  // New computation of U (which this time is not superluminal)
660  // as well as of gam_euler, ener_euler, etc...
661  // -----------------------------------
662 
663  hydro_euler() ;
664 
665 
666  //------------------------------------------------------
667  // First integral of motion
668  //------------------------------------------------------
669 
670  // Centrifugal potential :
671  if (relativistic) {
672  mlngamma = - log( gam_euler ) ;
673  }
674  else {
675  mlngamma = - 0.5 * uuu*uuu ;
676  }
677 
678  Tenseur mag(mp) ;
679  if (is_conduct()) {
680  mag = mu0*M_j(A_phi, a_j) ;}
681  else{
682  mag = mu0*M_j(omega*A_phi-A_t, a_j) ;}
683 
684  // Equatorial values of various potentials :
685  double nuf_b = nuf()(l_b, k_b, j_b, i_b) ;
686  double nuq_b = nuq()(l_b, k_b, j_b, i_b) ;
687  double mlngamma_b = mlngamma()(l_b, k_b, j_b, i_b) ;
688  double mag_b = mag()(l_b, k_b, j_b, i_b) ;
689 
690  // Central values of various potentials :
691  double nuf_c = nuf()(0,0,0,0) ;
692  double nuq_c = nuq()(0,0,0,0) ;
693  double mlngamma_c = 0 ;
694  double mag_c = mag()(0,0,0,0) ;
695 
696  // Scale factor to ensure that the enthalpy is equal to ent_b at
697  // the equator
698  double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
699  + nuq_c - nuq_b + mag_c - mag_b)
700  / ( nuf_b - nuf_c ) ;
701  alpha_r = sqrt(alpha_r2) ;
702  cout << "alpha_r = " << alpha_r << endl ;
703 
704  // Readjustment of nu :
705  // -------------------
706 
707  logn = alpha_r2 * nuf + nuq ;
708  double nu_c = logn()(0,0,0,0) ;
709 
710 
711  // First integral --> enthalpy in all space
712  //-----------------
713  ent = (ent_c + nu_c + mlngamma_c + mag_c) - logn - mlngamma
714  - mag ;
715 
716  // Test: is the enthalpy negative somewhere in the equatorial plane
717  // inside the star ? If yes, this means that the Keplerian velocity
718  // has been overstep.
719 
720  kepler = false ;
721  for (int l=0; l<nzet; l++) {
722  int imax = mg->get_nr(l) - 1 ;
723  if (l == l_b) imax-- ; // The surface point is skipped
724  for (int i=0; i<imax; i++) {
725  if ( ent()(l, 0, j_b, i) < 0. ) {
726  kepler = true ;
727  cout << "ent < 0 for l, i : " << l << " " << i
728  << " ent = " << ent()(l, 0, j_b, i) << endl ;
729  }
730  }
731  }
732 
733  if ( kepler ) {
734  cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ;
735  omega /= fact_omega ; // Omega is decreased
736  cout << "New rotation frequency : "
737  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
738  omega_trop_grand = true ;
739  }
740 
741  } // End of while ( kepler )
742 
743  if ( omega_trop_grand ) { // fact_omega is decreased for the
744  // next step
745  fact_omega = sqrt( fact_omega ) ;
746  cout << "**** New fact_omega : " << fact_omega << endl ;
747  }
748 
749  //----------------------------------------------------
750  // Adaptation of the mapping to the new enthalpy field
751  //----------------------------------------------------
752 
753  // Shall the adaptation be performed (cusp) ?
754  // ------------------------------------------
755 
756  double dent_eq = ent().dsdr()(l_b, k_b, j_b, i_b) ;
757  double dent_pole = ent().dsdr()(l_b, k_b, 0, i_b) ;
758  double rap_dent = fabs( dent_eq / dent_pole ) ;
759  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
760 
761  if ( rap_dent < thres_adapt ) {
762  adapt_flag = 0 ; // No adaptation of the mapping
763  cout << "******* FROZEN MAPPING *********" << endl ;
764  }
765  else{
766  adapt_flag = 1 ; // The adaptation of the mapping is to be
767  // performed
768  }
769 
770  mp_prev = mp_et ;
771 
772  mp.adapt(ent(), par_adapt) ;
773 
774  //----------------------------------------------------
775  // Computation of the enthalpy at the new grid points
776  //----------------------------------------------------
777 
778  mp_prev.homothetie(alpha_r) ;
779 
780  mp.reevaluate(&mp_prev, nzet+1, ent.set()) ;
781 
782  //----------------------------------------------------
783  // Equation of state
784  //----------------------------------------------------
785 
786  equation_of_state() ; // computes new values for nbar (n), ener (e)
787  // and press (p) from the new ent (H)
788 
789  //---------------------------------------------------------
790  // Matter source terms in the gravitational field equations
791  //---------------------------------------------------------
792 
793  //## Computation of tnphi and nphi from the Cartesian components
794  // of the shift for the test in hydro_euler():
795 
796  fait_nphi() ;
797 
798  hydro_euler() ; // computes new values for ener_euler (E),
799  // s_euler (S) and u_euler (U^i)
800 
801  if (relativistic) {
802 
803  //-------------------------------------------------------
804  // 2-D Poisson equation for tggg
805  //-------------------------------------------------------
806 
807  mp.poisson2d(source_tggg(), mp.cmp_zero(), par_poisson_tggg,
808  tggg.set()) ;
809 
810  //-------------------------------------------------------
811  // 2-D Poisson equation for dzeta
812  //-------------------------------------------------------
813 
814  mp.poisson2d(source_dzf(), source_dzq(), par_poisson_dzeta,
815  dzeta.set()) ;
816 
817  err_grv2 = lbda_grv2 - 1;
818  cout << "GRV2: " << err_grv2 << endl ;
819 
820  }
821  else {
822  err_grv2 = grv2() ;
823  }
824 
825 
826  //---------------------------------------
827  // Computation of the metric coefficients (except for N^phi)
828  //---------------------------------------
829 
830  // Relaxations on nu and dzeta :
831 
832  if (mer >= 10) {
833  logn = relax * logn + relax_prev * logn_prev ;
834 
835  dzeta = relax * dzeta + relax_prev * dzeta_prev ;
836  }
837 
838  // Update of the metric coefficients N, A, B and computation of K_ij :
839 
840  update_metric() ;
841 
842  //-----------------------
843  // Informations display
844  //-----------------------
845 
846  // partial_display(cout) ;
847  fichfreq << " " << omega / (2*M_PI) * f_unit ;
848  fichevol << " " << rap_dent ;
849  fichevol << " " << ray_pole() / ray_eq() ;
850  fichevol << " " << ent_c ;
851 
852  //-----------------------------------------
853  // Convergence towards a given baryon mass
854  //-----------------------------------------
855 
856  if (mer > mer_mass) {
857 
858  double xx ;
859  if (mbar_wanted > 0.) {
860  xx = mass_b() / mbar_wanted - 1. ;
861  cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx
862  << endl ;
863  }
864  else{
865  xx = mass_g() / fabs(mbar_wanted) - 1. ;
866  cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx
867  << endl ;
868  }
869  double xprog = ( mer > 2*mer_mass) ? 1. :
870  double(mer-mer_mass)/double(mer_mass) ;
871  xx *= xprog ;
872  double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
873  double fact = pow(ax, aexp_mass) ;
874  cout << " xprog, xx, ax, fact : " << xprog << " " <<
875  xx << " " << ax << " " << fact << endl ;
876 
877  if ( change_ent ) {
878  ent_c *= fact ;
879  }
880  else {
881  if (mer%4 == 0) omega *= fact ;
882  }
883  }
884 
885 
886  //------------------------------------------------------------
887  // Relative change in enthalpy with respect to previous step
888  //------------------------------------------------------------
889 
890  Tbl diff_ent_tbl = diffrel( ent(), ent_prev() ) ;
891  diff_ent = diff_ent_tbl(0) ;
892  for (int l=1; l<nzet; l++) {
893  diff_ent += diff_ent_tbl(l) ;
894  }
895  diff_ent /= nzet ;
896 
897  fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
898  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
899 
900  //------------------------------
901  // Recycling for the next step
902  //------------------------------
903 
904  ent_prev = ent ;
905  logn_prev = logn ;
906  dzeta_prev = dzeta ;
907 
908  fichconv << endl ;
909  fichfreq << endl ;
910  fichevol << endl ;
911  fichconv.flush() ;
912  fichfreq.flush() ;
913  fichevol.flush() ;
914 
915  } // End of main loop
916 
917  //=========================================================================
918  // End of iteration
919  //=========================================================================
920 
921  fichconv.close() ;
922  fichfreq.close() ;
923  fichevol.close() ;
924 
925 
926 }
927 }
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:116
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
virtual void MHD_comput()
Computes the electromagnetic part of the stress-energy tensor.
Cmp A_phi
-component of the electromagnetic potential 1-form divided by .
Definition: et_rot_mag.h:155
double a_j
Amplitude of the curent/charge function.
Definition: et_rot_mag.h:180
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,...
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
Definition: et_rot_mag.h:150
virtual double mass_g() const
Gravitational mass.
bool is_conduct() const
Tells if the star is made of conducting or isolating material.
Definition: et_rot_mag.h:241
virtual double grv2() const
Error on the virial identity GRV2.
Tenseur E_em
electromagnetic energy density in the Eulerian frame
Definition: et_rot_mag.h:161
Tenseur Jp_em
component of the electromagnetic momentum density 3-vector, as measured in the Eulerian frame.
Definition: et_rot_mag.h:167
double Q
In the case of a perfect conductor, the requated baryonic charge.
Definition: et_rot_mag.h:179
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 aexp_mass, Tbl &diff, const double Q0, const double a_j0, Cmp(*f_j)(const Cmp &x, const double), Cmp(*M_j)(const Cmp &x, const double))
Computes an equilibrium configuration.
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition: etoile.h:1625
Tenseur uuu
Norm of u_euler.
Definition: etoile.h:1518
double omega
Rotation angular velocity ([f_unit] )
Definition: etoile.h:1501
Tenseur & logn
Metric potential = logn_auto.
Definition: etoile.h:1521
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
Tenseur nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: etoile.h:1531
Tenseur khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition: etoile.h:1560
Tenseur tggg
Metric potential .
Definition: etoile.h:1537
Tenseur nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: etoile.h:1526
Cmp ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition: etoile.h:1608
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
virtual double mass_b() const
Baryon mass.
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
void update_metric()
Computes metric coefficients from known potentials.
Definition: et_rot_upmetr.C:69
Tenseur ak_car
Scalar .
Definition: etoile.h:1586
Tenseur & dzeta
Metric potential = beta_auto.
Definition: etoile.h:1534
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
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition: etoile_rot.C:781
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
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition: etoile.h:1567
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1507
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
Tenseur w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition: etoile.h:1550
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:566
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:474
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:471
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
Tenseur press
Fluid pressure.
Definition: etoile.h:461
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:437
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:465
Tenseur shift
Total shift vector.
Definition: etoile.h:512
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:468
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:457
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
double ray_pole() const
Coordinate radius at [r_unit].
Basic integer array class.
Definition: itbl.h:122
Radial mapping of rather general form.
Definition: map.h:2752
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:905
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
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.
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
const Cmp & cmp_zero() const
Returns the null Cmp defined on *this.
Definition: map.h:807
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
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition: param.C:453
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1142
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:522
Basic array class.
Definition: tbl.h:161
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition: tenseur.C:1548
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
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
Definition: tenseur_pde.C:118
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
void coef_i() const
Computes the physical value of *this.
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Tenseur 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...
Lorene prototypes.
Definition: app_hor.h:64
Standard electro-magnetic units.