184 double precision,
dimension(nsavehi) :: tsave_log, tsave_dat, tsave_slice, &
185 tsave_collapsed, tsave_custom
186 double precision :: dtsave_log, dtsave_dat, dtsave_slice, &
187 dtsave_collapsed, dtsave_custom
188 double precision :: tsavestart_log, tsavestart_dat, tsavestart_slice, &
189 tsavestart_collapsed, tsavestart_custom
190 double precision :: sizeuniformpart^D
191 double precision :: im_delta,im_nu,rka54,rka51,rkb54,rka55
192 double precision :: dx_vec(^ND)
193 integer :: ditsave_log, ditsave_dat, ditsave_slice, &
194 ditsave_collapsed, ditsave_custom
195 integer :: windex, ipower
196 integer :: i, j, k, ifile, io_state
197 integer :: iB, isave, iw, level, idim, islice
198 integer :: nx_vec(^ND), block_nx_vec(^ND)
199 integer :: my_unit, iostate
201 logical :: fileopen, file_exists
204 character(len=80) :: fmt_string
205 character(len=std_len) :: err_msg, restart_from_file_arg
206 character(len=std_len) :: basename_full, basename_prev, dummy_file
207 character(len=std_len),
dimension(:),
allocatable :: &
208 typeboundary_min^D, typeboundary_max^D
209 character(len=std_len),
allocatable :: limiter(:)
210 character(len=std_len),
allocatable :: gradient_limiter(:)
211 character(len=name_len) :: stretch_dim(ndim)
214 character(len=std_len) :: typesourcesplit
216 character(len=std_len),
allocatable :: flux_scheme(:)
218 character(len=std_len) :: typeboundspeed
220 character(len=std_len) :: time_stepper
222 character(len=std_len) :: time_integrator
224 character(len=std_len) :: typecurl
226 character(len=std_len) :: typeprolonglimit
232 character(len=std_len) :: typecourant
246 tsave_log, tsave_dat, tsave_slice, tsave_collapsed, tsave_custom, &
247 dtsave_log, dtsave_dat, dtsave_slice, dtsave_collapsed, dtsave_custom, &
248 ditsave_log, ditsave_dat, ditsave_slice, ditsave_collapsed, ditsave_custom,&
249 tsavestart_log, tsavestart_dat, tsavestart_slice, tsavestart_collapsed,&
255 namelist /methodlist/ time_stepper,time_integrator, &
314 allocate(typeboundary_min^d(nwfluxbc))
315 allocate(typeboundary_max^d(nwfluxbc))
378 typeprolonglimit =
'default'
416 tsave(isave,ifile) = bigdouble
417 itsave(isave,ifile) = biginteger
426 tsave_log = bigdouble
427 tsave_dat = bigdouble
428 tsave_slice = bigdouble
429 tsave_collapsed = bigdouble
430 tsave_custom = bigdouble
432 dtsave_log = bigdouble
433 dtsave_dat = bigdouble
434 dtsave_slice = bigdouble
435 dtsave_collapsed = bigdouble
436 dtsave_custom = bigdouble
438 ditsave_log = biginteger
439 ditsave_dat = biginteger
440 ditsave_slice = biginteger
441 ditsave_collapsed = biginteger
442 ditsave_custom = biginteger
444 tsavestart_log = bigdouble
445 tsavestart_dat = bigdouble
446 tsavestart_slice = bigdouble
447 tsavestart_collapsed = bigdouble
448 tsavestart_custom = bigdouble
469 typecourant =
'maxsum'
479 typeboundspeed =
'Einfeldt'
481 time_stepper =
'twostep'
482 time_integrator =
'default'
514 flux_scheme(level) =
'tvdlf'
516 limiter(level) =
'minmod'
517 gradient_limiter(level) =
'minmod'
522 typesourcesplit =
'sfs'
549 if(i==j.or.j==k.or.k==i)
then
551 else if(i+1==j.or.i-2==j)
then
568 inquire(file=trim(
par_files(i)), exist=file_exists)
570 if (.not. file_exists)
then
571 write(err_msg, *)
"The parameter file " // trim(
par_files(i)) // &
581 read(
unitpar, filelist,
end=101)
584 read(unitpar, savelist,
end=102)
587 read(unitpar, stoplist,
end=103)
590 read(unitpar, methodlist,
end=104)
593 read(unitpar, boundlist,
end=105)
596 read(unitpar, meshlist,
end=106)
599 read(unitpar, paramlist,
end=107)
602 read(unitpar, emissionlist,
end=108)
607 if (base_filename /= basename_prev) &
608 basename_full = trim(basename_full) // trim(base_filename)
609 basename_prev = base_filename
612 base_filename = basename_full
616 dummy_file = trim(base_filename)//
"DUMMY"
617 open(newunit=my_unit, file=trim(dummy_file), iostat=iostate)
618 if (iostate /= 0)
then
619 call mpistop(
"Can't write to output directory (" // &
620 trim(base_filename) //
")")
622 close(my_unit, status=
'delete')
626 if(source_split_usr) any_source_split=.true.
629 if(restart_from_file_arg /= undefined) &
630 restart_from_file=restart_from_file_arg
634 if(restart_from_file == undefined)
then
637 do index_latest_data = 9999, 0, -1
643 if(.not.file_exists) index_latest_data=-1
649 call mpi_bcast(index_latest_data, 1, mpi_integer, 0, icomm, ierrmpi)
651 if (resume_previous_run)
then
652 if (index_latest_data == -1)
then
653 if(mype==0)
write(*,*)
"No snapshots found to resume from, start a new run..."
656 write(restart_from_file,
"(a,i4.4,a)") trim(base_filename),index_latest_data,
".dat"
660 if (restart_from_file == undefined)
then
665 call mpistop(
"Please restart from a snapshot when firstprocess=T")
667 call mpistop(
'Change convert to .false. for a new run!')
670 if (small_pressure < 0.d0)
call mpistop(
"small_pressure should be positive.")
671 if (small_density < 0.d0)
call mpistop(
"small_density should be positive.")
673 if (small_temperature>0.d0) small_pressure=small_density*small_temperature
675 if(convert) autoconvert=.false.
677 where (tsave_log < bigdouble) tsave(:, 1) = tsave_log
678 where (tsave_dat < bigdouble) tsave(:, 2) = tsave_dat
679 where (tsave_slice < bigdouble) tsave(:, 3) = tsave_slice
680 where (tsave_collapsed < bigdouble) tsave(:, 4) = tsave_collapsed
681 where (tsave_custom < bigdouble) tsave(:, 5) = tsave_custom
683 if (dtsave_log < bigdouble) dtsave(1) = dtsave_log
684 if (dtsave_dat < bigdouble) dtsave(2) = dtsave_dat
685 if (dtsave_slice < bigdouble) dtsave(3) = dtsave_slice
686 if (dtsave_collapsed < bigdouble) dtsave(4) = dtsave_collapsed
687 if (dtsave_custom < bigdouble) dtsave(5) = dtsave_custom
689 if (tsavestart_log < bigdouble) tsavestart(1) = tsavestart_log
690 if (tsavestart_dat < bigdouble) tsavestart(2) = tsavestart_dat
691 if (tsavestart_slice < bigdouble) tsavestart(3) = tsavestart_slice
692 if (tsavestart_collapsed < bigdouble) tsavestart(4) = tsavestart_collapsed
693 if (tsavestart_custom < bigdouble) tsavestart(5) = tsavestart_custom
695 if (ditsave_log < bigdouble) ditsave(1) = ditsave_log
696 if (ditsave_dat < bigdouble) ditsave(2) = ditsave_dat
697 if (ditsave_slice < bigdouble) ditsave(3) = ditsave_slice
698 if (ditsave_collapsed < bigdouble) ditsave(4) = ditsave_collapsed
699 if (ditsave_custom < bigdouble) ditsave(5) = ditsave_custom
701 if (wall_time_max < bigdouble) wall_time_max=wall_time_max*3600.d0
704 write(unitterm, *)
''
705 write(unitterm, *)
'Output type | tsavestart | dtsave | ditsave | itsave(1) | tsave(1)'
706 write(fmt_string, *)
'(A12," | ",E9.3E2," | ",E9.3E2," | ",I6," | "'//&
711 if (mype == 0)
write(unitterm, fmt_string) trim(output_names(ifile)), &
712 tsavestart(ifile), dtsave(ifile), ditsave(ifile), itsave(1, ifile), tsave(1, ifile)
715 if (mype == 0)
write(unitterm, *)
''
718 if(slicedir(islice) > ndim) &
719 write(uniterr,*)
'Warning in read_par_files: ', &
720 'Slice ', islice,
' direction',slicedir(islice),
'larger than ndim=',ndim
721 if(slicedir(islice) < 1) &
722 write(uniterr,*)
'Warning in read_par_files: ', &
723 'Slice ', islice,
' direction',slicedir(islice),
'too small, should be [',1,ndim,
']'
726 if(it_max==biginteger .and. time_max==bigdouble.and.mype==0)
write(uniterr,*) &
727 'Warning in read_par_files: it_max or time_max not given!'
729 select case (typecourant)
731 type_courant=type_maxsum
733 type_courant=type_summax
735 type_courant=type_minimum
737 write(unitterm,*)
'Unknown typecourant=',typecourant
738 call mpistop(
"Error from read_par_files: no such typecourant!")
742 select case (flux_scheme(level))
744 flux_method(level)=fs_hll
746 flux_method(level)=fs_hllc
748 flux_method(level)=fs_hlld
750 flux_method(level)=fs_hllcd
752 flux_method(level)=fs_tvdlf
754 flux_method(level)=fs_tvdmu
756 flux_method(level)=fs_tvd
758 flux_method(level)=fs_cd
760 flux_method(level)=fs_cd4
762 flux_method(level)=fs_fd
764 flux_method(level)=fs_source
766 flux_method(level)=fs_nul
768 call mpistop(
"unkown or bad flux scheme")
770 if(flux_scheme(level)==
'tvd'.and.time_stepper/=
'onestep') &
771 call mpistop(
" tvd is onestep method, reset time_stepper='onestep'")
772 if(flux_scheme(level)==
'tvd')
then
773 if(mype==0.and.(.not.dimsplit))
write(unitterm,*) &
774 'Warning: setting dimsplit=T for tvd, as used for level=',level
777 if(flux_scheme(level)==
'hlld'.and.physics_type/=
'mhd' .and. physics_type/=
'twofl') &
778 call mpistop(
"Cannot use hlld flux if not using MHD or 2FL only charges physics!")
780 if(flux_scheme(level)==
'hllc'.and.physics_type==
'mf') &
781 call mpistop(
"Cannot use hllc flux if using magnetofriction physics!")
783 if(flux_scheme(level)==
'tvd'.and.physics_type==
'mf') &
784 call mpistop(
"Cannot use tvd flux if using magnetofriction physics!")
786 if(flux_scheme(level)==
'tvdmu'.and.physics_type==
'mf') &
787 call mpistop(
"Cannot use tvdmu flux if using magnetofriction physics!")
789 if (typepred1(level)==0)
then
790 select case (flux_scheme(level))
792 typepred1(level)=fs_cd
794 typepred1(level)=fs_cd4
796 typepred1(level)=fs_fd
797 case (
'tvdlf',
'tvdmu')
798 typepred1(level)=fs_hancock
800 typepred1(level)=fs_hll
802 typepred1(level)=fs_hllc
804 typepred1(level)=fs_hllcd
806 typepred1(level)=fs_hlld
807 case (
'nul',
'source',
'tvd')
808 typepred1(level)=fs_nul
810 call mpistop(
"No default predictor for this full step")
816 if(any(flux_scheme==
'fd')) need_global_cmax=.true.
817 if(any(limiter==
'schmid1')) need_global_a2max=.true.
820 select case (typecurl)
826 type_curl=stokesbased
828 write(unitterm,*)
"typecurl=",typecurl
829 call mpistop(
"unkown type of curl operator in read_par_files")
833 select case (time_stepper)
837 if (time_integrator==
'default')
then
838 time_integrator=
"Forward_Euler"
840 select case (time_integrator)
841 case (
"Forward_Euler")
842 t_integrator=forward_euler
844 t_integrator=imex_euler
848 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
849 call mpistop(
"unkown onestep time_integrator in read_par_files")
851 use_imex_scheme=(t_integrator==imex_euler.or.t_integrator==imex_sp)
855 if (time_integrator==
'default')
then
856 time_integrator=
"Predictor_Corrector"
858 select case (time_integrator)
859 case (
"Predictor_Corrector")
860 t_integrator=predictor_corrector
865 case (
"IMEX_Midpoint")
866 t_integrator=imex_midpoint
867 case (
"IMEX_Trapezoidal")
868 t_integrator=imex_trapezoidal
870 t_integrator=imex_222
872 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
873 call mpistop(
"unkown twostep time_integrator in read_par_files")
875 use_imex_scheme=(t_integrator==imex_midpoint.or.t_integrator==imex_trapezoidal&
876 .or.t_integrator==imex_222)
877 if (t_integrator==rk2_alf)
then
878 if(rk2_alfa<smalldouble.or.rk2_alfa>one)
call mpistop(
"set rk2_alfa within [0,1]")
880 rk_b2=1.0d0/(2.0d0*rk2_alfa)
886 if (time_integrator==
'default')
then
887 time_integrator=
'ssprk3'
889 select case (time_integrator)
895 t_integrator=imex_ars3
897 t_integrator=imex_232
899 t_integrator=imex_cb3a
901 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
902 call mpistop(
"unkown threestep time_integrator in read_par_files")
904 if(t_integrator==rk3_bt)
then
905 select case(rk3_switch)
935 call mpistop(
"Unknown rk3_switch")
938 rk3_b3=1.0d0-rk3_b1-rk3_b2
940 rk3_c3=rk3_a31+rk3_a32
942 if(t_integrator==ssprk3)
then
943 select case(ssprk_order)
946 rk_beta22=1.0d0/4.0d0
947 rk_beta33=2.0d0/3.0d0
948 rk_alfa21=3.0d0/4.0d0
949 rk_alfa31=1.0d0/3.0d0
953 rk_beta11=1.0d0/2.0d0
954 rk_beta22=1.0d0/2.0d0
955 rk_beta33=1.0d0/3.0d0
957 rk_alfa31=1.0d0/3.0d0
961 call mpistop(
"Unknown ssprk3_order")
963 rk_alfa22=1.0d0-rk_alfa21
964 rk_alfa33=1.0d0-rk_alfa31
966 if(t_integrator==imex_ars3)
then
967 ars_gamma=(3.0d0+dsqrt(3.0d0))/6.0d0
969 if(t_integrator==imex_232)
then
970 select case(imex_switch)
972 im_delta=1.0d0-1.0d0/dsqrt(2.0d0)
973 im_nu=(3.0d0+2.0d0*dsqrt(2.0d0))/6.0d0
974 imex_a21=2.0d0*im_delta
977 imex_b1=1.0d0/(2.0d0*dsqrt(2.0d0))
978 imex_b2=1.0d0/(2.0d0*dsqrt(2.0d0))
983 imex_a21=0.711664700366941d0
984 imex_a31=0.077338168947683d0
985 imex_a32=0.917273367886007d0
986 imex_b1=0.398930808264688d0
987 imex_b2=0.345755244189623d0
988 imex_ha21=0.353842865099275d0
989 imex_ha22=0.353842865099275d0
991 call mpistop(
"Unknown imex_siwtch")
994 imex_c3=imex_a31+imex_a32
995 imex_b3=1.0d0-imex_b1-imex_b2
997 if(t_integrator==imex_cb3a)
then
998 imex_c2 = 0.8925502329346865
1001 imex_c3 = imex_c2 / (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0)
1003 imex_b2 = (3.0d0*imex_c2 - 1.0d0) / (6.0d0*imex_c2**2)
1004 imex_b3 = (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) / (6.0d0*imex_c2**2)
1005 imex_a33 = (1.0d0/6.0d0 - imex_b2*imex_c2**2 - imex_b3*imex_c2*imex_c3) / (imex_b3*(imex_c3-imex_c2))
1006 imex_a32 = imex_c3 - imex_a33
1020 use_imex_scheme=(t_integrator==imex_ars3.or.t_integrator==imex_232.or.t_integrator==imex_cb3a)
1024 if (time_integrator==
'default')
then
1025 time_integrator=
"ssprk4"
1027 select case (time_integrator)
1033 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
1034 call mpistop(
"unkown fourstep time_integrator in read_par_files")
1036 if(t_integrator==ssprk4)
then
1037 select case(ssprk_order)
1039 rk_beta11=1.0d0/2.0d0
1040 rk_beta22=1.0d0/2.0d0
1041 rk_beta33=1.0d0/6.0d0
1042 rk_beta44=1.0d0/2.0d0
1044 rk_alfa31=2.0d0/3.0d0
1050 rk_beta11=1.0d0/3.0d0
1051 rk_beta22=1.0d0/3.0d0
1052 rk_beta33=1.0d0/3.0d0
1053 rk_beta44=1.0d0/4.0d0
1056 rk_alfa41=1.0d0/4.0d0
1061 call mpistop(
"Unknown ssprk_order")
1063 rk_alfa22=1.0d0-rk_alfa21
1064 rk_alfa33=1.0d0-rk_alfa31
1065 rk_alfa44=1.0d0-rk_alfa41
1070 if (time_integrator==
'default')
then
1071 time_integrator=
"ssprk5"
1073 select case (time_integrator)
1077 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
1078 call mpistop(
"unkown fivestep time_integrator in read_par_files")
1080 if(t_integrator==ssprk5)
then
1081 select case(ssprk_order)
1084 rk_beta11=0.391752226571890d0
1085 rk_beta22=0.368410593050371d0
1086 rk_beta33=0.251891774271694d0
1087 rk_beta44=0.544974750228521d0
1088 rk_beta54=0.063692468666290d0
1089 rk_beta55=0.226007483236906d0
1090 rk_alfa21=0.444370493651235d0
1091 rk_alfa31=0.620101851488403d0
1092 rk_alfa41=0.178079954393132d0
1093 rk_alfa53=0.517231671970585d0
1094 rk_alfa54=0.096059710526147d0
1096 rk_alfa22=0.555629506348765d0
1097 rk_alfa33=0.379898148511597d0
1098 rk_alfa44=0.821920045606868d0
1099 rk_alfa55=0.386708617503269d0
1100 rk_alfa22=1.0d0-rk_alfa21
1101 rk_alfa33=1.0d0-rk_alfa31
1102 rk_alfa44=1.0d0-rk_alfa41
1103 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1105 rk_c3=rk_alfa22*rk_c2+rk_beta22
1106 rk_c4=rk_alfa33*rk_c3+rk_beta33
1107 rk_c5=rk_alfa44*rk_c4+rk_beta44
1109 rk_beta11=0.39175222700392d0
1110 rk_beta22=0.36841059262959d0
1111 rk_beta33=0.25189177424738d0
1112 rk_beta44=0.54497475021237d0
1115 rk_alfa21=0.44437049406734d0
1116 rk_alfa31=0.62010185138540d0
1117 rk_alfa41=0.17807995410773d0
1118 rk_alfa53=0.51723167208978d0
1121 rk_alfa22=1.0d0-rk_alfa21
1122 rk_alfa33=1.0d0-rk_alfa31
1123 rk_alfa44=1.0d0-rk_alfa41
1125 rka51=0.00683325884039d0
1126 rka54=0.12759831133288d0
1127 rkb54=0.08460416338212d0
1128 rk_beta54=rkb54-rk_beta44*rka51/rk_alfa41
1129 rk_alfa54=rka54-rk_alfa44*rka51/rk_alfa41
1131 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1133 rk_c3=rk_alfa22*rk_c2+rk_beta22
1134 rk_c4=rk_alfa33*rk_c3+rk_beta33
1135 rk_c5=rk_alfa44*rk_c4+rk_beta44
1136 rk_beta55=1.0d0-rk_beta54-rk_alfa53*rk_c3-rk_alfa54*rk_c4-rk_alfa55*rk_c5
1138 call mpistop(
"Unknown ssprk_order")
1147 use_imex_scheme=.false.
1149 call mpistop(
"Unknown time_stepper in read_par_files")
1153 select case (stretch_dim(i))
1154 case (undefined,
'none')
1155 stretch_type(i) = stretch_none
1156 stretched_dim(i) = .false.
1157 case (
'uni',
'uniform')
1158 stretch_type(i) = stretch_uni
1159 stretched_dim(i) = .true.
1160 case (
'symm',
'symmetric')
1161 stretch_type(i) = stretch_symm
1162 stretched_dim(i) = .true.
1164 stretch_type(i) = stretch_none
1165 stretched_dim(i) = .false.
1166 if (mype == 0) print *,
'Got stretch_type = ', stretch_type(i)
1167 call mpistop(
'Unknown stretch type')
1172 if(typedimsplit ==
'default'.and. dimsplit) typedimsplit=
'xyyx'
1173 if(typedimsplit ==
'default'.and..not.dimsplit) typedimsplit=
'unsplit'
1174 dimsplit = typedimsplit /=
'unsplit'
1177 select case (typesourcesplit)
1179 sourcesplit=sourcesplit_sfs
1181 sourcesplit=sourcesplit_sf
1183 sourcesplit=sourcesplit_ssf
1185 sourcesplit=sourcesplit_ssfss
1187 write(unitterm,*)
'No such typesourcesplit=',typesourcesplit
1188 call mpistop(
"Error: Unknown typesourcesplit!")
1191 if(coordinate==-1)
then
1192 coordinate=cartesian
1194 write(*,*)
'Warning: coordinate system is not specified!'
1195 write(*,*)
'call set_coordinate_system in usr_init in mod_usr.t'
1196 write(*,*)
'Now use Cartesian coordinate'
1200 if(coordinate==cartesian)
then
1203 if(any(stretched_dim))
then
1204 coordinate=cartesian_stretched
1205 slab_uniform=.false.
1209 slab_uniform=.false.
1212 if(coordinate==spherical)
then
1214 if(mype==0)print *,
'Warning: spherical symmetry needs dimsplit=F, resetting'
1219 if (ndim==1) dimsplit=.false.
1222 select case(typeprolonglimit)
1237 allocate(type_limiter(nlevelshi))
1238 allocate(type_gradient_limiter(nlevelshi))
1240 do level=1,nlevelshi
1241 type_limiter(level) = limiter_type(limiter(level))
1242 type_gradient_limiter(level) = limiter_type(gradient_limiter(level))
1245 if (any(limiter(1:nlevelshi)==
'ppm')&
1246 .and.(flatsh.and.physics_type==
'rho'))
then
1247 call mpistop(
" PPM with flatsh=.true. can not be used with physics_type='rho'!")
1253 select case(typeboundary_min^d(iw))
1255 typeboundary(iw,2*^d-1)=bc_special
1257 typeboundary(iw,2*^d-1)=bc_cont
1259 typeboundary(iw,2*^d-1)=bc_symm
1261 typeboundary(iw,2*^d-1)=bc_asymm
1263 typeboundary(iw,2*^d-1)=bc_periodic
1265 typeboundary(iw,2*^d-1)=bc_aperiodic
1267 typeboundary(iw,2*^d-1)=bc_noinflow
1269 typeboundary(iw,2*^d-1)=12
1271 typeboundary(iw,2*^d-1)=bc_data
1273 typeboundary(iw,2*^d-1)=bc_icarus
1275 typeboundary(iw,2*^d-1)=bc_character
1277 write (unitterm,*)
"Undefined boundarytype found in read_par_files", &
1278 typeboundary_min^d(iw),
"for variable iw=",iw,
" and side iB=",2*^d-1
1282 select case(typeboundary_max^d(iw))
1284 typeboundary(iw,2*^d)=bc_special
1286 typeboundary(iw,2*^d)=bc_cont
1288 typeboundary(iw,2*^d)=bc_symm
1290 typeboundary(iw,2*^d)=bc_asymm
1292 typeboundary(iw,2*^d)=bc_periodic
1294 typeboundary(iw,2*^d)=bc_aperiodic
1296 typeboundary(iw,2*^d)=bc_noinflow
1298 typeboundary(iw,2*^d)=12
1300 typeboundary(iw,2*^d)=bc_data
1302 typeboundary(iw,2*^d-1)=bc_icarus
1303 case(
"bc_character")
1304 typeboundary(iw,2*^d)=bc_character
1306 write (unitterm,*)
"Undefined boundarytype found in read_par_files", &
1307 typeboundary_max^d(iw),
"for variable iw=",iw,
" and side iB=",2*^d
1313 if (nwfluxbc<nw_recon)
then
1314 do iw=nwfluxbc+1,nw_recon
1315 typeboundary(iw,:) = typeboundary(1, :)
1320 do iw=nw_recon+1, nw_recon+nwaux
1321 typeboundary(iw,:) = typeboundary(1, :)
1325 if (any(typeboundary == 0))
then
1326 call mpistop(
"Not all boundary conditions have been defined")
1330 periodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_periodic))
1331 aperiodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_aperiodic))
1332 if (periodb(idim).or.aperiodb(idim))
then
1334 if (typeboundary(iw,2*idim-1) .ne. typeboundary(iw,2*idim)) &
1335 call mpistop(
"Wrong counterpart in periodic boundary")
1337 if (typeboundary(iw,2*idim-1) /= bc_periodic .and. &
1338 typeboundary(iw,2*idim-1) /= bc_aperiodic)
then
1339 call mpistop(
"Each dimension should either have all "//&
1340 "or no variables periodic, some can be aperiodic")
1347 if(any(typeboundary(:,2*idim-1)==12))
then
1348 if(any(typeboundary(:,2*idim-1)/=12)) typeboundary(:,2*idim-1)=12
1349 if(phys_energy)
then
1354 typeboundary(:,2*idim-1)=bc_symm
1355 if(physics_type/=
'rho')
then
1356 select case(coordinate)
1358 typeboundary(phi_+1,2*idim-1)=bc_asymm
1359 if(physics_type==
'mhd') typeboundary(ndir+windex+phi_,2*idim-1)=bc_asymm
1361 typeboundary(3:ndir+1,2*idim-1)=bc_asymm
1362 if(physics_type==
'mhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim-1)=bc_asymm
1364 call mpistop(
'Pole is in cylindrical, polar, spherical coordinates!')
1368 if(any(typeboundary(:,2*idim)==12))
then
1369 if(any(typeboundary(:,2*idim)/=12)) typeboundary(:,2*idim)=12
1370 if(phys_energy)
then
1375 typeboundary(:,2*idim)=bc_symm
1376 if(physics_type/=
'rho')
then
1377 select case(coordinate)
1379 typeboundary(phi_+1,2*idim)=bc_asymm
1380 if(physics_type==
'mhd') typeboundary(ndir+windex+phi_,2*idim)=bc_asymm
1382 typeboundary(3:ndir+1,2*idim)=bc_asymm
1383 if(physics_type==
'mhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim)=bc_asymm
1385 call mpistop(
'Pole is in cylindrical, polar, spherical coordinates!')
1392 if(.not.phys_energy)
then
1397 if(any(limiter(1:nlevelshi)==
'mp5'))
then
1398 nghostcells=max(nghostcells,3)
1401 if(any(limiter(1:nlevelshi)==
'weno5'))
then
1402 nghostcells=max(nghostcells,3)
1405 if(any(limiter(1:nlevelshi)==
'weno5nm'))
then
1406 nghostcells=max(nghostcells,3)
1409 if(any(limiter(1:nlevelshi)==
'wenoz5'))
then
1410 nghostcells=max(nghostcells,3)
1413 if(any(limiter(1:nlevelshi)==
'wenoz5nm'))
then
1414 nghostcells=max(nghostcells,3)
1417 if(any(limiter(1:nlevelshi)==
'wenozp5'))
then
1418 nghostcells=max(nghostcells,3)
1421 if(any(limiter(1:nlevelshi)==
'wenozp5nm'))
then
1422 nghostcells=max(nghostcells,3)
1425 if(any(limiter(1:nlevelshi)==
'teno5ad'))
then
1426 nghostcells=max(nghostcells,3)
1429 if(any(limiter(1:nlevelshi)==
'weno5cu6'))
then
1430 nghostcells=max(nghostcells,3)
1433 if(any(limiter(1:nlevelshi)==
'ppm'))
then
1434 if(flatsh .or. flatcd)
then
1435 nghostcells=max(nghostcells,4)
1437 nghostcells=max(nghostcells,3)
1441 if(any(limiter(1:nlevelshi)==
'weno7'))
then
1442 nghostcells=max(nghostcells,4)
1445 if(any(limiter(1:nlevelshi)==
'mpweno7'))
then
1446 nghostcells=max(nghostcells,4)
1450 nghostcells = nghostcells + phys_wider_stencil
1453 if(stagger_grid .and. refine_max_level>1 .and. mod(nghostcells,2)/=0)
then
1454 nghostcells=nghostcells+1
1457 select case (coordinate)
1460 xprob^lim^de=xprob^lim^de*two*dpi;
1465 xprob^lim^d=xprob^lim^d*two*dpi;
1471 {ixghi^d = block_nx^d + 2*nghostcells\}
1472 {ixgshi^d = ixghi^d\}
1474 nx_vec = [{domain_nx^d|, }]
1475 block_nx_vec = [{block_nx^d|, }]
1477 if (any(nx_vec < 4) .or. any(mod(nx_vec, 2) == 1)) &
1478 call mpistop(
'Grid size (domain_nx^D) has to be even and >= 4')
1480 if (any(block_nx_vec < 4) .or. any(mod(block_nx_vec, 2) == 1)) &
1481 call mpistop(
'Block size (block_nx^D) has to be even and >= 4')
1483 {
if(mod(domain_nx^d,block_nx^d)/=0) &
1484 call mpistop(
'Grid (domain_nx^D) and block (block_nx^D) must be consistent') \}
1486 if(refine_max_level>nlevelshi.or.refine_max_level<1)
then
1487 write(unitterm,*)
'Error: refine_max_level',refine_max_level,
'>nlevelshi ',nlevelshi
1488 call mpistop(
"Reset nlevelshi and recompile!")
1491 if (any(stretched_dim))
then
1492 allocate(qstretch(0:nlevelshi,1:ndim),dxfirst(0:nlevelshi,1:ndim),&
1493 dxfirst_1mq(0:nlevelshi,1:ndim),dxmid(0:nlevelshi,1:ndim))
1494 allocate(nstretchedblocks(1:nlevelshi,1:ndim))
1495 qstretch(0:nlevelshi,1:ndim)=0.0d0
1496 dxfirst(0:nlevelshi,1:ndim)=0.0d0
1497 nstretchedblocks(1:nlevelshi,1:ndim)=0
1498 {
if (stretch_type(^d) == stretch_uni)
then
1500 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble)
then
1502 write(*,*)
'stretched grid needs finite qstretch_baselevel>1'
1503 write(*,*)
'will try default value for qstretch_baselevel in dimension', ^d
1505 if(xprobmin^d>smalldouble)
then
1506 qstretch_baselevel(^d)=(xprobmax^d/xprobmin^d)**(1.d0/dble(domain_nx^d))
1508 call mpistop(
"can not set qstretch_baselevel automatically")
1511 if(mod(block_nx^d,2)==1) &
1512 call mpistop(
"stretched grid needs even block size block_nxD")
1513 if(mod(domain_nx^d/block_nx^d,2)/=0) &
1514 call mpistop(
"number level 1 blocks in D must be even")
1515 qstretch(1,^d)=qstretch_baselevel(^d)
1516 dxfirst(1,^d)=(xprobmax^d-xprobmin^d) &
1517 *(1.0d0-qstretch(1,^d))/(1.0d0-qstretch(1,^d)**domain_nx^d)
1518 qstretch(0,^d)=qstretch(1,^d)**2
1519 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1520 if(refine_max_level>1)
then
1521 do ilev=2,refine_max_level
1522 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1523 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1524 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1529 {
if(stretch_type(^d) == stretch_uni)
then
1530 write(*,*)
'Stretched dimension ', ^d
1531 write(*,*)
'Using stretched grid with qs=',qstretch(0:refine_max_level,^d)
1532 write(*,*)
' and first cell sizes=',dxfirst(0:refine_max_level,^d)
1535 {
if(stretch_type(^d) == stretch_symm)
then
1537 write(*,*)
'will apply symmetric stretch in dimension', ^d
1539 if(mod(block_nx^d,2)==1) &
1540 call mpistop(
"stretched grid needs even block size block_nxD")
1542 if(nstretchedblocks_baselevel(^d)==0) &
1543 call mpistop(
"need finite even number of stretched blocks at baselevel")
1544 if(mod(nstretchedblocks_baselevel(^d),2)==1) &
1545 call mpistop(
"need even number of stretched blocks at baselevel")
1546 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) &
1547 call mpistop(
'stretched grid needs finite qstretch_baselevel>1')
1549 ipower=(nstretchedblocks_baselevel(^d)/2)*block_nx^d
1550 if(nstretchedblocks_baselevel(^d)==domain_nx^d/block_nx^d)
then
1551 xstretch^d=0.5d0*(xprobmax^d-xprobmin^d)
1553 xstretch^d=(xprobmax^d-xprobmin^d) &
1554 /(2.0d0+dble(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d) &
1555 *(1.0d0-qstretch_baselevel(^d))/(1.0d0-qstretch_baselevel(^d)**ipower))
1557 if(xstretch^d>(xprobmax^d-xprobmin^d)*0.5d0) &
1558 call mpistop(
" stretched grid part should not exceed full domain")
1559 dxfirst(1,^d)=xstretch^d*(1.0d0-qstretch_baselevel(^d)) &
1560 /(1.0d0-qstretch_baselevel(^d)**ipower)
1561 nstretchedblocks(1,^d)=nstretchedblocks_baselevel(^d)
1562 qstretch(1,^d)=qstretch_baselevel(^d)
1563 qstretch(0,^d)=qstretch(1,^d)**2
1564 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1565 dxmid(1,^d)=dxfirst(1,^d)
1566 dxmid(0,^d)=dxfirst(1,^d)*2.0d0
1567 if(refine_max_level>1)
then
1568 do ilev=2,refine_max_level
1569 nstretchedblocks(ilev,^d)=2*nstretchedblocks(ilev-1,^d)
1570 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1571 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1572 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1573 dxmid(ilev,^d)=dxmid(ilev-1,^d)/2.0d0
1577 sizeuniformpart^d=dxfirst(1,^d) &
1578 *(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
1580 print *,
'uniform part of size=',sizeuniformpart^d
1581 print *,
'setting of domain is then=',2*xstretch^d+sizeuniformpart^d
1582 print *,
'versus=',xprobmax^d-xprobmin^d
1584 if(dabs(xprobmax^d-xprobmin^d-2*xstretch^d-sizeuniformpart^d)>smalldouble)
then
1585 call mpistop(
'mismatch in domain size!')
1588 dxfirst_1mq(0:refine_max_level,1:ndim)=dxfirst(0:refine_max_level,1:ndim) &
1589 /(1.0d0-qstretch(0:refine_max_level,1:ndim))
1592 dx_vec = [{xprobmax^d-xprobmin^d|, }] / nx_vec
1595 write(c_ndim,
'(I1)') ^nd
1596 write(unitterm,
'(A30,' // c_ndim //
'(I0," "))') &
1597 ' Domain size (cells): ', nx_vec
1598 write(unitterm,
'(A30,' // c_ndim //
'(E9.3," "))') &
1599 ' Level one dx: ', dx_vec
1602 if (any(dx_vec < smalldouble)) &
1603 call mpistop(
"Incorrect domain size (too small grid spacing)")
1607 if(sum(w_refine_weight(:))==0) w_refine_weight(1) = 1.d0
1608 if(dabs(sum(w_refine_weight(:))-1.d0)>smalldouble)
then
1609 write(unitterm,*)
"Sum of all elements in w_refine_weight be 1.d0"
1610 call mpistop(
"Reset w_refine_weight so the sum is 1.d0")
1613 select case (typeboundspeed)
1618 case(
'cmaxleftright')
1621 call mpistop(
"set typeboundspeed='Einfieldt' or 'cmaxmean' or 'cmaxleftright'")
1624 if (mype==0)
write(unitterm,
'(A30)', advance=
'no')
'Refine estimation: '
1626 select case (refine_criterion)
1628 if (mype==0)
write(unitterm,
'(A)')
"user defined"
1630 if (mype==0)
write(unitterm,
'(A)')
"relative error"
1632 if (mype==0)
write(unitterm,
'(A)')
"Lohner's original scheme"
1634 if (mype==0)
write(unitterm,
'(A)')
"Lohner's scheme"
1636 call mpistop(
"Unknown error estimator, change refine_criterion")
1639 if (tfixgrid<bigdouble/2.0d0)
then
1640 if(mype==0)print*,
'Warning, at time=',tfixgrid,
'the grid will be fixed'
1642 if (itfixgrid<biginteger/2)
then
1643 if(mype==0)print*,
'Warning, at iteration=',itfixgrid,
'the grid will be fixed'
1645 if (ditregrid>1)
then
1646 if(mype==0)print*,
'Note, Grid is reconstructed once every',ditregrid,
'iterations'
1651 select case(slicedir(islice))
1653 if(slicecoord(islice)<xprobmin^d.or.slicecoord(islice)>xprobmax^d) &
1654 write(uniterr,*)
'Warning in read_par_files: ', &
1655 'Slice ', islice,
' coordinate',slicecoord(islice),
'out of bounds for dimension ',slicedir(islice)
1661 write(unitterm,
'(A30,A,A)')
'restart_from_file: ',
' ', trim(restart_from_file)
1662 write(unitterm,
'(A30,L1)')
'converting: ', convert
1663 write(unitterm,
'(A)')
''
1666 deallocate(flux_scheme)
1820 integer,
intent(in) :: fh
1821 integer(MPI_OFFSET_KIND),
intent(out) :: offset_tree
1822 integer(MPI_OFFSET_KIND),
intent(out) :: offset_block
1824 double precision :: rbuf(ndim)
1825 double precision,
allocatable :: params(:)
1826 integer :: i, version
1827 integer :: ibuf(ndim), iw
1828 integer :: er, n_par, tmp_int
1829 integer,
dimension(MPI_STATUS_SIZE) :: st
1830 logical :: periodic(ndim)
1831 character(len=name_len),
allocatable :: var_names(:), param_names(:)
1832 character(len=name_len) :: phys_name, geom_name
1835 call mpi_file_read(fh, version, 1, mpi_integer, st, er)
1837 call mpistop(
"Incompatible file version (maybe old format?)")
1841 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1842 offset_tree = ibuf(1)
1845 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1846 offset_block = ibuf(1)
1849 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1851 if (nw /= ibuf(1))
then
1852 write(*,*)
"nw=",nw,
" and nw found in restart file=",ibuf(1)
1853 write(*,*)
"Please be aware of changes in w at restart."
1858 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1859 if (ibuf(1) /=
ndir)
then
1860 write(*,*)
"ndir in restart file = ",ibuf(1)
1861 write(*,*)
"ndir = ",
ndir
1862 call mpistop(
"reset ndir to ndir in restart file")
1866 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1867 if (ibuf(1) /= ndim)
then
1868 write(*,*)
"ndim in restart file = ",ibuf(1)
1869 write(*,*)
"ndim = ",ndim
1870 call mpistop(
"reset ndim to ndim in restart file")
1874 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1876 write(*,*)
"number of levels in restart file = ",ibuf(1)
1878 call mpistop(
"refine_max_level < num. levels in restart file")
1882 call mpi_file_read(fh,
nleafs, 1, mpi_integer, st, er)
1885 call mpi_file_read(fh,
nparents, 1, mpi_integer, st, er)
1888 call mpi_file_read(fh,
it, 1, mpi_integer, st, er)
1891 call mpi_file_read(fh,
global_time, 1, mpi_double_precision, st, er)
1894 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1895 if (maxval(abs(rbuf(1:ndim) - [ xprobmin^
d ])) > 0)
then
1896 write(*,*)
"Error: xprobmin differs from restart data: ", rbuf(1:ndim)
1897 call mpistop(
"change xprobmin^D in par file")
1901 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1902 if (maxval(abs(rbuf(1:ndim) - [ xprobmax^
d ])) > 0)
then
1903 write(*,*)
"Error: xprobmax differs from restart data: ", rbuf(1:ndim)
1904 call mpistop(
"change xprobmax^D in par file")
1908 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
1909 if (any(ibuf(1:ndim) /= [
domain_nx^
d ]))
then
1910 write(*,*)
"Error: mesh size differs from restart data: ", ibuf(1:ndim)
1911 call mpistop(
"change domain_nx^D in par file")
1915 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
1916 if (any(ibuf(1:ndim) /= [
block_nx^
d ]))
then
1917 write(*,*)
"Error: block size differs from restart data:", ibuf(1:ndim)
1918 call mpistop(
"change block_nx^D in par file")
1922 if (version > 4)
then
1923 call mpi_file_read(fh, periodic, ndim, mpi_logical, st, er)
1924 if ({periodic(^
d) .and. .not.
periodb(^
d) .or. .not.periodic(^
d) .and.
periodb(^
d)| .or. }) &
1925 call mpistop(
"change in periodicity in par file")
1927 call mpi_file_read(fh, geom_name, name_len, mpi_character, st, er)
1930 write(*,*)
"type of coordinates in data is: ", geom_name
1931 call mpistop(
"select the correct coordinates in mod_usr.t file")
1934 call mpi_file_read(fh, stagger_mark_dat, 1, mpi_logical, st, er)
1936 write(*,*)
"Warning: stagger grid flag differs from restart data:", stagger_mark_dat
1942 if (version > 3)
then
1946 call mpi_file_read(fh, var_names(iw), name_len, mpi_character, st, er)
1950 call mpi_file_read(fh, phys_name, name_len, mpi_character, st, er)
1956 call mpi_file_read(fh, n_par, 1, mpi_integer, st, er)
1957 allocate(params(n_par))
1958 allocate(param_names(n_par))
1959 call mpi_file_read(fh, params, n_par, mpi_double_precision, st, er)
1960 call mpi_file_read(fh, param_names, name_len * n_par, mpi_character, st, er)
1963 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1968 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1971 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)