LORENE
gravastar_equil.C
1 /*
2  * Method Gravastar::equilibrium
3  *
4  * (see file star_h.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2010 Frederic Vincent
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 
30 
31 
32 
33 
34 // Headers C
35 #include <cmath>
36 
37 // Headers Lorene
38 #include "gravastar.h"
39 #include "param.h"
40 #include "tenseur.h"
41 
42 #include "graphique.h"
43 #include "utilitaires.h"
44 #include "unites.h"
45 
46 namespace Lorene {
47 void Gravastar::equilibrium(double omega0, double fact_omega,
48  int nzadapt, const Tbl& ent_limit,
49  const Itbl& icontrol, const Tbl& control,
50  Tbl& diff, Param*) {
51 
52  // Fundamental constants and units
53  // -------------------------------
54 
55  using namespace Unites ;
56 
57  // For the display
58  // ---------------
59  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
60  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
61 
62  // Grid parameters
63  // ---------------
64 
65  const Mg3d* mg = mp.get_mg() ;
66  int nz = mg->get_nzone() ; // total number of domains
67  int nzm1 = nz - 1 ;
68 
69  // The following is required to initialize mp_prev as a Map_et:
70  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
71 
72  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
73  assert(mg->get_type_t() == SYM) ;
74  int l_b = nzet - 1 ;
75  int i_b = mg->get_nr(l_b) - 1 ;
76  int j_b = mg->get_nt(l_b) - 1 ;
77  int k_b = 0 ;
78 
79  // Value of the enthalpy defining the surface of the star
80  double ent_b = ent_limit(nzet-1) ;
81 
82  // Parameters to control the iteration
83  // -----------------------------------
84 
85  int mer_max = icontrol(0) ;
86  int mer_rot = icontrol(1) ;
87  int mer_change_omega = icontrol(2) ;
88  int mer_fix_omega = icontrol(3) ;
89  int mermax_poisson = icontrol(4) ;
90  int delta_mer_kep = icontrol(5) ;
91 
92  // Protections:
93  if (mer_change_omega < mer_rot) {
94  cout << "Gravastar::equilibrium: mer_change_omega < mer_rot !" << endl ;
95  cout << " mer_change_omega = " << mer_change_omega << endl ;
96  cout << " mer_rot = " << mer_rot << endl ;
97  abort() ;
98  }
99  if (mer_fix_omega < mer_change_omega) {
100  cout << "Gravastar::equilibrium: mer_fix_omega < mer_change_omega !"
101  << endl ;
102  cout << " mer_fix_omega = " << mer_fix_omega << endl ;
103  cout << " mer_change_omega = " << mer_change_omega << endl ;
104  abort() ;
105  }
106 
107  /*
108  // In order to converge to a given baryon mass, shall the central
109  // enthalpy be varied or Omega ?
110  bool change_ent = true ;
111  if (mer_mass < 0) {
112  change_ent = false ;
113  mer_mass = abs(mer_mass) ;
114  }
115  */
116 
117  double precis = control(0) ;
118  double omega_ini = control(1) ;
119  double relax = control(2) ;
120  double relax_prev = double(1) - relax ;
121  double relax_poisson = control(3) ;
122  double thres_adapt = control(4) ;
123  double precis_adapt = control(5) ;
124 
125 
126  // Error indicators
127  // ----------------
128 
129  diff.set_etat_qcq() ;
130  double& diff_ent = diff.set(0) ;
131  double& diff_nuf = diff.set(1) ;
132  double& diff_nuq = diff.set(2) ;
133 // double& diff_dzeta = diff.set(3) ;
134 // double& diff_ggg = diff.set(4) ;
135  double& diff_shift_x = diff.set(5) ;
136  double& diff_shift_y = diff.set(6) ;
137  // double& vit_triax = diff.set(7) ;
138 
139  // Parameters for the function Map_et::adapt
140  // -----------------------------------------
141 
142  Param par_adapt ;
143  int nitermax = 100 ;
144  int niter ;
145  int adapt_flag = 1 ; // 1 = performs the full computation,
146  // 0 = performs only the rescaling by
147  // the factor alpha_r
148  int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
149  // isosurfaces
150  double alpha_r ;
151  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
152 
153  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
154  // locate zeros by the secant method
155  par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
156  // to the isosurfaces of ent is to be
157  // performed
158  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
159  // the enthalpy isosurface
160  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
161  // 0 = performs only the rescaling by
162  // the factor alpha_r
163  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
164  // (theta_*, phi_*)
165  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
166  // (theta_*, phi_*)
167 
168  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
169  // the secant method
170 
171  par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
172  // the determination of zeros by
173  // the secant method
174  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
175 
176  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
177  // distances will be multiplied
178 
179  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
180  // to define the isosurfaces.
181 
182  // Parameters for the function Map_et::poisson for nuf
183  // ----------------------------------------------------
184 
185  double precis_poisson = 1.e-16 ;
186 
187  Cmp cssjm1_nuf(ssjm1_nuf) ;
188  Cmp cssjm1_nuq(ssjm1_nuq) ;
189  Cmp cssjm1_tggg(ssjm1_tggg) ;
190  Cmp cssjm1_khi(ssjm1_khi) ;
191  Tenseur cssjm1_wshift(mp, 1, CON, mp.get_bvect_cart() ) ;
192  cssjm1_wshift.set_etat_qcq() ;
193  for (int i=1; i<=3; i++) {
194  cssjm1_wshift.set(i-1) = ssjm1_wshift(i) ;
195  }
196 
197  Tenseur cshift(mp, 1, CON, mp.get_bvect_cart() ) ;
198  cshift.set_etat_qcq() ;
199  for (int i=1; i<=3; i++) {
200  cshift.set(i-1) = -beta(i) ;
201  }
202 
203  Tenseur cw_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
204  cw_shift.set_etat_qcq() ;
205  for (int i=1; i<=3; i++) {
206  cw_shift.set(i-1) = w_shift(i) ;
207  }
208 
209  Tenseur ckhi_shift(mp) ;
210  ckhi_shift.set_etat_qcq() ;
211  ckhi_shift.set() = khi_shift ;
212 
213  Param par_poisson_nuf ;
214  par_poisson_nuf.add_int(mermax_poisson, 0) ; // maximum number of iterations
215  par_poisson_nuf.add_double(relax_poisson, 0) ; // relaxation parameter
216  par_poisson_nuf.add_double(precis_poisson, 1) ; // required precision
217  par_poisson_nuf.add_int_mod(niter, 0) ; // number of iterations actually used
218  par_poisson_nuf.add_cmp_mod( cssjm1_nuf ) ;
219 
220  Param par_poisson_nuq ;
221  par_poisson_nuq.add_int(mermax_poisson, 0) ; // maximum number of iterations
222  par_poisson_nuq.add_double(relax_poisson, 0) ; // relaxation parameter
223  par_poisson_nuq.add_double(precis_poisson, 1) ; // required precision
224  par_poisson_nuq.add_int_mod(niter, 0) ; // number of iterations actually used
225  par_poisson_nuq.add_cmp_mod( cssjm1_nuq ) ;
226 
227  Param par_poisson_tggg ;
228  par_poisson_tggg.add_int(mermax_poisson, 0) ; // maximum number of iterations
229  par_poisson_tggg.add_double(relax_poisson, 0) ; // relaxation parameter
230  par_poisson_tggg.add_double(precis_poisson, 1) ; // required precision
231  par_poisson_tggg.add_int_mod(niter, 0) ; // number of iterations actually used
232  par_poisson_tggg.add_cmp_mod( cssjm1_tggg ) ;
233  double lambda_tggg ;
234  par_poisson_tggg.add_double_mod( lambda_tggg ) ;
235 
236  Param par_poisson_dzeta ;
237  double lbda_grv2 ;
238  par_poisson_dzeta.add_double_mod( lbda_grv2 ) ;
239 
240  // Parameters for the function Scalar::poisson_vect
241  // -------------------------------------------------
242 
243  Param par_poisson_vect ;
244 
245  par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of iterations
246  par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
247  par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
248  par_poisson_vect.add_cmp_mod( cssjm1_khi ) ;
249  par_poisson_vect.add_tenseur_mod( cssjm1_wshift ) ;
250  par_poisson_vect.add_int_mod(niter, 0) ;
251 
252 
253  // Initializations
254  // ---------------
255 
256  // Initial angular velocity
257  omega = 0 ;
258 
259  double accrois_omega = (omega0 - omega_ini) /
260  double(mer_fix_omega - mer_change_omega) ;
261 
262 
263 
264  update_metric() ; // update of the metric coefficients
265 
266  //des_profile(logn,0.,4.,0.,0.,"log(N) avant");
267 
268  equation_of_state() ; // update of the density, pressure, etc...
269  //des_profile(logn,0.,4.,0.,0.,"log(N) apres");
270 
271  hydro_euler() ; // update of the hydro quantities relative to the
272  // Eulerian observer
273 
274  // Quantities at the previous step :
275 
276  ent.annule_domain(0);//A VERIFIER: j'annule de force enth dans le coeur (en toute rigueur je prendrais enth=p_coeur, mais j'ai peur que enth<0 pose pb)
277 
278  Map_et mp_prev = mp_et ;
279  Scalar ent_prev = ent ;
280  Scalar logn_prev = logn ;
281  Scalar dzeta_prev = dzeta ;
282 
283  /*des_profile(ent,0.,10.,0.,0.,"enthalpy init");
284  des_profile(press,0.,10.,0.,0.,"press init");
285  des_profile(ener,0.,10.,0.,0.,"ener init");
286  des_profile(uuu,0.,10.,0.,0.,"Uinit"); */
287 
288  // Creation of uninitialized tensors:
289  Scalar source_nuf(mp) ; // source term in the equation for nuf
290  Scalar source_nuq(mp) ; // source term in the equation for nuq
291  Scalar source_dzf(mp) ; // matter source term in the eq. for dzeta
292  Scalar source_dzq(mp) ; // quadratic source term in the eq. for dzeta
293  Scalar source_tggg(mp) ; // source term in the eq. for tggg
294  Vector source_shift(mp, CON, mp.get_bvect_cart()) ;
295  // source term for shift
296  Scalar mlngamma(mp) ; // centrifugal potential
297 
298 
299  ofstream fichconv("convergence.d") ; // Output file for diff_ent
300  fichconv << "# diff_ent GRV2 max_triax vit_triax" << endl ;
301 
302  ofstream fichfreq("frequency.d") ; // Output file for omega
303  fichfreq << "# f [Hz]" << endl ;
304 
305  ofstream fichevol("evolution.d") ; // Output file for various quantities
306  fichevol <<
307  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_cr"
308  << endl ;
309 
310  diff_ent = 1 ;
311  double err_grv2 = 1 ;
312  // double max_triax_prev = 0 ; // Triaxial amplitude at previous step
313 
314  //=========================================================================
315  // Start of iteration
316  //=========================================================================
317 
318  for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {
319 
320  cout << "-----------------------------------------------" << endl ;
321  cout << "step: " << mer << endl ;
322  cout << "diff_ent = " << display_bold << diff_ent << display_normal
323  << endl ;
324  cout << "err_grv2 = " << err_grv2 << endl ;
325  fichconv << mer ;
326  fichfreq << mer ;
327  fichevol << mer ;
328 
329  if (mer >= mer_rot) {
330 
331  if (mer < mer_change_omega) {
332  omega = omega_ini ;
333  }
334  else {
335  if (mer <= mer_fix_omega) {
336  omega = omega_ini + accrois_omega *
337  (mer - mer_change_omega) ;
338  }
339  }
340 
341  }
342 
343  //-----------------------------------------------
344  // Sources of the Poisson equations
345  //-----------------------------------------------
346 
347  // Source for nu
348  // -------------
349  Scalar bet = log(bbb) ;
350  bet.std_spectral_base() ;
351 
352  Vector d_logn = logn.derive_cov( mp.flat_met_spher() ) ;
353  Vector d_bet = bet.derive_cov( mp.flat_met_spher() ) ;
354 
355  if (relativistic) {
356  source_nuf = qpig * a_car *( ener_euler + s_euler ) ;
357  //des_profile(source_nuf,0.,10.,0.,0.,"source_nuf");
358  //cout << "source_nuf= " << source_nuf << endl;
359 
360  source_nuq = ak_car - d_logn(1)*(d_logn(1)+d_bet(1))
361  - d_logn(2)*(d_logn(2)+d_bet(2))
362  - d_logn(3)*(d_logn(3)+d_bet(3)) ;
363  // cout << "source_nuq= " << source_nuq << endl;
364  }
365  else {
366  source_nuf = qpig * nbar ;
367 
368  source_nuq = 0 ;
369  }
370  source_nuf.std_spectral_base() ;
371  source_nuq.std_spectral_base() ;
372 
373  //TEST!
374  //ener_euler.std_spectral_base() ;
375 
376  /*des_profile(s_euler,0.,10.,0.,0.,"s_euler");
377  //cout << "compa 1= " << s_euler << endl;
378  //cout << "compa 2= " << ener_euler << endl;
379  des_profile(ener_euler,0.,10.,0.,0.,"ener_euler");
380  des_profile(press,0.,10.,0.,0.,"press");
381  des_profile(ener,0.,10.,0.,0.,"ener");
382  des_profile(source_nuf,0.,10.,0.,0.,"source_nuf");
383  des_profile(source_nuq,0.,10.,0.,0.,"source_nuq");*/
384 
385  //des_profile(source_nuf,0.,10.,0.,0.,"source_nuf");
386 
387  // Source for dzeta
388  // ----------------
389  source_dzf = 2 * qpig * a_car * (press + (ener_euler+press) * uuu*uuu ) ;
390  source_dzf.std_spectral_base() ;
391 
392  source_dzq = 1.5 * ak_car
393  - d_logn(1)*d_logn(1) - d_logn(2)*d_logn(2) - d_logn(3)*d_logn(3) ;
394  source_dzq.std_spectral_base() ;
395 
396  // Source for tggg
397  // ---------------
398 
399  source_tggg = 4 * qpig * nn * a_car * bbb * press ;
400  source_tggg.std_spectral_base() ;
401 
402  source_tggg.mult_rsint() ;
403 
404 
405  // Source for shift
406  // ----------------
407 
408  // Matter term:
409  source_shift = (-4*qpig) * nn * a_car * (ener_euler + press)
410  * u_euler ;
411 
412  // Quadratic terms:
413  Vector vtmp = 6 * bet.derive_con( mp.flat_met_spher() )
414  - 2 * logn.derive_con( mp.flat_met_spher() ) ;
415  vtmp.change_triad(mp.get_bvect_cart()) ;
416 
417  Vector squad = nn * contract(tkij, 1, vtmp, 0) ;
418 
419  source_shift = source_shift + squad.up(0, mp.flat_met_cart() ) ;
420 
421  //----------------------------------------------
422  // Resolution of the Poisson equation for nuf
423  //----------------------------------------------
424 
425  source_nuf.poisson(par_poisson_nuf, nuf) ;
426 
427  cout << "$$$???$$$ Test of the Poisson equation for nuf :" << endl ;
428  Tbl err = source_nuf.test_poisson(nuf, cout, true) ;
429  diff_nuf = err(0, 0) ;
430 
431  //---------------------------------------
432  // Triaxial perturbation of nuf
433  //---------------------------------------
434  /*
435  if (mer == mer_triax) {
436 
437  if ( mg->get_np(0) == 1 ) {
438  cout <<
439  "Gravastar::equilibrium: np must be stricly greater than 1"
440  << endl << " to set a triaxial perturbation !" << endl ;
441  abort() ;
442  }
443 
444  const Coord& phi = mp.phi ;
445  const Coord& sint = mp.sint ;
446  Scalar perturb(mp) ;
447  perturb = 1 + ampli_triax * sint*sint * cos(2*phi) ;
448  nuf = nuf * perturb ;
449 
450  nuf.std_spectral_base() ; // set the bases for spectral expansions
451  // to be the standard ones for a
452  // scalar field
453 
454  }
455 
456  // Monitoring of the triaxial perturbation
457  // ---------------------------------------
458 
459  const Valeur& va_nuf = nuf.get_spectral_va() ;
460  va_nuf.coef() ; // Computes the spectral coefficients
461  double max_triax = 0 ;
462 
463  if ( mg->get_np(0) > 1 ) {
464 
465  for (int l=0; l<nz; l++) { // loop on the domains
466  for (int j=0; j<mg->get_nt(l); j++) {
467  for (int i=0; i<mg->get_nr(l); i++) {
468 
469  // Coefficient of cos(2 phi) :
470  double xcos2p = (*(va_nuf.c_cf))(l, 2, j, i) ;
471 
472  // Coefficient of sin(2 phi) :
473  double xsin2p = (*(va_nuf.c_cf))(l, 3, j, i) ;
474 
475  double xx = sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ;
476 
477  max_triax = ( xx > max_triax ) ? xx : max_triax ;
478  }
479  }
480  }
481 
482  }
483 
484  cout << "Triaxial part of nuf : " << max_triax << endl ;
485  */
486  if (relativistic) {
487 
488  //----------------------------------------------
489  // Resolution of the Poisson equation for nuq
490  //----------------------------------------------
491 
492  source_nuq.poisson(par_poisson_nuq, nuq) ;
493 
494  cout << "Test of the Poisson equation for nuq :" << endl ;
495  err = source_nuq.test_poisson(nuq, cout, true) ;
496  diff_nuq = err(0, 0) ;
497 
498  //---------------------------------------------------------
499  // Resolution of the vector Poisson equation for the shift
500  //---------------------------------------------------------
501 
502 
503  for (int i=1; i<=3; i++) {
504  if(source_shift(i).get_etat() != ETATZERO) {
505  if(source_shift(i).dz_nonzero()) {
506  assert( source_shift(i).get_dzpuis() == 4 ) ;
507  }
508  else{
509  (source_shift.set(i)).set_dzpuis(4) ;
510  }
511  }
512  }
513 
514  double lambda_shift = double(1) / double(3) ;
515 
516  if ( mg->get_np(0) == 1 ) {
517  lambda_shift = 0 ;
518  }
519 
520  Tenseur csource_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
521  csource_shift.set_etat_qcq() ;
522  for (int i=1; i<=3; i++) {
523  csource_shift.set(i-1) = source_shift(i) ;
524  }
525  csource_shift.set(2).set_etat_zero() ; //## bizarre...
526 
527  csource_shift.poisson_vect(lambda_shift, par_poisson_vect,
528  cshift, cw_shift, ckhi_shift) ;
529 
530  for (int i=1; i<=3; i++) {
531  beta.set(i) = - cshift(i-1) ;
532  beta.set(i).set_dzpuis(0) ; //## bizarre...
533  w_shift.set(i) = cw_shift(i-1) ;
534  }
535  khi_shift = ckhi_shift() ;
536 
537  cout << "Test of the Poisson equation for shift_x :" << endl ;
538  err = source_shift(1).test_poisson(-beta(1), cout, true) ;
539  diff_shift_x = err(0, 0) ;
540 
541  cout << "Test of the Poisson equation for shift_y :" << endl ;
542  err = source_shift(2).test_poisson(-beta(2), cout, true) ;
543  diff_shift_y = err(0, 0) ;
544 
545  // Computation of tnphi and nphi from the Cartesian components
546  // of the shift
547  // -----------------------------------------------------------
548 
549  fait_nphi() ;
550 
551  }
552 
553  //-----------------------------------------
554  // Determination of the fluid velociy U
555  //-----------------------------------------
556 
557  if (mer > mer_fix_omega + delta_mer_kep) {
558 
559  omega *= fact_omega ; // Increase of the angular velocity if
560  } // fact_omega != 1
561 
562  bool omega_trop_grand = false ;
563  bool kepler = true ;
564 
565  while ( kepler ) {
566 
567  // Possible decrease of Omega to ensure a velocity < c
568 
569  bool superlum = true ;
570 
571  while ( superlum ) {
572 
573  // New fluid velocity U :
574 
575  Scalar tmp = omega - nphi ;
576  tmp.annule_domain(nzm1) ;
577  tmp.std_spectral_base() ;
578 
579  tmp.mult_rsint() ; // Multiplication by r sin(theta)
580 
581  uuu = bbb / nn * tmp ;
582 
583  if (uuu.get_etat() == ETATQCQ) {
584  // Same basis as (Omega -N^phi) r sin(theta) :
585  (uuu.set_spectral_va()).set_base( tmp.get_spectral_va().get_base() ) ;
586  }
587 
588  // Is the new velocity larger than c in the equatorial plane ?
589 
590  superlum = false ;
591 
592  for (int l=0; l<nzet; l++) {
593  for (int i=0; i<mg->get_nr(l); i++) {
594 
595  double u1 = uuu.val_grid_point(l, 0, j_b, i) ;
596  if (u1 >= 1.) { // superluminal velocity
597  superlum = true ;
598  cout << "U > c for l, i : " << l << " " << i
599  << " U = " << u1 << endl ;
600  }
601  }
602  }
603  if ( superlum ) {
604  cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
605  omega /= fact_omega ; // Decrease of Omega
606  cout << "New rotation frequency : "
607  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
608  omega_trop_grand = true ;
609  }
610  } // end of while ( superlum )
611 
612 
613  // New computation of U (which this time is not superluminal)
614  // as well as of gam_euler, ener_euler, etc...
615  // -----------------------------------
616 
617  hydro_euler() ;
618  //des_profile(uuu,0.,10.,0.,0.,"U");
619 
620  //------------------------------------------------------
621  // First integral of motion
622  //------------------------------------------------------
623 
624  // Centrifugal potential :
625  if (relativistic) {
626  mlngamma = - log( gam_euler ) ;
627  }
628  else {
629  mlngamma = - 0.5 * uuu*uuu ;
630  }
631 
632  // Equatorial values of various potentials :
633  double nuf_b = nuf.val_grid_point(l_b, k_b, j_b, i_b) ;
634  double nuq_b = nuq.val_grid_point(l_b, k_b, j_b, i_b) ;
635  double mlngamma_b = mlngamma.val_grid_point(l_b, k_b, j_b, i_b) ;
636 
637  /* //Central values of various potentials :
638  double nuf_c = nuf.val_grid_point(0,0,0,0) ;
639  double nuq_c = nuq.val_grid_point(0,0,0,0) ;
640  double mlngamma_c = 0 ;*/
641 
642  //des_profile(nuf,0.,10.,1.5708,0.,"nuf");
643  des_profile(nuf,0.,10.,1.5708,0.,"nuf (debut step suivant)");
644 
645  //des_profile(nuq,0.,10.,0.,0.,"nuq");
646  //des_profile(logn,0.,10.,0.,0.,"nu");
647 
648  //Potentials at equatorial point of inner crust boundary
649  //***NB: a voir, valuers de l et i
650  //int l_cr = 0 ; //no : I want to be inside the crust, at r=rcrust^{+}
651  int l_cr = 1 ;
652  //int i_cr = mg->get_nr(l_cr) - 1 ;
653  int i_cr = 0 ;
654  int j_cr = mg->get_nt(l_cr) - 1 ;
655  int k_cr = 0 ;
656  double nuf_cr = nuf.val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
657  double nuq_cr = nuq.val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
658  double mlngamma_cr = mlngamma.val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
659 
660 
661  // Scale factor to ensure that the enthalpy is equal to ent_b at
662  // the equator
663  /*double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
664  + nuq_c - nuq_b) / ( nuf_b - nuf_c ) ;*/
665  /*int l_crp = 1 ; //r = r_{inner crust}^{+}
666  int i_crp = 0 ;
667  int j_crp = mg->get_nt(l_crp) - 1 ;
668  int k_crp = 0 ; */
669 
670  //double ent_cr=ent_limit(0);//enthalpy at inner crust boundary
671  //***NB: doit etre identique... :
672  double ent_cr=ent.val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
673  double alpha_r2 = ( ent_cr - ent_b + mlngamma_cr - mlngamma_b
674  + nuq_cr - nuq_b) / ( nuf_b - nuf_cr ) ;
675  alpha_r = sqrt(alpha_r2) ;
676  cout << "ent_cr,ent_b,mlngamma_cr,mlngamma_b,nuq_cr,nuq_b " << ent_cr << " " << ent_b << " " << mlngamma_cr << " " << mlngamma_b << " " << nuq_cr << " " << nuq_b << endl;
677  cout << "nuf_b, nuf_cr= " << nuf_b << " " << nuf_cr << endl;
678  cout << "num= " << ent_cr - ent_b + mlngamma_cr - mlngamma_b
679  + nuq_cr - nuq_b << endl;
680  cout << "deno= " << nuf_b - nuf_cr << endl;
681  cout << "alpha_r = " << alpha_r << endl ;
682  if (alpha_r != alpha_r){
683  cout << "alpha_r nan!" << endl;
684  abort();
685  }
686 
687  // Readjustment of nu :
688  // -------------------
689  //double nu_cr1 = logn.val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
690  logn = alpha_r2 * nuf + nuq ;
691  double nu_cr = logn.val_grid_point(l_cr,k_cr,j_cr,i_cr) ;
692 
693  // First integral --> enthalpy in all space
694  //-----------------
695  cout << "ent_cr, ent_crbis, nu_cr= " << ent_cr << " " << ent.val_grid_point(l_cr,k_cr,j_cr,i_cr) << " " << nu_cr << endl;
696  //des_profile(ent,0.,10.,0.,0.,"enthalpy avant cst");
697 
698  ent = (ent_cr + nu_cr + mlngamma_cr) - logn - mlngamma ;
699 
700  cout << "new ent_cr= " << ent.val_grid_point(l_cr,k_cr,j_cr,i_cr) << endl;
701 
702  //A VERIFIER: j'annule de force enth dans le coeur (en toute rigueur je prendrais enth=p_coeur, mais j'ai peur que enth<0 pose pb) et a l'ext de l'etoile
703 
704  ent.annule_domain(0); //***NB: si je le mets, des oscillations apparaissent dans ce domaine apres le remapping
705 
706  //ent.annule(nzet,nz-1); //***NB: si je le mets, erreur division par zero dans le remapping
707 
708  //des_profile(logn,0.,10.,0.,0.,"nu");
709  des_profile(ent,0.,10.,0.,0.,"enthalpy apres cst");
710  des_profile(ent,1.001,1.2,0.,0.,"enthalpy apres zoom");
711 
712  // Test: is the enthalpy negative somewhere in the equatorial plane
713  // inside the star ? If yes, this means that the Keplerian velocity
714  // has been overstep.
715 
716  kepler = false ;
717  for (int l=0; l<nzet; l++) {
718  int imax = mg->get_nr(l) - 1 ;
719  if (l == l_b) imax-- ; // The surface point is skipped
720  for (int i=0; i<imax; i++) {
721  if ( ent.val_grid_point(l, 0, j_b, i) < 0. ) {
722  //kepler = true ;
723  cout << "ent < 0 for l, i : " << l << " " << i
724  << " ent = " << ent.val_grid_point(l, 0, j_b, i) << endl ;
725  }
726  }
727  }
728 
729  if ( kepler ) {
730  cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ;
731  omega /= fact_omega ; // Omega is decreased
732  cout << "New rotation frequency : "
733  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
734  omega_trop_grand = true ;
735  }
736 
737  } // End of while ( kepler )
738 
739  if ( omega_trop_grand ) { // fact_omega is decreased for the
740  // next step
741  fact_omega = sqrt( fact_omega ) ;
742  cout << "**** New fact_omega : " << fact_omega << endl ;
743  }
744 
745  //----------------------------------------------------
746  // Adaptation of the mapping to the new enthalpy field
747  //----------------------------------------------------
748 
749  // Shall the adaptation be performed (cusp) ?
750  // ------------------------------------------
751 
752  double dent_eq = ent.dsdr().val_grid_point(l_b, k_b, j_b, i_b) ;
753  double dent_pole = ent.dsdr().val_grid_point(l_b, k_b, 0, i_b) ;
754  double rap_dent = fabs( dent_eq / dent_pole ) ;
755  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
756 
757  if ( rap_dent < thres_adapt ) {
758  adapt_flag = 0 ; // No adaptation of the mapping
759  cout << "******* FROZEN MAPPING *********" << endl ;
760  }
761  else{
762  adapt_flag = 1 ; // The adaptation of the mapping is to be
763  // performed
764  }
765 
766  mp_prev = mp_et ;
767 
768  //cout << "ent crust= " << ent << endl;
769 
770  //des_profile(ent,0.,10.,0.,0.,"enthalpy");
771 
772  //des_profile(logn,0.9,1.1,0.,0.,"log(N)");
773 
774  Cmp cent(ent) ;
775  mp.adapt(cent, par_adapt) ;
776 
777  //----------------------------------------------------
778  // Computation of the enthalpy at the new grid points
779  //----------------------------------------------------
780  //***NB***: pourquoi lorsque j'impose enth(domaine 0)=0 est-ce que le rayon du domaine 0 augmente quand meme? cf Remapping.png avec les courbes d'enth just apres le calcule de cst du mvt et juste apres remapping -> pourquoi la limite du domaine 0 est-elle la ou elle est dans la fig de droite??
781 
782  mp_prev.homothetie(alpha_r) ;
783 
784  // mp.reevaluate(&mp_prev, nzet+1, cent) ; //***NB: annule le champ (cent) pour les domaines [nzet+1,nz] -> pourquoi mettre nzet+1 et pas nzet???
785  mp.reevaluate(&mp_prev, nzet, cent) ;
786  ent = cent ;
787  des_profile(ent,0.,10.,0.,0.,"enthalpy apres mapping (used for eos)");
788 
789  //----------------------------------------------------
790  // Equation of state
791  //----------------------------------------------------
792 
793  equation_of_state() ; // computes new values for nbar (n), ener (e)
794  // and press (p) from the new ent (H)
795 
796  //---------------------------------------------------------
797  // Matter source terms in the gravitational field equations
798  //---------------------------------------------------------
799 
800  //## Computation of tnphi and nphi from the Cartesian components
801  // of the shift for the test in hydro_euler():
802 
803  fait_nphi() ;
804 
805  hydro_euler() ; // computes new values for ener_euler (E),
806  // s_euler (S) and u_euler (U^i)
807 
808  /*des_profile(press,0.,10.,0.,0.,"press");
809  des_profile(ener,0.,10.,0.,0.,"ener");
810  des_profile(uuu,0.,10.,0.,0.,"U (last fig for iteration)"); */
811 
812  if (relativistic) {
813 
814  //-------------------------------------------------------
815  // 2-D Poisson equation for tggg
816  //-------------------------------------------------------
817 
818  Cmp csource_tggg(source_tggg) ;
819  Cmp ctggg(tggg) ;
820  mp.poisson2d(csource_tggg, mp.cmp_zero(), par_poisson_tggg,
821  ctggg) ;
822  tggg = ctggg ;
823 
824 
825  //-------------------------------------------------------
826  // 2-D Poisson equation for dzeta
827  //-------------------------------------------------------
828 
829  Cmp csource_dzf(source_dzf) ;
830  Cmp csource_dzq(source_dzq) ;
831  Cmp cdzeta(dzeta) ;
832  mp.poisson2d(csource_dzf, csource_dzq, par_poisson_dzeta,
833  cdzeta) ;
834  dzeta = cdzeta ;
835 
836  err_grv2 = lbda_grv2 - 1;
837  cout << "GRV2: " << err_grv2 << endl ;
838 
839  }
840  else {
841  err_grv2 = grv2() ;
842  }
843 
844 
845  //---------------------------------------
846  // Computation of the metric coefficients (except for N^phi)
847  //---------------------------------------
848 
849  // Relaxations on nu and dzeta :
850 
851  if (mer >= 10) {
852  logn = relax * logn + relax_prev * logn_prev ;
853 
854  dzeta = relax * dzeta + relax_prev * dzeta_prev ;
855  }
856 
857  // Update of the metric coefficients N, A, B and computation of K_ij :
858 
859  update_metric() ;
860 
861  //-----------------------
862  // Informations display
863  //-----------------------
864 
865  partial_display(cout) ;
866  fichfreq << " " << omega / (2*M_PI) * f_unit ;
867  fichevol << " " << rap_dent ;
868  fichevol << " " << ray_pole() / ray_eq() ;
869  // fichevol << " " << ent_c ;
870  //fichevol << " " << ent_cr ; fait une erreur??
871 
872  //-----------------------------------------
873  // Convergence towards a given baryon mass
874  //-----------------------------------------
875  /*
876  if (mer > mer_mass) {
877 
878  double xx ;
879  if (mbar_wanted > 0.) {
880  xx = mass_b() / mbar_wanted - 1. ;
881  cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx
882  << endl ;
883  }
884  else{
885  xx = mass_g() / fabs(mbar_wanted) - 1. ;
886  cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx
887  << endl ;
888  }
889  double xprog = ( mer > 2*mer_mass) ? 1. :
890  double(mer-mer_mass)/double(mer_mass) ;
891  xx *= xprog ;
892  double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
893  double fact = pow(ax, aexp_mass) ;
894  cout << " xprog, xx, ax, fact : " << xprog << " " <<
895  xx << " " << ax << " " << fact << endl ;
896 
897  if ( change_ent ) {
898  ent_c *= fact ;
899  }
900  else {
901  if (mer%4 == 0) omega *= fact ;
902  }
903  }
904  */
905 
906  //------------------------------------------------------------
907  // Relative change in enthalpy with respect to previous step
908  //------------------------------------------------------------
909 
910  Tbl diff_ent_tbl = diffrel( ent, ent_prev ) ;
911  diff_ent = diff_ent_tbl(0) ;
912  for (int l=1; l<nzet; l++) {
913  diff_ent += diff_ent_tbl(l) ;
914  }
915  diff_ent /= nzet ;
916 
917  fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
918  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
919  /* fichconv << " " << log10( fabs(max_triax) + 1.e-16 ) ;
920 
921  vit_triax = 0 ;
922  if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
923  vit_triax = (max_triax - max_triax_prev) / max_triax_prev ;
924  }
925 
926  fichconv << " " << vit_triax ;*/
927 
928  //------------------------------
929  // Recycling for the next step
930  //------------------------------
931 
932  ent_prev = ent ;
933  logn_prev = logn ;
934  dzeta_prev = dzeta ;
935  // max_triax_prev = max_triax ;
936 
937  fichconv << endl ;
938  fichfreq << endl ;
939  fichevol << endl ;
940  fichconv.flush() ;
941  fichfreq.flush() ;
942  fichevol.flush() ;
943 
944  } // End of main loop
945 
946  //=========================================================================
947  // End of iteration
948  //=========================================================================
949 
950  ssjm1_nuf = cssjm1_nuf ;
951  ssjm1_nuq = cssjm1_nuq ;
952  ssjm1_tggg = cssjm1_tggg ;
953  ssjm1_khi = cssjm1_khi ;
954  for (int i=1; i<=3; i++) {
955  ssjm1_wshift.set(i) = cssjm1_wshift(i-1) ;
956  }
957 
958  fichconv.close() ;
959  fichfreq.close() ;
960  fichevol.close() ;
961 
962 }
963 }
void equation_of_state()
Allows to computes the proper baryon and energy density, as well as pressure from the enthalpy...
Definition: gravastar.C:67
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
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
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for ...
Definition: star_rot.h:237
Scalar dzeta
Metric potential .
Definition: star_rot.h:146
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Scalar a_car
Square of the metric factor A.
Definition: star_rot.h:116
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
Map & mp
Mapping associated with the star.
Definition: star.h:180
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: star_rot.h:138
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
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: star_rot.h:172
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: star_rot.h:106
Scalar bbb
Metric factor B.
Definition: star_rot.h:119
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
virtual double grv2() const
Error on the virial identity GRV2.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
void des_profile(const float *uutab, int nx, float xmin, float xmax, const char *nomx, const char *nomy, const char *title, const char *device=0x0, int nbound=0, float *xbound=0x0, bool logscale=false)
Basic routine for drawing a single profile with uniform x sampling.
Definition: des_profile.C:97
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
Basic integer array class.
Definition: itbl.h:122
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
void update_metric()
Computes metric coefficients from known potentials.
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: star_rot.h:143
double omega
Rotation angular velocity ([f_unit] )
Definition: star_rot.h:113
Scalar nbar
Baryon density in the fluid frame.
Definition: star.h:192
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Scalar ent
Log-enthalpy.
Definition: star.h:190
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Vector beta
Shift vector.
Definition: star.h:228
Tbl test_poisson(const Scalar &uu, ostream &ostr, bool detail=false) const
Checks if a Poisson equation with *this as a source has been correctly solved.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:334
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg ...
Definition: star_rot.h:220
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
void add_double_mod(double &x, int position=0)
Adds the address of a new modifiable double to the list.
Definition: param.C:456
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar nphi
Metric coefficient .
Definition: star_rot.h:125
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
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.
Scalar press
Fluid pressure.
Definition: star.h:194
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition: star_rot.h:210
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
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:928
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Sym_tensor tkij
Tensor related to the extrinsic curvature tensor by .
Definition: star_rot.h:179
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:281
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Definition: star_rot.C:613
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
Scalar nn
Lapse function N .
Definition: star.h:225
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: star_rot.h:228
Scalar uuu
Norm of u_euler.
Definition: star_rot.h:133
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
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:111
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition: star_rot.h:204
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
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
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
Scalar tggg
Metric potential .
Definition: star_rot.h:149
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift ...
Definition: star_rot.C:779
Scalar ak_car
Scalar .
Definition: star_rot.h:198
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
void equilibrium(double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, Tbl &diff, Param *=0x0)
Computes an equilibrium configuration.
Vector w_shift
Vector used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: star_rot.h:162
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:324
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