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
27char 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
142namespace Lorene {
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
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
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
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
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 ;
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) {
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) ;
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
439
440 if (alpha_r2_jk > alpha_r2) {
442 k_b = k ;
443 j_b = j ;
444 }
445
446 }
447 }
448
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
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
473 cout.precision(8) ;
474 cout << "pot" << norme(pot_ext) << endl ;
475
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) ;
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) ;
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) ;
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
591
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 ;
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
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) ;
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
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) ;
645
647 dcov_hdirac.inc_dzpuis() ;
649 dcov_hdirac_auto.inc_dzpuis() ;
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
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
684
685 source5 = - contract(hdirac, 0, dcov_logn_auto, 0) ;
686
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) ;
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
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 ;
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
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
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) ;
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
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 ;
796
797 //--------------------------------------------------------
798 // Vector Poisson equation for beta_auto
799 //--------------------------------------------------------
800
801 // Source
802 //--------
803
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
830 + 2./3. * hdirac * divbeta_auto
831 - contract(hdirac, 0, dcov_beta_auto, 1) ;
832
833 source_beta = source1_beta + source2_beta + source3_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 }
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
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
967
968
969 if (!conf_flat) {
970
971 //--------------------------------------------------------
972 // Poisson equation for hij
973 //--------------------------------------------------------
974
975
976 // Declaration of all sources
977 //---------------------------
978
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
996
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
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
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
1074
1075 // Source 4
1076 // ---------
1077
1079 source_4.inc_dzpuis() ;
1080 source_4 += - 2./3. * beta.divergence(flat) * tkij_auto ;
1081
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 *
1106
1107 tens_temp = nn / (3.*psi4) * hdirac.divergence(flat)*hij_auto ;
1108 tens_temp.inc_dzpuis() ;
1109
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,
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,
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
1149 contract(dcov_hdirac, 1, hij_auto, 1) ;
1150
1151 source_Rij += contract(hdirac, 0, dcov_hij_auto, 2) ;
1152
1154 1, 3) ;
1155
1157 gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
1158
1160 dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
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
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
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() ;
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() ;
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() ;
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() ;
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() ;
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() ;
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 ;
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
1365
1366 }
1367
1368 // End of relativistic equations
1369
1370
1371 //-------------------------------------------------
1372 // Relative change in enthalpy
1373 //-------------------------------------------------
1374
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
Class Connection.
Definition connection.h:113
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Radial mapping of rather general form.
Definition map.h:2752
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 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
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
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
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
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
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
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 .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
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].
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].
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].
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
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
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 ...
Tensor field of valence 1.
Definition vector.h:188
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
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
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
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.