29char init_data_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Isol_hor/init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $" ;
146#include "graphique.h"
149#include "utilitaires.h"
153void Isol_hor::init_data(
int bound_nn,
double lim_nn,
int bound_psi,
154 int bound_beta,
int solve_lapse,
int solve_psi,
155 int solve_shift,
double precis,
156 double relax_nn,
double relax_psi,
double relax_beta,
int niter) {
164 ofstream conv(
"resconv.d") ;
165 ofstream kss(
"kss.d") ;
166 conv <<
" # diff_nn diff_psi diff_beta " << endl ;
175 for (
int mer=0; mer<niter; mer++) {
184 double relax_init = 0.05 ;
185 double relax_speed = 0.005 ;
187 double corr = 1 - (1 - relax_init) *
exp (- relax_speed * mer) ;
193 cout <<
"nn = " << mer <<
" corr = " << corr << endl ;
197 cout <<
" relax_nn = " << relax_nn << endl ;
198 cout <<
" relax_psi = " << relax_psi << endl ;
199 cout <<
" relax_beta = " << relax_beta << endl ;
204 if (solve_lapse == 1) {
205 Valeur nn_bound (
mp.
get_mg()-> get_angu()) ;
211 nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
216 nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
221 nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
226 nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
231 nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
236 nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
241 nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
249 cout <<
"Unexpected type of boundary conditions for the lapse!"
251 <<
" bound_nn = " << bound_nn << endl ;
259 maxabs(nn_jp1.laplacian() - sou_nn,
260 "Absolute error in the resolution of the equation for N") ;
266 nn_jp1 = relax_nn * nn_jp1 + (1 - relax_nn) *
nn() ;
274 Scalar psi_jp1 (
mp) ;
275 if (solve_psi == 1) {
276 Valeur psi_bound (
mp.
get_mg()-> get_angu()) ;
282 psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
287 psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
292 psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
297 psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
302 psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
307 psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
311 cout <<
"Unexpected type of boundary conditions for psi!"
313 <<
" bound_psi = " << bound_psi << endl ;
321 maxabs(psi_jp1.laplacian() - sou_psi,
322 "Absolute error in the resolution of the equation for Psi") ;
324 psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) *
psi() ;
332 Vector beta_jp1(
beta()) ;
334 if (solve_shift == 1) {
340 source_vector = source_vector + source_reg ;
348 Valeur boundary_x (
mp.
get_mg()-> get_angu()) ;
349 Valeur boundary_y (
mp.
get_mg()-> get_angu()) ;
350 Valeur boundary_z (
mp.
get_mg()-> get_angu()) ;
352 switch (bound_beta) {
367 cout <<
"Unexpected type of boundary conditions for psi!"
369 <<
" bound_psi = " << bound_psi << endl ;
383 double precision = 1e-8 ;
384 poisson_vect_boundary(lambda, source_vector, beta_jp1, boundary_x,
385 boundary_y, boundary_z, 0, precision, 20) ;
439 source_vector.dec_dzpuis() ;
441 + lambda * beta_jp1.divergence(
ff)
442 .derive_con(
ff) - source_vector,
443 "Absolute error in the resolution of the equation for beta") ;
453 boost_vect.set(2) = 0. ;
454 boost_vect.set(3) = 0. ;
455 boost_vect.std_spectral_base() ;
457 beta_jp1 = beta_jp1 + boost_vect ;
462 boost_vect.set(2) = 0. ;
463 boost_vect.set(3) = 0. ;
464 boost_vect.std_spectral_base() ;
466 beta_jp1 = beta_jp1 + boost_vect ;
470 beta_jp1 = relax_beta * beta_jp1 + (1 - relax_beta) *
beta() ;
477 double diff_nn, diff_psi, diff_beta ;
481 if (solve_lapse == 1)
485 if (solve_shift == 1)
488 cout <<
"step = " << mer <<
" : diff_psi = " << diff_psi
489 <<
" diff_nn = " << diff_nn
490 <<
" diff_beta = " << diff_beta << endl ;
491 cout <<
"----------------------------------------------" << endl ;
492 if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
495 if (mer>0) {conv << mer <<
" " <<
log10(diff_nn) <<
" " <<
log10(diff_psi)
496 <<
" " <<
log10(diff_beta) << endl ; } ;
505 if (solve_lapse == 1)
507 if (solve_shift == 1)
510 if (solve_shift == 1)
517 gam().radial_vect(), 0, 1)) ;
518 double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
519 double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
521 Scalar aaLss (
pow(
psi(), 6) * kkss) ;
522 double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
523 double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
526 double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
527 double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
532 for (
int k=0 ; k<nnp ; k++)
533 for (
int j=0 ; j<nnt ; j++){
534 if (kkss.val_grid_point(1, k, j, 0) > max_kss)
535 max_kss = kkss.val_grid_point(1, k, j, 0) ;
536 if (kkss.val_grid_point(1, k, j, 0) < min_kss)
537 min_kss = kkss.val_grid_point(1, k, j, 0) ;
539 if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
540 max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
541 if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
542 min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
544 if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
545 max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
546 if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
547 min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
552 kss << mer <<
" " << max_kss <<
" " << min_kss <<
" " << max_aaLss <<
" " << min_aaLss
553 <<
" " << -1 * max_hh_tilde <<
" " << -1 * min_hh_tilde << endl ;
922void Isol_hor::init_data_spher(
int bound_nn,
double lim_nn,
int bound_psi,
923 int bound_beta,
int solve_lapse,
int solve_psi,
924 int solve_shift,
double precis,
925 double relax,
int niter) {
933 ofstream conv(
"resconv.d") ;
934 ofstream kss(
"kss.d") ;
935 conv <<
" # diff_nn diff_psi diff_beta " << endl ;
939 for (
int mer=0; mer<niter; mer++) {
954 Vector tem_vect (
beta() ) ;
959 Scalar dbdr (
b_tilde().dsdr() ) ;
969 Scalar psisr (
psi()) ;
971 psisr.inc_dzpuis(2) ;
977 Scalar source_psi_spher(
mp) ;
978 source_psi_spher = -1./12. *
psi4()*
psi()/(
nn() *
nn()) * (dbdr - bsr) * (dbdr - bsr) ;
982 Scalar source_nn_spher(
mp) ;
983 source_nn_spher = 2./3. *
psi4() /
nn() * (dbdr - bsr) * (dbdr - bsr)
988 Scalar source_btilde_spher(
mp) ;
990 Scalar tmp ( -1./3. * (dbdr + 2 * bsr).dsdr() ) ;
991 tmp.std_spectral_base() ;
994 source_btilde_spher = tmp + 2 * bsr2
997 Scalar source_btilde_trun(
mp) ;
999 source_btilde_trun = tmp +
1016 sourcebeta -= source_reg ;
1017 Scalar source_btilde (sourcebeta(1) ) ;
1021 Scalar mag_sou_psi ( source_psi_spher ) ;
1022 mag_sou_psi.dec_dzpuis(4) ;
1023 Scalar mag_sou_nn ( source_nn_spher ) ;
1024 mag_sou_nn.dec_dzpuis(4) ;
1025 Scalar mag_sou_btilde ( source_btilde_trun ) ;
1026 mag_sou_btilde.dec_dzpuis(4) ;
1028 Scalar diff_sou_psi ( source_psi_spher - sourcepsi) ;
1029 diff_sou_psi.dec_dzpuis(4) ;
1030 Scalar diff_sou_nn ( source_nn_spher - sourcenn) ;
1031 diff_sou_nn.dec_dzpuis(4) ;
1032 Scalar diff_sou_btilde ( source_btilde_trun - source_btilde) ;
1033 diff_sou_btilde.dec_dzpuis(4) ;
1054 bound_nn = 1 ; lim_nn = 1. ; bound_psi = 1 ; bound_beta = 1 ;
1056 double kappa_0 = 0.2 - 1. ;
1060 kappa.std_spectral_base() ;
1061 kappa.inc_dzpuis(2) ;
1068 Valeur psi_bound (
mp.
get_mg()-> get_angu()) ;
1069 Valeur nn_bound (
mp.
get_mg()-> get_angu()) ;
1070 Valeur btilde_bound (
mp.
get_mg()-> get_angu()) ;
1075 Scalar tmp_psi = -1./4. * (2 * psisr +
1076 2./3. *
psi4()/(
psi() *
nn()) * (dbdr - bsr) ) ;
1078 Scalar tmp_nn = kappa ;
1080 Scalar tmp_btilde =
nn() / (
psi() *
psi()) ;
1083 for (
int k=0 ; k<nnp ; k++)
1084 for (
int j=0 ; j<nnt ; j++){
1086 nn_bound.set(0, k, j, 0) = tmp_nn.val_grid_point(1, k, j, 0) ;
1087 btilde_bound.set(0, k, j, 0) = tmp_btilde.val_grid_point(1, k, j, 0) ;
1090 psi_bound.std_base_scal() ;
1091 nn_bound.std_base_scal() ;
1092 btilde_bound.std_base_scal() ;
1101 Scalar psi_jp1 (
mp) ;
1102 if (solve_psi == 1) {
1104 psi_jp1 = source_psi_spher.poisson_neumann(psi_bound, 0) + 1. ;
1107 maxabs(psi_jp1.laplacian() - source_psi_spher,
1108 "Absolute error in the resolution of the equation for Psi") ;
1110 psi_jp1 = relax * psi_jp1 + (1 - relax) *
psi() ;
1115 Scalar nn_jp1 (
mp) ;
1116 if (solve_lapse == 1) {
1118 nn_jp1 = source_nn_spher.poisson_dirichlet(nn_bound, 0) + 1. ;
1121 maxabs(nn_jp1.laplacian() - source_nn_spher,
1122 "Absolute error in the resolution of the equation for N") ;
1128 nn_jp1 = relax * nn_jp1 + (1 - relax) *
nn() ;
1134 Scalar btilde_jp1 (
mp) ;
1135 if (solve_shift == 1) {
1137 btilde_jp1 = source_btilde_spher.poisson_dirichlet(btilde_bound, 0) ;
1140 maxabs(btilde_jp1.laplacian() - source_btilde_spher,
1141 "Absolute error in the resolution of the equation for btilde") ;
1143 btilde_jp1 = relax * btilde_jp1 + (1 - relax) *
b_tilde() ;
1151 double diff_nn, diff_psi, diff_btilde ;
1154 diff_btilde = 1.e-16 ;
1155 if (solve_lapse == 1)
1159 if (solve_shift == 1)
1162 cout <<
"step = " << mer <<
" : diff_psi = " << diff_psi
1163 <<
" diff_nn = " << diff_nn
1164 <<
" diff_btilde = " << diff_btilde << endl ;
1165 cout <<
"----------------------------------------------" << endl ;
1166 if ((diff_psi<precis) && (diff_nn<precis) && (diff_btilde<precis))
1169 if (mer>0) {conv << mer <<
" " <<
log10(diff_nn) <<
" " <<
log10(diff_psi)
1170 <<
" " <<
log10(diff_btilde) << endl ; } ;
1179 if (solve_lapse == 1)
1181 if (solve_shift == 1)
1183 Vector beta_jp1 (btilde_jp1 *
tgam().radial_vect()) ;
1187 if (solve_shift == 1 || solve_lapse == 1)
1196 gam().radial_vect(), 0, 1)) ;
1197 double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1198 double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1200 Scalar aaLss (
pow(
psi(), 6) * kkss) ;
1201 double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1202 double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1205 double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1206 double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1210 for (
int k=0 ; k<nnp ; k++)
1211 for (
int j=0 ; j<nnt ; j++){
1212 if (kkss.val_grid_point(1, k, j, 0) > max_kss)
1213 max_kss = kkss.val_grid_point(1, k, j, 0) ;
1214 if (kkss.val_grid_point(1, k, j, 0) < min_kss)
1215 min_kss = kkss.val_grid_point(1, k, j, 0) ;
1217 if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
1218 max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1219 if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
1220 min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1222 if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
1223 max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1224 if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
1225 min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1230 kss << mer <<
" " << max_kss <<
" " << min_kss <<
" " << max_aaLss <<
" " << min_aaLss
1231 <<
" " << -1 * max_hh_tilde <<
" " << -1 * min_hh_tilde << endl ;
1241void Isol_hor::init_data_alt(
int,
double,
int,
1242 int,
int solve_lapse,
int solve_psi,
1243 int solve_shift,
double precis,
1244 double relax,
int niter) {
1252 ofstream conv(
"resconv.d") ;
1253 ofstream kss(
"kss.d") ;
1254 conv <<
" # diff_nn diff_psi diff_beta " << endl ;
1256 Scalar psi_j (
psi()) ;
1257 Scalar nn_j (
nn()) ;
1259 Scalar diffb ( btilde_j -
b_tilde() ) ;
1260 Scalar theta_j (
beta().divergence(
ff)) ;
1261 theta_j.dec_dzpuis(2) ;
1269 for (
int mer=0; mer<niter; mer++) {
1286 Scalar psisr (psi_j) ;
1288 psisr.inc_dzpuis(2) ;
1290 Scalar dchidr ( chi_j.dsdr() ) ;
1292 Scalar chisr (chi_j) ;
1294 chisr.inc_dzpuis(2) ;
1296 Scalar rdthetadr (theta_j.dsdr() ) ;
1297 rdthetadr.mult_r() ;
1298 rdthetadr.inc_dzpuis(2) ;
1300 Scalar theta_dz4 (theta_j) ;
1301 theta_dz4.inc_dzpuis(4) ;
1303 Scalar dbmb (dchidr - 2 * chisr) ;
1309 Scalar source_psi_spher(
mp) ;
1311 source_psi_spher = -1./12. * psi_j*psi_j*psi_j*psi_j*psi_j/(nn_j * nn_j)
1317 Scalar source_nn_spher(
mp) ;
1318 source_nn_spher = 2./3. * psi_j*psi_j*psi_j*psi_j/nn_j * dbmb *dbmb
1319 - 2 * psi_j.dsdr()/psi_j * nn_j.dsdr() ;
1325 double kappa_0 = 0.2 - 1. ;
1329 kappa.std_spectral_base() ;
1330 kappa.inc_dzpuis(2) ;
1335 Valeur psi_bound (
mp.
get_mg()-> get_angu()) ;
1336 Valeur nn_bound (
mp.
get_mg()-> get_angu()) ;
1343 Scalar tmp_psi = -1./4. * (2 * psisr +
1344 2./3. * psi_j*psi_j*psi_j/ nn_j * dbmb ) ;
1351 Scalar tmp_nn = kappa ;
1355 for (
int k=0 ; k<nnp ; k++)
1356 for (
int j=0 ; j<nnt ; j++){
1357 psi_bound.set(0, k, j, 0) = tmp_psi.val_grid_point(1, k, j, 0) ;
1358 nn_bound.set(0, k, j, 0) = tmp_nn.val_grid_point(1, k, j, 0) ;
1361 psi_bound.std_base_scal() ;
1362 nn_bound.std_base_scal() ;
1372 Scalar psi_jp1 (
mp) ;
1373 if (solve_psi == 1) {
1375 psi_jp1 = source_psi_spher.poisson_neumann(psi_bound, 0) + 1. ;
1378 maxabs(psi_jp1.laplacian() - source_psi_spher,
1379 "Absolute error in the resolution of the equation for Psi") ;
1381 psi_jp1 = relax * psi_jp1 + (1 - relax) * psi_j ;
1386 Scalar nn_jp1 (
mp) ;
1387 if (solve_lapse == 1) {
1389 nn_jp1 = source_nn_spher.poisson_dirichlet(nn_bound, 0) + 1. ;
1392 maxabs(nn_jp1.laplacian() - source_nn_spher,
1393 "Absolute error in the resolution of the equation for N") ;
1399 nn_jp1 = relax * nn_jp1 + (1 - relax) * nn_j ;
1406 Scalar theta_jp1 (
mp) ;
1407 Scalar chi_jp1 (
mp) ;
1409 if (solve_shift == 1) {
1412 Scalar theta_i(theta_j) ;
1413 Scalar chi_i(chi_j) ;
1417 for (
int i=0 ; i<niter ; i++) {
1429 Scalar source_theta_spher(
mp) ;
1430 source_theta_spher = (dbmb * (nn_j.dsdr()/nn_j - 6 * psi_j.dsdr()/psi_j)).dsdr() ;
1431 source_theta_spher.dec_dzpuis() ;
1435 Scalar source_chi_spher(
mp) ;
1436 source_chi_spher = 4./3. * (dchidr - 2 * chisr) * ( nn_j.dsdr()/nn_j - 6 * psi_j.dsdr()/psi_j )
1437 - 1./3. * rdthetadr + 2 * theta_dz4 ;
1440 Valeur theta_bound (
mp.
get_mg()-> get_angu()) ;
1441 Valeur chi_bound (
mp.
get_mg()-> get_angu()) ;
1447 Scalar tmp_theta = dchidr ;
1448 tmp_theta.dec_dzpuis(2) ;
1449 tmp_theta += nn_j/(psi_j*psi_j) ;
1453 Scalar tmp_chi = nn_j/(psi_j*psi_j) ;
1456 for (
int k=0 ; k<nnp ; k++)
1457 for (
int j=0 ; j<nnt ; j++){
1458 theta_bound.set(0, k, j, 0) = tmp_theta.val_grid_point(1, k, j, 0) ;
1459 chi_bound.set(0, k, j, 0) = tmp_chi.val_grid_point(1, k, j, 0) ;
1461 theta_bound.std_base_scal() ;
1462 chi_bound.std_base_scal() ;
1465 Scalar theta_ip1(
mp) ;
1466 Scalar chi_ip1(
mp) ;
1468 theta_ip1 = source_theta_spher.poisson_dirichlet(theta_bound, 0) ;
1469 chi_ip1 = source_chi_spher.poisson_dirichlet(chi_bound, 0) ;
1472 maxabs(theta_ip1.laplacian() - source_theta_spher,
1473 "Absolute error in the resolution of the equation for Theta") ;
1474 maxabs(chi_ip1.laplacian() - source_chi_spher,
1475 "Absolute error in the resolution of the equation for chi") ;
1478 theta_ip1 = relax * theta_ip1 + (1 - relax) * theta_i ;
1479 chi_ip1 = relax * chi_ip1 + (1 - relax) * chi_i ;
1482 double diff_theta_int, diff_chi_int, int_precis ;
1483 diff_theta_int = 1.e-16 ;
1484 diff_chi_int = 1.e-16 ;
1485 int_precis = 1.e-3 ;
1487 diff_theta_int =
max(
diffrel(theta_ip1, theta_i) ) ;
1488 diff_chi_int =
max(
diffrel(chi_ip1, chi_i) ) ;
1491 cout <<
"internal step = " << i
1492 <<
" diff_theta_int = " << diff_theta_int
1493 <<
" diff_chi_int = " << diff_chi_int << endl ;
1494 cout <<
"----------------------------------------------" << endl ;
1495 if ((diff_theta_int<int_precis) && (diff_chi_int<int_precis))
1497 theta_jp1 = theta_ip1 ;
1502 theta_i = theta_ip1 ;
1512 double diff_nn, diff_psi, diff_theta, diff_chi ;
1515 diff_theta = 1.e-16 ;
1518 if (solve_lapse == 1)
1522 if (solve_shift == 1)
1523 diff_theta =
max(
diffrel(theta_jp1, theta_j) ) ;
1524 if (solve_shift == 1)
1528 cout <<
"step = " << mer <<
" : diff_psi = " << diff_psi
1529 <<
" diff_nn = " << diff_nn
1530 <<
" diff_theta = " << diff_theta
1531 <<
" diff_chi = " << diff_chi << endl ;
1532 cout <<
"----------------------------------------------" << endl ;
1533 if ((diff_psi<precis) && (diff_nn<precis) && (diff_theta<precis) && (diff_chi<precis))
1536 if (mer>0) {conv << mer <<
" " <<
log10(diff_nn) <<
" " <<
log10(diff_psi)
1537 <<
" " <<
log10(diff_theta)
1538 <<
" " <<
log10(diff_chi) << endl ; } ;
1548 if (solve_lapse == 1)
1551 if (solve_shift == 1)
1553 theta_j = theta_jp1 ;
1556 Vector beta_jp1 (chi_jp1 *
tgam().radial_vect()) ;
1560 if (solve_shift == 1 || solve_lapse == 1)
1569 gam().radial_vect(), 0, 1)) ;
1570 double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1571 double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1573 Scalar aaLss (
pow(
psi(), 6) * kkss) ;
1574 double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1575 double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1578 double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1579 double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1583 for (
int k=0 ; k<nnp ; k++)
1584 for (
int j=0 ; j<nnt ; j++){
1585 if (kkss.val_grid_point(1, k, j, 0) > max_kss)
1586 max_kss = kkss.val_grid_point(1, k, j, 0) ;
1587 if (kkss.val_grid_point(1, k, j, 0) < min_kss)
1588 min_kss = kkss.val_grid_point(1, k, j, 0) ;
1590 if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
1591 max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1592 if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
1593 min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1595 if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
1596 max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1597 if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
1598 min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1603 kss << mer <<
" " << max_kss <<
" " << min_kss <<
" " << max_aaLss <<
" " << min_aaLss
1604 <<
" " << -1 * max_hh_tilde <<
" " << -1 * min_hh_tilde << endl ;
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution)
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial)
const Vector source_beta() const
Source for .
double omega
Angular velocity in LORENE's units.
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif)
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial)
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
Metric met_gamt
3 metric tilde
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook's boundary condition.
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
const Scalar source_psi() const
Source for .
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
const Valeur boundary_nn_Dir_lapl(int mer=1) const
Dirichlet boundary condition for N fixing the divergence of the connection form .
double boost_z
Boost velocity in z-direction.
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
const Valeur boundary_beta_z() const
Component z of boundary value of .
const Valeur boundary_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
Map_af & mp
Affine mapping.
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial)
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif) .
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
double boost_x
Boost velocity in x-direction.
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial)
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution)
const Scalar source_nn() const
Source for N.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
const Scalar & dsdr() const
Returns of *this .
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 const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
const Scalar & psi4() const
Factor at the current time step (jtime ).
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
int jtime
Time step index of the latest slice.
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Metric & gam() const
Induced metric at the current time step (jtime )
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Evolution_std< double > the_time
Time label of each slice.
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Cmp log10(const Cmp &)
Basis 10 logarithm.
Cmp exp(const Cmp &)
Exponential.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
void arrete(int a=0)
Setting a stop point in a code.
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Standard units of space, time and mass.