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 
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 }
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
void annule_hard()
Sets the Cmp to zero in a hard way.
Definition: cmp.C:338
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
Class Connection.
Definition: connection.h:113
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
Basic integer array class.
Definition: itbl.h:122
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
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition: map_et.C:905
Coord ya
Absolute y coordinate.
Definition: map.h:731
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Coord r
r coordinate centered on the grid
Definition: map.h:718
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
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.
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Coord za
Absolute z coordinate.
Definition: map.h:732
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
Coord xa
Absolute x coordinate.
Definition: map.h:730
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.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
Multi-domain grid.
Definition: grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Parameter storage.
Definition: param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition: param.C:1004
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition: param.C:1142
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition: param.C:522
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition: scalar_manip.C:151
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
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
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
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
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
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
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
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...
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: star.h:512
Scalar lnq_auto
Scalar field generated principally by the star.
Definition: star.h:543
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi.
Definition: star.h:634
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Definition: star.h:538
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
Definition: star.h:557
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Scalar ssjm1_h11
Effective source at the previous step for the resolution of the Poisson equation for h00_auto.
Definition: star.h:641
Scalar ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto.
Definition: star.h:623
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Definition: star.h:535
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:527
bool irrotational
true for an irrotational star, false for a corotating one
Definition: star.h:491
Scalar kcar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Definition: star.h:618
Scalar logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:532
Scalar ssjm1_h21
Effective source at the previous step for the resolution of the Poisson equation for h10_auto.
Definition: star.h:646
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Sym_tensor tkij_comp
Part of the extrinsic curvature tensor generated by beta_comp.
Definition: star.h:606
Scalar ssjm1_h33
Effective source at the previous step for the resolution of the Poisson equation for h22_auto.
Definition: star.h:666
Scalar kcar_auto
Part of the scalar generated by beta_auto, i.e.
Definition: star.h:612
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:570
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition: star.h:581
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition: star.h:588
Scalar psi4
Conformal factor .
Definition: star.h:552
Scalar ssjm1_h31
Effective source at the previous step for the resolution of the Poisson equation for h20_auto.
Definition: star.h:651
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Definition: star.h:600
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Definition: star.h:555
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition: star.h:562
Scalar ssjm1_h22
Effective source at the previous step for the resolution of the Poisson equation for h11_auto.
Definition: star.h:656
Scalar ssjm1_h32
Effective source at the previous step for the resolution of the Poisson equation for h21_auto.
Definition: star.h:661
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.
Metric gtilde
Conformal metric .
Definition: star.h:565
Scalar pot_centri
Centrifugal potential.
Definition: star.h:521
Scalar ssjm1_lnq
Effective source at the previous step for the resolution of the Poisson equation for lnq_auto.
Definition: star.h:628
bool conf_flat
true if the 3-metric is conformally flat, false for a more general metric.
Definition: star.h:681
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Scalar nn
Lapse function N .
Definition: star.h:225
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:462
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition: star.h:212
Scalar press
Fluid pressure.
Definition: star.h:194
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition: star_global.C:186
Scalar ent
Log-enthalpy.
Definition: star.h:190
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
Map & mp
Mapping associated with the star.
Definition: star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Vector beta
Shift vector.
Definition: star.h:228
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
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
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Basic array class.
Definition: tbl.h:161
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void poisson_vect_oohara(double lambda, Param &par, Tenseur &shift, Tenseur &scal) const
Solves the vectorial Poisson equation .
Definition: tenseur_pde.C:218
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Tensor handling.
Definition: tensor.h:288
void smooth(int nzet, Valeur &uuva) const
Changes the function *this as a smooth one when there exists a discontinuity between the nucleus and ...
Definition: valeur_smooth.C:77
Tensor field of valence 1.
Definition: vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
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
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
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
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...
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:926
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.