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