183 double precision,
dimension(nsavehi) :: tsave_log, tsave_dat, tsave_slice, &
184 tsave_collapsed, tsave_custom
185 double precision :: dtsave_log, dtsave_dat, dtsave_slice, &
186 dtsave_collapsed, dtsave_custom
187 double precision :: tsavestart_log, tsavestart_dat, tsavestart_slice, &
188 tsavestart_collapsed, tsavestart_custom
189 double precision :: sizeuniformpart^D
190 double precision :: im_delta,im_nu,rka54,rka51,rkb54,rka55
191 double precision :: dx_vec(^ND)
192 integer :: ditsave_log, ditsave_dat, ditsave_slice, &
193 ditsave_collapsed, ditsave_custom
194 integer :: windex, ipower
195 integer :: i, j, k, ifile, io_state
196 integer :: iB, isave, iw, level, idim, islice
197 integer :: nx_vec(^ND), block_nx_vec(^ND)
198 integer :: my_unit, iostate
200 logical :: fileopen, file_exists
203 character(len=80) :: fmt_string
204 character(len=std_len) :: err_msg, restart_from_file_arg
205 character(len=std_len) :: basename_full, basename_prev, dummy_file
206 character(len=std_len),
dimension(:),
allocatable :: &
207 typeboundary_min^D, typeboundary_max^D
208 character(len=std_len),
allocatable :: limiter(:)
209 character(len=std_len),
allocatable :: gradient_limiter(:)
210 character(len=name_len) :: stretch_dim(ndim)
213 character(len=std_len) :: typesourcesplit
215 character(len=std_len),
allocatable :: flux_scheme(:)
217 character(len=std_len) :: typeboundspeed
219 character(len=std_len) :: time_stepper
221 character(len=std_len) :: time_integrator
223 character(len=std_len) :: typecurl
225 character(len=std_len) :: typeprolonglimit
231 character(len=std_len) :: typecourant
245 tsave_log, tsave_dat, tsave_slice, tsave_collapsed, tsave_custom, &
246 dtsave_log, dtsave_dat, dtsave_slice, dtsave_collapsed, dtsave_custom, &
247 ditsave_log, ditsave_dat, ditsave_slice, ditsave_collapsed, ditsave_custom,&
248 tsavestart_log, tsavestart_dat, tsavestart_slice, tsavestart_collapsed,&
254 namelist /methodlist/ time_stepper,time_integrator, &
313 allocate(typeboundary_min^d(nwfluxbc))
314 allocate(typeboundary_max^d(nwfluxbc))
377 typeprolonglimit =
'default'
415 tsave(isave,ifile) = bigdouble
416 itsave(isave,ifile) = biginteger
425 tsave_log = bigdouble
426 tsave_dat = bigdouble
427 tsave_slice = bigdouble
428 tsave_collapsed = bigdouble
429 tsave_custom = bigdouble
431 dtsave_log = bigdouble
432 dtsave_dat = bigdouble
433 dtsave_slice = bigdouble
434 dtsave_collapsed = bigdouble
435 dtsave_custom = bigdouble
437 ditsave_log = biginteger
438 ditsave_dat = biginteger
439 ditsave_slice = biginteger
440 ditsave_collapsed = biginteger
441 ditsave_custom = biginteger
443 tsavestart_log = bigdouble
444 tsavestart_dat = bigdouble
445 tsavestart_slice = bigdouble
446 tsavestart_collapsed = bigdouble
447 tsavestart_custom = bigdouble
468 typecourant =
'maxsum'
478 typeboundspeed =
'Einfeldt'
480 time_stepper =
'twostep'
481 time_integrator =
'default'
513 flux_scheme(level) =
'tvdlf'
515 limiter(level) =
'minmod'
516 gradient_limiter(level) =
'minmod'
521 typesourcesplit =
'sfs'
548 if(i==j.or.j==k.or.k==i)
then
550 else if(i+1==j.or.i-2==j)
then
567 inquire(file=trim(
par_files(i)), exist=file_exists)
569 if (.not. file_exists)
then
570 write(err_msg, *)
"The parameter file " // trim(
par_files(i)) // &
580 read(
unitpar, filelist,
end=101)
583 read(unitpar, savelist,
end=102)
586 read(unitpar, stoplist,
end=103)
589 read(unitpar, methodlist,
end=104)
592 read(unitpar, boundlist,
end=105)
595 read(unitpar, meshlist,
end=106)
598 read(unitpar, paramlist,
end=107)
601 read(unitpar, emissionlist,
end=108)
606 if (base_filename /= basename_prev) &
607 basename_full = trim(basename_full) // trim(base_filename)
608 basename_prev = base_filename
611 base_filename = basename_full
615 dummy_file = trim(base_filename)//
"DUMMY"
616 open(newunit=my_unit, file=trim(dummy_file), iostat=iostate)
617 if (iostate /= 0)
then
618 call mpistop(
"Can't write to output directory (" // &
619 trim(base_filename) //
")")
621 close(my_unit, status=
'delete')
625 if(source_split_usr) any_source_split=.true.
628 if(restart_from_file_arg /= undefined) &
629 restart_from_file=restart_from_file_arg
633 if(restart_from_file == undefined)
then
636 do index_latest_data = 9999, 0, -1
642 if(.not.file_exists) index_latest_data=-1
648 call mpi_bcast(index_latest_data, 1, mpi_integer, 0, icomm, ierrmpi)
650 if (resume_previous_run)
then
651 if (index_latest_data == -1)
then
652 if(mype==0)
write(*,*)
"No snapshots found to resume from, start a new run..."
655 write(restart_from_file,
"(a,i4.4,a)") trim(base_filename),index_latest_data,
".dat"
659 if (restart_from_file == undefined)
then
664 call mpistop(
"Please restart from a snapshot when firstprocess=T")
666 call mpistop(
'Change convert to .false. for a new run!')
669 if (small_pressure < 0.d0)
call mpistop(
"small_pressure should be positive.")
670 if (small_density < 0.d0)
call mpistop(
"small_density should be positive.")
672 if (small_temperature>0.d0) small_pressure=small_density*small_temperature
674 if(convert) autoconvert=.false.
676 where (tsave_log < bigdouble) tsave(:, 1) = tsave_log
677 where (tsave_dat < bigdouble) tsave(:, 2) = tsave_dat
678 where (tsave_slice < bigdouble) tsave(:, 3) = tsave_slice
679 where (tsave_collapsed < bigdouble) tsave(:, 4) = tsave_collapsed
680 where (tsave_custom < bigdouble) tsave(:, 5) = tsave_custom
682 if (dtsave_log < bigdouble) dtsave(1) = dtsave_log
683 if (dtsave_dat < bigdouble) dtsave(2) = dtsave_dat
684 if (dtsave_slice < bigdouble) dtsave(3) = dtsave_slice
685 if (dtsave_collapsed < bigdouble) dtsave(4) = dtsave_collapsed
686 if (dtsave_custom < bigdouble) dtsave(5) = dtsave_custom
688 if (tsavestart_log < bigdouble) tsavestart(1) = tsavestart_log
689 if (tsavestart_dat < bigdouble) tsavestart(2) = tsavestart_dat
690 if (tsavestart_slice < bigdouble) tsavestart(3) = tsavestart_slice
691 if (tsavestart_collapsed < bigdouble) tsavestart(4) = tsavestart_collapsed
692 if (tsavestart_custom < bigdouble) tsavestart(5) = tsavestart_custom
694 if (ditsave_log < bigdouble) ditsave(1) = ditsave_log
695 if (ditsave_dat < bigdouble) ditsave(2) = ditsave_dat
696 if (ditsave_slice < bigdouble) ditsave(3) = ditsave_slice
697 if (ditsave_collapsed < bigdouble) ditsave(4) = ditsave_collapsed
698 if (ditsave_custom < bigdouble) ditsave(5) = ditsave_custom
700 if (wall_time_max < bigdouble) wall_time_max=wall_time_max*3600.d0
703 write(unitterm, *)
''
704 write(unitterm, *)
'Output type | tsavestart | dtsave | ditsave | itsave(1) | tsave(1)'
705 write(fmt_string, *)
'(A12," | ",E9.3E2," | ",E9.3E2," | ",I6," | "'//&
710 if (mype == 0)
write(unitterm, fmt_string) trim(output_names(ifile)), &
711 tsavestart(ifile), dtsave(ifile), ditsave(ifile), itsave(1, ifile), tsave(1, ifile)
714 if (mype == 0)
write(unitterm, *)
''
717 if(slicedir(islice) > ndim) &
718 write(uniterr,*)
'Warning in read_par_files: ', &
719 'Slice ', islice,
' direction',slicedir(islice),
'larger than ndim=',ndim
720 if(slicedir(islice) < 1) &
721 write(uniterr,*)
'Warning in read_par_files: ', &
722 'Slice ', islice,
' direction',slicedir(islice),
'too small, should be [',1,ndim,
']'
725 if(it_max==biginteger .and. time_max==bigdouble.and.mype==0)
write(uniterr,*) &
726 'Warning in read_par_files: it_max or time_max not given!'
728 select case (typecourant)
730 type_courant=type_maxsum
732 type_courant=type_summax
734 type_courant=type_minimum
736 write(unitterm,*)
'Unknown typecourant=',typecourant
737 call mpistop(
"Error from read_par_files: no such typecourant!")
741 select case (flux_scheme(level))
743 flux_method(level)=fs_hll
745 flux_method(level)=fs_hllc
747 flux_method(level)=fs_hlld
749 flux_method(level)=fs_hllcd
751 flux_method(level)=fs_tvdlf
753 flux_method(level)=fs_tvdmu
755 flux_method(level)=fs_tvd
757 flux_method(level)=fs_cd
759 flux_method(level)=fs_cd4
761 flux_method(level)=fs_fd
763 flux_method(level)=fs_source
765 flux_method(level)=fs_nul
767 call mpistop(
"unkown or bad flux scheme")
769 if(flux_scheme(level)==
'tvd'.and.time_stepper/=
'onestep') &
770 call mpistop(
" tvd is onestep method, reset time_stepper='onestep'")
771 if(flux_scheme(level)==
'tvd')
then
772 if(mype==0.and.(.not.dimsplit))
write(unitterm,*) &
773 'Warning: setting dimsplit=T for tvd, as used for level=',level
776 if(flux_scheme(level)==
'hlld'.and.physics_type/=
'mhd' .and. physics_type/=
'twofl') &
777 call mpistop(
"Cannot use hlld flux if not using MHD or 2FL only charges physics!")
779 if(flux_scheme(level)==
'hllc'.and.physics_type==
'mf') &
780 call mpistop(
"Cannot use hllc flux if using magnetofriction physics!")
782 if(flux_scheme(level)==
'tvd'.and.physics_type==
'mf') &
783 call mpistop(
"Cannot use tvd flux if using magnetofriction physics!")
785 if(flux_scheme(level)==
'tvdmu'.and.physics_type==
'mf') &
786 call mpistop(
"Cannot use tvdmu flux if using magnetofriction physics!")
788 if (typepred1(level)==0)
then
789 select case (flux_scheme(level))
791 typepred1(level)=fs_cd
793 typepred1(level)=fs_cd4
795 typepred1(level)=fs_fd
796 case (
'tvdlf',
'tvdmu')
797 typepred1(level)=fs_hancock
799 typepred1(level)=fs_hll
801 typepred1(level)=fs_hllc
803 typepred1(level)=fs_hllcd
805 typepred1(level)=fs_hlld
806 case (
'nul',
'source',
'tvd')
807 typepred1(level)=fs_nul
809 call mpistop(
"No default predictor for this full step")
815 if(any(flux_scheme==
'fd')) need_global_cmax=.true.
816 if(any(limiter==
'schmid1')) need_global_a2max=.true.
819 select case (typecurl)
825 type_curl=stokesbased
827 write(unitterm,*)
"typecurl=",typecurl
828 call mpistop(
"unkown type of curl operator in read_par_files")
832 select case (time_stepper)
836 if (time_integrator==
'default')
then
837 time_integrator=
"Forward_Euler"
839 select case (time_integrator)
840 case (
"Forward_Euler")
841 t_integrator=forward_euler
843 t_integrator=imex_euler
847 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
848 call mpistop(
"unkown onestep time_integrator in read_par_files")
850 use_imex_scheme=(t_integrator==imex_euler.or.t_integrator==imex_sp)
854 if (time_integrator==
'default')
then
855 time_integrator=
"Predictor_Corrector"
857 select case (time_integrator)
858 case (
"Predictor_Corrector")
859 t_integrator=predictor_corrector
864 case (
"IMEX_Midpoint")
865 t_integrator=imex_midpoint
866 case (
"IMEX_Trapezoidal")
867 t_integrator=imex_trapezoidal
869 t_integrator=imex_222
871 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
872 call mpistop(
"unkown twostep time_integrator in read_par_files")
874 use_imex_scheme=(t_integrator==imex_midpoint.or.t_integrator==imex_trapezoidal&
875 .or.t_integrator==imex_222)
876 if (t_integrator==rk2_alf)
then
877 if(rk2_alfa<smalldouble.or.rk2_alfa>one)
call mpistop(
"set rk2_alfa within [0,1]")
879 rk_b2=1.0d0/(2.0d0*rk2_alfa)
885 if (time_integrator==
'default')
then
886 time_integrator=
'ssprk3'
888 select case (time_integrator)
894 t_integrator=imex_ars3
896 t_integrator=imex_232
898 t_integrator=imex_cb3a
900 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
901 call mpistop(
"unkown threestep time_integrator in read_par_files")
903 if(t_integrator==rk3_bt)
then
904 select case(rk3_switch)
934 call mpistop(
"Unknown rk3_switch")
937 rk3_b3=1.0d0-rk3_b1-rk3_b2
939 rk3_c3=rk3_a31+rk3_a32
941 if(t_integrator==ssprk3)
then
942 select case(ssprk_order)
945 rk_beta22=1.0d0/4.0d0
946 rk_beta33=2.0d0/3.0d0
947 rk_alfa21=3.0d0/4.0d0
948 rk_alfa31=1.0d0/3.0d0
952 rk_beta11=1.0d0/2.0d0
953 rk_beta22=1.0d0/2.0d0
954 rk_beta33=1.0d0/3.0d0
956 rk_alfa31=1.0d0/3.0d0
960 call mpistop(
"Unknown ssprk3_order")
962 rk_alfa22=1.0d0-rk_alfa21
963 rk_alfa33=1.0d0-rk_alfa31
965 if(t_integrator==imex_ars3)
then
966 ars_gamma=(3.0d0+dsqrt(3.0d0))/6.0d0
968 if(t_integrator==imex_232)
then
969 select case(imex_switch)
971 im_delta=1.0d0-1.0d0/dsqrt(2.0d0)
972 im_nu=(3.0d0+2.0d0*dsqrt(2.0d0))/6.0d0
973 imex_a21=2.0d0*im_delta
976 imex_b1=1.0d0/(2.0d0*dsqrt(2.0d0))
977 imex_b2=1.0d0/(2.0d0*dsqrt(2.0d0))
982 imex_a21=0.711664700366941d0
983 imex_a31=0.077338168947683d0
984 imex_a32=0.917273367886007d0
985 imex_b1=0.398930808264688d0
986 imex_b2=0.345755244189623d0
987 imex_ha21=0.353842865099275d0
988 imex_ha22=0.353842865099275d0
990 call mpistop(
"Unknown imex_siwtch")
993 imex_c3=imex_a31+imex_a32
994 imex_b3=1.0d0-imex_b1-imex_b2
996 if(t_integrator==imex_cb3a)
then
997 imex_c2 = 0.8925502329346865
1000 imex_c3 = imex_c2 / (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0)
1002 imex_b2 = (3.0d0*imex_c2 - 1.0d0) / (6.0d0*imex_c2**2)
1003 imex_b3 = (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) / (6.0d0*imex_c2**2)
1004 imex_a33 = (1.0d0/6.0d0 - imex_b2*imex_c2**2 - imex_b3*imex_c2*imex_c3) / (imex_b3*(imex_c3-imex_c2))
1005 imex_a32 = imex_c3 - imex_a33
1019 use_imex_scheme=(t_integrator==imex_ars3.or.t_integrator==imex_232.or.t_integrator==imex_cb3a)
1023 if (time_integrator==
'default')
then
1024 time_integrator=
"ssprk4"
1026 select case (time_integrator)
1032 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
1033 call mpistop(
"unkown fourstep time_integrator in read_par_files")
1035 if(t_integrator==ssprk4)
then
1036 select case(ssprk_order)
1038 rk_beta11=1.0d0/2.0d0
1039 rk_beta22=1.0d0/2.0d0
1040 rk_beta33=1.0d0/6.0d0
1041 rk_beta44=1.0d0/2.0d0
1043 rk_alfa31=2.0d0/3.0d0
1049 rk_beta11=1.0d0/3.0d0
1050 rk_beta22=1.0d0/3.0d0
1051 rk_beta33=1.0d0/3.0d0
1052 rk_beta44=1.0d0/4.0d0
1055 rk_alfa41=1.0d0/4.0d0
1060 call mpistop(
"Unknown ssprk_order")
1062 rk_alfa22=1.0d0-rk_alfa21
1063 rk_alfa33=1.0d0-rk_alfa31
1064 rk_alfa44=1.0d0-rk_alfa41
1069 if (time_integrator==
'default')
then
1070 time_integrator=
"ssprk5"
1072 select case (time_integrator)
1076 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
1077 call mpistop(
"unkown fivestep time_integrator in read_par_files")
1079 if(t_integrator==ssprk5)
then
1080 select case(ssprk_order)
1083 rk_beta11=0.391752226571890d0
1084 rk_beta22=0.368410593050371d0
1085 rk_beta33=0.251891774271694d0
1086 rk_beta44=0.544974750228521d0
1087 rk_beta54=0.063692468666290d0
1088 rk_beta55=0.226007483236906d0
1089 rk_alfa21=0.444370493651235d0
1090 rk_alfa31=0.620101851488403d0
1091 rk_alfa41=0.178079954393132d0
1092 rk_alfa53=0.517231671970585d0
1093 rk_alfa54=0.096059710526147d0
1095 rk_alfa22=0.555629506348765d0
1096 rk_alfa33=0.379898148511597d0
1097 rk_alfa44=0.821920045606868d0
1098 rk_alfa55=0.386708617503269d0
1099 rk_alfa22=1.0d0-rk_alfa21
1100 rk_alfa33=1.0d0-rk_alfa31
1101 rk_alfa44=1.0d0-rk_alfa41
1102 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1104 rk_c3=rk_alfa22*rk_c2+rk_beta22
1105 rk_c4=rk_alfa33*rk_c3+rk_beta33
1106 rk_c5=rk_alfa44*rk_c4+rk_beta44
1108 rk_beta11=0.39175222700392d0
1109 rk_beta22=0.36841059262959d0
1110 rk_beta33=0.25189177424738d0
1111 rk_beta44=0.54497475021237d0
1114 rk_alfa21=0.44437049406734d0
1115 rk_alfa31=0.62010185138540d0
1116 rk_alfa41=0.17807995410773d0
1117 rk_alfa53=0.51723167208978d0
1120 rk_alfa22=1.0d0-rk_alfa21
1121 rk_alfa33=1.0d0-rk_alfa31
1122 rk_alfa44=1.0d0-rk_alfa41
1124 rka51=0.00683325884039d0
1125 rka54=0.12759831133288d0
1126 rkb54=0.08460416338212d0
1127 rk_beta54=rkb54-rk_beta44*rka51/rk_alfa41
1128 rk_alfa54=rka54-rk_alfa44*rka51/rk_alfa41
1130 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1132 rk_c3=rk_alfa22*rk_c2+rk_beta22
1133 rk_c4=rk_alfa33*rk_c3+rk_beta33
1134 rk_c5=rk_alfa44*rk_c4+rk_beta44
1135 rk_beta55=1.0d0-rk_beta54-rk_alfa53*rk_c3-rk_alfa54*rk_c4-rk_alfa55*rk_c5
1137 call mpistop(
"Unknown ssprk_order")
1146 use_imex_scheme=.false.
1148 call mpistop(
"Unknown time_stepper in read_par_files")
1152 select case (stretch_dim(i))
1153 case (undefined,
'none')
1154 stretch_type(i) = stretch_none
1155 stretched_dim(i) = .false.
1156 case (
'uni',
'uniform')
1157 stretch_type(i) = stretch_uni
1158 stretched_dim(i) = .true.
1159 case (
'symm',
'symmetric')
1160 stretch_type(i) = stretch_symm
1161 stretched_dim(i) = .true.
1163 stretch_type(i) = stretch_none
1164 stretched_dim(i) = .false.
1165 if (mype == 0) print *,
'Got stretch_type = ', stretch_type(i)
1166 call mpistop(
'Unknown stretch type')
1171 if(typedimsplit ==
'default'.and. dimsplit) typedimsplit=
'xyyx'
1172 if(typedimsplit ==
'default'.and..not.dimsplit) typedimsplit=
'unsplit'
1173 dimsplit = typedimsplit /=
'unsplit'
1176 select case (typesourcesplit)
1178 sourcesplit=sourcesplit_sfs
1180 sourcesplit=sourcesplit_sf
1182 sourcesplit=sourcesplit_ssf
1184 sourcesplit=sourcesplit_ssfss
1186 write(unitterm,*)
'No such typesourcesplit=',typesourcesplit
1187 call mpistop(
"Error: Unknown typesourcesplit!")
1190 if(coordinate==-1)
then
1191 coordinate=cartesian
1193 write(*,*)
'Warning: coordinate system is not specified!'
1194 write(*,*)
'call set_coordinate_system in usr_init in mod_usr.t'
1195 write(*,*)
'Now use Cartesian coordinate'
1199 if(coordinate==cartesian)
then
1202 if(any(stretched_dim))
then
1203 coordinate=cartesian_stretched
1204 slab_uniform=.false.
1208 slab_uniform=.false.
1211 if(coordinate==spherical)
then
1213 if(mype==0)print *,
'Warning: spherical symmetry needs dimsplit=F, resetting'
1218 if (ndim==1) dimsplit=.false.
1221 select case(typeprolonglimit)
1236 allocate(type_limiter(nlevelshi))
1237 allocate(type_gradient_limiter(nlevelshi))
1239 do level=1,nlevelshi
1240 type_limiter(level) = limiter_type(limiter(level))
1241 type_gradient_limiter(level) = limiter_type(gradient_limiter(level))
1244 if (any(limiter(1:nlevelshi)==
'ppm')&
1245 .and.(flatsh.and.physics_type==
'rho'))
then
1246 call mpistop(
" PPM with flatsh=.true. can not be used with physics_type='rho'!")
1252 select case(typeboundary_min^d(iw))
1254 typeboundary(iw,2*^d-1)=bc_special
1256 typeboundary(iw,2*^d-1)=bc_cont
1258 typeboundary(iw,2*^d-1)=bc_symm
1260 typeboundary(iw,2*^d-1)=bc_asymm
1262 typeboundary(iw,2*^d-1)=bc_periodic
1264 typeboundary(iw,2*^d-1)=bc_aperiodic
1266 typeboundary(iw,2*^d-1)=bc_noinflow
1268 typeboundary(iw,2*^d-1)=12
1270 typeboundary(iw,2*^d-1)=bc_data
1272 typeboundary(iw,2*^d-1)=bc_icarus
1274 typeboundary(iw,2*^d-1)=bc_character
1276 write (unitterm,*)
"Undefined boundarytype found in read_par_files", &
1277 typeboundary_min^d(iw),
"for variable iw=",iw,
" and side iB=",2*^d-1
1281 select case(typeboundary_max^d(iw))
1283 typeboundary(iw,2*^d)=bc_special
1285 typeboundary(iw,2*^d)=bc_cont
1287 typeboundary(iw,2*^d)=bc_symm
1289 typeboundary(iw,2*^d)=bc_asymm
1291 typeboundary(iw,2*^d)=bc_periodic
1293 typeboundary(iw,2*^d)=bc_aperiodic
1295 typeboundary(iw,2*^d)=bc_noinflow
1297 typeboundary(iw,2*^d)=12
1299 typeboundary(iw,2*^d)=bc_data
1301 typeboundary(iw,2*^d-1)=bc_icarus
1302 case(
"bc_character")
1303 typeboundary(iw,2*^d)=bc_character
1305 write (unitterm,*)
"Undefined boundarytype found in read_par_files", &
1306 typeboundary_max^d(iw),
"for variable iw=",iw,
" and side iB=",2*^d
1312 if (nwfluxbc<nwflux)
then
1313 do iw=nwfluxbc+1,nwflux
1314 typeboundary(iw,:) = typeboundary(1, :)
1319 do iw=nwflux+1, nwflux+nwaux
1320 typeboundary(iw,:) = typeboundary(1, :)
1324 if (any(typeboundary == 0))
then
1325 call mpistop(
"Not all boundary conditions have been defined")
1329 periodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_periodic))
1330 aperiodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_aperiodic))
1331 if (periodb(idim).or.aperiodb(idim))
then
1333 if (typeboundary(iw,2*idim-1) .ne. typeboundary(iw,2*idim)) &
1334 call mpistop(
"Wrong counterpart in periodic boundary")
1336 if (typeboundary(iw,2*idim-1) /= bc_periodic .and. &
1337 typeboundary(iw,2*idim-1) /= bc_aperiodic)
then
1338 call mpistop(
"Each dimension should either have all "//&
1339 "or no variables periodic, some can be aperiodic")
1346 if(any(typeboundary(:,2*idim-1)==12))
then
1347 if(any(typeboundary(:,2*idim-1)/=12)) typeboundary(:,2*idim-1)=12
1348 select case(physics_type)
1349 case (
'rho',
'ard',
'rd',
'nonlinear',
'ffhd')
1351 typeboundary(:,2*idim-1)=bc_symm
1352 if(mype==0) print *,
'symmetric minimal pole'
1353 case (
'hd',
'rhd',
'srhd',
'mhd',
'rmhd')
1354 typeboundary(:,2*idim-1)=bc_symm
1356 if(phys_energy)
then
1361 select case(coordinate)
1363 typeboundary(phi_+1,2*idim-1)=bc_asymm
1364 if(physics_type==
'mhd'.or.physics_type==
'rmhd') typeboundary(ndir+windex+phi_,2*idim-1)=bc_asymm
1366 typeboundary(3:ndir+1,2*idim-1)=bc_asymm
1367 if(physics_type==
'mhd'.or.physics_type==
'rmhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim-1)=bc_asymm
1369 call mpistop(
'Pole is in cylindrical, polar, spherical coordinates!')
1372 call mpistop(
'Pole treatment for twofl or mf not implemented yet')
1374 call mpistop(
'unknown physics type for setting minimal pole boundary treatment')
1377 if(any(typeboundary(:,2*idim)==12))
then
1378 if(any(typeboundary(:,2*idim)/=12)) typeboundary(:,2*idim)=12
1379 select case(physics_type)
1380 case (
'rho',
'ard',
'rd',
'nonlinear',
'ffhd')
1382 typeboundary(:,2*idim)=bc_symm
1383 if(mype==0) print *,
'symmetric maximal pole'
1384 case (
'hd',
'rhd',
'srhd',
'mhd',
'rmhd')
1385 typeboundary(:,2*idim)=bc_symm
1387 if(phys_energy)
then
1392 select case(coordinate)
1394 typeboundary(phi_+1,2*idim)=bc_asymm
1395 if(physics_type==
'mhd'.or.physics_type==
'rmhd') typeboundary(ndir+windex+phi_,2*idim)=bc_asymm
1397 typeboundary(3:ndir+1,2*idim)=bc_asymm
1398 if(physics_type==
'mhd'.or.physics_type==
'rmhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim)=bc_asymm
1400 call mpistop(
'Pole is in cylindrical, polar, spherical coordinates!')
1403 call mpistop(
'Pole treatment for twofl or mf not implemented yet')
1405 call mpistop(
'unknown physics type for setting maximal pole boundary treatment')
1411 if(.not.phys_energy)
then
1416 if(any(limiter(1:nlevelshi)==
'mp5'))
then
1417 nghostcells=max(nghostcells,3)
1420 if(any(limiter(1:nlevelshi)==
'weno5'))
then
1421 nghostcells=max(nghostcells,3)
1424 if(any(limiter(1:nlevelshi)==
'weno5nm'))
then
1425 nghostcells=max(nghostcells,3)
1428 if(any(limiter(1:nlevelshi)==
'wenoz5'))
then
1429 nghostcells=max(nghostcells,3)
1432 if(any(limiter(1:nlevelshi)==
'wenoz5nm'))
then
1433 nghostcells=max(nghostcells,3)
1436 if(any(limiter(1:nlevelshi)==
'wenozp5'))
then
1437 nghostcells=max(nghostcells,3)
1440 if(any(limiter(1:nlevelshi)==
'wenozp5nm'))
then
1441 nghostcells=max(nghostcells,3)
1444 if(any(limiter(1:nlevelshi)==
'teno5ad'))
then
1445 nghostcells=max(nghostcells,3)
1448 if(any(limiter(1:nlevelshi)==
'weno5cu6'))
then
1449 nghostcells=max(nghostcells,3)
1452 if(any(limiter(1:nlevelshi)==
'ppm'))
then
1453 if(flatsh .or. flatcd)
then
1454 nghostcells=max(nghostcells,4)
1456 nghostcells=max(nghostcells,3)
1460 if(any(limiter(1:nlevelshi)==
'weno7'))
then
1461 nghostcells=max(nghostcells,4)
1464 if(any(limiter(1:nlevelshi)==
'mpweno7'))
then
1465 nghostcells=max(nghostcells,4)
1469 nghostcells = nghostcells + phys_wider_stencil
1472 if(stagger_grid .and. refine_max_level>1 .and. mod(nghostcells,2)/=0)
then
1473 nghostcells=nghostcells+1
1476 select case (coordinate)
1479 xprob^lim^de=xprob^lim^de*two*dpi;
1484 xprob^lim^d=xprob^lim^d*two*dpi;
1490 {ixghi^d = block_nx^d + 2*nghostcells\}
1491 {ixgshi^d = ixghi^d\}
1493 nx_vec = [{domain_nx^d|, }]
1494 block_nx_vec = [{block_nx^d|, }]
1496 if (any(nx_vec < 4) .or. any(mod(nx_vec, 2) == 1)) &
1497 call mpistop(
'Grid size (domain_nx^D) has to be even and >= 4')
1499 if (any(block_nx_vec < 4) .or. any(mod(block_nx_vec, 2) == 1)) &
1500 call mpistop(
'Block size (block_nx^D) has to be even and >= 4')
1502 {
if(mod(domain_nx^d,block_nx^d)/=0) &
1503 call mpistop(
'Grid (domain_nx^D) and block (block_nx^D) must be consistent') \}
1505 if(refine_max_level>nlevelshi.or.refine_max_level<1)
then
1506 write(unitterm,*)
'Error: refine_max_level',refine_max_level,
'>nlevelshi ',nlevelshi
1507 call mpistop(
"Reset nlevelshi and recompile!")
1510 if (any(stretched_dim))
then
1511 allocate(qstretch(0:nlevelshi,1:ndim),dxfirst(0:nlevelshi,1:ndim),&
1512 dxfirst_1mq(0:nlevelshi,1:ndim),dxmid(0:nlevelshi,1:ndim))
1513 allocate(nstretchedblocks(1:nlevelshi,1:ndim))
1514 qstretch(0:nlevelshi,1:ndim)=0.0d0
1515 dxfirst(0:nlevelshi,1:ndim)=0.0d0
1516 nstretchedblocks(1:nlevelshi,1:ndim)=0
1517 {
if (stretch_type(^d) == stretch_uni)
then
1519 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble)
then
1521 write(*,*)
'stretched grid needs finite qstretch_baselevel>1'
1522 write(*,*)
'will try default value for qstretch_baselevel in dimension', ^d
1524 if(xprobmin^d>smalldouble)
then
1525 qstretch_baselevel(^d)=(xprobmax^d/xprobmin^d)**(1.d0/dble(domain_nx^d))
1527 call mpistop(
"can not set qstretch_baselevel automatically")
1530 if(mod(block_nx^d,2)==1) &
1531 call mpistop(
"stretched grid needs even block size block_nxD")
1532 if(mod(domain_nx^d/block_nx^d,2)/=0) &
1533 call mpistop(
"number level 1 blocks in D must be even")
1534 qstretch(1,^d)=qstretch_baselevel(^d)
1535 dxfirst(1,^d)=(xprobmax^d-xprobmin^d) &
1536 *(1.0d0-qstretch(1,^d))/(1.0d0-qstretch(1,^d)**domain_nx^d)
1537 qstretch(0,^d)=qstretch(1,^d)**2
1538 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1539 if(refine_max_level>1)
then
1540 do ilev=2,refine_max_level
1541 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1542 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1543 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1548 {
if(stretch_type(^d) == stretch_uni)
then
1549 write(*,*)
'Stretched dimension ', ^d
1550 write(*,*)
'Using stretched grid with qs=',qstretch(0:refine_max_level,^d)
1551 write(*,*)
' and first cell sizes=',dxfirst(0:refine_max_level,^d)
1554 {
if(stretch_type(^d) == stretch_symm)
then
1556 write(*,*)
'will apply symmetric stretch in dimension', ^d
1558 if(mod(block_nx^d,2)==1) &
1559 call mpistop(
"stretched grid needs even block size block_nxD")
1561 if(nstretchedblocks_baselevel(^d)==0) &
1562 call mpistop(
"need finite even number of stretched blocks at baselevel")
1563 if(mod(nstretchedblocks_baselevel(^d),2)==1) &
1564 call mpistop(
"need even number of stretched blocks at baselevel")
1565 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) &
1566 call mpistop(
'stretched grid needs finite qstretch_baselevel>1')
1568 ipower=(nstretchedblocks_baselevel(^d)/2)*block_nx^d
1569 if(nstretchedblocks_baselevel(^d)==domain_nx^d/block_nx^d)
then
1570 xstretch^d=0.5d0*(xprobmax^d-xprobmin^d)
1572 xstretch^d=(xprobmax^d-xprobmin^d) &
1573 /(2.0d0+dble(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d) &
1574 *(1.0d0-qstretch_baselevel(^d))/(1.0d0-qstretch_baselevel(^d)**ipower))
1576 if(xstretch^d>(xprobmax^d-xprobmin^d)*0.5d0) &
1577 call mpistop(
" stretched grid part should not exceed full domain")
1578 dxfirst(1,^d)=xstretch^d*(1.0d0-qstretch_baselevel(^d)) &
1579 /(1.0d0-qstretch_baselevel(^d)**ipower)
1580 nstretchedblocks(1,^d)=nstretchedblocks_baselevel(^d)
1581 qstretch(1,^d)=qstretch_baselevel(^d)
1582 qstretch(0,^d)=qstretch(1,^d)**2
1583 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1584 dxmid(1,^d)=dxfirst(1,^d)
1585 dxmid(0,^d)=dxfirst(1,^d)*2.0d0
1586 if(refine_max_level>1)
then
1587 do ilev=2,refine_max_level
1588 nstretchedblocks(ilev,^d)=2*nstretchedblocks(ilev-1,^d)
1589 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1590 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1591 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1592 dxmid(ilev,^d)=dxmid(ilev-1,^d)/2.0d0
1596 sizeuniformpart^d=dxfirst(1,^d) &
1597 *(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
1599 print *,
'uniform part of size=',sizeuniformpart^d
1600 print *,
'setting of domain is then=',2*xstretch^d+sizeuniformpart^d
1601 print *,
'versus=',xprobmax^d-xprobmin^d
1603 if(dabs(xprobmax^d-xprobmin^d-2*xstretch^d-sizeuniformpart^d)>smalldouble)
then
1604 call mpistop(
'mismatch in domain size!')
1607 dxfirst_1mq(0:refine_max_level,1:ndim)=dxfirst(0:refine_max_level,1:ndim) &
1608 /(1.0d0-qstretch(0:refine_max_level,1:ndim))
1611 dx_vec = [{xprobmax^d-xprobmin^d|, }] / nx_vec
1614 write(c_ndim,
'(I1)') ^nd
1615 write(unitterm,
'(A30,' // c_ndim //
'(I0," "))') &
1616 ' Domain size (cells): ', nx_vec
1617 write(unitterm,
'(A30,' // c_ndim //
'(E9.3," "))') &
1618 ' Level one dx: ', dx_vec
1621 if (any(dx_vec < smalldouble)) &
1622 call mpistop(
"Incorrect domain size (too small grid spacing)")
1626 if(sum(w_refine_weight(:))==0) w_refine_weight(1) = 1.d0
1627 if(dabs(sum(w_refine_weight(:))-1.d0)>smalldouble)
then
1628 write(unitterm,*)
"Sum of all elements in w_refine_weight be 1.d0"
1629 call mpistop(
"Reset w_refine_weight so the sum is 1.d0")
1632 select case (typeboundspeed)
1637 case(
'cmaxleftright')
1640 call mpistop(
"set typeboundspeed='Einfeldt' or 'cmaxmean' or 'cmaxleftright'")
1643 if (mype==0)
write(unitterm,
'(A30)', advance=
'no')
'Refine estimation: '
1645 select case (refine_criterion)
1647 if (mype==0)
write(unitterm,
'(A)')
"user defined"
1649 if (mype==0)
write(unitterm,
'(A)')
"relative error"
1651 if (mype==0)
write(unitterm,
'(A)')
"Lohner's original scheme"
1653 if (mype==0)
write(unitterm,
'(A)')
"Lohner's scheme"
1655 call mpistop(
"Unknown error estimator, change refine_criterion")
1658 if (tfixgrid<bigdouble/2.0d0)
then
1659 if(mype==0)print*,
'Warning, at time=',tfixgrid,
'the grid will be fixed'
1661 if (itfixgrid<biginteger/2)
then
1662 if(mype==0)print*,
'Warning, at iteration=',itfixgrid,
'the grid will be fixed'
1664 if (ditregrid>1)
then
1665 if(mype==0)print*,
'Note, Grid is reconstructed once every',ditregrid,
'iterations'
1670 select case(slicedir(islice))
1672 if(slicecoord(islice)<xprobmin^d.or.slicecoord(islice)>xprobmax^d) &
1673 write(uniterr,*)
'Warning in read_par_files: ', &
1674 'Slice ', islice,
' coordinate',slicecoord(islice),
'out of bounds for dimension ',slicedir(islice)
1680 write(unitterm,
'(A30,A,A)')
'restart_from_file: ',
' ', trim(restart_from_file)
1681 write(unitterm,
'(A30,L1)')
'converting: ', convert
1682 write(unitterm,
'(A)')
''
1685 deallocate(flux_scheme)
1837 integer,
intent(in) :: fh
1838 integer(MPI_OFFSET_KIND),
intent(out) :: offset_tree
1839 integer(MPI_OFFSET_KIND),
intent(out) :: offset_block
1841 double precision :: rbuf(ndim)
1842 double precision,
allocatable :: params(:)
1843 integer :: i, version
1844 integer :: ibuf(ndim), iw
1845 integer :: er, n_par, tmp_int
1846 integer,
dimension(MPI_STATUS_SIZE) :: st
1847 logical :: periodic(ndim)
1848 character(len=name_len),
allocatable :: var_names(:), param_names(:)
1849 character(len=name_len) :: phys_name, geom_name
1852 call mpi_file_read(fh, version, 1, mpi_integer, st, er)
1854 call mpistop(
"Incompatible file version (maybe old format?)")
1858 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1859 offset_tree = ibuf(1)
1862 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1863 offset_block = ibuf(1)
1866 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1868 if (nw /= ibuf(1))
then
1869 write(*,*)
"nw=",nw,
" and nw found in restart file=",ibuf(1)
1870 write(*,*)
"Please be aware of changes in w at restart."
1875 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1876 if (ibuf(1) /=
ndir)
then
1877 write(*,*)
"ndir in restart file = ",ibuf(1)
1878 write(*,*)
"ndir = ",
ndir
1879 call mpistop(
"reset ndir to ndir in restart file")
1883 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1884 if (ibuf(1) /= ndim)
then
1885 write(*,*)
"ndim in restart file = ",ibuf(1)
1886 write(*,*)
"ndim = ",ndim
1887 call mpistop(
"reset ndim to ndim in restart file")
1891 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1893 write(*,*)
"number of levels in restart file = ",ibuf(1)
1895 call mpistop(
"refine_max_level < num. levels in restart file")
1899 call mpi_file_read(fh,
nleafs, 1, mpi_integer, st, er)
1902 call mpi_file_read(fh,
nparents, 1, mpi_integer, st, er)
1905 call mpi_file_read(fh,
it, 1, mpi_integer, st, er)
1908 call mpi_file_read(fh,
global_time, 1, mpi_double_precision, st, er)
1911 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1912 if (maxval(abs(rbuf(1:ndim) - [ xprobmin^
d ])) > 0)
then
1913 write(*,*)
"Error: xprobmin differs from restart data: ", rbuf(1:ndim)
1914 call mpistop(
"change xprobmin^D in par file")
1918 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1919 if (maxval(abs(rbuf(1:ndim) - [ xprobmax^
d ])) > 0)
then
1920 write(*,*)
"Error: xprobmax differs from restart data: ", rbuf(1:ndim)
1921 call mpistop(
"change xprobmax^D in par file")
1925 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
1926 if (any(ibuf(1:ndim) /= [
domain_nx^
d ]))
then
1927 write(*,*)
"Error: mesh size differs from restart data: ", ibuf(1:ndim)
1928 call mpistop(
"change domain_nx^D in par file")
1932 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
1933 if (any(ibuf(1:ndim) /= [
block_nx^
d ]))
then
1934 write(*,*)
"Error: block size differs from restart data:", ibuf(1:ndim)
1935 call mpistop(
"change block_nx^D in par file")
1939 if (version > 4)
then
1940 call mpi_file_read(fh, periodic, ndim, mpi_logical, st, er)
1941 if ({periodic(^
d) .and. .not.
periodb(^
d) .or. .not.periodic(^
d) .and.
periodb(^
d)| .or. }) &
1942 call mpistop(
"change in periodicity in par file")
1944 call mpi_file_read(fh, geom_name, name_len, mpi_character, st, er)
1947 write(*,*)
"type of coordinates in data is: ", geom_name
1948 call mpistop(
"select the correct coordinates in mod_usr.t file")
1951 call mpi_file_read(fh, stagger_mark_dat, 1, mpi_logical, st, er)
1953 write(*,*)
"Warning: stagger grid flag differs from restart data:", stagger_mark_dat
1959 if (version > 3)
then
1963 call mpi_file_read(fh, var_names(iw), name_len, mpi_character, st, er)
1967 call mpi_file_read(fh, phys_name, name_len, mpi_character, st, er)
1973 call mpi_file_read(fh, n_par, 1, mpi_integer, st, er)
1974 allocate(params(n_par))
1975 allocate(param_names(n_par))
1976 call mpi_file_read(fh, params, n_par, mpi_double_precision, st, er)
1977 call mpi_file_read(fh, param_names, name_len * n_par, mpi_character, st, er)
1980 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1985 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1988 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)