LORENE
boson_star_equil.C
1 /*
2  * Method Boson_star::equilibrium
3  *
4  * (see file boson_star.h for documentation).
5  */
6 
7 /*
8  * Copyright (c) 2012 Claire Some, Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 /*
32  * $Id: boson_star_equil.C,v 1.7 2016/12/05 16:17:49 j_novak Exp $
33  * $Log: boson_star_equil.C,v $
34  * Revision 1.7 2016/12/05 16:17:49 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.6 2014/10/13 08:52:49 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.5 2014/10/06 15:13:04 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.4 2013/04/03 12:10:13 e_gourgoulhon
44  * Added member kk to Compobj; suppressed tkij
45  *
46  * Revision 1.3 2012/12/03 15:27:30 c_some
47  * Small changes
48  *
49  * Revision 1.2 2012/11/23 15:43:05 c_some
50  * Small changes
51  *
52  * Revision 1.1 2012/11/22 16:04:12 c_some
53  * New class Boson_star
54  *
55  *
56  * $Header: /cvsroot/Lorene/C++/Source/Compobj/boson_star_equil.C,v 1.7 2016/12/05 16:17:49 j_novak Exp $
57  *
58  */
59 
60 // Headers C
61 #include <cmath>
62 
63 // Headers Lorene
64 #include "boson_star.h"
65 #include "param.h"
66 #include "tenseur.h"
67 
68 #include "graphique.h"
69 #include "utilitaires.h"
70 #include "unites.h"
71 
72 namespace Lorene {
73 void Boson_star::equilibrium(double, double,
74  int nzadapt, const Tbl& phi_limit, const Itbl& icontrol,
75  const Tbl& control, Tbl& diff, Param*) {
76 
77  // Fundamental constants and units
78  // -------------------------------
79 
80  using namespace Unites ;
81 
82  // For the display
83  // ---------------
84  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
85  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
86 
87  // Grid parameters
88  // ---------------
89 
90  const Mg3d* mg = mp.get_mg() ;
91  int nz = mg->get_nzone() ; // total number of domains
92 
93  // The following is required to initialize mp_prev as a Map_et:
94  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
95 
96  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
97 
98  int nzet = nzadapt ; //## to be checked
99 
100  assert(mg->get_type_t() == SYM) ;
101  int l_b = nzet - 1 ;
102  int j_b = mg->get_nt(l_b) - 1 ;
103  int k_b = 0 ;
104 
105  // Value of the enthalpy defining the surface of the star
106  // double ent_b = phi_limit(nzet-1) ;
107 
108  // Parameters to control the iteration
109  // -----------------------------------
110 
111  int mer_max = icontrol(0) ;
112 // int mer_rot = icontrol(1) ;
113 // int mer_change_omega = icontrol(2) ;
114 // int mer_fix_omega = icontrol(3) ;
115 //## int mer_mass = icontrol(4) ;
116  int mermax_poisson = icontrol(5) ;
117  int mer_triax = icontrol(6) ;
118 //## int delta_mer_kep = icontrol(7) ;
119 
120 
121  double precis = control(0) ;
122 // double omega_ini = control(1) ;
123  double relax = control(2) ;
124  double relax_prev = double(1) - relax ;
125  double relax_poisson = control(3) ;
126 // double thres_adapt = control(4) ;
127  double ampli_triax = control(5) ;
128  double precis_adapt = control(6) ;
129 
130 
131  // Error indicators
132  // ----------------
133 
134  diff.set_etat_qcq() ;
135  double& diff_phi = diff.set(0) ;
136  double& diff_nuf = diff.set(1) ;
137  double& diff_nuq = diff.set(2) ;
138 // double& diff_dzeta = diff.set(3) ;
139 // double& diff_ggg = diff.set(4) ;
140  double& diff_shift_x = diff.set(5) ;
141  double& diff_shift_y = diff.set(6) ;
142  double& vit_triax = diff.set(7) ;
143 
144  // Parameters for the function Map_et::adapt
145  // -----------------------------------------
146 
147  Param par_adapt ;
148  int nitermax = 100 ;
149  int niter ;
150  int adapt_flag = 1 ; // 1 = performs the full computation,
151  // 0 = performs only the rescaling by
152  // the factor alpha_r
153  int nz_search = nzet + 1 ; // Number of domains for searching the enthalpy
154  // isosurfaces
155  double alpha_r ;
156  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
157 
158  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
159  // locate zeros by the secant method
160  par_adapt.add_int(nzadapt, 1) ; // number of domains where the adjustment
161  // to the isosurfaces of ent is to be
162  // performed
163  par_adapt.add_int(nz_search, 2) ; // number of domains to search for
164  // the enthalpy isosurface
165  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
166  // 0 = performs only the rescaling by
167  // the factor alpha_r
168  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
169  // (theta_*, phi_*)
170  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
171  // (theta_*, phi_*)
172 
173  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
174  // the secant method
175 
176  par_adapt.add_double(precis_adapt, 0) ; // required absolute precision in
177  // the determination of zeros by
178  // the secant method
179  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping, 0 = contracting mapping
180 
181  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
182  // distances will be multiplied
183 
184  par_adapt.add_tbl(phi_limit, 0) ; // array of values of the field Phi
185  // to define the isosurfaces.
186 
187  // Parameters for the function Map_et::poisson for nuf
188  // ----------------------------------------------------
189 
190  double precis_poisson = 1.e-16 ;
191 
192  // Preparations
193  // ------------
194 
195  // Cartesian components of the shift vector are required
197 
198 
199  Cmp cssjm1_nuf(ssjm1_nuf) ;
200  Cmp cssjm1_nuq(ssjm1_nuq) ;
201  Cmp cssjm1_tggg(ssjm1_tggg) ;
202  Cmp cssjm1_khi(ssjm1_khi) ;
203  Tenseur cssjm1_wshift(mp, 1, CON, mp.get_bvect_cart() ) ;
204  cssjm1_wshift.set_etat_qcq() ;
205  for (int i=1; i<=3; i++) {
206  cssjm1_wshift.set(i-1) = ssjm1_wshift(i) ;
207  }
208 
209  Tenseur cshift(mp, 1, CON, mp.get_bvect_cart() ) ;
210  cshift.set_etat_qcq() ;
211  for (int i=1; i<=3; i++) {
212  cshift.set(i-1) = -beta(i) ;
213  }
214 
215  Tenseur cw_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
216  cw_shift.set_etat_qcq() ;
217  for (int i=1; i<=3; i++) {
218  cw_shift.set(i-1) = w_shift(i) ;
219  }
220 
221  Tenseur ckhi_shift(mp) ;
222  ckhi_shift.set_etat_qcq() ;
223  ckhi_shift.set() = khi_shift ;
224 
225  Param par_poisson_nuf ;
226  par_poisson_nuf.add_int(mermax_poisson, 0) ; // maximum number of iterations
227  par_poisson_nuf.add_double(relax_poisson, 0) ; // relaxation parameter
228  par_poisson_nuf.add_double(precis_poisson, 1) ; // required precision
229  par_poisson_nuf.add_int_mod(niter, 0) ; // number of iterations actually used
230  par_poisson_nuf.add_cmp_mod( cssjm1_nuf ) ;
231 
232  Param par_poisson_nuq ;
233  par_poisson_nuq.add_int(mermax_poisson, 0) ; // maximum number of iterations
234  par_poisson_nuq.add_double(relax_poisson, 0) ; // relaxation parameter
235  par_poisson_nuq.add_double(precis_poisson, 1) ; // required precision
236  par_poisson_nuq.add_int_mod(niter, 0) ; // number of iterations actually used
237  par_poisson_nuq.add_cmp_mod( cssjm1_nuq ) ;
238 
239  Param par_poisson_tggg ;
240  par_poisson_tggg.add_int(mermax_poisson, 0) ; // maximum number of iterations
241  par_poisson_tggg.add_double(relax_poisson, 0) ; // relaxation parameter
242  par_poisson_tggg.add_double(precis_poisson, 1) ; // required precision
243  par_poisson_tggg.add_int_mod(niter, 0) ; // number of iterations actually used
244  par_poisson_tggg.add_cmp_mod( cssjm1_tggg ) ;
245  double lambda_tggg ;
246  par_poisson_tggg.add_double_mod( lambda_tggg ) ;
247 
248  Param par_poisson_dzeta ;
249  double lbda_grv2 ;
250  par_poisson_dzeta.add_double_mod( lbda_grv2 ) ;
251 
252  // Parameters for the function Scalar::poisson_vect
253  // -------------------------------------------------
254 
255  Param par_poisson_vect ;
256 
257  par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of iterations
258  par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
259  par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
260  par_poisson_vect.add_cmp_mod( cssjm1_khi ) ;
261  par_poisson_vect.add_tenseur_mod( cssjm1_wshift ) ;
262  par_poisson_vect.add_int_mod(niter, 0) ;
263 
264 
265  // Initializations
266  // ---------------
267 
268  //## Spherical components of the shift vector are restored
270  update_metric() ;
271  //## Back to Cartesian components
273 
274 
275  // Quantities at the previous step :
276  Map_et mp_prev = mp_et ;
277  Scalar rphi_prev = rphi ;
278  Scalar logn_prev = logn ;
279  Scalar dzeta_prev = dzeta ;
280 
281  // Creation of uninitialized tensors:
282  Scalar source_nuf(mp) ; // source term in the equation for nuf
283  Scalar source_nuq(mp) ; // source term in the equation for nuq
284  Scalar source_dzf(mp) ; // matter source term in the eq. for dzeta
285  Scalar source_dzq(mp) ; // quadratic source term in the eq. for dzeta
286  Scalar source_tggg(mp) ; // source term in the eq. for tggg
287  Vector source_shift(mp, CON, mp.get_bvect_cart()) ;
288  // source term for shift
289 
290 
291  ofstream fichconv("convergence.d") ; // Output file for diff_phi
292  fichconv << "# diff_phi GRV2 max_triax vit_triax" << endl ;
293 
294 
295  ofstream fichevol("evolution.d") ; // Output file for various quantities
296  fichevol <<
297  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq rphi_c"
298  << endl ;
299 
300  diff_phi = 1 ;
301  double err_grv2 = 1 ;
302  double max_triax_prev = 0 ; // Triaxial amplitude at previous step
303 
304  //=========================================================================
305  // Start of iteration
306  //=========================================================================
307 
308  for(int mer=0 ; (diff_phi > precis) && (mer<mer_max) ; mer++ ) {
309 
310  cout << "-----------------------------------------------" << endl ;
311  cout << "step: " << mer << endl ;
312  cout << "diff_phi = " << display_bold << diff_phi << display_normal
313  << endl ;
314  cout << "err_grv2 = " << err_grv2 << endl ;
315  fichconv << mer ;
316  fichevol << mer ;
317 
318 
319  //-----------------------------------------------
320  // Sources of the Poisson equations
321  //-----------------------------------------------
322 
323  // Source for nu
324  // -------------
325  Scalar bet = log(bbb) ;
326  bet.std_spectral_base() ;
327 
328  Vector d_logn = logn.derive_cov( mp.flat_met_spher() ) ;
329  Vector d_bet = bet.derive_cov( mp.flat_met_spher() ) ;
330 
331  Scalar s_euler = stress_euler.trace(gamma) ;
332  source_nuf = qpig * a_car *( ener_euler + s_euler ) ;
333 
334  source_nuq = ak_car - d_logn(1)*(d_logn(1)+d_bet(1))
335  - d_logn(2)*(d_logn(2)+d_bet(2))
336  - d_logn(3)*(d_logn(3)+d_bet(3)) ;
337 
338  source_nuf.std_spectral_base() ;
339  source_nuq.std_spectral_base() ;
340 
341  // Source for dzeta
342  // ----------------
343  source_dzf = 2 * qpig * a_car * b_car * stress_euler(3,3) ;
344  source_dzf.std_spectral_base() ;
345 
346  source_dzq = 1.5 * ak_car
347  - d_logn(1)*d_logn(1) - d_logn(2)*d_logn(2) - d_logn(3)*d_logn(3) ;
348  source_dzq.std_spectral_base() ;
349 
350  // Source for tggg
351  // ---------------
352 
353  source_tggg = 2 * qpig * nn * a_car * bbb * (s_euler - b_car * stress_euler(3,3)) ;
354  source_tggg.std_spectral_base() ;
355 
356  source_tggg.mult_rsint() ;
357 
358 
359  // Source for shift
360  // ----------------
361 
362  // Matter term:
363  Vector mom_euler_cart = mom_euler ;
364  mom_euler_cart.change_triad(mp.get_bvect_cart()) ;
365  source_shift = (-4*qpig) * nn * a_car * mom_euler_cart ;
366 
367  // Quadratic terms:
368  Vector vtmp = 6 * bet.derive_con( mp.flat_met_spher() )
369  - 2 * logn.derive_con( mp.flat_met_spher() ) ;
370 
371  Vector squad = nn * contract(kk, 1, vtmp, 0) / b_car ;
372  squad.change_triad(mp.get_bvect_cart()) ;
373 
374  source_shift = source_shift + squad.up(0, mp.flat_met_cart() ) ;
375 
376  //----------------------------------------------
377  // Resolution of the Poisson equation for nuf
378  //----------------------------------------------
379 
380  source_nuf.poisson(par_poisson_nuf, nuf) ;
381 
382  cout << "Test of the Poisson equation for nuf :" << endl ;
383  Tbl err = source_nuf.test_poisson(nuf, cout, true) ;
384  diff_nuf = err(0, 0) ;
385 
386  //---------------------------------------
387  // Triaxial perturbation of nuf
388  //---------------------------------------
389 
390  if (mer == mer_triax) {
391 
392  if ( mg->get_np(0) == 1 ) {
393  cout <<
394  "Boson_star::equilibrium: np must be stricly greater than 1"
395  << endl << " to set a triaxial perturbation !" << endl ;
396  abort() ;
397  }
398 
399  const Coord& phi = mp.phi ;
400  const Coord& sint = mp.sint ;
401  Scalar perturb(mp) ;
402  perturb = 1 + ampli_triax * sint*sint * cos(2*phi) ;
403  nuf = nuf * perturb ;
404 
405  nuf.std_spectral_base() ; // set the bases for spectral expansions
406  // to be the standard ones for a
407  // scalar field
408 
409  }
410 
411  // Monitoring of the triaxial perturbation
412  // ---------------------------------------
413 
414  const Valeur& va_nuf = nuf.get_spectral_va() ;
415  va_nuf.coef() ; // Computes the spectral coefficients
416  double max_triax = 0 ;
417 
418  if ( mg->get_np(0) > 1 ) {
419 
420  for (int l=0; l<nz; l++) { // loop on the domains
421  for (int j=0; j<mg->get_nt(l); j++) {
422  for (int i=0; i<mg->get_nr(l); i++) {
423 
424  // Coefficient of cos(2 phi) :
425  double xcos2p = (*(va_nuf.c_cf))(l, 2, j, i) ;
426 
427  // Coefficient of sin(2 phi) :
428  double xsin2p = (*(va_nuf.c_cf))(l, 3, j, i) ;
429 
430  double xx = sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ;
431 
432  max_triax = ( xx > max_triax ) ? xx : max_triax ;
433  }
434  }
435  }
436 
437  }
438 
439  cout << "Triaxial part of nuf : " << max_triax << endl ;
440 
441  //----------------------------------------------
442  // Resolution of the Poisson equation for nuq
443  //----------------------------------------------
444 
445  source_nuq.poisson(par_poisson_nuq, nuq) ;
446 
447  cout << "Test of the Poisson equation for nuq :" << endl ;
448  err = source_nuq.test_poisson(nuq, cout, true) ;
449  diff_nuq = err(0, 0) ;
450 
451  //---------------------------------------------------------
452  // Resolution of the vector Poisson equation for the shift
453  //---------------------------------------------------------
454 
455 
456  for (int i=1; i<=3; i++) {
457  if(source_shift(i).get_etat() != ETATZERO) {
458  if(source_shift(i).dz_nonzero()) {
459  assert( source_shift(i).get_dzpuis() == 4 ) ;
460  }
461  else{
462  (source_shift.set(i)).set_dzpuis(4) ;
463  }
464  }
465  }
466 
467  double lambda_shift = double(1) / double(3) ;
468 
469  if ( mg->get_np(0) == 1 ) {
470  lambda_shift = 0 ;
471  }
472 
473  Tenseur csource_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
474  csource_shift.set_etat_qcq() ;
475  for (int i=1; i<=3; i++) {
476  csource_shift.set(i-1) = source_shift(i) ;
477  }
478  csource_shift.set(2).set_etat_zero() ; //## bizarre...
479 
480  csource_shift.poisson_vect(lambda_shift, par_poisson_vect,
481  cshift, cw_shift, ckhi_shift) ;
482 
483  for (int i=1; i<=3; i++) {
484  beta.set(i) = - cshift(i-1) ;
485  beta.set(i).set_dzpuis(0) ; //## bizarre...
486  w_shift.set(i) = cw_shift(i-1) ;
487  }
488  khi_shift = ckhi_shift() ;
489 
490  cout << "Test of the Poisson equation for shift_x :" << endl ;
491  err = source_shift(1).test_poisson(-beta(1), cout, true) ;
492  diff_shift_x = err(0, 0) ;
493 
494  cout << "Test of the Poisson equation for shift_y :" << endl ;
495  err = source_shift(2).test_poisson(-beta(2), cout, true) ;
496  diff_shift_y = err(0, 0) ;
497 
498  // Computation of tnphi and nphi from the Cartesian components
499  // of the shift
500  // -----------------------------------------------------------
501 
502  fait_nphi() ;
503 
504 
505 
506  //----------------------------------------------------
507  // Adaptation of the mapping to the new enthalpy field
508  //----------------------------------------------------
509 
510  // Shall the adaptation be performed ?
511  // ---------------------------------
512 
513  adapt_flag = 0 ; // No adaptation of the mapping
514 
515  mp_prev = mp_et ;
516 
517 
518  //---------------------------------------------------------
519  // Matter source terms in the gravitational field equations
520  //---------------------------------------------------------
521 
522  //## Computation of tnphi and nphi from the Cartesian components
523  // of the shift for the test in hydro_euler():
524 
525  fait_nphi() ;
526 
527  // Update of the energy-momentum tensor
528  update_ener_mom() ;
529 
530 
531  //-------------------------------------------------------
532  // 2-D Poisson equation for tggg
533  //-------------------------------------------------------
534 
535  Cmp csource_tggg(source_tggg) ;
536  Cmp ctggg(tggg) ;
537  mp.poisson2d(csource_tggg, mp.cmp_zero(), par_poisson_tggg,
538  ctggg) ;
539  tggg = ctggg ;
540 
541 
542  //-------------------------------------------------------
543  // 2-D Poisson equation for dzeta
544  //-------------------------------------------------------
545 
546  Cmp csource_dzf(source_dzf) ;
547  Cmp csource_dzq(source_dzq) ;
548  Cmp cdzeta(dzeta) ;
549  mp.poisson2d(csource_dzf, csource_dzq, par_poisson_dzeta,
550  cdzeta) ;
551  dzeta = cdzeta ;
552 
553  err_grv2 = lbda_grv2 - 1;
554  cout << "GRV2: " << err_grv2 << endl ;
555 
556 
557  //---------------------------------------
558  // Computation of the metric coefficients (except for N^phi)
559  //---------------------------------------
560 
561  // Relaxations on nu and dzeta :
562 
563  if (mer >= 10) {
564  logn = relax * logn + relax_prev * logn_prev ;
565 
566  dzeta = relax * dzeta + relax_prev * dzeta_prev ;
567  }
568 
569  // Update of the metric coefficients N, A, B and computation of K_ij :
570 
571  //## Spherical components of the shift vector are restored
573  update_metric() ;
574  //## Back to Cartesian components
576 
577  //-----------------------
578  // Informations display
579  //-----------------------
580 
581  cout << *this << endl ;
582 
583 
584  //------------------------------------------------------------
585  // Relative change in Phi with respect to previous step
586  //------------------------------------------------------------
587 
588  Tbl diff_phi_tbl = diffrel( rphi, rphi_prev ) ;
589  diff_phi = diff_phi_tbl(0) ;
590  for (int l=1; l<nzet; l++) {
591  diff_phi += diff_phi_tbl(l) ;
592  }
593  diff_phi /= nzet ;
594 
595  fichconv << " " << log10( fabs(diff_phi) + 1.e-16 ) ;
596  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
597  fichconv << " " << log10( fabs(max_triax) + 1.e-16 ) ;
598 
599  vit_triax = 0 ;
600  if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
601  vit_triax = (max_triax - max_triax_prev) / max_triax_prev ;
602  }
603 
604  fichconv << " " << vit_triax ;
605 
606  //------------------------------
607  // Recycling for the next step
608  //------------------------------
609 
610  rphi_prev = rphi ;
611  logn_prev = logn ;
612  dzeta_prev = dzeta ;
613  max_triax_prev = max_triax ;
614 
615  fichconv << endl ;
616  fichevol << endl ;
617  fichconv.flush() ;
618  fichevol.flush() ;
619 
620  } // End of main loop
621 
622  //=========================================================================
623  // End of iteration
624  //=========================================================================
625 
626  ssjm1_nuf = cssjm1_nuf ;
627  ssjm1_nuq = cssjm1_nuq ;
628  ssjm1_tggg = cssjm1_tggg ;
629  ssjm1_khi = cssjm1_khi ;
630  for (int i=1; i<=3; i++) {
631  ssjm1_wshift.set(i) = cssjm1_wshift(i-1) ;
632  }
633 
634  // Spherical components of the shift vector are restored
636 
637  fichconv.close() ;
638  fichevol.close() ;
639 
640 }
641 }
Scalar logn
Logarithm of the lapse N .
Definition: compobj.h:498
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: compobj.h:542
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
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 mom_euler
Total 3-momentum density in the Eulerian frame.
Definition: compobj.h:150
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 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
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Scalar ak_car
Scalar .
Definition: compobj.h:318
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
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
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: compobj.h:513
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
void update_ener_mom()
Computes the 3+1 components of the energy-momentum tensor (E, P_i and S_{ij}) from the values of the ...
Definition: boson_star.C:214
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
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition: compobj.h:554
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: compobj.h:572
Sym_tensor kk
Extrinsic curvature tensor .
Definition: compobj.h:156
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
Scalar a_car
Square of the metric factor A.
Definition: compobj.h:290
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
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 nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: compobj.h:508
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
Coord phi
coordinate centered on the grid
Definition: map.h:738
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:364
Scalar b_car
Square of the metric factor B.
Definition: compobj.h:296
Coord sint
Definition: map.h:739
Scalar tggg
Metric potential .
Definition: compobj.h:519
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 bbb
Metric factor B.
Definition: compobj.h:293
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg ...
Definition: compobj.h:564
Vector beta
Shift vector .
Definition: compobj.h:141
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
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for ...
Definition: compobj.h:581
Scalar rphi
Real part of the scalar field Phi.
Definition: boson_star.h:72
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Scalar dzeta
Metric potential .
Definition: compobj.h:516
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
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
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
Scalar ener_euler
Total energy density E in the Eulerian frame.
Definition: compobj.h:147
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition: compobj.h:548
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
Scalar nn
Lapse function N .
Definition: compobj.h:138
virtual void equilibrium(double rphi_c, double iphi_c, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, Tbl &diff, Param *=0x0)
Solves the equation satisfied by the scalar field.
Metric gamma
3-metric
Definition: compobj.h:144
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
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
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
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
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift ...
Definition: star_QI.C:389
void update_metric()
Computes metric coefficients from known potentials.
Definition: star_QI.C:413
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Vector w_shift
Vector used in the decomposition of shift , following Shibata&#39;s prescription [Prog.
Definition: compobj.h:532
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:135
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
Definition: compobj.h:153