17 integer,
private :: itag
20 double precision,
private :: rk2_alfa
26 character(len=*),
parameter ::
fmt_r =
'es16.8'
27 character(len=*),
parameter ::
fmt_r2 =
'es10.2'
28 character(len=*),
parameter ::
fmt_i =
'i8'
44 integer :: len, stat, n, i, ipars
45 integer,
parameter :: max_files = 20
46 integer :: n_par_files
47 character(len=max_files*std_len) :: all_par_files
48 character(len=std_len) :: tmp_files(max_files), arg
49 logical :: unknown_arg, help, morepars
52 print *,
'-----------------------------------------------------------------------------'
53 print *,
'-----------------------------------------------------------------------------'
54 print *,
'| __ __ ____ ___ _ __ __ ______ ___ ____ |'
55 print *,
'| | \/ | _ \_ _| / \ | \/ | _ \ \ / / \ / ___| |'
56 print *,
'| | |\/| | |_) | |_____ / _ \ | |\/| | |_) \ \ / / _ \| | |'
57 print *,
'| | | | | __/| |_____/ ___ \| | | | _ < \ V / ___ \ |___ |'
58 print *,
'| |_| |_|_| |___| /_/ \_\_| |_|_| \_\ \_/_/ \_\____| |'
59 print *,
'-----------------------------------------------------------------------------'
60 print *,
'-----------------------------------------------------------------------------'
66 all_par_files=
"amrvac.par"
79 CALL get_command_argument(i, arg)
80 IF (len_trim(arg) == 0)
EXIT
84 CALL get_command_argument(i, arg)
86 all_par_files=trim(arg)
89 do while (ipars<max_files.and.morepars)
90 CALL get_command_argument(i+ipars,arg)
92 if (index(trim(arg),
"-")==1.or.len_trim(arg)==0)
then
96 all_par_files=trim(all_par_files)//
" "//trim(arg)
105 CALL get_command_argument(i, arg)
110 CALL get_command_argument(i, arg)
113 case(
"-collapsenext")
115 CALL get_command_argument(i, arg)
118 case(
"-snapshotnext")
120 CALL get_command_argument(i, arg)
128 if(
mype==0)print *,
'convert specified: convert=T'
129 case(
"--help",
"-help")
140 if (unknown_arg)
then
141 print*,
"======================================="
142 print*,
"Error: Command line argument ' ",trim(arg),
" ' not recognized"
143 print*,
"======================================="
150 print *,
'Usage example:'
151 print *,
'mpirun -np 4 ./amrvac -i file.par [file2.par ...]'
152 print *,
' (later .par files override earlier ones)'
154 print *,
'Optional arguments:'
155 print *,
'-convert Convert snapshot files'
156 print *,
'-if file0001.dat Use this snapshot to restart from'
157 print *,
' (you can modify e.g. output names)'
158 print *,
'-resume Automatically resume previous run'
159 print *,
' (useful for long runs on HPC systems)'
160 print *,
'-snapshotnext N Manual index for next snapshot'
161 print *,
'-slicenext N Manual index for next slice output'
162 print *,
'-collapsenext N Manual index for next collapsed output'
171 tmp_files, n_par_files)
188 logical :: fileopen, file_exists
189 integer :: i, j, k, ifile, io_state
190 integer :: iB, isave, iw, level, idim, islice
191 integer :: nx_vec(^ND), block_nx_vec(^ND)
192 integer :: my_unit, iostate
194 double precision :: dx_vec(^ND)
197 character(len=80) :: fmt_string
198 character(len=std_len) :: err_msg, restart_from_file_arg
199 character(len=std_len) :: basename_full, basename_prev, dummy_file
200 character(len=std_len),
dimension(:),
allocatable :: &
201 typeboundary_min^D, typeboundary_max^D
202 character(len=std_len),
allocatable :: limiter(:)
203 character(len=std_len),
allocatable :: gradient_limiter(:)
204 character(len=name_len) :: stretch_dim(ndim)
207 character(len=std_len) :: typesourcesplit
209 character(len=std_len),
allocatable :: flux_scheme(:)
211 character(len=std_len) :: typeboundspeed
213 character(len=std_len) :: time_stepper
215 character(len=std_len) :: time_integrator
217 character(len=std_len) :: typecurl
219 character(len=std_len) :: typeprolonglimit
225 character(len=std_len) :: typecourant
227 double precision,
dimension(nsavehi) :: tsave_log, tsave_dat, tsave_slice, &
228 tsave_collapsed, tsave_custom
229 double precision :: dtsave_log, dtsave_dat, dtsave_slice, &
230 dtsave_collapsed, dtsave_custom
231 integer :: ditsave_log, ditsave_dat, ditsave_slice, &
232 ditsave_collapsed, ditsave_custom
233 double precision :: tsavestart_log, tsavestart_dat, tsavestart_slice, &
234 tsavestart_collapsed, tsavestart_custom
235 integer :: windex, ipower
236 double precision :: sizeuniformpart^D
237 double precision :: im_delta,im_nu,rka54,rka51,rkb54,rka55
250 tsave_log, tsave_dat, tsave_slice, tsave_collapsed, tsave_custom, &
251 dtsave_log, dtsave_dat, dtsave_slice, dtsave_collapsed, dtsave_custom, &
252 ditsave_log, ditsave_dat, ditsave_slice, ditsave_collapsed, ditsave_custom,&
253 tsavestart_log, tsavestart_dat, tsavestart_slice, tsavestart_collapsed,&
259 namelist /methodlist/ time_stepper,time_integrator, &
315 allocate(typeboundary_min^d(nwfluxbc))
316 allocate(typeboundary_max^d(nwfluxbc))
379 typeprolonglimit =
'default'
413 tsave(isave,ifile) = bigdouble
414 itsave(isave,ifile) = biginteger
423 tsave_log = bigdouble
424 tsave_dat = bigdouble
425 tsave_slice = bigdouble
426 tsave_collapsed = bigdouble
427 tsave_custom = bigdouble
429 dtsave_log = bigdouble
430 dtsave_dat = bigdouble
431 dtsave_slice = bigdouble
432 dtsave_collapsed = bigdouble
433 dtsave_custom = bigdouble
435 ditsave_log = biginteger
436 ditsave_dat = biginteger
437 ditsave_slice = biginteger
438 ditsave_collapsed = biginteger
439 ditsave_custom = biginteger
441 tsavestart_log = bigdouble
442 tsavestart_dat = bigdouble
443 tsavestart_slice = bigdouble
444 tsavestart_collapsed = bigdouble
445 tsavestart_custom = bigdouble
465 typecourant =
'maxsum'
475 typeboundspeed =
'Einfeldt'
477 time_stepper =
'twostep'
478 time_integrator =
'default'
507 flux_scheme(level) =
'tvdlf'
509 limiter(level) =
'minmod'
510 gradient_limiter(level) =
'minmod'
515 typesourcesplit =
'sfs'
542 if(i==j.or.j==k.or.k==i)
then
544 else if(i+1==j.or.i-2==j)
then
561 inquire(file=trim(
par_files(i)), exist=file_exists)
563 if (.not. file_exists)
then
564 write(err_msg, *)
"The parameter file " // trim(
par_files(i)) // &
574 read(
unitpar, filelist,
end=101)
577 read(unitpar, savelist,
end=102)
580 read(unitpar, stoplist,
end=103)
583 read(unitpar, methodlist,
end=104)
586 read(unitpar, boundlist,
end=105)
589 read(unitpar, meshlist,
end=106)
592 read(unitpar, paramlist,
end=107)
595 read(unitpar, emissionlist,
end=108)
600 if (base_filename /= basename_prev) &
601 basename_full = trim(basename_full) // trim(base_filename)
602 basename_prev = base_filename
605 base_filename = basename_full
609 dummy_file = trim(base_filename)//
"DUMMY"
610 open(newunit=my_unit, file=trim(dummy_file), iostat=iostate)
611 if (iostate /= 0)
then
612 call mpistop(
"Can't write to output directory (" // &
613 trim(base_filename) //
")")
615 close(my_unit, status=
'delete')
620 if(restart_from_file_arg /= undefined) &
621 restart_from_file=restart_from_file_arg
625 if(restart_from_file == undefined)
then
626 do index_latest_data = -1, 9998
631 if (index_latest_data == -1)
then
634 do index_latest_data = 9999, 0, -1
643 call mpi_bcast(index_latest_data, 1, mpi_integer, 0, icomm, ierrmpi)
645 if (resume_previous_run)
then
646 if (index_latest_data == -1)
then
647 if(mype==0)
write(*,*)
"No snapshots found to resume from, start a new run..."
650 write(restart_from_file,
"(a,i4.4,a)") trim(base_filename),index_latest_data,
".dat"
654 if (restart_from_file == undefined)
then
659 call mpistop(
"Please restart from a snapshot when firstprocess=T")
661 call mpistop(
'Change convert to .false. for a new run!')
664 if (small_pressure < 0.d0)
call mpistop(
"small_pressure should be positive.")
665 if (small_density < 0.d0)
call mpistop(
"small_density should be positive.")
667 if (small_temperature>0.d0) small_pressure=small_density*small_temperature
669 if(convert) autoconvert=.false.
671 where (tsave_log < bigdouble) tsave(:, 1) = tsave_log
672 where (tsave_dat < bigdouble) tsave(:, 2) = tsave_dat
673 where (tsave_slice < bigdouble) tsave(:, 3) = tsave_slice
674 where (tsave_collapsed < bigdouble) tsave(:, 4) = tsave_collapsed
675 where (tsave_custom < bigdouble) tsave(:, 5) = tsave_custom
677 if (dtsave_log < bigdouble) dtsave(1) = dtsave_log
678 if (dtsave_dat < bigdouble) dtsave(2) = dtsave_dat
679 if (dtsave_slice < bigdouble) dtsave(3) = dtsave_slice
680 if (dtsave_collapsed < bigdouble) dtsave(4) = dtsave_collapsed
681 if (dtsave_custom < bigdouble) dtsave(5) = dtsave_custom
683 if (tsavestart_log < bigdouble) tsavestart(1) = tsavestart_log
684 if (tsavestart_dat < bigdouble) tsavestart(2) = tsavestart_dat
685 if (tsavestart_slice < bigdouble) tsavestart(3) = tsavestart_slice
686 if (tsavestart_collapsed < bigdouble) tsavestart(4) = tsavestart_collapsed
687 if (tsavestart_custom < bigdouble) tsavestart(5) = tsavestart_custom
689 if (ditsave_log < bigdouble) ditsave(1) = ditsave_log
690 if (ditsave_dat < bigdouble) ditsave(2) = ditsave_dat
691 if (ditsave_slice < bigdouble) ditsave(3) = ditsave_slice
692 if (ditsave_collapsed < bigdouble) ditsave(4) = ditsave_collapsed
693 if (ditsave_custom < bigdouble) ditsave(5) = ditsave_custom
695 if (wall_time_max < bigdouble) wall_time_max=wall_time_max*3600.d0
698 write(unitterm, *)
''
699 write(unitterm, *)
'Output type | tsavestart | dtsave | ditsave | itsave(1) | tsave(1)'
700 write(fmt_string, *)
'(A12," | ",E9.3E2," | ",E9.3E2," | ",I6," | "'//&
705 if (mype == 0)
write(unitterm, fmt_string) trim(output_names(ifile)), &
706 tsavestart(ifile), dtsave(ifile), ditsave(ifile), itsave(1, ifile), tsave(1, ifile)
709 if (mype == 0)
write(unitterm, *)
''
712 if(slicedir(islice) > ndim) &
713 write(uniterr,*)
'Warning in read_par_files: ', &
714 'Slice ', islice,
' direction',slicedir(islice),
'larger than ndim=',ndim
715 if(slicedir(islice) < 1) &
716 write(uniterr,*)
'Warning in read_par_files: ', &
717 'Slice ', islice,
' direction',slicedir(islice),
'too small, should be [',1,ndim,
']'
720 if(it_max==biginteger .and. time_max==bigdouble.and.mype==0)
write(uniterr,*) &
721 'Warning in read_par_files: it_max or time_max not given!'
723 select case (typecourant)
725 type_courant=type_maxsum
727 type_courant=type_summax
729 type_courant=type_minimum
731 write(unitterm,*)
'Unknown typecourant=',typecourant
732 call mpistop(
"Error from read_par_files: no such typecourant!")
736 select case (flux_scheme(level))
738 flux_method(level)=fs_hll
740 flux_method(level)=fs_hllc
742 flux_method(level)=fs_hlld
744 flux_method(level)=fs_hllcd
746 flux_method(level)=fs_tvdlf
748 flux_method(level)=fs_tvdmu
750 flux_method(level)=fs_tvd
752 flux_method(level)=fs_cd
754 flux_method(level)=fs_cd4
756 flux_method(level)=fs_fd
758 flux_method(level)=fs_source
760 flux_method(level)=fs_nul
762 call mpistop(
"unkown or bad flux scheme")
764 if(flux_scheme(level)==
'tvd'.and.time_stepper/=
'onestep') &
765 call mpistop(
" tvd is onestep method, reset time_stepper='onestep'")
766 if(flux_scheme(level)==
'tvd')
then
767 if(mype==0.and.(.not.dimsplit))
write(unitterm,*) &
768 'Warning: setting dimsplit=T for tvd, as used for level=',level
771 if(flux_scheme(level)==
'hlld'.and.physics_type/=
'mhd' .and. physics_type/=
'twofl') &
772 call mpistop(
"Cannot use hlld flux if not using MHD or 2FL only charges physics!")
774 if(flux_scheme(level)==
'hllc'.and.physics_type==
'mf') &
775 call mpistop(
"Cannot use hllc flux if using magnetofriction physics!")
777 if(flux_scheme(level)==
'tvd'.and.physics_type==
'mf') &
778 call mpistop(
"Cannot use tvd flux if using magnetofriction physics!")
780 if(flux_scheme(level)==
'tvdmu'.and.physics_type==
'mf') &
781 call mpistop(
"Cannot use tvdmu flux if using magnetofriction physics!")
783 if (typepred1(level)==0)
then
784 select case (flux_scheme(level))
786 typepred1(level)=fs_cd
788 typepred1(level)=fs_cd4
790 typepred1(level)=fs_fd
791 case (
'tvdlf',
'tvdmu')
792 typepred1(level)=fs_hancock
794 typepred1(level)=fs_hll
796 typepred1(level)=fs_hllc
798 typepred1(level)=fs_hllcd
800 typepred1(level)=fs_hlld
801 case (
'nul',
'source',
'tvd')
802 typepred1(level)=fs_nul
804 call mpistop(
"No default predictor for this full step")
810 if(any(flux_scheme==
'fd')) need_global_cmax=.true.
811 if(any(limiter==
'schmid1')) need_global_a2max=.true.
814 select case (typecurl)
820 type_curl=stokesbased
822 write(unitterm,*)
"typecurl=",typecurl
823 call mpistop(
"unkown type of curl operator in read_par_files")
827 select case (time_stepper)
831 if (time_integrator==
'default')
then
832 time_integrator=
"Forward_Euler"
834 select case (time_integrator)
835 case (
"Forward_Euler")
836 t_integrator=forward_euler
838 t_integrator=imex_euler
842 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
843 call mpistop(
"unkown onestep time_integrator in read_par_files")
845 use_imex_scheme=(t_integrator==imex_euler.or.t_integrator==imex_sp)
849 if (time_integrator==
'default')
then
850 time_integrator=
"Predictor_Corrector"
852 select case (time_integrator)
853 case (
"Predictor_Corrector")
854 t_integrator=predictor_corrector
859 case (
"IMEX_Midpoint")
860 t_integrator=imex_midpoint
861 case (
"IMEX_Trapezoidal")
862 t_integrator=imex_trapezoidal
864 t_integrator=imex_222
866 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
867 call mpistop(
"unkown twostep time_integrator in read_par_files")
869 use_imex_scheme=(t_integrator==imex_midpoint.or.t_integrator==imex_trapezoidal&
870 .or.t_integrator==imex_222)
871 if (t_integrator==rk2_alf)
then
872 if(rk2_alfa<smalldouble.or.rk2_alfa>one)
call mpistop(
"set rk2_alfa within [0,1]")
874 rk_b2=1.0d0/(2.0d0*rk2_alfa)
880 if (time_integrator==
'default')
then
881 time_integrator=
'ssprk3'
883 select case (time_integrator)
889 t_integrator=imex_ars3
891 t_integrator=imex_232
893 t_integrator=imex_cb3a
895 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
896 call mpistop(
"unkown threestep time_integrator in read_par_files")
898 if(t_integrator==rk3_bt)
then
899 select case(rk3_switch)
929 call mpistop(
"Unknown rk3_switch")
932 rk3_b3=1.0d0-rk3_b1-rk3_b2
934 rk3_c3=rk3_a31+rk3_a32
936 if(t_integrator==ssprk3)
then
937 select case(ssprk_order)
940 rk_beta22=1.0d0/4.0d0
941 rk_beta33=2.0d0/3.0d0
942 rk_alfa21=3.0d0/4.0d0
943 rk_alfa31=1.0d0/3.0d0
947 rk_beta11=1.0d0/2.0d0
948 rk_beta22=1.0d0/2.0d0
949 rk_beta33=1.0d0/3.0d0
951 rk_alfa31=1.0d0/3.0d0
955 call mpistop(
"Unknown ssprk3_order")
957 rk_alfa22=1.0d0-rk_alfa21
958 rk_alfa33=1.0d0-rk_alfa31
960 if(t_integrator==imex_ars3)
then
961 ars_gamma=(3.0d0+dsqrt(3.0d0))/6.0d0
963 if(t_integrator==imex_232)
then
964 select case(imex_switch)
966 im_delta=1.0d0-1.0d0/dsqrt(2.0d0)
967 im_nu=(3.0d0+2.0d0*dsqrt(2.0d0))/6.0d0
968 imex_a21=2.0d0*im_delta
971 imex_b1=1.0d0/(2.0d0*dsqrt(2.0d0))
972 imex_b2=1.0d0/(2.0d0*dsqrt(2.0d0))
977 imex_a21=0.711664700366941d0
978 imex_a31=0.077338168947683d0
979 imex_a32=0.917273367886007d0
980 imex_b1=0.398930808264688d0
981 imex_b2=0.345755244189623d0
982 imex_ha21=0.353842865099275d0
983 imex_ha22=0.353842865099275d0
985 call mpistop(
"Unknown imex_siwtch")
988 imex_c3=imex_a31+imex_a32
989 imex_b3=1.0d0-imex_b1-imex_b2
991 if(t_integrator==imex_cb3a)
then
992 imex_c2 = 0.8925502329346865
995 imex_c3 = imex_c2 / (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0)
997 imex_b2 = (3.0d0*imex_c2 - 1.0d0) / (6.0d0*imex_c2**2)
998 imex_b3 = (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) / (6.0d0*imex_c2**2)
999 imex_a33 = (1.0d0/6.0d0 - imex_b2*imex_c2**2 - imex_b3*imex_c2*imex_c3) / (imex_b3*(imex_c3-imex_c2))
1000 imex_a32 = imex_c3 - imex_a33
1014 use_imex_scheme=(t_integrator==imex_ars3.or.t_integrator==imex_232.or.t_integrator==imex_cb3a)
1018 if (time_integrator==
'default')
then
1019 time_integrator=
"ssprk4"
1021 select case (time_integrator)
1027 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
1028 call mpistop(
"unkown fourstep time_integrator in read_par_files")
1030 if(t_integrator==ssprk4)
then
1031 select case(ssprk_order)
1033 rk_beta11=1.0d0/2.0d0
1034 rk_beta22=1.0d0/2.0d0
1035 rk_beta33=1.0d0/6.0d0
1036 rk_beta44=1.0d0/2.0d0
1038 rk_alfa31=2.0d0/3.0d0
1044 rk_beta11=1.0d0/3.0d0
1045 rk_beta22=1.0d0/3.0d0
1046 rk_beta33=1.0d0/3.0d0
1047 rk_beta44=1.0d0/4.0d0
1050 rk_alfa41=1.0d0/4.0d0
1055 call mpistop(
"Unknown ssprk_order")
1057 rk_alfa22=1.0d0-rk_alfa21
1058 rk_alfa33=1.0d0-rk_alfa31
1059 rk_alfa44=1.0d0-rk_alfa41
1064 if (time_integrator==
'default')
then
1065 time_integrator=
"ssprk5"
1067 select case (time_integrator)
1071 write(unitterm,*)
"time_integrator=",time_integrator,
"time_stepper=",time_stepper
1072 call mpistop(
"unkown fivestep time_integrator in read_par_files")
1074 if(t_integrator==ssprk5)
then
1075 select case(ssprk_order)
1078 rk_beta11=0.391752226571890d0
1079 rk_beta22=0.368410593050371d0
1080 rk_beta33=0.251891774271694d0
1081 rk_beta44=0.544974750228521d0
1082 rk_beta54=0.063692468666290d0
1083 rk_beta55=0.226007483236906d0
1084 rk_alfa21=0.444370493651235d0
1085 rk_alfa31=0.620101851488403d0
1086 rk_alfa41=0.178079954393132d0
1087 rk_alfa53=0.517231671970585d0
1088 rk_alfa54=0.096059710526147d0
1090 rk_alfa22=0.555629506348765d0
1091 rk_alfa33=0.379898148511597d0
1092 rk_alfa44=0.821920045606868d0
1093 rk_alfa55=0.386708617503269d0
1094 rk_alfa22=1.0d0-rk_alfa21
1095 rk_alfa33=1.0d0-rk_alfa31
1096 rk_alfa44=1.0d0-rk_alfa41
1097 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1099 rk_c3=rk_alfa22*rk_c2+rk_beta22
1100 rk_c4=rk_alfa33*rk_c3+rk_beta33
1101 rk_c5=rk_alfa44*rk_c4+rk_beta44
1103 rk_beta11=0.39175222700392d0
1104 rk_beta22=0.36841059262959d0
1105 rk_beta33=0.25189177424738d0
1106 rk_beta44=0.54497475021237d0
1109 rk_alfa21=0.44437049406734d0
1110 rk_alfa31=0.62010185138540d0
1111 rk_alfa41=0.17807995410773d0
1112 rk_alfa53=0.51723167208978d0
1115 rk_alfa22=1.0d0-rk_alfa21
1116 rk_alfa33=1.0d0-rk_alfa31
1117 rk_alfa44=1.0d0-rk_alfa41
1119 rka51=0.00683325884039d0
1120 rka54=0.12759831133288d0
1121 rkb54=0.08460416338212d0
1122 rk_beta54=rkb54-rk_beta44*rka51/rk_alfa41
1123 rk_alfa54=rka54-rk_alfa44*rka51/rk_alfa41
1125 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1127 rk_c3=rk_alfa22*rk_c2+rk_beta22
1128 rk_c4=rk_alfa33*rk_c3+rk_beta33
1129 rk_c5=rk_alfa44*rk_c4+rk_beta44
1130 rk_beta55=1.0d0-rk_beta54-rk_alfa53*rk_c3-rk_alfa54*rk_c4-rk_alfa55*rk_c5
1132 call mpistop(
"Unknown ssprk_order")
1141 use_imex_scheme=.false.
1143 call mpistop(
"Unknown time_stepper in read_par_files")
1147 select case (stretch_dim(i))
1148 case (undefined,
'none')
1149 stretch_type(i) = stretch_none
1150 stretched_dim(i) = .false.
1151 case (
'uni',
'uniform')
1152 stretch_type(i) = stretch_uni
1153 stretched_dim(i) = .true.
1154 case (
'symm',
'symmetric')
1155 stretch_type(i) = stretch_symm
1156 stretched_dim(i) = .true.
1158 stretch_type(i) = stretch_none
1159 stretched_dim(i) = .false.
1160 if (mype == 0) print *,
'Got stretch_type = ', stretch_type(i)
1161 call mpistop(
'Unknown stretch type')
1166 if(typedimsplit ==
'default'.and. dimsplit) typedimsplit=
'xyyx'
1167 if(typedimsplit ==
'default'.and..not.dimsplit) typedimsplit=
'unsplit'
1168 dimsplit = typedimsplit /=
'unsplit'
1171 select case (typesourcesplit)
1173 sourcesplit=sourcesplit_sfs
1175 sourcesplit=sourcesplit_sf
1177 sourcesplit=sourcesplit_ssf
1179 sourcesplit=sourcesplit_ssfss
1181 write(unitterm,*)
'No such typesourcesplit=',typesourcesplit
1182 call mpistop(
"Error: Unknown typesourcesplit!")
1185 if(coordinate==-1)
then
1186 coordinate=cartesian
1188 write(*,*)
'Warning: coordinate system is not specified!'
1189 write(*,*)
'call set_coordinate_system in usr_init in mod_usr.t'
1190 write(*,*)
'Now use Cartesian coordinate'
1194 if(coordinate==cartesian)
then
1197 if(any(stretched_dim))
then
1198 coordinate=cartesian_stretched
1199 slab_uniform=.false.
1203 slab_uniform=.false.
1206 if(coordinate==spherical)
then
1208 if(mype==0)print *,
'Warning: spherical symmetry needs dimsplit=F, resetting'
1213 if (ndim==1) dimsplit=.false.
1216 select case(typeprolonglimit)
1231 allocate(type_limiter(nlevelshi))
1232 allocate(type_gradient_limiter(nlevelshi))
1234 do level=1,nlevelshi
1235 type_limiter(level) = limiter_type(limiter(level))
1236 type_gradient_limiter(level) = limiter_type(gradient_limiter(level))
1239 if (any(limiter(1:nlevelshi)==
'ppm')&
1240 .and.(flatsh.and.physics_type==
'rho'))
then
1241 call mpistop(
" PPM with flatsh=.true. can not be used with physics_type='rho'!")
1247 select case(typeboundary_min^d(iw))
1249 typeboundary(iw,2*^d-1)=bc_special
1251 typeboundary(iw,2*^d-1)=bc_cont
1253 typeboundary(iw,2*^d-1)=bc_symm
1255 typeboundary(iw,2*^d-1)=bc_asymm
1257 typeboundary(iw,2*^d-1)=bc_periodic
1259 typeboundary(iw,2*^d-1)=bc_aperiodic
1261 typeboundary(iw,2*^d-1)=bc_noinflow
1263 typeboundary(iw,2*^d-1)=12
1265 typeboundary(iw,2*^d-1)=bc_data
1267 typeboundary(iw,2*^d-1)=bc_character
1269 write (unitterm,*)
"Undefined boundarytype found in read_par_files", &
1270 typeboundary_min^d(iw),
"for variable iw=",iw,
" and side iB=",2*^d-1
1274 select case(typeboundary_max^d(iw))
1276 typeboundary(iw,2*^d)=bc_special
1278 typeboundary(iw,2*^d)=bc_cont
1280 typeboundary(iw,2*^d)=bc_symm
1282 typeboundary(iw,2*^d)=bc_asymm
1284 typeboundary(iw,2*^d)=bc_periodic
1286 typeboundary(iw,2*^d)=bc_aperiodic
1288 typeboundary(iw,2*^d)=bc_noinflow
1290 typeboundary(iw,2*^d)=12
1292 typeboundary(iw,2*^d)=bc_data
1293 case(
"bc_character")
1294 typeboundary(iw,2*^d)=bc_character
1296 write (unitterm,*)
"Undefined boundarytype found in read_par_files", &
1297 typeboundary_max^d(iw),
"for variable iw=",iw,
" and side iB=",2*^d
1303 if (nwfluxbc<nwflux)
then
1304 do iw = nwfluxbc+1, nwflux
1305 typeboundary(iw,:) = typeboundary(1, :)
1310 do iw = nwflux+1, nwflux+nwaux
1311 typeboundary(iw,:) = typeboundary(1, :)
1315 if (any(typeboundary == 0))
then
1316 call mpistop(
"Not all boundary conditions have been defined")
1320 periodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_periodic))
1321 aperiodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_aperiodic))
1322 if (periodb(idim).or.aperiodb(idim))
then
1324 if (typeboundary(iw,2*idim-1) .ne. typeboundary(iw,2*idim)) &
1325 call mpistop(
"Wrong counterpart in periodic boundary")
1327 if (typeboundary(iw,2*idim-1) /= bc_periodic .and. &
1328 typeboundary(iw,2*idim-1) /= bc_aperiodic)
then
1329 call mpistop(
"Each dimension should either have all "//&
1330 "or no variables periodic, some can be aperiodic")
1337 if(any(typeboundary(:,2*idim-1)==12))
then
1338 if(any(typeboundary(:,2*idim-1)/=12)) typeboundary(:,2*idim-1)=12
1339 if(phys_energy)
then
1344 typeboundary(:,2*idim-1)=bc_symm
1345 if(physics_type/=
'rho')
then
1346 select case(coordinate)
1348 typeboundary(phi_+1,2*idim-1)=bc_asymm
1349 if(physics_type==
'mhd') typeboundary(ndir+windex+phi_,2*idim-1)=bc_asymm
1351 typeboundary(3:ndir+1,2*idim-1)=bc_asymm
1352 if(physics_type==
'mhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim-1)=bc_asymm
1354 call mpistop(
'Pole is in cylindrical, polar, spherical coordinates!')
1358 if(any(typeboundary(:,2*idim)==12))
then
1359 if(any(typeboundary(:,2*idim)/=12)) typeboundary(:,2*idim)=12
1360 if(phys_energy)
then
1365 typeboundary(:,2*idim)=bc_symm
1366 if(physics_type/=
'rho')
then
1367 select case(coordinate)
1369 typeboundary(phi_+1,2*idim)=bc_asymm
1370 if(physics_type==
'mhd') typeboundary(ndir+windex+phi_,2*idim)=bc_asymm
1372 typeboundary(3:ndir+1,2*idim)=bc_asymm
1373 if(physics_type==
'mhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim)=bc_asymm
1375 call mpistop(
'Pole is in cylindrical, polar, spherical coordinates!')
1382 if(any(limiter(1:nlevelshi)==
'mp5'))
then
1383 nghostcells=max(nghostcells,3)
1386 if(any(limiter(1:nlevelshi)==
'weno5'))
then
1387 nghostcells=max(nghostcells,3)
1390 if(any(limiter(1:nlevelshi)==
'weno5nm'))
then
1391 nghostcells=max(nghostcells,3)
1394 if(any(limiter(1:nlevelshi)==
'wenoz5'))
then
1395 nghostcells=max(nghostcells,3)
1398 if(any(limiter(1:nlevelshi)==
'wenoz5nm'))
then
1399 nghostcells=max(nghostcells,3)
1402 if(any(limiter(1:nlevelshi)==
'wenozp5'))
then
1403 nghostcells=max(nghostcells,3)
1406 if(any(limiter(1:nlevelshi)==
'wenozp5nm'))
then
1407 nghostcells=max(nghostcells,3)
1410 if(any(limiter(1:nlevelshi)==
'teno5ad'))
then
1411 nghostcells=max(nghostcells,3)
1414 if(any(limiter(1:nlevelshi)==
'weno5cu6'))
then
1415 nghostcells=max(nghostcells,3)
1418 if(any(limiter(1:nlevelshi)==
'ppm'))
then
1419 nghostcells=max(nghostcells,4)
1422 if(any(limiter(1:nlevelshi)==
'weno7'))
then
1423 nghostcells=max(nghostcells,4)
1426 if(any(limiter(1:nlevelshi)==
'mpweno7'))
then
1427 nghostcells=max(nghostcells,4)
1431 nghostcells = nghostcells + phys_wider_stencil
1434 if(stagger_grid .and. refine_max_level>1 .and. mod(nghostcells,2)/=0)
then
1435 nghostcells=nghostcells+1
1438 select case (coordinate)
1441 xprob^lim^de=xprob^lim^de*two*dpi;
1446 xprob^lim^d=xprob^lim^d*two*dpi;
1452 {ixghi^d = block_nx^d + 2*nghostcells\}
1453 {ixgshi^d = ixghi^d\}
1455 nx_vec = [{domain_nx^d|, }]
1456 block_nx_vec = [{block_nx^d|, }]
1458 if (any(nx_vec < 4) .or. any(mod(nx_vec, 2) == 1)) &
1459 call mpistop(
'Grid size (domain_nx^D) has to be even and >= 4')
1461 if (any(block_nx_vec < 4) .or. any(mod(block_nx_vec, 2) == 1)) &
1462 call mpistop(
'Block size (block_nx^D) has to be even and >= 4')
1464 {
if(mod(domain_nx^d,block_nx^d)/=0) &
1465 call mpistop(
'Grid (domain_nx^D) and block (block_nx^D) must be consistent') \}
1467 if(refine_max_level>nlevelshi.or.refine_max_level<1)
then
1468 write(unitterm,*)
'Error: refine_max_level',refine_max_level,
'>nlevelshi ',nlevelshi
1469 call mpistop(
"Reset nlevelshi and recompile!")
1472 if (any(stretched_dim))
then
1473 allocate(qstretch(0:nlevelshi,1:ndim),dxfirst(0:nlevelshi,1:ndim),&
1474 dxfirst_1mq(0:nlevelshi,1:ndim),dxmid(0:nlevelshi,1:ndim))
1475 allocate(nstretchedblocks(1:nlevelshi,1:ndim))
1476 qstretch(0:nlevelshi,1:ndim)=0.0d0
1477 dxfirst(0:nlevelshi,1:ndim)=0.0d0
1478 nstretchedblocks(1:nlevelshi,1:ndim)=0
1479 {
if (stretch_type(^d) == stretch_uni)
then
1481 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble)
then
1483 write(*,*)
'stretched grid needs finite qstretch_baselevel>1'
1484 write(*,*)
'will try default value for qstretch_baselevel in dimension', ^d
1486 if(xprobmin^d>smalldouble)
then
1487 qstretch_baselevel(^d)=(xprobmax^d/xprobmin^d)**(1.d0/dble(domain_nx^d))
1489 call mpistop(
"can not set qstretch_baselevel automatically")
1492 if(mod(block_nx^d,2)==1) &
1493 call mpistop(
"stretched grid needs even block size block_nxD")
1494 if(mod(domain_nx^d/block_nx^d,2)/=0) &
1495 call mpistop(
"number level 1 blocks in D must be even")
1496 qstretch(1,^d)=qstretch_baselevel(^d)
1497 dxfirst(1,^d)=(xprobmax^d-xprobmin^d) &
1498 *(1.0d0-qstretch(1,^d))/(1.0d0-qstretch(1,^d)**domain_nx^d)
1499 qstretch(0,^d)=qstretch(1,^d)**2
1500 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1501 if(refine_max_level>1)
then
1502 do ilev=2,refine_max_level
1503 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1504 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1505 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1510 {
if(stretch_type(^d) == stretch_uni)
then
1511 write(*,*)
'Stretched dimension ', ^d
1512 write(*,*)
'Using stretched grid with qs=',qstretch(0:refine_max_level,^d)
1513 write(*,*)
' and first cell sizes=',dxfirst(0:refine_max_level,^d)
1516 {
if(stretch_type(^d) == stretch_symm)
then
1518 write(*,*)
'will apply symmetric stretch in dimension', ^d
1520 if(mod(block_nx^d,2)==1) &
1521 call mpistop(
"stretched grid needs even block size block_nxD")
1523 if(nstretchedblocks_baselevel(^d)==0) &
1524 call mpistop(
"need finite even number of stretched blocks at baselevel")
1525 if(mod(nstretchedblocks_baselevel(^d),2)==1) &
1526 call mpistop(
"need even number of stretched blocks at baselevel")
1527 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) &
1528 call mpistop(
'stretched grid needs finite qstretch_baselevel>1')
1530 ipower=(nstretchedblocks_baselevel(^d)/2)*block_nx^d
1531 if(nstretchedblocks_baselevel(^d)==domain_nx^d/block_nx^d)
then
1532 xstretch^d=0.5d0*(xprobmax^d-xprobmin^d)
1534 xstretch^d=(xprobmax^d-xprobmin^d) &
1535 /(2.0d0+dble(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d) &
1536 *(1.0d0-qstretch_baselevel(^d))/(1.0d0-qstretch_baselevel(^d)**ipower))
1538 if(xstretch^d>(xprobmax^d-xprobmin^d)*0.5d0) &
1539 call mpistop(
" stretched grid part should not exceed full domain")
1540 dxfirst(1,^d)=xstretch^d*(1.0d0-qstretch_baselevel(^d)) &
1541 /(1.0d0-qstretch_baselevel(^d)**ipower)
1542 nstretchedblocks(1,^d)=nstretchedblocks_baselevel(^d)
1543 qstretch(1,^d)=qstretch_baselevel(^d)
1544 qstretch(0,^d)=qstretch(1,^d)**2
1545 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1546 dxmid(1,^d)=dxfirst(1,^d)
1547 dxmid(0,^d)=dxfirst(1,^d)*2.0d0
1548 if(refine_max_level>1)
then
1549 do ilev=2,refine_max_level
1550 nstretchedblocks(ilev,^d)=2*nstretchedblocks(ilev-1,^d)
1551 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1552 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1553 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1554 dxmid(ilev,^d)=dxmid(ilev-1,^d)/2.0d0
1558 sizeuniformpart^d=dxfirst(1,^d) &
1559 *(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
1561 print *,
'uniform part of size=',sizeuniformpart^d
1562 print *,
'setting of domain is then=',2*xstretch^d+sizeuniformpart^d
1563 print *,
'versus=',xprobmax^d-xprobmin^d
1565 if(dabs(xprobmax^d-xprobmin^d-2*xstretch^d-sizeuniformpart^d)>smalldouble)
then
1566 call mpistop(
'mismatch in domain size!')
1569 dxfirst_1mq(0:refine_max_level,1:ndim)=dxfirst(0:refine_max_level,1:ndim) &
1570 /(1.0d0-qstretch(0:refine_max_level,1:ndim))
1573 dx_vec = [{xprobmax^d-xprobmin^d|, }] / nx_vec
1576 write(c_ndim,
'(I1)') ^nd
1577 write(unitterm,
'(A30,' // c_ndim //
'(I0," "))') &
1578 ' Domain size (cells): ', nx_vec
1579 write(unitterm,
'(A30,' // c_ndim //
'(E9.3," "))') &
1580 ' Level one dx: ', dx_vec
1583 if (any(dx_vec < smalldouble)) &
1584 call mpistop(
"Incorrect domain size (too small grid spacing)")
1588 if(sum(w_refine_weight(:))==0) w_refine_weight(1) = 1.d0
1589 if(dabs(sum(w_refine_weight(:))-1.d0)>smalldouble)
then
1590 write(unitterm,*)
"Sum of all elements in w_refine_weight be 1.d0"
1591 call mpistop(
"Reset w_refine_weight so the sum is 1.d0")
1594 select case (typeboundspeed)
1599 case(
'cmaxleftright')
1602 call mpistop(
"set typeboundspeed='Einfieldt' or 'cmaxmean' or 'cmaxleftright'")
1605 if (mype==0)
write(unitterm,
'(A30)', advance=
'no')
'Refine estimation: '
1607 select case (refine_criterion)
1609 if (mype==0)
write(unitterm,
'(A)')
"user defined"
1611 if (mype==0)
write(unitterm,
'(A)')
"relative error"
1613 if (mype==0)
write(unitterm,
'(A)')
"Lohner's original scheme"
1615 if (mype==0)
write(unitterm,
'(A)')
"Lohner's scheme"
1617 call mpistop(
"Unknown error estimator, change refine_criterion")
1620 if (tfixgrid<bigdouble/2.0d0)
then
1621 if(mype==0)print*,
'Warning, at time=',tfixgrid,
'the grid will be fixed'
1623 if (itfixgrid<biginteger/2)
then
1624 if(mype==0)print*,
'Warning, at iteration=',itfixgrid,
'the grid will be fixed'
1626 if (ditregrid>1)
then
1627 if(mype==0)print*,
'Note, Grid is reconstructed once every',ditregrid,
'iterations'
1632 select case(slicedir(islice))
1634 if(slicecoord(islice)<xprobmin^d.or.slicecoord(islice)>xprobmax^d) &
1635 write(uniterr,*)
'Warning in read_par_files: ', &
1636 'Slice ', islice,
' coordinate',slicecoord(islice),
'out of bounds for dimension ',slicedir(islice)
1642 write(unitterm,
'(A30,A,A)')
'restart_from_file: ',
' ', trim(restart_from_file)
1643 write(unitterm,
'(A30,L1)')
'converting: ', convert
1644 write(unitterm,
'(A)')
''
1647 deallocate(flux_scheme)
1654 character(len=*),
intent(in) :: line
1656 character(len=*),
intent(in) :: delims
1658 integer,
intent(in) :: n_max
1660 integer,
intent(inout) :: n_found
1662 character(len=*),
intent(inout) :: fields(n_max)
1663 logical,
intent(out),
optional :: fully_read
1665 integer :: ixs_start(n_max)
1666 integer :: ixs_end(n_max)
1667 integer :: ix, ix_prev
1672 do while (n_found < n_max)
1674 ix = verify(line(ix_prev+1:), delims)
1677 n_found = n_found + 1
1678 ixs_start(n_found) = ix_prev + ix
1681 ix = scan(line(ixs_start(n_found)+1:), delims) - 1
1684 ixs_end(n_found) = len(line)
1686 ixs_end(n_found) = ixs_start(n_found) + ix
1689 fields(n_found) = line(ixs_start(n_found):ixs_end(n_found))
1690 ix_prev = ixs_end(n_found)
1693 if (
present(fully_read))
then
1694 ix = verify(line(ix_prev+1:), delims)
1695 fully_read = (ix == 0)
1728 case (
'regression_test')
1732 call mpistop(
"usr_print_log not defined")
1737 call mpistop(
"Error in SaveFile: Unknown typefilelog")
1744 write(*,*)
'No save method is defined for ifile=',ifile
1756 integer,
intent(out) :: fh
1757 integer,
intent(in) :: ix
1758 character(len=*),
intent(in) :: extension
1759 character(len=*),
intent(in),
optional :: suffix
1760 character(len=std_len) :: filename
1763 if (ix >= 10000)
then
1764 call mpistop(
"Number of output files is limited to 10000 (0...9999)")
1767 if (
present(suffix))
then
1769 trim(suffix), ix, extension
1771 write(filename,
"(a,i4.4,a)") trim(
base_filename), ix, extension
1778 amode = ior(mpi_mode_create, mpi_mode_wronly)
1779 call mpi_file_open(mpi_comm_self,filename,amode, &
1783 print *,
"Error, cannot create file ", trim(filename)
1792 integer,
intent(in) :: ix
1793 character(len=std_len) :: filename
1795 write(filename,
"(a,i4.4,a)") trim(
base_filename), ix,
".dat"
1800 character(len=*),
intent(in) :: filename
1804 i = len_trim(filename) - 7
1818 integer,
intent(in) :: fh
1819 integer(kind=MPI_OFFSET_KIND),
intent(in) :: offset_tree
1820 integer(kind=MPI_OFFSET_KIND),
intent(in) :: offset_block
1821 character(len=*),
intent(in) :: dataset_names(:)
1822 integer,
intent(in) :: nw_vars
1823 integer,
dimension(MPI_STATUS_SIZE) :: st
1826 character(len=name_len) :: dname
1829 call mpi_file_write(fh, int(offset_tree), 1, mpi_integer, st, er)
1830 call mpi_file_write(fh, int(offset_block), 1, mpi_integer, st, er)
1831 call mpi_file_write(fh, nw_vars, 1, mpi_integer, st, er)
1832 call mpi_file_write(fh,
ndir, 1, mpi_integer, st, er)
1833 call mpi_file_write(fh,
ndim, 1, mpi_integer, st, er)
1834 call mpi_file_write(fh,
levmax, 1, mpi_integer, st, er)
1835 call mpi_file_write(fh,
nleafs, 1, mpi_integer, st, er)
1836 call mpi_file_write(fh,
nparents, 1, mpi_integer, st, er)
1837 call mpi_file_write(fh,
it, 1, mpi_integer, st, er)
1840 call mpi_file_write(fh,
global_time, 1, mpi_double_precision, st, er)
1842 call mpi_file_write(fh, [ xprobmin^
d ],
ndim, &
1843 mpi_double_precision, st, er)
1844 call mpi_file_write(fh, [ xprobmax^
d ],
ndim, &
1845 mpi_double_precision, st, er)
1847 mpi_integer, st, er)
1849 mpi_integer, st, er)
1852 call mpi_file_write(fh,
periodb,
ndim, mpi_logical, st, er)
1856 name_len, mpi_character, st, er)
1859 call mpi_file_write(fh,
stagger_grid, 1, mpi_logical, st, er)
1864 dname = trim(adjustl((dataset_names(iw))))
1865 call mpi_file_write(fh, dname, name_len, mpi_character, st, er)
1869 call mpi_file_write(fh,
physics_type, name_len, mpi_character, st, er)
1880 call mpi_file_write(fh,
snapshotnext, 1, mpi_integer, st, er)
1882 call mpi_file_write(fh,
snapshotnext+1, 1, mpi_integer, st, er)
1884 call mpi_file_write(fh,
slicenext, 1, mpi_integer, st, er)
1885 call mpi_file_write(fh,
collapsenext, 1, mpi_integer, st, er)
1898 integer,
intent(in) :: fh
1899 integer(kind=MPI_OFFSET_KIND),
intent(in) :: offset_tree
1900 integer(kind=MPI_OFFSET_KIND),
intent(in) :: offset_block
1913 integer,
intent(in) :: fh
1914 integer(MPI_OFFSET_KIND),
intent(out) :: offset_tree
1915 integer(MPI_OFFSET_KIND),
intent(out) :: offset_block
1916 integer :: i, version
1917 integer :: ibuf(ndim), iw
1918 double precision :: rbuf(ndim)
1919 integer,
dimension(MPI_STATUS_SIZE) :: st
1920 character(len=name_len),
allocatable :: var_names(:), param_names(:)
1921 double precision,
allocatable :: params(:)
1922 character(len=name_len) :: phys_name, geom_name
1923 integer :: er, n_par, tmp_int
1924 logical :: stagger_mark_dat, periodic(ndim)
1927 call mpi_file_read(fh, version, 1, mpi_integer, st, er)
1929 call mpistop(
"Incompatible file version (maybe old format?)")
1933 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1934 offset_tree = ibuf(1)
1937 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1938 offset_block = ibuf(1)
1941 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1943 if (nw /= ibuf(1))
then
1944 write(*,*)
"nw=",nw,
" and nw found in restart file=",ibuf(1)
1945 write(*,*)
"Please be aware of changes in w at restart."
1950 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1951 if (ibuf(1) /=
ndir)
then
1952 write(*,*)
"ndir in restart file = ",ibuf(1)
1953 write(*,*)
"ndir = ",
ndir
1954 call mpistop(
"reset ndir to ndir in restart file")
1958 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1959 if (ibuf(1) /= ndim)
then
1960 write(*,*)
"ndim in restart file = ",ibuf(1)
1961 write(*,*)
"ndim = ",ndim
1962 call mpistop(
"reset ndim to ndim in restart file")
1966 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1968 write(*,*)
"number of levels in restart file = ",ibuf(1)
1970 call mpistop(
"refine_max_level < num. levels in restart file")
1974 call mpi_file_read(fh,
nleafs, 1, mpi_integer, st, er)
1977 call mpi_file_read(fh,
nparents, 1, mpi_integer, st, er)
1980 call mpi_file_read(fh,
it, 1, mpi_integer, st, er)
1983 call mpi_file_read(fh,
global_time, 1, mpi_double_precision, st, er)
1986 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1987 if (maxval(abs(rbuf(1:ndim) - [ xprobmin^
d ])) > 0)
then
1988 write(*,*)
"Error: xprobmin differs from restart data: ", rbuf(1:ndim)
1989 call mpistop(
"change xprobmin^D in par file")
1993 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1994 if (maxval(abs(rbuf(1:ndim) - [ xprobmax^
d ])) > 0)
then
1995 write(*,*)
"Error: xprobmax differs from restart data: ", rbuf(1:ndim)
1996 call mpistop(
"change xprobmax^D in par file")
2000 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
2001 if (any(ibuf(1:ndim) /= [
domain_nx^
d ]))
then
2002 write(*,*)
"Error: mesh size differs from restart data: ", ibuf(1:ndim)
2003 call mpistop(
"change domain_nx^D in par file")
2007 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
2008 if (any(ibuf(1:ndim) /= [
block_nx^
d ]))
then
2009 write(*,*)
"Error: block size differs from restart data:", ibuf(1:ndim)
2010 call mpistop(
"change block_nx^D in par file")
2014 if (version > 4)
then
2015 call mpi_file_read(fh, periodic, ndim, mpi_logical, st, er)
2016 if ({periodic(^
d) .and. .not.
periodb(^
d) .or. .not.periodic(^
d) .and.
periodb(^
d)| .or. }) &
2017 call mpistop(
"change in periodicity in par file")
2019 call mpi_file_read(fh, geom_name, name_len, mpi_character, st, er)
2022 write(*,*)
"type of coordinates in data is: ", geom_name
2023 call mpistop(
"select the correct coordinates in mod_usr.t file")
2026 call mpi_file_read(fh, stagger_mark_dat, 1, mpi_logical, st, er)
2028 write(*,*)
"Error: stagger grid flag differs from restart data:", stagger_mark_dat
2029 call mpistop(
"change parameter to use stagger grid")
2034 if (version > 3)
then
2038 call mpi_file_read(fh, var_names(iw), name_len, mpi_character, st, er)
2042 call mpi_file_read(fh, phys_name, name_len, mpi_character, st, er)
2048 call mpi_file_read(fh, n_par, 1, mpi_integer, st, er)
2049 allocate(params(n_par))
2050 allocate(param_names(n_par))
2051 call mpi_file_read(fh, params, n_par, mpi_double_precision, st, er)
2052 call mpi_file_read(fh, param_names, name_len * n_par, mpi_character, st, er)
2055 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
2060 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
2063 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
2081 integer,
intent(in) :: ixo^l
2083 count_ix = product([ ixomax^d ] - [ ixomin^d ] + 1)
2091 integer,
intent(in) :: igrid
2094 integer,
intent(out) :: n_ghost(2*
ndim)
2096 integer,
intent(out) :: ixo^
l
2098 integer,
intent(out) :: n_values
2107 if(ps(igrid)%is_physical_boundary(2*idim-1)) n_ghost(idim)=
nghostcells
2109 if(ps(igrid)%is_physical_boundary(2*idim)) n_ghost(
ndim+idim)=
nghostcells
2113 {ixomin^
d = ixmlo^
d - n_ghost(^
d)\}
2114 {ixomax^
d = ixmhi^
d + n_ghost(
ndim+^
d)\}
2125 integer :: file_handle, igrid, Morton_no, iwrite
2126 integer :: ipe, ix_buffer(2*ndim+1), n_values
2127 integer :: ixO^L, n_ghost(2*ndim)
2128 integer :: ixOs^L,n_values_stagger
2129 integer :: iorecvstatus(MPI_STATUS_SIZE)
2130 integer :: ioastatus(MPI_STATUS_SIZE)
2131 integer :: igrecvstatus(MPI_STATUS_SIZE)
2132 integer :: istatus(MPI_STATUS_SIZE)
2134 integer(kind=MPI_OFFSET_KIND) :: offset_tree_info
2135 integer(kind=MPI_OFFSET_KIND) :: offset_block_data
2136 integer(kind=MPI_OFFSET_KIND) :: offset_offsets
2137 double precision,
allocatable :: w_buffer(:)
2139 integer,
allocatable :: block_ig(:, :)
2140 integer,
allocatable :: block_lvl(:)
2141 integer(kind=MPI_OFFSET_KIND),
allocatable :: block_offset(:)
2148 n_values = n_values +
count_ix(ixgs^
ll) * nws
2150 allocate(w_buffer(n_values))
2153 allocate(block_ig(ndim,
nleafs))
2154 allocate(block_lvl(
nleafs))
2155 allocate(block_offset(
nleafs+1))
2162 offset_tree_info = -1
2163 offset_block_data = -1
2167 call mpi_file_get_position(file_handle, offset_tree_info,
ierrmpi)
2174 igrid =
sfc(1, morton_no)
2175 ipe =
sfc(2, morton_no)
2178 block_ig(:, morton_no) = [ pnode%ig^
d ]
2179 block_lvl(morton_no) = pnode%level
2180 block_offset(morton_no) = 0
2183 call mpi_file_write(file_handle, block_lvl,
size(block_lvl), &
2184 mpi_integer, istatus,
ierrmpi)
2186 call mpi_file_write(file_handle, block_ig,
size(block_ig), &
2187 mpi_integer, istatus,
ierrmpi)
2190 call mpi_file_get_position(file_handle, offset_offsets,
ierrmpi)
2191 call mpi_file_write(file_handle, block_offset(1:
nleafs),
nleafs, &
2194 call mpi_file_get_position(file_handle, offset_block_data,
ierrmpi)
2197 if (offset_block_data - offset_tree_info /= &
2199 nleafs * ((1+ndim) * size_int + 2 * size_int))
then
2201 print *,
"Warning: MPI_OFFSET type /= 8 bytes"
2202 print *,
"This *could* cause problems when reading .dat files"
2206 block_offset(1) = offset_block_data
2216 w_buffer(1:n_values) = pack(ps(igrid)%w(ixo^s, 1:nw), .true.)
2217 {ixosmin^
d = ixomin^
d -1\}
2218 {ixosmax^
d = ixomax^
d \}
2219 n_values_stagger=
count_ix(ixos^l)*nws
2220 w_buffer(n_values+1:n_values+n_values_stagger) = pack(ps(igrid)%ws(ixos^s, 1:nws), .true.)
2221 n_values=n_values+n_values_stagger
2223 w_buffer(1:n_values) = pack(ps(igrid)%w(ixo^s, 1:nw), .true.)
2225 ix_buffer(1) = n_values
2226 ix_buffer(2:) = n_ghost
2229 call mpi_send(ix_buffer, 2*ndim+1, &
2231 call mpi_send(w_buffer, n_values, &
2235 call mpi_file_write(file_handle, ix_buffer(2:), &
2236 2*ndim, mpi_integer, istatus,
ierrmpi)
2237 call mpi_file_write(file_handle, w_buffer, &
2238 n_values, mpi_double_precision, istatus,
ierrmpi)
2241 block_offset(iwrite+1) = block_offset(iwrite) + &
2242 int(n_values, mpi_offset_kind) * size_double + &
2254 call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, ipe, itag,
icomm,&
2256 n_values = ix_buffer(1)
2258 call mpi_recv(w_buffer, n_values, mpi_double_precision,&
2261 call mpi_file_write(file_handle, ix_buffer(2:), &
2262 2*ndim, mpi_integer, istatus,
ierrmpi)
2263 call mpi_file_write(file_handle, w_buffer, &
2264 n_values, mpi_double_precision, istatus,
ierrmpi)
2267 block_offset(iwrite+1) = block_offset(iwrite) + &
2268 int(n_values, mpi_offset_kind) * size_double + &
2274 call mpi_file_seek(file_handle, offset_offsets, mpi_seek_set,
ierrmpi)
2275 call mpi_file_write(file_handle, block_offset(1:
nleafs),
nleafs, &
2279 call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set,
ierrmpi)
2283 call mpi_file_close(file_handle,
ierrmpi)
2298 integer :: ix_buffer(2*ndim+1), n_values, n_values_stagger
2299 integer :: ixO^L, ixOs^L
2300 integer :: file_handle, amode, igrid, Morton_no, iread
2301 integer :: istatus(MPI_STATUS_SIZE)
2302 integer :: iorecvstatus(MPI_STATUS_SIZE)
2303 integer :: ipe,inrecv,nrecv, file_version
2305 integer(MPI_OFFSET_KIND) :: offset_tree_info
2306 integer(MPI_OFFSET_KIND) :: offset_block_data
2307 double precision,
allocatable :: w_buffer(:)
2308 double precision,
dimension(:^D&,:),
allocatable :: w
2309 double precision :: ws(ixGs^T,1:ndim)
2316 mpi_info_null,file_handle,
ierrmpi)
2317 call mpi_file_read(file_handle, file_version, 1, mpi_integer, &
2321 call mpi_bcast(file_version,1,mpi_integer,0,
icomm,
ierrmpi)
2324 if (
mype == 0) print *,
"Unknown version, trying old snapshot reader..."
2325 call mpi_file_close(file_handle,
ierrmpi)
2339 else if (
mype == 0)
then
2340 call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set,
ierrmpi)
2359 n_values = n_values +
count_ix(ixgs^
ll) * nws
2361 allocate(w_buffer(n_values))
2367 call mpi_file_seek(file_handle, offset_tree_info, &
2379 call mpi_file_seek(file_handle, offset_block_data, mpi_seek_set,
ierrmpi)
2387 call mpi_file_read(file_handle,ix_buffer(1:2*ndim), 2*ndim, &
2391 {ixomin^
d = ixmlo^
d - ix_buffer(^
d)\}
2392 {ixomax^
d = ixmhi^
d + ix_buffer(ndim+^
d)\}
2395 {ixosmin^
d = ixomin^
d - 1\}
2396 {ixosmax^
d = ixomax^
d\}
2397 n_values_stagger = n_values
2398 n_values = n_values +
count_ix(ixos^l) * nws
2401 call mpi_file_read(file_handle, w_buffer, n_values, &
2402 mpi_double_precision, istatus,
ierrmpi)
2404 if (
mype == ipe)
then
2408 w(ixo^s, 1:
nw_found) = reshape(w_buffer(1:n_values_stagger), &
2410 ps(igrid)%ws(ixos^s,1:nws)=reshape(w_buffer(n_values_stagger+1:n_values), &
2411 shape(ws(ixos^s, 1:nws)))
2413 w(ixo^s, 1:
nw_found) = reshape(w_buffer(1:n_values), &
2426 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2429 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2432 call mpi_send([ ixo^l, n_values ], 2*ndim+1, &
2434 call mpi_send(w_buffer, n_values, &
2440 call mpi_file_close(file_handle,
ierrmpi)
2449 call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, 0, itag,
icomm,&
2451 {ixomin^
d = ix_buffer(^
d)\}
2452 {ixomax^
d = ix_buffer(ndim+^
d)\}
2453 n_values = ix_buffer(2*ndim+1)
2455 call mpi_recv(w_buffer, n_values, mpi_double_precision,&
2460 {ixosmin^
d = ixomin^
d - 1\}
2461 {ixosmax^
d = ixomax^
d\}
2462 w(ixo^s, 1:
nw_found) = reshape(w_buffer(1:n_values_stagger), &
2464 ps(igrid)%ws(ixos^s,1:nws)=reshape(w_buffer(n_values_stagger+1:n_values), &
2465 shape(ws(ixos^s, 1:nws)))
2467 w(ixo^s, 1:
nw_found) = reshape(w_buffer(1:n_values), &
2480 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2483 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2496 double precision :: wio(ixG^T,1:nw)
2497 integer :: fh, igrid, Morton_no, iread
2498 integer :: levmaxini, ndimini, ndirini
2499 integer :: nwini, neqparini, nxini^D
2500 integer(kind=MPI_OFFSET_KIND) :: offset
2501 integer :: istatus(MPI_STATUS_SIZE)
2502 integer,
allocatable :: iorecvstatus(:,:)
2503 integer :: ipe,inrecv,nrecv
2504 integer :: sendini(7+^ND)
2505 character(len=80) :: filename
2507 double precision :: eqpar_dummy(100)
2511 mpi_mode_rdonly,mpi_info_null,fh,
ierrmpi)
2513 offset=-int(7*size_int+size_double,kind=mpi_offset_kind)
2514 call mpi_file_seek(fh,offset,mpi_seek_end,
ierrmpi)
2518 call mpi_file_read(fh,levmaxini,1,mpi_integer,istatus,
ierrmpi)
2519 call mpi_file_read(fh,ndimini,1,mpi_integer,istatus,
ierrmpi)
2520 call mpi_file_read(fh,ndirini,1,mpi_integer,istatus,
ierrmpi)
2521 call mpi_file_read(fh,nwini,1,mpi_integer,istatus,
ierrmpi)
2522 call mpi_file_read(fh,neqparini,1,mpi_integer,istatus,
ierrmpi)
2523 call mpi_file_read(fh,
it,1,mpi_integer,istatus,
ierrmpi)
2528 write(*,*)
"number of levels in restart file = ",levmaxini
2530 call mpistop(
"refine_max_level < number of levels in restart file")
2532 if (ndimini/=
ndim)
then
2533 write(*,*)
"ndim in restart file = ",ndimini
2534 write(*,*)
"ndim = ",
ndim
2535 call mpistop(
"reset ndim to ndim in restart file")
2537 if (ndirini/=
ndir)
then
2538 write(*,*)
"ndir in restart file = ",ndirini
2539 write(*,*)
"ndir = ",
ndir
2540 call mpistop(
"reset ndir to ndir in restart file")
2543 write(*,*)
"nw=",nw,
" and nw in restart file=",nwini
2544 call mpistop(
"currently, changing nw at restart is not allowed")
2547 offset=offset-int(ndimini*size_int+neqparini*size_double,kind=mpi_offset_kind)
2548 call mpi_file_seek(fh,offset,mpi_seek_end,
ierrmpi)
2550 {
call mpi_file_read(fh,nxini^d,1,mpi_integer,istatus,
ierrmpi)\}
2552 write(*,*)
"Error: reset resolution to ",nxini^d+2*
nghostcells
2553 call mpistop(
"change with setamrvac")
2556 call mpi_file_read(fh,eqpar_dummy,neqparini, &
2557 mpi_double_precision,istatus,
ierrmpi)
2563 sendini=(/
nleafs,levmaxini,ndimini,ndirini,nwini,neqparini,
it ,^d&nxini^d /)
2565 call mpi_bcast(sendini,7+^nd,mpi_integer,0,
icomm,
ierrmpi)
2566 nleafs=sendini(1);levmaxini=sendini(2);ndimini=sendini(3);
2567 ndirini=sendini(4);nwini=sendini(5);
2568 neqparini=sendini(6);
it=sendini(7);
2569 nxini^d=sendini(7+^d);
2576 int(
nleafs,kind=mpi_offset_kind)
2577 call mpi_file_seek(fh,offset,mpi_seek_set,
ierrmpi)
2594 *int(morton_no-1,kind=mpi_offset_kind)
2595 call mpi_file_read_at(fh,offset,ps(igrid)%w,1,
type_block_io, &
2604 *int(morton_no-1,kind=mpi_offset_kind)
2611 call mpi_file_close(fh,
ierrmpi)
2614 allocate(iorecvstatus(mpi_status_size,nrecv))
2621 iorecvstatus(:,inrecv),
ierrmpi)
2623 deallocate(iorecvstatus)
2638 integer :: i, iw, level
2639 double precision :: wmean(1:nw), total_volume
2640 double precision :: volume_coverage(refine_max_level)
2641 integer :: nx^D, nc, ncells, dit
2642 double precision :: dtTimeLast, now, cellupdatesPerSecond
2643 double precision :: activeBlocksPerCore, wctPerCodeTime, timeToFinish
2644 character(len=40) :: fmt_string
2645 character(len=80) :: filename
2646 character(len=2048) :: line
2647 logical,
save :: opened = .false.
2648 integer :: amode, istatus(MPI_STATUS_SIZE)
2649 integer,
parameter :: my_unit = 20
2660 nx^d=ixmhi^d-ixmlo^d+1;
2670 cellupdatespersecond = dble(ncells) * dble(
nstep) * &
2671 dble(dit) / (dttimelast * dble(
npe))
2677 wctpercodetime = dttimelast / max(dit *
dt, epsilon(1.0d0))
2683 if (.not. opened)
then
2689 open(unit=my_unit,file=trim(filename),status=
'replace')
2690 close(my_unit, status=
'delete')
2693 amode = ior(mpi_mode_create,mpi_mode_wronly)
2694 amode = ior(amode,mpi_mode_append)
2696 call mpi_file_open(mpi_comm_self, filename, amode, &
2702 line =
"it global_time dt"
2704 i = len_trim(line) + 2
2705 write(line(i:),
"(a,a)") trim(cons_wnames(level)),
" "
2709 do level = 1, refine_max_level
2710 i = len_trim(line) + 2
2711 write(line(i:),
"(a,i0)")
"c", level
2715 do level=1,refine_max_level
2716 i = len_trim(line) + 2
2717 write(line(i:),
"(a,i0)")
"n", level
2721 line = trim(line) //
" | Xload Xmemory 'Cell_Updates /second/core'"
2722 line = trim(line) //
" 'Active_Blocks/Core' 'Wct Per Code Time [s]'"
2723 line = trim(line) //
" 'TimeToFinish [hrs]'"
2727 call mpi_file_write(
log_fh, trim(line) // new_line(
'a'), &
2728 len_trim(line)+1, mpi_character, istatus,
ierrmpi)
2734 fmt_string =
'(' //
fmt_i //
',2' //
fmt_r //
')'
2736 i = len_trim(line) + 2
2738 write(fmt_string,
'(a,i0,a)')
'(', nw,
fmt_r //
')'
2739 write(line(i:), fmt_string) wmean(1:nw)
2740 i = len_trim(line) + 2
2742 write(fmt_string,
'(a,i0,a)')
'(', refine_max_level,
fmt_r //
')'
2743 write(line(i:), fmt_string) volume_coverage(1:refine_max_level)
2744 i = len_trim(line) + 2
2746 write(fmt_string,
'(a,i0,a)')
'(', refine_max_level,
fmt_i //
')'
2747 write(line(i:), fmt_string)
nleafs_level(1:refine_max_level)
2748 i = len_trim(line) + 2
2750 fmt_string =
'(a,6' //
fmt_r2 //
')'
2751 write(line(i:), fmt_string)
'| ',
xload,
xmemory, cellupdatespersecond, &
2752 activeblockspercore, wctpercodetime, timetofinish
2754 call mpi_file_write(
log_fh, trim(line) // new_line(
'a') , &
2755 len_trim(line)+1, mpi_character, istatus,
ierrmpi)
2765 integer,
parameter :: n_modes = 2
2766 integer,
parameter :: my_unit = 123
2767 character(len=40) :: fmt_string
2768 logical,
save :: file_open = .false.
2770 double precision :: modes(nw, n_modes), volume
2772 do power = 1, n_modes
2777 if (.not. file_open)
then
2781 write(my_unit, *)
"# time mean(w) mean(w**2)"
2784 write(fmt_string,
"(a,i0,a)")
"(", nw * n_modes + 1,
fmt_r //
")"
2795 integer,
intent(in) :: power
2796 double precision,
intent(out) :: mode(nw)
2797 double precision,
intent(out) :: volume
2798 integer :: iigrid, igrid, iw
2799 double precision :: wsum(nw+1)
2800 double precision :: dsum_recv(1:nw+1)
2805 do iigrid = 1, igridstail
2806 igrid = igrids(iigrid)
2809 wsum(nw+1) = wsum(nw+1) + sum(ps(igrid)%dvolume(
ixm^t))
2813 wsum(iw) = wsum(iw) + &
2814 sum(ps(igrid)%dvolume(
ixm^t)*ps(igrid)%w(
ixm^t,iw)**power)
2819 call mpi_allreduce(wsum, dsum_recv, nw+1, mpi_double_precision, &
2823 volume = dsum_recv(nw+1)
2824 mode = dsum_recv(1:nw) / volume
2833 double precision,
intent(out) :: vol_cov(1:refine_max_level)
2834 double precision :: dsum_recv(1:refine_max_level)
2835 integer :: iigrid, igrid, iw, level
2838 vol_cov(1:refine_max_level)=zero
2840 do iigrid = 1, igridstail
2841 igrid = igrids(iigrid);
2843 vol_cov(level) = vol_cov(level)+ &
2848 call mpi_allreduce(vol_cov, dsum_recv, refine_max_level, mpi_double_precision, &
2852 vol_cov = dsum_recv / sum(dsum_recv)
2860 pure function func(w_vec, w_size)
result(val)
2861 integer,
intent(in) :: w_size
2862 double precision,
intent(in) :: w_vec(w_size)
2863 double precision :: val
2866 double precision,
intent(out) :: f_avg
2867 double precision,
intent(out) :: volume
2868 integer :: iigrid, igrid, i^D
2869 double precision :: wsum(2)
2870 double precision :: dsum_recv(2)
2875 do iigrid = 1, igridstail
2876 igrid = igrids(iigrid)
2879 wsum(2) = wsum(2) + sum(ps(igrid)%dvolume(
ixm^t))
2882 {
do i^d = ixmlo^d, ixmhi^d\}
2883 wsum(1) = wsum(1) + ps(igrid)%dvolume(i^d) * &
2884 func(ps(igrid)%w(i^d, :), nw)
2889 call mpi_allreduce(wsum, dsum_recv, 2, mpi_double_precision, &
2890 mpi_sum, icomm, ierrmpi)
2893 volume = dsum_recv(2)
2894 f_avg = dsum_recv(1) / volume
2902 double precision,
intent(out) :: wmax(nw)
2904 integer :: iigrid, igrid, iw
2905 double precision :: wmax_mype(nw),wmax_recv(nw)
2907 wmax_mype(1:nw) = -bigdouble
2910 do iigrid = 1, igridstail
2911 igrid = igrids(iigrid)
2913 wmax_mype(iw)=max(wmax_mype(iw),maxval(ps(igrid)%w(
ixm^t,iw)))
2918 call mpi_allreduce(wmax_mype, wmax_recv, nw, mpi_double_precision, &
2921 wmax(1:nw)=wmax_recv(1:nw)
2929 double precision,
intent(out) :: wmin(nw)
2931 integer :: iigrid, igrid, iw
2932 double precision :: wmin_mype(nw),wmin_recv(nw)
2934 wmin_mype(1:nw) = bigdouble
2937 do iigrid = 1, igridstail
2938 igrid = igrids(iigrid)
2940 wmin_mype(iw)=min(wmin_mype(iw),minval(ps(igrid)%w(
ixm^t,iw)))
2945 call mpi_allreduce(wmin_mype, wmin_recv, nw, mpi_double_precision, &
2948 wmin(1:nw)=wmin_recv(1:nw)
subroutine alloc_node(igrid)
allocate arrays on igrid node
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine generate_plotfile
subroutine read_forest(file_handle)
subroutine write_forest(file_handle)
Collapses D-dimensional output to D-1 view by line-of-sight integration.
subroutine write_collapsed
Module with basic grid data structures.
integer, dimension(:), allocatable, save nleafs_level
How many leaves are present per refinement level.
integer, dimension(:), allocatable, save sfc_to_igrid
Go from a Morton number to an igrid index (for a single processor)
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
integer, save nleafs
Number of leaf block.
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
integer, dimension(:,:), allocatable, save sfc
Array to go from a Morton number to an igrid and processor index. Sfc(1:3, MN) contains [igrid,...
integer, save nparents
Number of parent blocks.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Module with geometry-related routines (e.g., divergence, curl)
This module contains definitions of global parameters and variables and some generic functions/subrou...
character(len=std_len), dimension(:), allocatable typeentropy
Which type of entropy fix to use with Riemann-type solvers.
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
type(state), pointer block
Block pointer for using one block and its previous state.
double precision xload
Stores the memory and load imbalance, used in printlog.
integer nstep
How many sub-steps the time integrator takes.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
integer it_max
Stop the simulation after this many time steps have been taken.
logical internalboundary
if there is an internal boundary
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
double precision tvdlfeps
integer spectrum_wl
wave length for spectrum
logical nocartesian
IO switches for conversion.
integer, dimension(:), allocatable typepred1
The spatial discretization for the predictor step when using a two step PC method.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
logical reset_it
If true, reset iteration count to 0.
double precision small_pressure
integer ixghi
Upper index of grid block arrays.
character(len=std_len) geometry_name
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
logical source_split_usr
Use split or unsplit way to add user's source terms, default: unsplit.
integer domain_nx
number of cells for each dimension in level-one mesh
integer, parameter unitpar
file handle for IO
character(len=std_len) filename_spectrum
Base file name for synthetic EUV spectrum output.
logical resume_previous_run
If true, restart a previous run from the latest snapshot.
double precision global_time
The global simulation time.
integer, dimension(nsavehi, nfile) itsave
Save output of type N on iterations itsave(:, N)
double precision time_max
End time for the simulation.
integer, dimension(3, 3) kr
Kronecker delta tensor.
logical, dimension(:), allocatable logflag
double precision time_init
Start time for the simulation.
logical stretch_uncentered
If true, adjust mod_geometry routines to account for grid stretching (but the flux computation will n...
logical firstprocess
If true, call initonegrid_usr upon restarting.
integer snapshotini
Resume from the snapshot with this index.
double precision small_temperature
error handling
double precision xprob
minimum and maximum domain boundaries for each dimension
integer it
Number of time steps taken.
character(len=std_len) filename_euv
Base file name for synthetic EUV emission output.
logical, dimension(:), allocatable loglimit
double precision, dimension(:), allocatable dg
extent of grid blocks in domain per dimension, in array over levels
integer it_init
initial iteration count
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
integer ditregrid
Reconstruct the AMR grid once every ditregrid iteration(s)
logical saveprim
If true, convert from conservative to primitive variables in output.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer itfixgrid
Fix the AMR grid after this many time steps.
integer, parameter filecollapse_
Constant indicating collapsed output.
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
double precision, dimension(:), allocatable amr_wavefilter
refinement: lohner estimate wavefilter setting
integer, parameter rpxmin
integer, parameter nlevelshi
The maximum number of levels in the grid refinement.
double precision location_slit
location of the slit
logical save_physical_boundary
True for save physical boundary cells in dat files.
logical stagger_grid
True for using stagger grid.
double precision time_convert_factor
Conversion factor for time unit.
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
logical reset_time
If true, reset iteration count and global_time to original values, and start writing snapshots at ind...
integer, parameter nsavehi
Maximum number of saves that can be defined by tsave or itsave.
integer mype
The rank of the current MPI task.
double precision dtpar
If dtpar is positive, it sets the timestep dt, otherwise courantpar is used to limit the time step ba...
character(len=std_len) typediv
integer block_nx
number of cells for each dimension in grid block excluding ghostcells
integer type_block_io
MPI type for IO: block excluding ghost cells.
double precision, dimension(nfile) tsavestart
Start of read out (not counting specified read outs)
integer, dimension(nfile) ditsave
Repeatedly save output of type N when ditsave(N) time steps have passed.
integer, parameter plevel_
logical, dimension(ndim) collapse
If collapse(DIM) is true, generate output integrated over DIM.
integer, dimension(:), allocatable, parameter d
character(len=std_len) usr_filename
User parameter file.
logical ghost_copy
whether copy values instead of interpolation in ghost cells of finer blocks
double precision length_convert_factor
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision imex222_lambda
IMEX-222(lambda) one-parameter family of schemes.
double precision courantpar
The Courant (CFL) number used for the simulation.
double precision wall_time_max
Ending wall time (in hours) for the simulation.
integer ixm
the mesh range (within a block with ghost cells)
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
integer, dimension(:), allocatable flux_method
Which flux scheme of spatial discretization to use (per grid level)
character(len=std_len) collapse_type
integer slowsteps
If > 1, then in the first slowsteps-1 time steps dt is reduced by a factor .
integer ssprk_order
SSPRK choice of methods (both threestep and fourstep, Shu-Osher 2N* implementation) also fivestep SSP...
double precision image_rotate
rotation of image
character(len=std_len) resolution_sxr
resolution of the SXR image
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
integer npe
The number of MPI tasks.
integer nwauxio
Number of auxiliary variables that are only included in output.
double precision, dimension(ndim) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
double precision, dimension(ndim, 2) writespshift
integer imex_switch
IMEX_232 choice and parameters.
double precision time_between_print
to monitor timeintegration loop at given wall-clock time intervals
integer, parameter unitterm
Unit for standard output.
double precision, dimension(nfile) dtsave
Repeatedly save output of type N when dtsave(N) simulation time has passed.
logical, dimension(:), allocatable w_write
integer iprob
problem switch allowing different setups in same usr_mod.t
character(len=std_len) restart_from_file
If not 'unavailable', resume from snapshot with this base file name.
logical, dimension(ndim) periodb
True for dimensions with periodic boundaries.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer, parameter filelog_
Constant indicating log output.
logical final_dt_reduction
If true, allow final dt reduction for matching time_max on output.
logical, dimension(:), allocatable writelevel
integer, parameter fileout_
Constant indicating regular output.
double precision los_theta
direction of the line of sight (LOS)
character(len=std_len) typetvd
Which type of TVD method to use.
integer nbufferx
Number of cells as buffer zone.
double precision, dimension(:), allocatable entropycoef
double precision, dimension(:,:), allocatable dx
double precision tfixgrid
Fix the AMR grid after this time.
integer, parameter unitsnapshot
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(:), allocatable w_refine_weight
Weights of variables used to calculate error for mesh refinement.
character(len=std_len) typedimsplit
character(len=std_len) typeaverage
character(len= *), parameter undefined
double precision, dimension(nsavehi, nfile) tsave
Save output of type N on times tsave(:, N)
double precision spectrum_window_max
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
character(len=std_len) resolution_euv
resolution of the EUV image
integer wavelength
wavelength for output
integer collapselevel
The level at which to produce line-integrated / collapsed output.
logical reset_grid
If true, rebuild the AMR grid upon restarting.
integer, parameter rpxmax
character(len=std_len) resolution_spectrum
resolution of the spectrum
double precision, dimension(:), allocatable refine_threshold
Error tolerance for refinement decision.
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
logical big_image
big image
double precision small_density
double precision spectrum_window_min
spectral window
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
integer refine_max_level
Maximal number of AMR levels.
integer, parameter fileslice_
Constant indicating slice output.
double precision, dimension(:), allocatable derefine_ratio
integer, parameter nfile
Number of output methods.
integer max_blocks
The maximum number of grid blocks in a processor.
integer, parameter fileanalysis_
Constant indicating analysis output (see Writing a custom analysis subroutine)
integer rk3_switch
RK3 Butcher table.
character(len=std_len) typefilelog
Which type of log to write: 'normal', 'special', 'regression_test'.
integer direction_slit
direction of the slit
logical pass_wall_time
If true, wall time is up, modify snapshotnext for later overwrite.
double precision, dimension(1:3) x_origin
where the is the origin (X=0,Y=0) of image
logical final_dt_exit
Force timeloop exit when final dt < dtmin.
integer, dimension(nfile) isaveit
integer, dimension(:,:), allocatable node
integer, dimension(nfile) isavet
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
integer log_fh
MPI file handle for logfile.
Module with slope/flux limiters.
double precision cada3_radius
radius of the asymptotic region [0.001, 10], larger means more accurate in smooth region but more ove...
double precision schmid_rad
Module containing all the particle routines.
This module defines the procedures of a physics module. It contains function pointers for the various...
procedure(sub_write_info), pointer phys_write_info
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
character(len=name_len) physics_type
String describing the physics type of the simulation.
logical phys_energy
Solve energy equation or not.
Writes D-1 slice, can do so in various formats, depending on slice_type.
integer slicenext
the file number of slices
character(len=std_len) slice_type
choose data type of slice: vtu, vtuCC, dat, or csv
integer, dimension(nslicemax) slicedir
The slice direction for each slice.
integer nslices
Number of slices to output.
double precision, dimension(nslicemax) slicecoord
Slice coordinates, see Slice output.
Module for handling problematic values in simulations, such as negative pressures.
integer, public small_values_daverage
Average over this many cells in each direction.
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Module for handling split source terms (split from the fluxes)
double precision timelast
Module with all the methods that users can customize in AMRVAC.
procedure(p_no_args), pointer usr_print_log
procedure(p_no_args), pointer usr_write_analysis
procedure(transform_w), pointer usr_transform_w
The data structure that contains information about a tree node/grid block.