27 char star_bin_equilibrium_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_equilibrium.C,v 1.29 2014/10/13 08:53:38 j_novak Exp $" ;
132 #include "metrique.h"
135 #include "graphique.h"
136 #include "utilitaires.h"
144 int mermax_poisson,
double relax_poisson,
145 double relax_potvit,
double thres_adapt,
146 Tbl& diff,
double om) {
163 int i_b = mg->
get_nr(l_b) - 1 ;
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) ;
200 int nz_search =
nzet ;
203 double precis_secant = 1.e-14 ;
205 double reg_map = 1. ;
210 par_adapt.
add_int(nitermax, 0) ;
215 par_adapt.
add_int(nz_search, 2) ;
217 par_adapt.
add_int(adapt_flag, 3) ;
237 par_adapt.
add_tbl(ent_limit, 0) ;
261 double precis_poisson = 1.e-16 ;
268 par_logn.
add_int(mermax_poisson, 0) ;
279 par_lnq.
add_int(mermax_poisson, 0) ;
290 par_beta2.
add_int(mermax_poisson, 0) ;
298 for (
int i=0; i<3; i++) {
299 ssjm1wbeta.
set(i) =
Cmp(ssjm1_wbeta(i+1)) ;
310 par_h11.
add_int(mermax_poisson, 0) ;
321 par_h21.
add_int(mermax_poisson, 0) ;
332 par_h31.
add_int(mermax_poisson, 0) ;
343 par_h22.
add_int(mermax_poisson, 0) ;
354 par_h32.
add_int(mermax_poisson, 0) ;
365 par_h33.
add_int(mermax_poisson, 0) ;
395 for(
int mer=0 ; mer<mermax ; mer++ ) {
397 cout <<
"-----------------------------------------------" << endl ;
398 cout <<
"step: " << mer << endl ;
399 cout <<
"diff_ent = " << diff_ent << endl ;
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++) {
436 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b) /
438 ( logn_auto_b - logn_auto_c ) ;
440 if (alpha_r2_jk > alpha_r2) {
441 alpha_r2 = alpha_r2_jk ;
449 alpha_r =
sqrt(alpha_r2) ;
451 cout <<
"k_b, j_b, alpha_r: " << k_b <<
" " << j_b <<
" "
472 ent = (ent_c + logn_auto_c + pot_ext_c) -
logn_auto - pot_ext ;
474 cout <<
"pot" <<
norme(pot_ext) << endl ;
487 double rap_dent = fabs( dent_eq / dent_pole ) ;
488 cout <<
"| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
490 if ( rap_dent < thres_adapt ) {
492 cout <<
"******* FROZEN MAPPING *********" << endl ;
500 for (
int l=0; l<
nzet; l++) {
503 ent_limit.
set(
nzet-1) = ent_b ;
516 double ray_eqq =
ray_eq() ;
518 double new_rr_out_2 = (separation - ray_eqq) * 0.95 ;
519 double new_rr_out_3 = (separation + ray_eqq_pi) * 1.05 ;
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.) ;
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) ;
530 for (
int dd=4; dd<=nz-2; dd++) {
532 mp.
val_r(dd, 1., M_PI/2, 0.)) ;
537 cout <<
"too small number of domains" << endl ;
550 double ent_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++) {
556 if (xx > ent_s_max) {
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 ;
641 double lambda_dirac = 0. ;
680 0, dcov_logn_auto, 0,
true) ;
682 source4 = -
contract(
hij, 0, 1, dcovdcov_logn_auto +
685 source5 = -
contract(hdirac, 0, dcov_logn_auto, 0) ;
687 source_tot = source1 + source2 + source3 + source4 + source5 ;
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 ;
709 cout <<
"logn_auto" << endl <<
norme(
logn_auto/(nr*nt*np)) << endl ;
716 "Relative error in the resolution of the equation for logn_auto : "
718 for (
int l=0; l<nz; l++) {
719 cout << tdiff_logn(l) <<
" " ;
722 diff_logn =
max(
abs(tdiff_logn)) ;
735 source3 = -
contract(dcon_lnq, 0, dcov_lnq_auto, 0,
true) ;
738 dcov_phi_auto + dcov_logn_auto, 0,
true) ;
741 dcov_hij_auto, 0, 1, dcov_gtilde, 0, 1), 0, 1) ;
744 dcov_hij_auto, 0, 1, dcov_gtilde, 0, 2), 0, 1) ;
746 source7 = -
contract(
hij, 0, 1, dcovdcov_lnq_auto + dcov_lnq_auto *
749 source8 = - 0.25 *
contract(dcov_hdirac_auto, 0, 1)
750 -
contract(hdirac, 0, dcov_lnq_auto, 0) ;
752 source_tot = source1 + source2 + source3 + source4 + source5 +
753 source6 + source7 + source8 ;
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 ;
782 cout <<
"lnq_auto" << endl <<
norme(
lnq_auto/(nr*nt*np)) << endl ;
789 "Relative error in the resolution of the equation for lnq : "
791 for (
int l=0; l<nz; l++) {
792 cout << tdiff_lnq(l) <<
" " ;
795 diff_lnq =
max(
abs(tdiff_lnq)) ;
812 source1_beta = (4.*qpig) *
nn %
psi4
821 source4_beta = -
contract(
hij, 0, 1, dcovdcov_beta_auto, 1, 2) ;
824 dcovdcov_beta_auto, 0, 1), 0) ;
830 + 2./3. * hdirac * divbeta_auto
831 -
contract(hdirac, 0, dcov_beta_auto, 1) ;
833 source_beta = source1_beta + source2_beta + source3_beta
834 + source4_beta + source5_beta + source6_beta + source7_beta ;
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 ;
875 for (
int i=1; i<=3; i++) {
876 if (source_beta(i).get_etat() != ETATZERO)
880 for (
int i=1; i<=3; i++) {
881 if(source_beta(i).dz_nonzero()) {
882 assert( source_beta(i).get_dzpuis() == 4 ) ;
885 (source_beta.
set(i)).set_dzpuis(4) ;
889 double lambda = double(1) / double(3) ;
893 for (
int i=0; i<3; i++) {
894 source_p.
set(i) =
Cmp(source_beta(i+1)) ;
899 for (
int i=0; i<3 ;i++){
900 vect_auxi.
set(i) = 0. ;
909 for (
int i=0; i<3 ;i++)
918 for (
int i=1; i<=3; i++)
922 for (
int i=0; i<3; i++){
923 ssjm1_wbeta.
set(i+1) = ssjm1wbeta(i) ;
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)) ;
946 "Relative error in the resolution of the equation for beta_auto : "
948 cout <<
"x component : " ;
949 for (
int l=0; l<nz; l++) {
950 cout << tdiff_beta_x(l) <<
" " ;
953 cout <<
"y component : " ;
954 for (
int l=0; l<nz; l++) {
955 cout << tdiff_beta_y(l) <<
" " ;
958 cout <<
"z component : " ;
959 for (
int l=0; l<nz; l++) {
960 cout << tdiff_beta_z(l) <<
" " ;
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)) ;
995 source_1 =
contract(dcon_hij, 1, dcov_lnq_auto, 0) ;
997 source_2 = -
contract(dcon_hij, 2, dcov_lnq_auto, 0)
1008 double r0 =
mp.
val_r(nz-2, 1, 0, 0) ;
1009 double sigma = 1.*r0 ;
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() ;
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. ;
1059 omdsdp.
set(1) = - om * yya * ff ;
1060 omdsdp.
set(2) = om * xxa * ff ;
1064 omdsdp.
set(1) = om * yya * ff ;
1065 omdsdp.
set(2) = - om * xxa * ff ;
1082 source_5 = dcon_hdirac_auto ;
1098 source_Sij += -
nn / (3.*
psi4) * gtilde_con *
1100 dcov_gtilde, 0, 1), 0, 1)
1102 dcov_gtilde, 0, 2), 0, 1)) ;
1104 source_Sij += - 8.*
nn / (3.*
psi4) * gtilde_con *
1110 source_Sij += tens_temp ;
1112 source_Sij += - 8./(3.*
psi4) *
contract(dcov_phi_auto, 0,
1119 - 0.33333333333333333 *
s_euler * gtilde_con ) ;
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) ;
1126 dcov_hij_auto, 2), 1, dcov_qq, 0) -
1128 hij, 1), 1, dcov_qq, 0) ;
1131 dcov_hij_auto, 2), 0, dcov_qq, 0) ;
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)
1151 source_Rij +=
contract(hdirac, 0, dcov_hij_auto, 2) ;
1157 gtilde_con, 0, dcov_hij_auto, 2), 0, dcov_hij, 2), 1, 3) ;
1160 dcov_hij_auto, 1), 2, gtilde_con, 1), 0, dcov_hij, 1), 0, 3) +
1162 dcov_hij_auto, 1), 0, dcov_hij, 1), 0, 3), 0, gtilde_con, 1) ;
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) ;
1167 source_Rij = source_Rij * 0.5 ;
1169 for(
int i=1; i<=3; i++)
1170 for(
int j=1; j<=i; j++) {
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) ;
1179 source3 = 2.*
psi4/
nn * (source_3a(i,j) + source_3a(j,i) +
1182 source_tot_hij = source_tot_hij + source3 ;
1186 cout <<
"source_mat" << endl
1188 - 0.33333333333333333 *
s_euler * gtilde_con ))
1189 (i,j))/(nr*nt*np) << endl ;
1190 cout <<
"max source_mat" << endl
1192 - 0.33333333333333333 *
s_euler * gtilde_con ))
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 ;
1242 source_tot_hij.poisson(par_h11,
hij_auto.
set(1,1)) ;
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) <<
" " ;
1253 diff_h11 =
max(
abs(tdiff_h11)) ;
1258 source_tot_hij.poisson(par_h21,
hij_auto.
set(2,1)) ;
1262 Tbl tdiff_h21 =
diffrel(laplac, source_tot_hij) ;
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) <<
" " ;
1271 diff_h21 =
max(
abs(tdiff_h21)) ;
1276 source_tot_hij.poisson(par_h31,
hij_auto.
set(3,1)) ;
1280 Tbl tdiff_h31 =
diffrel(laplac, source_tot_hij) ;
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) <<
" " ;
1289 diff_h31 =
max(
abs(tdiff_h31)) ;
1294 source_tot_hij.poisson(par_h22,
hij_auto.
set(2,2)) ;
1298 Tbl tdiff_h22 =
diffrel(laplac, source_tot_hij) ;
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) <<
" " ;
1307 diff_h22 =
max(
abs(tdiff_h22)) ;
1312 source_tot_hij.poisson(par_h32,
hij_auto.
set(3,2)) ;
1316 Tbl tdiff_h32 =
diffrel(laplac, source_tot_hij) ;
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) <<
" " ;
1325 diff_h32 =
max(
abs(tdiff_h32)) ;
1330 source_tot_hij.poisson(par_h33,
hij_auto.
set(3,3)) ;
1334 Tbl tdiff_h33 =
diffrel(laplac, source_tot_hij) ;
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) <<
" " ;
1343 diff_h33 =
max(
abs(tdiff_h33)) ;
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++){
1376 diff_ent = diff_ent_tbl(0) ;
1377 for (
int l=1; l<
nzet; l++) {
1378 diff_ent += diff_ent_tbl(l) ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void annule_hard()
Sets the Cmp to zero in a hard way.
Tbl & set(int l)
Read/write of the value in a given domain.
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Basic integer array class.
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Radial mapping of rather general form.
virtual void homothetie(double lambda)
Sets a new radial scale.
Coord ya
Absolute y coordinate.
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Coord r
r coordinate centered on the grid
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
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.
double get_ori_x() const
Returns the x coordinate of the origin.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord xa
Absolute x coordinate.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
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.
Tensor field of valence 0 (or component of a tensorial field).
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.
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.
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 .
void annule_hard()
Sets the Scalar to zero in a hard way.
Valeur & set_spectral_va()
Returns va (read/write version)
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.
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.
Scalar lnq_auto
Scalar field generated principally by the star.
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for khi.
Vector dcon_logn
Contravariant derivative of the total logarithm of the lapse.
Vector dcon_phi
Contravariant derivative of the logarithm of the conformal factor.
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.
Scalar ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto.
Vector dcov_logn
Covariant derivative of the total logarithm of the lapse.
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
bool irrotational
true for an irrotational star, false for a corotating one
Scalar kcar_comp
Part of the scalar generated by beta_auto and beta_comp, i.e.
Scalar logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Scalar ssjm1_h21
Effective source at the previous step for the resolution of the Poisson equation for h10_auto.
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.
Scalar ssjm1_h33
Effective source at the previous step for the resolution of the Poisson equation for h22_auto.
Scalar kcar_auto
Part of the scalar generated by beta_auto, i.e.
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Scalar psi4
Conformal factor .
Scalar ssjm1_h31
Effective source at the previous step for the resolution of the Poisson equation for h20_auto.
Sym_tensor tkij_auto
Part of the extrinsic curvature tensor generated by beta_auto.
Vector dcov_phi
Covariant derivative of the logarithm of the conformal factor.
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Scalar ssjm1_h22
Effective source at the previous step for the resolution of the Poisson equation for h11_auto.
Scalar ssjm1_h32
Effective source at the previous step for the resolution of the Poisson equation for h21_auto.
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 .
Scalar pot_centri
Centrifugal potential.
Scalar ssjm1_lnq
Effective source at the previous step for the resolution of the Poisson equation for lnq_auto.
bool conf_flat
true if the 3-metric is conformally flat, false for a more general metric.
Scalar logn
Logarithm of the lapse N .
Scalar nn
Lapse function N .
Scalar ener_euler
Total energy density in the Eulerian frame.
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
double ray_eq() const
Coordinate radius at , [r_unit].
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Scalar press
Fluid pressure.
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Map & mp
Mapping associated with the star.
int nzet
Number of domains of *mp occupied by the star.
double ray_pole() const
Coordinate radius at [r_unit].
Class intended to describe valence-2 symmetric tensors.
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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).
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.
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.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Scalar & set(int)
Read/write access to a component.
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 norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
Cmp abs(const Cmp &)
Absolute value.
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...
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
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_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 .
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Standard units of space, time and mass.