30char et_bin_bhns_extr_eq_ylm_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_eq_ylm.C,v 1.7 2014/10/13 08:52:54 j_novak Exp $" ;
68#include "et_bin_bhns_extr.h"
75#include "utilitaires.h"
112 int i_b = mg->
get_nr(l_b) - 1 ;
122 double& diff_ent = diff.
set(0) ;
123 double& diff_vel_pot = diff.
set(1) ;
124 double& diff_logn = diff.
set(2) ;
125 double& diff_beta = diff.
set(3) ;
126 double& diff_shift_x = diff.
set(4) ;
127 double& diff_shift_y = diff.
set(5) ;
128 double& diff_shift_z = diff.
set(6) ;
139 int nz_search =
nzet ;
141 double precis_secant = 1.e-14 ;
143 double reg_map = 1. ;
147 par_adapt.
add_int(nitermax, 0) ;
152 par_adapt.
add_int(nz_search, 2) ;
154 par_adapt.
add_int(adapt_flag, 3) ;
170 par_adapt.
add_tbl(ent_limit, 0) ;
176 double precis_poisson = 1.e-16 ;
180 par_poisson1.
add_int(mermax_poisson, 0) ;
192 par_poisson2.
add_int(mermax_poisson, 0) ;
202 Param par_poisson_vect ;
204 par_poisson_vect.
add_int(mermax_poisson, 0) ;
206 par_poisson_vect.
add_double(relax_poisson, 0) ;
207 par_poisson_vect.
add_double(precis_poisson, 1) ;
230 cout <<
"We don't have those ylm programmed yet!!!!" << endl ;
234 int nylm = (l_max+1) * (l_max+1) ;
236 Cmp** ylmvec =
new Cmp* [nylm] ;
242 for(
int mer=0 ; mer<mermax ; mer++ ) {
244 cout <<
"-----------------------------------------------" << endl ;
245 cout <<
"step: " << mer << endl ;
246 cout <<
"diff_ent = " << diff_ent << endl ;
255 precis_poisson, relax_potvit) ;
266 double logn_auto_c =
logn_auto()(0, 0, 0, 0) ;
267 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
273 double alpha_r2 = 0 ;
274 for (
int k=0; k<mg->
get_np(l_b); k++) {
275 for (
int j=0; j<mg->
get_nt(l_b); j++) {
277 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
278 double logn_auto_b =
logn_auto()(l_b, k, j, i_b) ;
281 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
282 / ( logn_auto_b - logn_auto_c ) ;
284 if (alpha_r2_jk > alpha_r2) {
285 alpha_r2 = alpha_r2_jk ;
293 alpha_r =
sqrt(alpha_r2) ;
295 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" "
317 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
324 double dentdx =
ent().dsdx()(0, 0, 0, 0) ;
325 double dentdy =
ent().dsdy()(0, 0, 0, 0) ;
327 cout <<
"dH/dx|_center = " << dentdx << endl ;
328 cout <<
"dH/dy|_center = " << dentdy << endl ;
330 double dec_fact = 1. ;
334 func_in.
set() = 1. - dec_fact * (dentdx/ent_c) *
mp.
x
335 - dec_fact * (dentdy/ent_c) *
mp.
y ;
347 ent.
set() =
ent() * (func_in() + func_ex()) ;
351 double dentdx_new =
ent().dsdx()(0, 0, 0, 0) ;
352 double dentdy_new =
ent().dsdy()(0, 0, 0, 0) ;
353 cout <<
"dH/dx|_new = " << dentdx_new << endl ;
354 cout <<
"dH/dy|_new = " << dentdy_new << endl ;
363 double dent_eq =
ent().dsdr().val_point(
ray_eq_pi(),M_PI/2.,M_PI) ;
364 double dent_pole =
ent().dsdr().val_point(
ray_pole(),0.,0.) ;
365 double rap_dent = fabs( dent_eq / dent_pole ) ;
366 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
368 if ( rap_dent < thres_adapt ) {
370 cout <<
"******* FROZEN MAPPING *********" << endl ;
378 for (
int l=0; l<
nzet; l++) {
379 ent_limit.
set(l) =
ent()(l, k_b, j_b, i_b) ;
382 ent_limit.
set(
nzet-1) = ent_b ;
395 double fact_resize = 1. / alpha_r ;
396 for (
int l=
nzet; l<nz-1; l++) {
397 mp_et.
resize(l, fact_resize) ;
401 double ent_s_max = -1 ;
404 for (
int k=0; k<mg->
get_np(l_b); k++) {
405 for (
int j=0; j<mg->
get_nt(l_b); j++) {
406 double xx = fabs(
ent()(l_b, k, j, i_b) ) ;
407 if (xx > ent_s_max) {
414 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1"
415 <<
" and nzet : " << endl ;
416 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max
417 <<
" and j = " << j_s_max << endl ;
426 for (
int l=0; l<=l_max; l++) {
427 for (
int m=0; m<= l; m++) {
429 int index = l*l + 2*m ;
430 if(m > 0) index -= 1 ;
431 if(index != oldindex+1) {
432 cout <<
"error!! " << l <<
" " << m <<
" " << index
433 <<
" " << oldindex << endl ;
438 if ((l+m) % 2 == 0) {
455 if(oldindex+1 != nylm) {
456 cout <<
"ERROR!!! " << oldindex <<
" "<< nylm << endl ;
490 r_bh.
set() =
pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
495 xx_cov.
set(0) = xx + sepa ;
501 xsr_cov = xx_cov / r_bh ;
505 msr = ggrav * mass / r_bh ;
509 lapse_bh = 1. /
sqrt( 1.+2.*msr ) ;
513 lapse_bh2 = 1. / (1.+2.*msr) ;
529 for (
int i=0; i<3; i++)
530 for (
int j=0; j<3; j++)
531 lltkij.
set() += xsr_cov(i) * xsr_cov(j) *
tkij_auto(i, j) ;
565 + 2.*lapse_bh2 % msr % (ldnu % ldbeta + ldldnu)
566 + lapse_bh2 % lapse_bh2 % msr % (2.*(ldnu + 4.*msr % ldnu)
568 - (4.*
a_car % lapse_bh2 % lapse_bh2 % msr / 3. /
nnn / r_bh)
569 * (2.+3.*msr) * (3.+4.*msr) * lltkij
570 + (2.*
a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
571 /
nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
572 + (4.*
pow(lapse_bh2, 3.) % msr % msr / 3. / r_bh / r_bh)
575 - 3.*(
a_car%lapse_bh/
nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
579 cout <<
"The computation of BH-NS binary systems"
580 <<
" should be relativistic !!!" << endl ;
589 double* intvec_logn =
new double[nylm] ;
595 if(nu_int[0] != -1.0) {
596 for (
int i=0; i<nylm; i++) {
597 intvec_logn[i] *= relax_ylm ;
598 intvec_logn[i] += (1.0 - relax_ylm) * nu_int[i] ;
602 source().poisson_ylm(par_poisson1,
logn_auto.
set(), nylm,
605 for (
int i=0; i<nylm; i++) {
606 nu_int[i] = intvec_logn[i] ;
609 delete [] intvec_logn ;
622 "Relative error in the resolution of the equation for logn_auto : "
625 for (
int l=0; l<nz; l++) {
626 cout << tdiff_logn(l) <<
" " ;
629 diff_logn =
max(
abs(tdiff_logn)) ;
643 + lapse_bh2 % msr % (ldnu%ldnu + ldbeta%ldbeta + 2.*ldldbeta)
644 + lapse_bh2 % lapse_bh2 % msr % (2.*(1.+4.*msr) * ldbeta
646 - (
a_car % lapse_bh2 % lapse_bh2 % msr /
nnn / r_bh)
647 * (2.+3.*msr) * (3.+4.*msr) * lltkij
648 + (2.*
a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
649 /
nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
650 + (2.*
pow(lapse_bh2, 3.) % msr % msr / r_bh / r_bh)
653 - 2.*(
a_car%lapse_bh/
nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
660 double* intvec_beta =
new double[nylm] ;
665 if(beta_int[0] != -1.0) {
666 for (
int i=0; i<nylm; i++) {
667 intvec_beta[i] *= relax_ylm ;
668 intvec_beta[i] += (1.0 - relax_ylm) * beta_int[i] ;
675 for (
int i=0; i<nylm; i++) {
676 beta_int[i] = intvec_beta[i] ;
679 delete [] intvec_beta ;
686 cout <<
"Relative error in the resolution of the equation for "
687 <<
"beta_auto : " << endl ;
688 for (
int l=0; l<nz; l++) {
689 cout << tdiff_beta(l) <<
" " ;
692 diff_beta =
max(
abs(tdiff_beta)) ;
703 xx_con.
set(0) = xx + sepa ;
709 xsr_con = xx_con / r_bh ;
751 eldn.
set(0) = ldn_cov(0) ;
752 eldn.
set(1) = ldn_cov(1) ;
753 eldn.
set(2) = ldn_cov(2) ;
757 Tenseur ldivn = xsr_con % divn ;
778 evtmp.
set(0) = vtmp(0) ;
779 evtmp.
set(1) = vtmp(1) ;
780 evtmp.
set(2) = vtmp(2) ;
794 - 2.*
nnn % lapse_bh2 * msr / r_bh
796 + 2.*lapse_bh2 * msr * (3.*ldldn + ldldivn) / 3.
797 - lapse_bh2 * msr / r_bh
798 * (4.*ldivn - lapse_bh2 % (3.*ldn_con + 8.*msr * ldn_con)
799 - (eldn + 2.*lapse_bh2*(9.+11.*msr)*lldn_cov%xsr_con) / 3.)
800 - 2.*lapse_bh2 % lapse_bh2 * msr / r_bh / r_bh
802 - lapse_bh2 * (12.+51.*msr+46.*msr*msr) * ln % xsr_con )
804 + 8.*
pow(lapse_bh2, 4.) * msr / r_bh / r_bh
805 % (lapse_ns - 1.) * (2.+10.*msr+9.*msr*msr) * xsr_con / 3.
806 + 2.*
pow(lapse_bh2, 3.) * msr / r_bh * (2.+3.*msr)
807 * ( (1.+2.*msr) * evtmp - (3.+2.*msr) * lvtmp * xsr_con) / 3. ;
816 for (
int i=0; i<3; i++) {
817 for (
int l=0; l<nz; l++) {
818 if (source_shift(i).
get_etat() != ETATZERO)
839 double lambda_shift = double(1) / double(3) ;
841 double* intvec_shift =
new double[4*nylm] ;
842 for (
int i=0; i<3; i++) {
843 double* intvec =
new double[nylm];
846 for (
int j=0; j<nylm; j++) {
847 intvec_shift[i*nylm+j] = intvec[j] ;
849 if(shift_int[i*nylm] != -1.0) {
850 for (
int j=0; j<nylm; j++) {
851 intvec_shift[i*nylm+j] *= relax_ylm ;
852 intvec_shift[i*nylm+j] += (1.0 - relax_ylm)
853 * shift_int[i*nylm+j] ;
859 Cmp soscal (source_scal()) ;
860 double* intvec2 =
new double[nylm] ;
863 for (
int j=0; j<nylm; j++) {
864 intvec_shift[3*nylm+j] = intvec2[j] ;
866 if(shift_int[3*nylm] != -1.0) {
867 for (
int i=0; i<nylm; i++) {
868 intvec_shift[3*nylm+i] *= relax_ylm ;
869 intvec_shift[3*nylm+i] += (1.0 - relax_ylm)
870 *shift_int[3*nylm+i] ;
876 source_shift.poisson_vect_ylm(lambda_shift, par_poisson_vect,
880 for (
int i=0; i<4*nylm; i++) {
881 shift_int[i] = intvec_shift[i] ;
884 delete[] intvec_shift ;
900 for (
int i=0; i<3; i++) {
902 + lambda_shift * graddivn(i) ;
905 Tbl tdiff_shift_x =
diffrel(lap_shift(0), source_shift(0)) ;
906 Tbl tdiff_shift_y =
diffrel(lap_shift(1), source_shift(1)) ;
907 Tbl tdiff_shift_z =
diffrel(lap_shift(2), source_shift(2)) ;
909 cout <<
"Relative error in the resolution of the equation for "
910 <<
"shift_auto : " << endl ;
911 cout <<
"x component : " ;
912 for (
int l=0; l<nz; l++) {
913 cout << tdiff_shift_x(l) <<
" " ;
916 cout <<
"y component : " ;
917 for (
int l=0; l<nz; l++) {
918 cout << tdiff_shift_y(l) <<
" " ;
921 cout <<
"z component : " ;
922 for (
int l=0; l<nz; l++) {
923 cout << tdiff_shift_z(l) <<
" " ;
927 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
928 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
929 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
946 diff_ent = diff_ent_tbl(0) ;
947 for (
int l=1; l<
nzet; l++) {
948 diff_ent += diff_ent_tbl(l) ;
954 for (
int i=0; i<nylm; i++) {
977 double relax_poisson,
1001 int l_b =
nzet - 1 ;
1002 int i_b = mg->
get_nr(l_b) - 1 ;
1012 double& diff_ent = diff.
set(0) ;
1013 double& diff_vel_pot = diff.
set(1) ;
1014 double& diff_logn = diff.
set(2) ;
1015 double& diff_beta = diff.
set(3) ;
1016 double& diff_shift_x = diff.
set(4) ;
1017 double& diff_shift_y = diff.
set(5) ;
1018 double& diff_shift_z = diff.
set(6) ;
1024 int nitermax = 100 ;
1026 int adapt_flag = 1 ;
1029 int nz_search =
nzet ;
1031 double precis_secant = 1.e-14 ;
1033 double reg_map = 1. ;
1037 par_adapt.
add_int(nitermax, 0) ;
1042 par_adapt.
add_int(nz_search, 2) ;
1044 par_adapt.
add_int(adapt_flag, 3) ;
1060 par_adapt.
add_tbl(ent_limit, 0) ;
1066 double precis_poisson = 1.e-16 ;
1068 Param par_poisson1 ;
1070 par_poisson1.
add_int(mermax_poisson, 0) ;
1080 Param par_poisson2 ;
1082 par_poisson2.
add_int(mermax_poisson, 0) ;
1092 Param par_poisson_vect ;
1094 par_poisson_vect.
add_int(mermax_poisson, 0) ;
1096 par_poisson_vect.
add_double(relax_poisson, 0) ;
1097 par_poisson_vect.
add_double(precis_poisson, 1) ;
1120 cout <<
"We don't have those ylm programmed yet!!!!" << endl ;
1124 int nylm = (l_max+1) * (l_max+1) ;
1126 Cmp** ylmvec =
new Cmp* [nylm] ;
1132 for(
int mer=0 ; mer<mermax ; mer++ ) {
1134 cout <<
"-----------------------------------------------" << endl ;
1135 cout <<
"step: " << mer << endl ;
1136 cout <<
"diff_ent = " << diff_ent << endl ;
1145 precis_poisson, relax_potvit) ;
1156 double logn_auto_c =
logn_auto()(0, 0, 0, 0) ;
1157 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
1163 double alpha_r2 = 0 ;
1164 for (
int k=0; k<mg->
get_np(l_b); k++) {
1165 for (
int j=0; j<mg->
get_nt(l_b); j++) {
1167 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
1168 double logn_auto_b =
logn_auto()(l_b, k, j, i_b) ;
1171 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
1172 / ( logn_auto_b - logn_auto_c ) ;
1174 if (alpha_r2_jk > alpha_r2) {
1175 alpha_r2 = alpha_r2_jk ;
1183 alpha_r =
sqrt(alpha_r2) ;
1185 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" "
1186 << alpha_r << endl ;
1207 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
1213 double dentdx =
ent().dsdx()(0, 0, 0, 0) ;
1214 double dentdy =
ent().dsdy()(0, 0, 0, 0) ;
1216 cout <<
"dH/dx|_center = " << dentdx << endl ;
1217 cout <<
"dH/dy|_center = " << dentdy << endl ;
1219 double dec_fact = 1. ;
1223 func_in.
set() = 1. - dec_fact * (dentdx/ent_c) *
mp.
x ;
1229 func_ex.
set() = 1. ;
1235 ent.
set() =
ent() * (func_in() + func_ex()) ;
1239 double dentdx_new =
ent().dsdx()(0, 0, 0, 0) ;
1241 cout <<
"dH/dx|_new = " << dentdx_new << endl ;
1250 double dent_eq =
ent().dsdr().val_point(
ray_eq_pi(),M_PI/2.,M_PI) ;
1251 double dent_pole =
ent().dsdr().val_point(
ray_pole(),0.,0.) ;
1252 double rap_dent = fabs( dent_eq / dent_pole ) ;
1253 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
1255 if ( rap_dent < thres_adapt ) {
1257 cout <<
"******* FROZEN MAPPING *********" << endl ;
1265 for (
int l=0; l<
nzet; l++) {
1266 ent_limit.
set(l) =
ent()(l, k_b, j_b, i_b) ;
1269 ent_limit.
set(
nzet-1) = ent_b ;
1282 double fact_resize = 1. / alpha_r ;
1283 for (
int l=
nzet; l<nz-1; l++) {
1284 mp_et.
resize(l, fact_resize) ;
1288 double ent_s_max = -1 ;
1291 for (
int k=0; k<mg->
get_np(l_b); k++) {
1292 for (
int j=0; j<mg->
get_nt(l_b); j++) {
1293 double xx = fabs(
ent()(l_b, k, j, i_b) ) ;
1294 if (xx > ent_s_max) {
1301 cout <<
"Max. abs(enthalpy) at the boundary between domains nzet-1"
1302 <<
" and nzet : " << endl ;
1303 cout <<
" " << ent_s_max <<
" reached for k = " << k_s_max
1304 <<
" and j = " << j_s_max << endl ;
1313 for (
int l=0; l<=l_max; l++) {
1314 for (
int m=0; m<= l; m++) {
1316 int index = l*l + 2*m ;
1317 if(m > 0) index -= 1 ;
1318 if(index != oldindex+1) {
1319 cout <<
"error!! " << l <<
" " << m <<
" " << index
1320 <<
" " << oldindex << endl ;
1325 if ((l+m) % 2 == 0) {
1333 if((l+m) % 2 == 0) {
1342 if(oldindex+1 != nylm) {
1343 cout <<
"ERROR!!! " << oldindex <<
" "<< nylm << endl ;
1377 r_bh.
set() =
pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
1382 xx_cov.
set(0) = xx + sepa ;
1383 xx_cov.
set(1) = yy ;
1384 xx_cov.
set(2) = zz ;
1388 msr = ggrav * mass / r_bh ;
1392 tmp = 1. / ( 1. - 0.25*msr*msr ) ;
1414 - 0.5 * tmp % msr % msr % xdnu / r_bh / r_bh
1415 - tmp % msr % xdbeta / r_bh / r_bh ;
1419 cout <<
"The computation of BH-NS binary systems"
1420 <<
" should be relativistic !!!" << endl ;
1429 double* intvec_logn =
new double[nylm] ;
1435 if(nu_int[0] != -1.0) {
1436 for (
int i=0; i<nylm; i++) {
1437 intvec_logn[i] *= relax_ylm ;
1438 intvec_logn[i] += (1.0 - relax_ylm) * nu_int[i] ;
1442 source().poisson_ylm(par_poisson1,
logn_auto.
set(), nylm,
1445 for (
int i=0; i<nylm; i++) {
1446 nu_int[i] = intvec_logn[i] ;
1449 delete [] intvec_logn ;
1462 "Relative error in the resolution of the equation for logn_auto : "
1465 for (
int l=0; l<nz; l++) {
1466 cout << tdiff_logn(l) <<
" " ;
1469 diff_logn =
max(
abs(tdiff_logn)) ;
1483 - tmp % msr % xdnu / r_bh / r_bh
1484 - 0.5 * tmp % msr %msr % xdbeta / r_bh / r_bh ;
1491 double* intvec_beta =
new double[nylm] ;
1496 if(beta_int[0] != -1.0) {
1497 for (
int i=0; i<nylm; i++) {
1498 intvec_beta[i] *= relax_ylm ;
1499 intvec_beta[i] += (1.0 - relax_ylm) * beta_int[i] ;
1504 nylm, intvec_beta) ;
1506 for (
int i=0; i<nylm; i++) {
1507 beta_int[i] = intvec_beta[i] ;
1510 delete [] intvec_beta ;
1517 cout <<
"Relative error in the resolution of the equation for "
1518 <<
"beta_auto : " << endl ;
1519 for (
int l=0; l<nz; l++) {
1520 cout << tdiff_beta(l) <<
" " ;
1523 diff_beta =
max(
abs(tdiff_beta)) ;
1534 bhtmp = tmp % msr % (3.*msr-8.) % xx_cov / r_bh / r_bh ;
1553 for (
int i=0; i<3; i++) {
1554 for (
int l=0; l<nz; l++) {
1555 if (source_shift(i).get_etat() != ETATZERO)
1576 double lambda_shift = double(1) / double(3) ;
1578 double* intvec_shift =
new double[4*nylm] ;
1579 for (
int i=0; i<3; i++) {
1580 double* intvec =
new double[nylm];
1583 for (
int j=0; j<nylm; j++) {
1584 intvec_shift[i*nylm+j] = intvec[j] ;
1586 if(shift_int[i*nylm] != -1.0) {
1587 for (
int j=0; j<nylm; j++) {
1588 intvec_shift[i*nylm+j] *= relax_ylm ;
1589 intvec_shift[i*nylm+j] += (1.0 - relax_ylm)
1590 * shift_int[i*nylm+j] ;
1596 Cmp soscal (source_scal()) ;
1597 double* intvec2 =
new double[nylm] ;
1600 for (
int j=0; j<nylm; j++) {
1601 intvec_shift[3*nylm+j] = intvec2[j] ;
1603 if(shift_int[3*nylm] != -1.0) {
1604 for (
int i=0; i<nylm; i++) {
1605 intvec_shift[3*nylm+i] *= relax_ylm ;
1606 intvec_shift[3*nylm+i] += (1.0 - relax_ylm)
1607 *shift_int[3*nylm+i] ;
1613 source_shift.poisson_vect_ylm(lambda_shift, par_poisson_vect,
1617 for (
int i=0; i<4*nylm; i++) {
1618 shift_int[i] = intvec_shift[i] ;
1621 delete[] intvec_shift ;
1637 for (
int i=0; i<3; i++) {
1639 + lambda_shift * graddivn(i) ;
1642 Tbl tdiff_shift_x =
diffrel(lap_shift(0), source_shift(0)) ;
1643 Tbl tdiff_shift_y =
diffrel(lap_shift(1), source_shift(1)) ;
1644 Tbl tdiff_shift_z =
diffrel(lap_shift(2), source_shift(2)) ;
1646 cout <<
"Relative error in the resolution of the equation for "
1647 <<
"shift_auto : " << endl ;
1648 cout <<
"x component : " ;
1649 for (
int l=0; l<nz; l++) {
1650 cout << tdiff_shift_x(l) <<
" " ;
1653 cout <<
"y component : " ;
1654 for (
int l=0; l<nz; l++) {
1655 cout << tdiff_shift_y(l) <<
" " ;
1658 cout <<
"z component : " ;
1659 for (
int l=0; l<nz; l++) {
1660 cout << tdiff_shift_z(l) <<
" " ;
1664 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
1665 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
1666 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
1683 diff_ent = diff_ent_tbl(0) ;
1684 for (
int l=1; l<
nzet; l++) {
1685 diff_ent += diff_ent_tbl(l) ;
1691 for (
int i=0; i<nylm; i++) {
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
void annule(int l)
Sets the Cmp to zero in a given domain.
Active physical coordinates and mapping derivatives.
void equil_bhns_extr_ylm_ks(double ent_c, const double &mass, const double &sepa, double *nu_int, double *beta_int, double *shift_int, int mermax, int mermax_poisson, double relax_poisson, double relax_ylm, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the Kerr-Schild background metric u...
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
void get_ylm(int nylm, Cmp **ylmvec) const
Constructs spherical harmonics.
void get_integrals(int nylm, double *intvec, Cmp **ylmvec, Cmp source) const
Computes multipole moments.
double velocity_pot_extr(const double &mass, const double &sepa, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
void hydro_euler_extr(const double &mass, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
void equil_bhns_extr_ylm_cf(double ent_c, const double &mass, const double &sepa, double *nu_int, double *beta_int, double *shift_int, int mermax, int mermax_poisson, double relax_poisson, double relax_ylm, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the conformally flat background met...
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
bool irrotational
true for an irrotational star, false for a corotating one
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Tenseur pot_centri
Centrifugal potential.
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
int nzet
Number of domains of *mp occupied by the star.
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Map & mp
Mapping associated with the star.
Tenseur press
Fluid pressure.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Tenseur ener_euler
Total energy density in the Eulerian frame.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Tenseur a_car
Total conformal factor .
double ray_pole() const
Coordinate radius at [r_unit].
Radial mapping of rather general form.
void resize_extr(double lambda)
Rescales the outer boundary of the outermost domain in the case of non-compactified external domain.
virtual void homothetie(double lambda)
Sets a new radial scale.
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord y
y coordinate centered on the grid
virtual void reevaluate(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 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 z
z coordinate centered on the grid
Coord x
x coordinate centered on the grid
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
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.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
void set_std_base()
Set the standard spectal basis of decomposition for each component.
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
int get_etat() const
Returns the logical state.
Cmp sqrt(const Cmp &)
Square root.
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 .
Cmp abs(const Cmp &)
Absolute value.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Standard units of space, time and mass.