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) ;