LORENE
star_bin_equilibrium.C
1 
2 /*
3  * Method of class Star_bin to compute an equilibrium configuration
4  *
5  * (see file star.h for documentation).
6  */
7 /*
8  * Copyright (c) 2004 Francois Limousin
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 version 2
14  * as published by the Free Software Foundation.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 char star_bin_equilibrium_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $" ;
28 
29 /*
30  * $Id: star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $
31  * $Log: star_bin_equilibrium.C,v $
32  * Revision 1.29 2014/10/13 08:53:38 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.28 2014/10/06 15:13:16 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.27 2006/08/01 14:26:01 f_limousin
39  * Display
40  *
41  * Revision 1.26 2006/05/31 09:26:04 f_limousin
42  * Modif. of the size of the different domains
43  *
44  * Revision 1.25 2006/04/11 14:24:44 f_limousin
45  * New version of the code : improvement of the computation of some
46  * critical sources, estimation of the dirac gauge, helical symmetry...
47  *
48  * Revision 1.24 2005/11/03 13:27:09 f_limousin
49  * Final version for the letter.
50  *
51  * Revision 1.23 2005/09/14 12:48:02 f_limousin
52  * Comment graphical outputs.
53  *
54  * Revision 1.22 2005/09/14 12:30:52 f_limousin
55  * Saving of fields lnq and logn in class Star.
56  *
57  * Revision 1.21 2005/09/13 19:38:31 f_limousin
58  * Reintroduction of the resolution of the equations in cartesian coordinates.
59  *
60  * Revision 1.20 2005/04/08 12:36:44 f_limousin
61  * Just to avoid warnings...
62  *
63  * Revision 1.19 2005/02/24 16:27:21 f_limousin
64  * Add mermax_poisson and relax_poisson in the parameters of the function.
65  *
66  * Revision 1.18 2005/02/24 16:04:13 f_limousin
67  * Change the name of some variables (for instance dcov_logn --> dlogn).
68  * Improve the resolution of the tensorial poisson equation for hh.
69  *
70  * Revision 1.17 2005/02/18 13:14:18 j_novak
71  * Changing of malloc/free to new/delete + suppression of some unused variables
72  * (trying to avoid compilation warnings).
73  *
74  * Revision 1.16 2005/02/17 17:32:53 f_limousin
75  * Change the name of some quantities to be consistent with other classes
76  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
77  *
78  * Revision 1.15 2005/02/11 18:13:47 f_limousin
79  * Important modification : all the poisson equations for the metric
80  * quantities are now solved on an affine mapping.
81  *
82  * Revision 1.14 2004/12/17 16:23:19 f_limousin
83  * Modif. comments.
84  *
85  * Revision 1.13 2004/06/22 12:49:12 f_limousin
86  * Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
87  *
88  * Revision 1.12 2004/05/27 12:41:00 p_grandclement
89  * correction of some shadowed variables
90  *
91  * Revision 1.11 2004/05/25 14:18:00 f_limousin
92  * Include filters
93  *
94  * Revision 1.10 2004/05/10 10:26:22 f_limousin
95  * Minor changes to avoid warnings in the compilation of Lorene
96  *
97  * Revision 1.9 2004/04/08 16:32:48 f_limousin
98  * The new variable is ln(Q) instead of Q=psi^2*N. It improves the
99  * convergence of the code.
100  *
101  * Revision 1.8 2004/03/25 10:29:26 j_novak
102  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
103  *
104  * Revision 1.7 2004/03/23 09:56:09 f_limousin
105  * Many minor changes
106  *
107  * Revision 1.6 2004/02/27 21:16:32 e_gourgoulhon
108  * Function contract_desal replaced by function contract with
109  * argument desaliasing set to true.
110  *
111  * Revision 1.5 2004/02/27 09:51:51 f_limousin
112  * Many minor changes.
113  *
114  * Revision 1.4 2004/02/21 17:05:13 e_gourgoulhon
115  * Method Scalar::point renamed Scalar::val_grid_point.
116  * Method Scalar::set_point renamed Scalar::set_grid_point.
117  *
118  * Revision 1.3 2004/01/20 15:17:48 f_limousin
119  * First version
120  *
121  *
122  * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $
123  *
124  */
125 
126 // C headers
127 #include <cmath>
128 
129 // Lorene headers
130 #include "cmp.h"
131 #include "tenseur.h"
132 #include "metrique.h"
133 #include "star.h"
134 #include "param.h"
135 #include "graphique.h"
136 #include "utilitaires.h"
137 #include "tensor.h"
138 #include "nbr_spx.h"
139 #include "unites.h"
140 
141 
142 namespace Lorene {
143 void Star_bin::equilibrium(double ent_c, int mermax, int mermax_potvit,
144  int mermax_poisson, double relax_poisson,
145  double relax_potvit, double thres_adapt,
146  Tbl& diff, double om) {
147 
148  // Fundamental constants and units
149  // -------------------------------
150  using namespace Unites ;
151 
152  // Initializations
153  // ---------------
154 
155  const Mg3d* mg = mp.get_mg() ;
156  int nz = mg->get_nzone() ; // total number of domains
157 
158  // The following is required to initialize mp_prev as a Map_et:
159  Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
160 
161  // Domain and radial indices of points at the surface of the star:
162  int l_b = nzet - 1 ;
163  int i_b = mg->get_nr(l_b) - 1 ;
164  int k_b ;
165  int j_b ;
166 
167  // Value of the enthalpy defining the surface of the star
168  double ent_b = 0 ;
169 
170  // Error indicators
171  // ----------------
172 
173  double& diff_ent = diff.set(0) ;
174  double& diff_vel_pot = diff.set(1) ;
175  double& diff_logn = diff.set(2) ;
176  double& diff_lnq = diff.set(3) ;
177  double& diff_beta_x = diff.set(4) ;
178  double& diff_beta_y = diff.set(5) ;
179  double& diff_beta_z = diff.set(6) ;
180  double& diff_h11 = diff.set(7) ;
181  double& diff_h21 = diff.set(8) ;
182  double& diff_h31 = diff.set(9) ;
183  double& diff_h22 = diff.set(10) ;
184  double& diff_h32 = diff.set(11) ;
185  double& diff_h33 = diff.set(12) ;
186 
187 
188 
189  // Parameters for te function Map_et::adapt
190  // -----------------------------------------
191 
192  Param par_adapt ;
193  int nitermax = 100 ;
194  int niter ;
195  int adapt_flag = 1 ; // 1 = performs the full computation,
196  // 0 = performs only the rescaling by
197  // the factor alpha_r
198  //## int nz_search = nzet + 1 ; // Number of domains for searching the
199  // enthalpy
200  int nz_search = nzet ; // Number of domains for searching the enthalpy
201  // isosurfaces
202 
203  double precis_secant = 1.e-14 ;
204  double alpha_r ;
205  double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
206 
207  Tbl ent_limit(nz) ;
208 
209 
210  par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
211  // locate zeros by the secant method
212  par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
213  // to the isosurfaces of ent is to be
214  // performed
215  par_adapt.add_int(nz_search, 2) ; // number of domains to search
216  // the enthalpy isosurface
217  par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
218  // 0 = performs only the rescaling by
219  // the factor alpha_r
220  par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
221  // (theta_*, phi_*)
222  par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
223  // (theta_*, phi_*)
224 
225  par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used in
226  // the secant method
227 
228  par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
229  // the determination of zeros by
230  // the secant method
231  par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping,
232  // 0 = contracting mapping
233 
234  par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
235  // distances will be multiplied
236 
237  par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
238  // to define the isosurfaces.
239 
240 
241  Cmp ssjm1logn (ssjm1_logn) ;
242  Cmp ssjm1lnq (ssjm1_lnq) ;
243  Cmp ssjm1h11 (ssjm1_h11) ;
244  Cmp ssjm1h21 (ssjm1_h21) ;
245  Cmp ssjm1h31 (ssjm1_h31) ;
246  Cmp ssjm1h22 (ssjm1_h22) ;
247  Cmp ssjm1h32 (ssjm1_h32) ;
248  Cmp ssjm1h33 (ssjm1_h33) ;
249 
250 
251  ssjm1logn.set_etat_qcq() ;
252  ssjm1lnq.set_etat_qcq() ;
253  ssjm1h11.set_etat_qcq() ;
254  ssjm1h21.set_etat_qcq() ;
255  ssjm1h31.set_etat_qcq() ;
256  ssjm1h22.set_etat_qcq() ;
257  ssjm1h32.set_etat_qcq() ;
258  ssjm1h33.set_etat_qcq() ;
259 
260 
261  double precis_poisson = 1.e-16 ;
262 
263  // Parameters for the function Scalar::poisson for logn_auto
264  // ---------------------------------------------------------------
265 
266  Param par_logn ;
267 
268  par_logn.add_int(mermax_poisson, 0) ; // maximum number of iterations
269  par_logn.add_double(relax_poisson, 0) ; // relaxation parameter
270  par_logn.add_double(precis_poisson, 1) ; // required precision
271  par_logn.add_int_mod(niter, 0) ; // number of iterations actually used
272  par_logn.add_cmp_mod( ssjm1logn ) ;
273 
274  // Parameters for the function Scalar::poisson for lnq_auto
275  // ---------------------------------------------------------------
276 
277  Param par_lnq ;
278 
279  par_lnq.add_int(mermax_poisson, 0) ; // maximum number of iterations
280  par_lnq.add_double(relax_poisson, 0) ; // relaxation parameter
281  par_lnq.add_double(precis_poisson, 1) ; // required precision
282  par_lnq.add_int_mod(niter, 0) ; // number of iterations actually used -
283  par_lnq.add_cmp_mod( ssjm1lnq ) ;
284 
285  // Parameters for the function Vector::poisson for beta method 2
286  // ---------------------------------------------------------------
287 
288  Param par_beta2 ;
289 
290  par_beta2.add_int(mermax_poisson, 0) ; // maximum number of iterations
291  par_beta2.add_double(relax_poisson, 0) ; // relaxation parameter
292  par_beta2.add_double(precis_poisson, 1) ; // required precision
293  par_beta2.add_int_mod(niter, 0) ; // number of iterations actually used
294 
295  Cmp ssjm1khi (ssjm1_khi) ;
296  Tenseur ssjm1wbeta(mp, 1, CON, mp.get_bvect_cart()) ;
297  ssjm1wbeta.set_etat_qcq() ;
298  for (int i=0; i<3; i++) {
299  ssjm1wbeta.set(i) = Cmp(ssjm1_wbeta(i+1)) ;
300  }
301 
302  par_beta2.add_cmp_mod(ssjm1khi) ;
303  par_beta2.add_tenseur_mod(ssjm1wbeta) ;
304 
305  // Parameters for the function Scalar::poisson for h11_auto
306  // -------------------------------------------------------------
307 
308  Param par_h11 ;
309 
310  par_h11.add_int(mermax_poisson, 0) ; // maximum number of iterations
311  par_h11.add_double(relax_poisson, 0) ; // relaxation parameter
312  par_h11.add_double(precis_poisson, 1) ; // required precision
313  par_h11.add_int_mod(niter, 0) ; // number of iterations actually used
314  par_h11.add_cmp_mod( ssjm1h11 ) ;
315 
316  // Parameters for the function Scalar::poisson for h21_auto
317  // -------------------------------------------------------------
318 
319  Param par_h21 ;
320 
321  par_h21.add_int(mermax_poisson, 0) ; // maximum number of iterations
322  par_h21.add_double(relax_poisson, 0) ; // relaxation parameter
323  par_h21.add_double(precis_poisson, 1) ; // required precision
324  par_h21.add_int_mod(niter, 0) ; // number of iterations actually used
325  par_h21.add_cmp_mod( ssjm1h21 ) ;
326 
327  // Parameters for the function Scalar::poisson for h31_auto
328  // -------------------------------------------------------------
329 
330  Param par_h31 ;
331 
332  par_h31.add_int(mermax_poisson, 0) ; // maximum number of iterations
333  par_h31.add_double(relax_poisson, 0) ; // relaxation parameter
334  par_h31.add_double(precis_poisson, 1) ; // required precision
335  par_h31.add_int_mod(niter, 0) ; // number of iterations actually used
336  par_h31.add_cmp_mod( ssjm1h31 ) ;
337 
338  // Parameters for the function Scalar::poisson for h22_auto
339  // -------------------------------------------------------------
340 
341  Param par_h22 ;
342 
343  par_h22.add_int(mermax_poisson, 0) ; // maximum number of iterations
344  par_h22.add_double(relax_poisson, 0) ; // relaxation parameter
345  par_h22.add_double(precis_poisson, 1) ; // required precision
346  par_h22.add_int_mod(niter, 0) ; // number of iterations actually used
347  par_h22.add_cmp_mod( ssjm1h22 ) ;
348 
349  // Parameters for the function Scalar::poisson for h32_auto
350  // -------------------------------------------------------------
351 
352  Param par_h32 ;
353 
354  par_h32.add_int(mermax_poisson, 0) ; // maximum number of iterations
355  par_h32.add_double(relax_poisson, 0) ; // relaxation parameter
356  par_h32.add_double(precis_poisson, 1) ; // required precision
357  par_h32.add_int_mod(niter, 0) ; // number of iterations actually used
358  par_h32.add_cmp_mod( ssjm1h32 ) ;
359 
360  // Parameters for the function Scalar::poisson for h33_auto
361  // -------------------------------------------------------------
362 
363  Param par_h33 ;
364 
365  par_h33.add_int(mermax_poisson, 0) ; // maximum number of iterations
366  par_h33.add_double(relax_poisson, 0) ; // relaxation parameter
367  par_h33.add_double(precis_poisson, 1) ; // required precision
368  par_h33.add_int_mod(niter, 0) ; // number of iterations actually used
369  par_h33.add_cmp_mod( ssjm1h33 ) ;
370 
371 
372  // External potential
373  // See Eq (99) from Gourgoulhon et al. (2001)
374  // ------------------
375 
376  cout << "logn_comp" << norme(logn_comp) << endl ;
377  cout << "pot_centri" << norme(pot_centri) << endl ;
378  cout << "loggam" << norme(loggam) << endl ;
379  Scalar pot_ext = logn_comp + pot_centri + loggam ;
380 
381  Scalar ent_jm1 = ent ; // Enthalpy at previous step
382 
383  Scalar source_tot(mp) ; // source term in the equation for hij_auto,
384  // logn_auto and beta_auto
385 
386  Vector source_beta(mp, CON, mp.get_bvect_cart()) ; // source term
387  // in the equation for beta_auto
388 
389 
390 
391  //=========================================================================
392  // Start of iteration
393  //=========================================================================
394 
395  for(int mer=0 ; mer<mermax ; mer++ ) {
396 
397  cout << "-----------------------------------------------" << endl ;
398  cout << "step: " << mer << endl ;
399  cout << "diff_ent = " << diff_ent << endl ;
400 
401  //-----------------------------------------------------
402  // Resolution of the elliptic equation for the velocity
403  // scalar potential
404  //-----------------------------------------------------
405 
406  if (irrotational) {
407  diff_vel_pot = velocity_potential(mermax_potvit, precis_poisson,
408  relax_potvit) ;
409 
410  }
411 
412  diff_vel_pot = 0. ; // to avoid the warning
413 
414  //-----------------------------------------------------
415  // Computation of the new radial scale
416  //--------------------------------------------------
417 
418  // alpha_r (r = alpha_r r') is determined so that the enthalpy
419  // takes the requested value ent_b at the stellar surface
420 
421  // Values at the center of the star:
422  double logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
423  double pot_ext_c = pot_ext.val_grid_point(0, 0, 0, 0) ;
424 
425  // Search for the reference point (theta_*, phi_*) [notation of
426  // Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
427  // at the surface of the star
428  // ------------------------------------------------------------
429  double alpha_r2 = 0 ;
430  for (int k=0; k<mg->get_np(l_b); k++) {
431  for (int j=0; j<mg->get_nt(l_b); j++) {
432 
433  double pot_ext_b = pot_ext.val_grid_point(l_b, k, j, i_b) ;
434  double logn_auto_b = logn_auto.val_grid_point(l_b, k, j, i_b) ;
435  // See Eq (100) from Gourgoulhon et al. (2001)
436  double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
437 
438  ( logn_auto_b - logn_auto_c ) ;
439 
440  if (alpha_r2_jk > alpha_r2) {
441  alpha_r2 = alpha_r2_jk ;
442  k_b = k ;
443  j_b = j ;
444  }
445 
446  }
447  }
448 
449  alpha_r = sqrt(alpha_r2) ;
450 
451  cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
452  << alpha_r << endl ;
453 
454  // New value of logn_auto
455  // ----------------------
456 
457  logn_auto = alpha_r2 * logn_auto ;
458  logn_auto_c = logn_auto.val_grid_point(0, 0, 0, 0) ;
459 
460  //------------------------------------------------------------
461  // Change the values of the inner points of the second domain
462  // by those of the outer points of the first domain
463  //------------------------------------------------------------
464 
465  logn_auto.set_spectral_va().smooth(nzet, logn_auto.set_spectral_va()) ;
466 
467  //------------------------------------------
468  // First integral --> enthalpy in all space
469  // See Eq (98) from Gourgoulhon et al. (2001)
470  //-------------------------------------------
471 
472  ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
473  cout.precision(8) ;
474  cout << "pot" << norme(pot_ext) << endl ;
475 
476  (ent.set_spectral_va()).smooth(nzet, ent.set_spectral_va()) ;
477 
478  //----------------------------------------------------
479  // Adaptation of the mapping to the new enthalpy field
480  //----------------------------------------------------
481 
482  // Shall the adaptation be performed (cusp) ?
483  // ------------------------------------------
484 
485  double dent_eq = ent.dsdr().val_point(ray_eq(),M_PI/2.,0.) ;
486  double dent_pole = ent.dsdr().val_point(ray_pole(),0.,0.) ;
487  double rap_dent = fabs( dent_eq / dent_pole ) ;
488  cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
489 
490  if ( rap_dent < thres_adapt ) {
491  adapt_flag = 0 ; // No adaptation of the mapping
492  cout << "******* FROZEN MAPPING *********" << endl ;
493  }
494  else{
495  adapt_flag = 1 ; // The adaptation of the mapping is to be
496  // performed
497  }
498 
499  ent_limit.set_etat_qcq() ;
500  for (int l=0; l<nzet; l++) { // loop on domains inside the star
501  ent_limit.set(l) = ent.val_grid_point(l, k_b, j_b, i_b) ;
502  }
503  ent_limit.set(nzet-1) = ent_b ;
504 
505  Map_et mp_prev = mp_et ;
506 
507  Cmp ent_cmp(ent) ;
508  mp.adapt(ent_cmp, par_adapt) ;
509  ent = ent_cmp ;
510 
511  // Readjustment of the external boundary of domain l=nzet
512  // to keep a fixed ratio with respect to star's surface
513 
514  if (nz>= 5) {
515  double separation = 2. * fabs(mp.get_ori_x()) ;
516  double ray_eqq = ray_eq() ;
517  double ray_eqq_pi = ray_eq_pi() ;
518  double new_rr_out_2 = (separation - ray_eqq) * 0.95 ;
519  double new_rr_out_3 = (separation + ray_eqq_pi) * 1.05 ;
520 
521  double rr_in_1 = mp.val_r(1,-1., M_PI/2, 0.) ;
522  double rr_out_1 = mp.val_r(1, 1., M_PI/2, 0.) ;
523  double rr_out_2 = mp.val_r(2, 1., M_PI/2, 0.) ;
524  double rr_out_3 = mp.val_r(3, 1., M_PI/2, 0.) ;
525 
526  mp.resize(1, 0.5*(new_rr_out_2 + rr_in_1) / rr_out_1) ;
527  mp.resize(2, new_rr_out_2 / rr_out_2) ;
528  mp.resize(3, new_rr_out_3 / rr_out_3) ;
529 
530  for (int dd=4; dd<=nz-2; dd++) {
531  mp.resize(dd, new_rr_out_3 * pow(4., dd-3) /
532  mp.val_r(dd, 1., M_PI/2, 0.)) ;
533  }
534 
535  }
536  else {
537  cout << "too small number of domains" << endl ;
538  }
539 
540  //----------------------------------------------------
541  // Computation of the enthalpy at the new grid points
542  //----------------------------------------------------
543 
544  mp_prev.homothetie(alpha_r) ;
545 
546  Cmp ent_cmp2 (ent) ;
547  mp.reevaluate_symy(&mp_prev, nzet+1, ent_cmp2) ;
548  ent = ent_cmp2 ;
549 
550  double ent_s_max = -1 ;
551  int k_s_max = -1 ;
552  int j_s_max = -1 ;
553  for (int k=0; k<mg->get_np(l_b); k++) {
554  for (int j=0; j<mg->get_nt(l_b); j++) {
555  double xx = fabs( ent.val_grid_point(l_b, k, j, i_b) ) ;
556  if (xx > ent_s_max) {
557  ent_s_max = xx ;
558  k_s_max = k ;
559  j_s_max = j ;
560  }
561  }
562  }
563  cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
564  << " and nzet : " << endl ;
565  cout << " " << ent_s_max << " reached for k = " << k_s_max <<
566  " and j = " << j_s_max << endl ;
567 
568  //----------------------------------------------------
569  // Equation of state
570  //----------------------------------------------------
571 
572  equation_of_state() ; // computes new values for nbar (n), ener (e)
573  // and press (p) from the new ent (H)
574 
575  //---------------------------------------------------------
576  // Matter source terms in the gravitational field equations
577  //---------------------------------------------------------
578 
579  hydro_euler() ; // computes new values for ener_euler (E),
580  // s_euler (S) and u_euler (U^i)
581 
582 
583  // -------------------------------
584  // AUXILIARY QUANTITIES
585  // -------------------------------
586 
587  // Derivatives of N and logN
588  //--------------------------
589 
590  const Vector dcov_logn_auto = logn_auto.derive_cov(flat) ;
591 
592  Tensor dcovdcov_logn_auto = (logn_auto.derive_cov(flat))
593  .derive_cov(flat) ;
594  dcovdcov_logn_auto.inc_dzpuis() ;
595 
596  // Derivatives of lnq, phi and Q
597  //-------------------------------
598 
599  const Scalar phi (0.5 * (lnq - logn)) ;
600  const Scalar phi_auto (0.5 * (lnq_auto - logn_auto)) ;
601 
602  const Vector dcov_phi_auto = phi_auto.derive_cov(flat) ;
603 
604  const Vector dcov_lnq = 2*dcov_phi + dcov_logn ;
605  const Vector dcon_lnq = 2*dcon_phi + dcon_logn ;
606  const Vector dcov_lnq_auto = lnq_auto.derive_cov(flat) ;
607  Tensor dcovdcov_lnq_auto = dcov_lnq_auto.derive_cov(flat) ;
608  dcovdcov_lnq_auto.inc_dzpuis() ;
609 
610  Scalar qq = exp(lnq) ;
611  qq.std_spectral_base() ;
612  const Vector& dcov_qq = qq.derive_cov(flat) ;
613 
614  const Scalar& divbeta_auto = beta_auto.divergence(flat) ;
615  const Tensor& dcov_beta_auto = beta_auto.derive_cov(flat) ;
616  Tensor dcovdcov_beta_auto = beta_auto.derive_cov(flat)
617  .derive_cov(flat) ;
618  dcovdcov_beta_auto.inc_dzpuis() ;
619 
620 
621  // Derivatives of hij, gtilde...
622  //------------------------------
623 
624  Scalar psi2 (pow(psi4, 0.5)) ;
625  psi2.std_spectral_base() ;
626 
627  const Tensor& dcov_hij = hij.derive_cov(flat) ;
628  const Tensor& dcon_hij = hij.derive_con(flat) ;
629  const Tensor& dcov_hij_auto = hij_auto.derive_cov(flat) ;
630 
631  const Sym_tensor& gtilde_cov = gtilde.cov() ;
632  const Sym_tensor& gtilde_con = gtilde.con() ;
633  const Tensor& dcov_gtilde = gtilde_cov.derive_cov(flat) ;
634 
635  Connection gamijk (gtilde, flat) ;
636  const Tensor& deltaijk = gamijk.get_delta() ;
637 
638  // H^i and its derivatives ( = O in Dirac gauge)
639  // ---------------------------------------------
640 
641  double lambda_dirac = 0. ;
642 
643  const Vector hdirac = lambda_dirac * hij.divergence(flat) ;
644  const Vector hdirac_auto = lambda_dirac * hij_auto.divergence(flat) ;
645 
646  Tensor dcov_hdirac = lambda_dirac * hdirac.derive_cov(flat) ;
647  dcov_hdirac.inc_dzpuis() ;
648  Tensor dcov_hdirac_auto = lambda_dirac * hdirac_auto.derive_cov(flat) ;
649  dcov_hdirac_auto.inc_dzpuis() ;
650  Tensor dcon_hdirac_auto = lambda_dirac * hdirac_auto.derive_con(flat) ;
651  dcon_hdirac_auto.inc_dzpuis() ;
652 
653 
654  //--------------------------------------------------------
655  // Poisson equation for logn_auto (nu_auto)
656  //--------------------------------------------------------
657 
658  // Source
659  //--------
660 
661  int nr = mp.get_mg()->get_nr(0) ;
662  int nt = mp.get_mg()->get_nt(0) ;
663  int np = mp.get_mg()->get_np(0) ;
664 
665  Scalar source1(mp) ;
666  Scalar source2(mp) ;
667  Scalar source3(mp) ;
668  Scalar source4(mp) ;
669  Scalar source5(mp) ;
670  Scalar source6(mp) ;
671  Scalar source7(mp) ;
672  Scalar source8(mp) ;
673 
674  source1 = qpig * psi4 % (ener_euler + s_euler) ;
675 
676  source2 = psi4 % (kcar_auto + kcar_comp) ;
677 
678  source3 = - contract(dcov_logn_auto, 0, dcon_logn, 0, true)
679  - 2. * contract(contract(gtilde_con, 0, dcov_phi, 0),
680  0, dcov_logn_auto, 0, true) ;
681 
682  source4 = - contract(hij, 0, 1, dcovdcov_logn_auto +
683  dcov_logn_auto*dcov_logn, 0, 1) ;
684 
685  source5 = - contract(hdirac, 0, dcov_logn_auto, 0) ;
686 
687  source_tot = source1 + source2 + source3 + source4 + source5 ;
688 
689 
690  cout << "moyenne de la source 1 pour logn_auto" << endl ;
691  cout << norme(source1/(nr*nt*np)) << endl ;
692  cout << "moyenne de la source 2 pour logn_auto" << endl ;
693  cout << norme(source2/(nr*nt*np)) << endl ;
694  cout << "moyenne de la source 3 pour logn_auto" << endl ;
695  cout << norme(source3/(nr*nt*np)) << endl ;
696  cout << "moyenne de la source 4 pour logn_auto" << endl ;
697  cout << norme(source4/(nr*nt*np)) << endl ;
698  cout << "moyenne de la source 5 pour logn_auto" << endl ;
699  cout << norme(source5/(nr*nt*np)) << endl ;
700  cout << "moyenne de la source pour logn_auto" << endl ;
701  cout << norme(source_tot/(nr*nt*np)) << endl ;
702 
703  // Resolution of the Poisson equation
704  // ----------------------------------
705 
706  source_tot.poisson(par_logn, logn_auto) ;
707  ssjm1_logn = ssjm1logn ;
708 
709  cout << "logn_auto" << endl << norme(logn_auto/(nr*nt*np)) << endl ;
710 
711  // Check: has the Poisson equation been correctly solved ?
712  // -----------------------------------------------------
713 
714  Tbl tdiff_logn = diffrel(logn_auto.laplacian(), source_tot) ;
715  cout <<
716  "Relative error in the resolution of the equation for logn_auto : "
717  << endl ;
718  for (int l=0; l<nz; l++) {
719  cout << tdiff_logn(l) << " " ;
720  }
721  cout << endl ;
722  diff_logn = max(abs(tdiff_logn)) ;
723 
724  //--------------------------------------------------------
725  // Poisson equation for lnq_auto
726  //--------------------------------------------------------
727 
728  // Source
729  //--------
730 
731  source1 = qpig * psi4 % s_euler ;
732 
733  source2 = 0.75 * psi4 % (kcar_auto + kcar_comp) ;
734 
735  source3 = - contract(dcon_lnq, 0, dcov_lnq_auto, 0, true) ;
736 
737  source4 = 2. * contract(contract(gtilde_con, 0, dcov_phi, 0), 0,
738  dcov_phi_auto + dcov_logn_auto, 0, true) ;
739 
740  source5 = 0.0625 * contract(gtilde_con, 0, 1, contract(
741  dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
742 
743  source6 = - 0.125 * contract(gtilde_con, 0, 1, contract(
744  dcov_hij_auto, 0, 1, dcov_gtilde, 0, 2), 0, 1) ;
745 
746  source7 = - contract(hij, 0, 1, dcovdcov_lnq_auto + dcov_lnq_auto *
747  dcov_lnq, 0, 1) ;
748 
749  source8 = - 0.25 * contract(dcov_hdirac_auto, 0, 1)
750  - contract(hdirac, 0, dcov_lnq_auto, 0) ;
751 
752  source_tot = source1 + source2 + source3 + source4 + source5 +
753  source6 + source7 + source8 ;
754 
755 
756  cout << "moyenne de la source 1 pour lnq_auto" << endl ;
757  cout << norme(source1/(nr*nt*np)) << endl ;
758  cout << "moyenne de la source 2 pour lnq_auto" << endl ;
759  cout << norme(source2/(nr*nt*np)) << endl ;
760  cout << "moyenne de la source 3 pour lnq_auto" << endl ;
761  cout << norme(source3/(nr*nt*np)) << endl ;
762  cout << "moyenne de la source 4 pour lnq_auto" << endl ;
763  cout << norme(source4/(nr*nt*np)) << endl ;
764  cout << "moyenne de la source 5 pour lnq_auto" << endl ;
765  cout << norme(source5/(nr*nt*np)) << endl ;
766  cout << "moyenne de la source 6 pour lnq_auto" << endl ;
767  cout << norme(source6/(nr*nt*np)) << endl ;
768  cout << "moyenne de la source 7 pour lnq_auto" << endl ;
769  cout << norme(source7/(nr*nt*np)) << endl ;
770  cout << "moyenne de la source 8 pour lnq_auto" << endl ;
771  cout << norme(source8/(nr*nt*np)) << endl ;
772  cout << "moyenne de la source pour lnq_auto" << endl ;
773  cout << norme(source_tot/(nr*nt*np)) << endl ;
774 
775 
776  // Resolution of the Poisson equation
777  // ----------------------------------
778 
779  source_tot.poisson(par_lnq, lnq_auto) ;
780  ssjm1_lnq = ssjm1lnq ;
781 
782  cout << "lnq_auto" << endl << norme(lnq_auto/(nr*nt*np)) << endl ;
783 
784  // Check: has the Poisson equation been correctly solved
785  // -----------------------------------------------------
786 
787  Tbl tdiff_lnq = diffrel(lnq_auto.laplacian(), source_tot) ;
788  cout <<
789  "Relative error in the resolution of the equation for lnq : "
790  << endl ;
791  for (int l=0; l<nz; l++) {
792  cout << tdiff_lnq(l) << " " ;
793  }
794  cout << endl ;
795  diff_lnq = max(abs(tdiff_lnq)) ;
796 
797  //--------------------------------------------------------
798  // Vector Poisson equation for beta_auto
799  //--------------------------------------------------------
800 
801  // Source
802  //--------
803 
804  Vector source1_beta(mp, CON, mp.get_bvect_cart()) ;
805  Vector source2_beta(mp, CON, mp.get_bvect_cart()) ;
806  Vector source3_beta(mp, CON, mp.get_bvect_cart()) ;
807  Vector source4_beta(mp, CON, mp.get_bvect_cart()) ;
808  Vector source5_beta(mp, CON, mp.get_bvect_cart()) ;
809  Vector source6_beta(mp, CON, mp.get_bvect_cart()) ;
810  Vector source7_beta(mp, CON, mp.get_bvect_cart()) ;
811 
812  source1_beta = (4.*qpig) * nn % psi4
813  %(ener_euler + press) * u_euler ;
814 
815  source2_beta = 2. * nn * contract(tkij_auto, 1,
816  dcov_logn - 6 * dcov_phi, 0) ;
817 
818  source3_beta = - 2. * nn * contract(tkij_auto, 0, 1, deltaijk,
819  1, 2) ;
820 
821  source4_beta = - contract(hij, 0, 1, dcovdcov_beta_auto, 1, 2) ;
822 
823  source5_beta = - 0.3333333333333333 * contract(hij, 1, contract(
824  dcovdcov_beta_auto, 0, 1), 0) ;
825 
826  source6_beta.set_etat_zero() ; //hdirac_auto.derive_lie(omdsdp) ;
827  source6_beta.inc_dzpuis() ;
828 
829  source7_beta = contract(beta, 0, dcov_hdirac_auto, 1) ;
830  + 2./3. * hdirac * divbeta_auto
831  - contract(hdirac, 0, dcov_beta_auto, 1) ;
832 
833  source_beta = source1_beta + source2_beta + source3_beta
834  + source4_beta + source5_beta + source6_beta + source7_beta ;
835 
836 
837  cout << "moyenne de la source 1 pour beta_auto" << endl ;
838  cout << norme(source1_beta(1)/(nr*nt*np)) << endl ;
839  cout << norme(source1_beta(2)/(nr*nt*np)) << endl ;
840  cout << norme(source1_beta(3)/(nr*nt*np)) << endl ;
841  cout << "moyenne de la source 2 pour beta_auto" << endl ;
842  cout << norme(source2_beta(1)/(nr*nt*np)) << endl ;
843  cout << norme(source2_beta(2)/(nr*nt*np)) << endl ;
844  cout << norme(source2_beta(3)/(nr*nt*np)) << endl ;
845  cout << "moyenne de la source 3 pour beta_auto" << endl ;
846  cout << norme(source3_beta(1)/(nr*nt*np)) << endl ;
847  cout << norme(source3_beta(2)/(nr*nt*np)) << endl ;
848  cout << norme(source3_beta(3)/(nr*nt*np)) << endl ;
849  cout << "moyenne de la source 4 pour beta_auto" << endl ;
850  cout << norme(source4_beta(1)/(nr*nt*np)) << endl ;
851  cout << norme(source4_beta(2)/(nr*nt*np)) << endl ;
852  cout << norme(source4_beta(3)/(nr*nt*np)) << endl ;
853  cout << "moyenne de la source 5 pour beta_auto" << endl ;
854  cout << norme(source5_beta(1)/(nr*nt*np)) << endl ;
855  cout << norme(source5_beta(2)/(nr*nt*np)) << endl ;
856  cout << norme(source5_beta(3)/(nr*nt*np)) << endl ;
857  cout << "moyenne de la source 6 pour beta_auto" << endl ;
858  cout << norme(source6_beta(1)/(nr*nt*np)) << endl ;
859  cout << norme(source6_beta(2)/(nr*nt*np)) << endl ;
860  cout << norme(source6_beta(3)/(nr*nt*np)) << endl ;
861  cout << "moyenne de la source 7 pour beta_auto" << endl ;
862  cout << norme(source7_beta(1)/(nr*nt*np)) << endl ;
863  cout << norme(source7_beta(2)/(nr*nt*np)) << endl ;
864  cout << norme(source7_beta(3)/(nr*nt*np)) << endl ;
865  cout << "moyenne de la source pour beta_auto" << endl ;
866  cout << norme(source_beta(1)/(nr*nt*np)) << endl ;
867  cout << norme(source_beta(2)/(nr*nt*np)) << endl ;
868  cout << norme(source_beta(3)/(nr*nt*np)) << endl ;
869 
870  // Resolution of the Poisson equation
871  // ----------------------------------
872 
873  // Filter for the source of beta vector
874 
875  for (int i=1; i<=3; i++) {
876  if (source_beta(i).get_etat() != ETATZERO)
877  source_beta.set(i).filtre(4) ;
878  }
879 
880  for (int i=1; i<=3; i++) {
881  if(source_beta(i).dz_nonzero()) {
882  assert( source_beta(i).get_dzpuis() == 4 ) ;
883  }
884  else{
885  (source_beta.set(i)).set_dzpuis(4) ;
886  }
887  }
888 
889  double lambda = double(1) / double(3) ;
890 
891  Tenseur source_p(mp, 1, CON, mp.get_bvect_cart() ) ;
892  source_p.set_etat_qcq() ;
893  for (int i=0; i<3; i++) {
894  source_p.set(i) = Cmp(source_beta(i+1)) ;
895  }
896 
897  Tenseur vect_auxi (mp, 1, CON, mp.get_bvect_cart()) ;
898  vect_auxi.set_etat_qcq() ;
899  for (int i=0; i<3 ;i++){
900  vect_auxi.set(i) = 0. ;
901  }
902  Tenseur scal_auxi (mp) ;
903  scal_auxi.set_etat_qcq() ;
904  scal_auxi.set().annule_hard() ;
905  scal_auxi.set_std_base() ;
906 
907  Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart() ) ;
908  resu_p.set_etat_qcq() ;
909  for (int i=0; i<3 ;i++)
910  resu_p.set(i).annule_hard() ;
911  resu_p.set_std_base() ;
912 
913  //source_p.poisson_vect(lambda, par_beta2, resu_p, vect_auxi,
914  // scal_auxi) ;
915 
916  source_p.poisson_vect_oohara(lambda, par_beta2, resu_p, scal_auxi) ;
917 
918  for (int i=1; i<=3; i++)
919  beta_auto.set(i) = resu_p(i-1) ;
920 
921  ssjm1_khi = ssjm1khi ;
922  for (int i=0; i<3; i++){
923  ssjm1_wbeta.set(i+1) = ssjm1wbeta(i) ;
924  }
925 
926  cout << "beta_auto_x" << endl << norme(beta_auto(1)/(nr*nt*np))
927  << endl ;
928  cout << "beta_auto_y" << endl << norme(beta_auto(2)/(nr*nt*np))
929  << endl ;
930  cout << "beta_auto_z" << endl << norme(beta_auto(3)/(nr*nt*np))
931  << endl ;
932 
933 
934  // Check: has the equation for beta_auto been correctly solved ?
935  // --------------------------------------------------------------
936 
937  Vector lap_beta = (beta_auto.derive_con(flat)).divergence(flat)
938  + lambda* beta_auto.divergence(flat).derive_con(flat) ;
939 
940  source_beta.dec_dzpuis() ;
941  Tbl tdiff_beta_x = diffrel(lap_beta(1), source_beta(1)) ;
942  Tbl tdiff_beta_y = diffrel(lap_beta(2), source_beta(2)) ;
943  Tbl tdiff_beta_z = diffrel(lap_beta(3), source_beta(3)) ;
944 
945  cout <<
946  "Relative error in the resolution of the equation for beta_auto : "
947  << endl ;
948  cout << "x component : " ;
949  for (int l=0; l<nz; l++) {
950  cout << tdiff_beta_x(l) << " " ;
951  }
952  cout << endl ;
953  cout << "y component : " ;
954  for (int l=0; l<nz; l++) {
955  cout << tdiff_beta_y(l) << " " ;
956  }
957  cout << endl ;
958  cout << "z component : " ;
959  for (int l=0; l<nz; l++) {
960  cout << tdiff_beta_z(l) << " " ;
961  }
962  cout << endl ;
963 
964  diff_beta_x = max(abs(tdiff_beta_x)) ;
965  diff_beta_y = max(abs(tdiff_beta_y)) ;
966  diff_beta_z = max(abs(tdiff_beta_z)) ;
967 
968 
969  if (!conf_flat) {
970 
971  //--------------------------------------------------------
972  // Poisson equation for hij
973  //--------------------------------------------------------
974 
975 
976  // Declaration of all sources
977  //---------------------------
978 
979  Scalar source_tot_hij(mp) ;
980  Tensor source_Sij(mp, 2, CON, mp.get_bvect_cart()) ;
981  Tensor source_Rij(mp, 2, CON, mp.get_bvect_cart()) ;
982  Tensor tens_temp(mp, 2, CON, mp.get_bvect_cart()) ;
983 
984  Tensor source_1 (mp, 2, CON, mp.get_bvect_cart()) ;
985  Tensor source_2 (mp, 2, CON, mp.get_bvect_cart()) ;
986  Tensor source_3a (mp, 2, CON, mp.get_bvect_cart()) ;
987  Tensor source_3b (mp, 2, CON, mp.get_bvect_cart()) ;
988  Tensor source_4 (mp, 2, CON, mp.get_bvect_cart()) ;
989  Tensor source_5 (mp, 2, CON, mp.get_bvect_cart()) ;
990  Tensor source_6 (mp, 2, CON, mp.get_bvect_cart()) ;
991 
992  // Sources
993  //-----------
994 
995  source_1 = contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
996 
997  source_2 = - contract(dcon_hij, 2, dcov_lnq_auto, 0)
998  - 2./3. * contract(hdirac, 0, dcov_lnq_auto, 0) * flat.con() ;
999 
1000  // Lie derivative of A^{ij}
1001  // --------------------------
1002 
1003  Scalar decouple_logn = (logn_auto - 1.e-8)/(logn - 2.e-8) ;
1004 
1005  // Function exp(-(r-r_0)^2/sigma^2)
1006  // --------------------------------
1007 
1008  double r0 = mp.val_r(nz-2, 1, 0, 0) ;
1009  double sigma = 1.*r0 ;
1010 
1011  Scalar rr (mp) ;
1012  rr = mp.r ;
1013 
1014  Scalar ff (mp) ;
1015  ff = exp( -(rr - r0)*(rr - r0)/sigma/sigma ) ;
1016  for (int ii=0; ii<nz-1; ii++)
1017  ff.set_domain(ii) = 1. ;
1018  ff.set_outer_boundary(nz-1, 0) ;
1019  ff.std_spectral_base() ;
1020 
1021  // ff.annule_domain(nz-1) ;
1022  //des_profile(ff, 0, 20, 0, 0) ;
1023 
1024  // Construction of Omega d/dphi
1025  // ----------------------------
1026 
1027  // Construction of D_k \Phi^i
1028  Itbl type (2) ;
1029  type.set(0) = CON ;
1030  type.set(1) = COV ;
1031 
1032  Tensor dcov_omdsdphi (mp, 2, type, mp.get_bvect_cart()) ;
1033  dcov_omdsdphi.set(1,1) = 0. ;
1034  dcov_omdsdphi.set(2,1) = om * ff ;
1035  dcov_omdsdphi.set(3,1) = 0. ;
1036  dcov_omdsdphi.set(2,2) = 0. ;
1037  dcov_omdsdphi.set(3,2) = 0. ;
1038  dcov_omdsdphi.set(3,3) = 0. ;
1039  dcov_omdsdphi.set(1,2) = -om * ff ;
1040  dcov_omdsdphi.set(1,3) = 0. ;
1041  dcov_omdsdphi.set(2,3) = 0. ;
1042  dcov_omdsdphi.std_spectral_base() ;
1043 
1044  source_3a = contract(tkij_auto, 0, dcov_omdsdphi, 1) ;
1045  source_3a.inc_dzpuis() ;
1046 
1047  // Source 3b
1048  // ------------
1049 
1050  Vector omdsdp (mp, CON, mp.get_bvect_cart()) ;
1051  Scalar yya (mp) ;
1052  yya = mp.ya ;
1053  Scalar xxa (mp) ;
1054  xxa = mp.xa ;
1055  Scalar zza (mp) ;
1056  zza = mp.za ;
1057 
1058  if (fabs(mp.get_rot_phi()) < 1e-10){
1059  omdsdp.set(1) = - om * yya * ff ;
1060  omdsdp.set(2) = om * xxa * ff ;
1061  omdsdp.set(3).annule_hard() ;
1062  }
1063  else{
1064  omdsdp.set(1) = om * yya * ff ;
1065  omdsdp.set(2) = - om * xxa * ff ;
1066  omdsdp.set(3).annule_hard() ;
1067  }
1068 
1069  omdsdp.set(1).set_outer_boundary(nz-1, 0) ;
1070  omdsdp.set(2).set_outer_boundary(nz-1, 0) ;
1071  omdsdp.std_spectral_base() ;
1072 
1073  source_3b = - contract(omdsdp, 0, tkij_auto.derive_cov(flat), 2) ;
1074 
1075  // Source 4
1076  // ---------
1077 
1078  source_4 = - tkij_auto.derive_lie(beta) ;
1079  source_4.inc_dzpuis() ;
1080  source_4 += - 2./3. * beta.divergence(flat) * tkij_auto ;
1081 
1082  source_5 = dcon_hdirac_auto ;
1083 
1084  source_6 = - 2./3. * hdirac_auto.divergence(flat) * flat.con() ;
1085  source_6.inc_dzpuis() ;
1086 
1087  // Source terms for Sij
1088  //---------------------
1089 
1090  source_Sij = 8. * nn / psi4 * phi_auto.derive_con(gtilde) *
1091  contract(gtilde_con, 0, dcov_phi, 0) ;
1092 
1093  source_Sij += 4. / psi4 * phi_auto.derive_con(gtilde) *
1094  nn * contract(gtilde_con, 0, dcov_logn, 0) +
1095  4. / psi4 * nn * contract(gtilde_con, 0, dcov_logn, 0) *
1096  phi_auto.derive_con(gtilde) ;
1097 
1098  source_Sij += - nn / (3.*psi4) * gtilde_con *
1099  ( 0.25 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
1100  dcov_gtilde, 0, 1), 0, 1)
1101  - 0.5 * contract(gtilde_con, 0, 1, contract(dcov_hij_auto, 0, 1,
1102  dcov_gtilde, 0, 2), 0, 1)) ;
1103 
1104  source_Sij += - 8.*nn / (3.*psi4) * gtilde_con *
1105  contract(dcov_phi_auto, 0, contract(gtilde_con, 0, dcov_phi, 0), 0) ;
1106 
1107  tens_temp = nn / (3.*psi4) * hdirac.divergence(flat)*hij_auto ;
1108  tens_temp.inc_dzpuis() ;
1109 
1110  source_Sij += tens_temp ;
1111 
1112  source_Sij += - 8./(3.*psi4) * contract(dcov_phi_auto, 0,
1113  nn*contract(gtilde_con, 0, dcov_logn, 0), 0) * gtilde_con ;
1114 
1115  source_Sij += 2.*nn* contract(gtilde_cov, 0, 1, tkij_auto *
1116  (tkij_auto+tkij_comp), 1, 3) ;
1117 
1118  source_Sij += - 2. * qpig * nn * ( psi4 * stress_euler
1119  - 0.33333333333333333 * s_euler * gtilde_con ) ;
1120 
1121  source_Sij += - 1./(psi4*psi2) * contract(gtilde_con, 1,
1122  contract(gtilde_con, 1, qq*dcovdcov_lnq_auto +
1123  qq*dcov_lnq_auto*dcov_lnq, 0), 1) ;
1124 
1125  source_Sij += - 0.5/(psi4*psi2) * contract(contract(hij, 1,
1126  dcov_hij_auto, 2), 1, dcov_qq, 0) -
1127  0.5/(psi4*psi2) * contract(contract(dcov_hij_auto, 2,
1128  hij, 1), 1, dcov_qq, 0) ;
1129 
1130  source_Sij += 0.5/(psi4*psi2) * contract(contract(hij, 0,
1131  dcov_hij_auto, 2), 0, dcov_qq, 0) ;
1132 
1133  source_Sij += 1./(3.*psi4*psi2)*contract(gtilde_con, 0, 1,
1134  qq*dcovdcov_lnq_auto + qq*dcov_lnq_auto*dcov_lnq, 0, 1)
1135  *gtilde_con ;
1136 
1137  source_Sij += 1./(3.*psi4*psi2) * contract(hdirac, 0,
1138  dcov_qq, 0)*hij_auto ;
1139 
1140  // Source terms for Rij
1141  //---------------------
1142 
1143  source_Rij = contract(hij, 0, 1, dcov_hij_auto.derive_cov(flat),
1144  2, 3) ;
1145  source_Rij.inc_dzpuis() ;
1146 
1147 
1148  source_Rij += - contract(hij_auto, 1, dcov_hdirac, 1) -
1149  contract(dcov_hdirac, 1, hij_auto, 1) ;
1150 
1151  source_Rij += contract(hdirac, 0, dcov_hij_auto, 2) ;
1152 
1153  source_Rij += - contract(contract(dcov_hij_auto, 1, dcov_hij, 2),
1154  1, 3) ;
1155 
1156  source_Rij += - contract(gtilde_cov, 0, 1, contract(contract(
1157  gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
1158 
1159  source_Rij += contract(contract(contract(contract(gtilde_cov, 0,
1160  dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
1161  contract(contract(contract(contract(gtilde_cov, 0,
1162  dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
1163 
1164  source_Rij += 0.5 * contract(gtilde_con*gtilde_con, 1, 3,
1165  contract(dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
1166 
1167  source_Rij = source_Rij * 0.5 ;
1168 
1169  for(int i=1; i<=3; i++)
1170  for(int j=1; j<=i; j++) {
1171 
1172  source_tot_hij = source_1(i,j) + source_1(j,i)
1173  + source_2(i,j) + 2.*psi4/nn * (
1174  source_4(i,j) - source_Sij(i,j))
1175  - 2.* source_Rij(i,j) +
1176  source_5(i,j) + source_5(j,i) + source_6(i,j) ;
1177  source_tot_hij.dec_dzpuis() ;
1178 
1179  source3 = 2.*psi4/nn * (source_3a(i,j) + source_3a(j,i) +
1180  source_3b(i,j)) ;
1181 
1182  source_tot_hij = source_tot_hij + source3 ;
1183 
1184  //source_tot_hij.inc_dzpuis() ;
1185 
1186  cout << "source_mat" << endl
1187  << norme((- 2. * qpig * nn * ( psi4 * stress_euler
1188  - 0.33333333333333333 * s_euler * gtilde_con ))
1189  (i,j))/(nr*nt*np) << endl ;
1190  cout << "max source_mat" << endl
1191  << max((- 2. * qpig * nn * ( psi4 * stress_euler
1192  - 0.33333333333333333 * s_euler * gtilde_con ))
1193  (i,j)) << endl ;
1194 
1195  cout << "source1" << endl
1196  << norme(source_1(i,j)/(nr*nt*np)) << endl ;
1197  cout << "max source1" << endl
1198  << max(source_1(i,j)) << endl ;
1199  cout << "source2" << endl
1200  << norme(source_2(i,j)/(nr*nt*np)) << endl ;
1201  cout << "max source2" << endl
1202  << max(source_2(i,j)) << endl ;
1203  cout << "source3a" << endl
1204  << norme(source_3a(i,j)/(nr*nt*np)) << endl ;
1205  cout << "max source3a" << endl
1206  << max(source_3a(i,j)) << endl ;
1207  cout << "source3b" << endl
1208  << norme(source_3b(i,j)/(nr*nt*np)) << endl ;
1209  cout << "max source3b" << endl
1210  << max(source_3b(i,j)) << endl ;
1211  cout << "source4" << endl
1212  << norme(source_4(i,j)/(nr*nt*np)) << endl ;
1213  cout << "max source4" << endl
1214  << max(source_4(i,j)) << endl ;
1215  cout << "source5" << endl
1216  << norme(source_5(i,j)/(nr*nt*np)) << endl ;
1217  cout << "max source5" << endl
1218  << max(source_5(i,j)) << endl ;
1219  cout << "source6" << endl
1220  << norme(source_6(i,j)/(nr*nt*np)) << endl ;
1221  cout << "max source6" << endl
1222  << max(source_6(i,j)) << endl ;
1223  cout << "source_Rij" << endl
1224  << norme(source_Rij(i,j)/(nr*nt*np)) << endl ;
1225  cout << "max source_Rij" << endl
1226  << max(source_Rij(i,j)) << endl ;
1227  cout << "source_Sij" << endl
1228  << norme(source_Sij(i,j)/(nr*nt*np)) << endl ;
1229  cout << "max source_Sij" << endl
1230  << max(source_Sij(i,j)) << endl ;
1231  cout << "source_tot" << endl
1232  << norme(source_tot_hij/(nr*nt*np)) << endl ;
1233  cout << "max source_tot" << endl
1234  << max(source_tot_hij) << endl ;
1235 
1236  // Resolution of the Poisson equations and
1237  // Check: has the Poisson equation been correctly solved ?
1238  // -----------------------------------------------------
1239 
1240  if(i==1 && j==1) {
1241 
1242  source_tot_hij.poisson(par_h11, hij_auto.set(1,1)) ;
1243 
1244  Scalar laplac (hij_auto(1,1).laplacian()) ;
1245  laplac.dec_dzpuis() ;
1246  Tbl tdiff_h11 = diffrel(laplac, source_tot_hij) ;
1247  cout << "Relative error in the resolution of the equation for "
1248  << "h11_auto : " << endl ;
1249  for (int l=0; l<nz; l++) {
1250  cout << tdiff_h11(l) << " " ;
1251  }
1252  cout << endl ;
1253  diff_h11 = max(abs(tdiff_h11)) ;
1254  }
1255 
1256  if(i==2 && j==1) {
1257 
1258  source_tot_hij.poisson(par_h21, hij_auto.set(2,1)) ;
1259 
1260  Scalar laplac (hij_auto(2,1).laplacian()) ;
1261  laplac.dec_dzpuis() ;
1262  Tbl tdiff_h21 = diffrel(laplac, source_tot_hij) ;
1263 
1264  cout <<
1265  "Relative error in the resolution of the equation for "
1266  << "h21_auto : " << endl ;
1267  for (int l=0; l<nz; l++) {
1268  cout << tdiff_h21(l) << " " ;
1269  }
1270  cout << endl ;
1271  diff_h21 = max(abs(tdiff_h21)) ;
1272  }
1273 
1274  if(i==3 && j==1) {
1275 
1276  source_tot_hij.poisson(par_h31, hij_auto.set(3,1)) ;
1277 
1278  Scalar laplac (hij_auto(3,1).laplacian()) ;
1279  laplac.dec_dzpuis() ;
1280  Tbl tdiff_h31 = diffrel(laplac, source_tot_hij) ;
1281 
1282  cout <<
1283  "Relative error in the resolution of the equation for "
1284  << "h31_auto : " << endl ;
1285  for (int l=0; l<nz; l++) {
1286  cout << tdiff_h31(l) << " " ;
1287  }
1288  cout << endl ;
1289  diff_h31 = max(abs(tdiff_h31)) ;
1290  }
1291 
1292  if(i==2 && j==2) {
1293 
1294  source_tot_hij.poisson(par_h22, hij_auto.set(2,2)) ;
1295 
1296  Scalar laplac (hij_auto(2,2).laplacian()) ;
1297  laplac.dec_dzpuis() ;
1298  Tbl tdiff_h22 = diffrel(laplac, source_tot_hij) ;
1299 
1300  cout <<
1301  "Relative error in the resolution of the equation for "
1302  << "h22_auto : " << endl ;
1303  for (int l=0; l<nz; l++) {
1304  cout << tdiff_h22(l) << " " ;
1305  }
1306  cout << endl ;
1307  diff_h22 = max(abs(tdiff_h22)) ;
1308  }
1309 
1310  if(i==3 && j==2) {
1311 
1312  source_tot_hij.poisson(par_h32, hij_auto.set(3,2)) ;
1313 
1314  Scalar laplac (hij_auto(3,2).laplacian()) ;
1315  laplac.dec_dzpuis() ;
1316  Tbl tdiff_h32 = diffrel(laplac, source_tot_hij) ;
1317 
1318  cout <<
1319  "Relative error in the resolution of the equation for "
1320  << "h32_auto : " << endl ;
1321  for (int l=0; l<nz; l++) {
1322  cout << tdiff_h32(l) << " " ;
1323  }
1324  cout << endl ;
1325  diff_h32 = max(abs(tdiff_h32)) ;
1326  }
1327 
1328  if(i==3 && j==3) {
1329 
1330  source_tot_hij.poisson(par_h33, hij_auto.set(3,3)) ;
1331 
1332  Scalar laplac (hij_auto(3,3).laplacian()) ;
1333  laplac.dec_dzpuis() ;
1334  Tbl tdiff_h33 = diffrel(laplac, source_tot_hij) ;
1335 
1336  cout <<
1337  "Relative error in the resolution of the equation for "
1338  << "h33_auto : " << endl ;
1339  for (int l=0; l<nz; l++) {
1340  cout << tdiff_h33(l) << " " ;
1341  }
1342  cout << endl ;
1343  diff_h33 = max(abs(tdiff_h33)) ;
1344  }
1345 
1346  }
1347 
1348  cout << "Tenseur hij auto in cartesian coordinates" << endl ;
1349  for (int i=1; i<=3; i++)
1350  for (int j=1; j<=i; j++) {
1351  cout << " Comp. " << i << " " << j << " : " ;
1352  for (int l=0; l<nz; l++){
1353  cout << norme(hij_auto(i,j)/(nr*nt*np))(l) << " " ;
1354  }
1355  cout << endl ;
1356  }
1357  cout << endl ;
1358 
1359  ssjm1_h11 = ssjm1h11 ;
1360  ssjm1_h21 = ssjm1h21 ;
1361  ssjm1_h31 = ssjm1h31 ;
1362  ssjm1_h22 = ssjm1h22 ;
1363  ssjm1_h32 = ssjm1h32 ;
1364  ssjm1_h33 = ssjm1h33 ;
1365 
1366  }
1367 
1368  // End of relativistic equations
1369 
1370 
1371  //-------------------------------------------------
1372  // Relative change in enthalpy
1373  //-------------------------------------------------
1374 
1375  Tbl diff_ent_tbl = diffrel( ent, ent_jm1 ) ;
1376  diff_ent = diff_ent_tbl(0) ;
1377  for (int l=1; l<nzet; l++) {
1378  diff_ent += diff_ent_tbl(l) ;
1379  }
1380  diff_ent /= nzet ;
1381 
1382 
1383  ent_jm1 = ent ;
1384 
1385 
1386  } // End of main loop
1387 
1388  //=========================================================================
1389  // End of iteration
1390  //=========================================================================
1391 
1392 }
1393 
1394 
1395 
1396 
1397 
1398 
1399 
1400 }
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Definition: star.h:555
Coord xa
Absolute x coordinate.
Definition: map.h:730
bool irrotational
true for an irrotational star, false for a corotating one
Definition: star.h:491
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1142
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Radial mapping of rather general form.
Definition: map.h:2752
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
Map & mp
Mapping associated with the star.
Definition: star.h:180
Scalar ssjm1_lnq
Effective source at the previous step for the resolution of the Poisson equation for lnq_auto...
Definition: star.h:628
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition: star.h:588
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Scalar ssjm1_h31
Effective source at the previous step for the resolution of the Poisson equation for h20_auto...
Definition: star.h:651
Lorene prototypes.
Definition: app_hor.h:64
Scalar ssjm1_h11
Effective source at the previous step for the resolution of the Poisson equation for h00_auto...
Definition: star.h:641
Standard units of space, time and mass.
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition: star.h:212
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Scalar ssjm1_h33
Effective source at the previous step for the resolution of the Poisson equation for h22_auto...
Definition: star.h:666
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Scalar logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:532
Basic integer array class.
Definition: itbl.h:122
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: star.h:512
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:784
Scalar kcar_auto
Part of the scalar generated by beta_auto, i.e.
Definition: star.h:612
void annule_hard()
Sets the Cmp to zero in a hard way.
Definition: cmp.C:338
Scalar ent
Log-enthalpy.
Definition: star.h:190
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Vector beta
Shift vector.
Definition: star.h:228
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:134
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 logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:527
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
void inc_dzpuis()
dzpuis += 1 ;
Definition: tenseur.C:1117
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition: connection.h:271
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Scalar ssjm1_h21
Effective source at the previous step for the resolution of the Poisson equation for h10_auto...
Definition: star.h:646
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:637
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:808
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition: star.h:581
Scalar pot_centri
Centrifugal potential.
Definition: star.h:521
Scalar press
Fluid pressure.
Definition: star.h:194
Class Connection.
Definition: connection.h:113
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) ...
Definition: star.h:562
Scalar lnq_auto
Scalar field generated principally by the star.
Definition: star.h:543
void equilibrium(double ent_c, int mermax, int mermax_potvit, int mermax_poisson, double relax_poisson, double relax_potvit, double thres_adapt, Tbl &diff, double om)
Computes an equilibrium configuration.
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition: scalar_manip.C:315
Parameter storage.
Definition: param.h:125
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:522
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:905
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Sym_tensor tkij_comp
Part of the extrinsic curvature tensor generated by beta_comp.
Definition: star.h:606
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition: star_global.C:186
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Tensor handling.
Definition: tensor.h:288
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Definition: star.h:535
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition: scalar.C:890
Metric gtilde
Conformal metric .
Definition: star.h:565
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:816
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Definition: star.h:538
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
Coord ya
Absolute y coordinate.
Definition: map.h:731
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
Multi-domain grid.
Definition: grilles.h:273
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:721
virtual void reevaluate_symy(const Map *mp_prev, int nzet, Cmp &uu) const =0
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Scalar nn
Lapse function N .
Definition: star.h:225
Scalar ssjm1_h22
Effective source at the previous step for the resolution of the Poisson equation for h11_auto...
Definition: star.h:656
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition: star.h:600
Coord za
Absolute z coordinate.
Definition: map.h:732
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
bool conf_flat
true if the 3-metric is conformally flat, false for a more general metric.
Definition: star.h:681
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
Scalar ssjm1_h32
Effective source at the previous step for the resolution of the Poisson equation for h21_auto...
Definition: star.h:661
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:570
Scalar kcar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Definition: star.h:618
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi...
Definition: star.h:634
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition: scalar_deriv.C:390
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: sym_tensor.C:360
Scalar psi4
Conformal factor .
Definition: star.h:552
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
Definition: star.h:557
Coord r
r coordinate centered on the grid
Definition: map.h:718
Scalar ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto...
Definition: star.h:623