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 {
143void 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
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
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
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
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
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
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
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.
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
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
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.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
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].
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
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 .
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 ...
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.