MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_input_output.t
Go to the documentation of this file.
1!> Module for reading input and writing output
3 use mod_comm_lib, only: mpistop
4
5 implicit none
6 public
7
8
9 !> coefficient for rk2_alfa
10 double precision, private :: rk2_alfa
11
12 !> List of compatible versions
13 integer, parameter :: compatible_versions(3) = [3, 4, 5]
14
15 !> number of w found in dat files
16 integer :: nw_found
17
18 !> tag for MPI message
19 integer, private :: itag
20
21 !> whether staggered field is in dat
22 logical, private :: stagger_mark_dat=.false.
23
24 ! Formats used in output
25 character(len=*), parameter :: fmt_r = 'es16.8' ! Default precision
26 character(len=*), parameter :: fmt_r2 = 'es10.2' ! Two digits
27 character(len=*), parameter :: fmt_i = 'i8' ! Integer format
28
29 ! public methods
30 public :: snapshot_write_header
31
32contains
33
34 !> Read the command line arguments passed to amrvac
35 subroutine read_arguments()
37 use mod_slice, only: slicenext
38
39 integer :: len, stat, n, i, ipars
40 integer, parameter :: max_files = 20 ! Maximum number of par files
41 integer :: n_par_files
42 logical :: unknown_arg, help, morepars
43 character(len=max_files*std_len) :: all_par_files
44 character(len=std_len) :: tmp_files(max_files), arg
45
46 if (mype == 0) then
47 print *, '-----------------------------------------------------------------------------'
48 print *, '-----------------------------------------------------------------------------'
49 print *, '| __ __ ____ ___ _ __ __ ______ ___ ____ |'
50 print *, '| | \/ | _ \_ _| / \ | \/ | _ \ \ / / \ / ___| |'
51 print *, '| | |\/| | |_) | |_____ / _ \ | |\/| | |_) \ \ / / _ \| | |'
52 print *, '| | | | | __/| |_____/ ___ \| | | | _ < \ V / ___ \ |___ |'
53 print *, '| |_| |_|_| |___| /_/ \_\_| |_|_| \_\ \_/_/ \_\____| |'
54 print *, '-----------------------------------------------------------------------------'
55 print *, '-----------------------------------------------------------------------------'
56 end if
57
58 ! =============== Fortran 2003 command line reading ================
59
60 ! Default command line arguments
61 all_par_files="amrvac.par"
64 slicenext=-1
66 help=.false.
67 convert=.false.
69
70 ! Argument 0 is program name, so we start from one
71 i = 1
72 unknown_arg=.false.
73 DO
74 CALL get_command_argument(i, arg)
75 IF (len_trim(arg) == 0) EXIT
76 select case(arg)
77 case("-i")
78 i = i+1
79 CALL get_command_argument(i, arg)
80 !if(mype==0)print *,'found argument=',arg
81 all_par_files=trim(arg)
82 morepars=.true.
83 ipars=1
84 do while (ipars<max_files.and.morepars)
85 CALL get_command_argument(i+ipars,arg)
86 !if(mype==0)print *,'found argument=',arg
87 if (index(trim(arg),"-")==1.or.len_trim(arg)==0) then
88 morepars=.false.
89 else
90 ipars=ipars+1
91 all_par_files=trim(all_par_files)//" "//trim(arg)
92 !if(mype==0)print *,'now all_par_files=',TRIM(all_par_files)
93 endif
94 end do
95 !if(mype==0)print *,'-i identified ',ipars, ' par-file arguments'
96 i=i+ipars-1
97 !if(mype==0)print *,'-i arguments passed to all_par_files=',TRIM(all_par_files)
98 case("-if")
99 i = i+1
100 CALL get_command_argument(i, arg)
101 restart_from_file=trim(arg)
102 !if(mype==0)print *,'-if has argument=',TRIM(arg),' passed to restart_from_file=',TRIM(restart_from_file)
103 case("-slicenext")
104 i = i+1
105 CALL get_command_argument(i, arg)
106 read(arg,*,iostat=stat) slicenext
107 !if(mype==0)print *,'-slicenext has argument=',arg,' passed to slicenext=',slicenext
108 case("-collapsenext")
109 i = i+1
110 CALL get_command_argument(i, arg)
111 read(arg,*,iostat=stat) collapsenext
112 !if(mype==0)print *,'-collapsenext has argument=',arg,' passed to collapsenext=',collapsenext
113 case("-snapshotnext")
114 i = i+1
115 CALL get_command_argument(i, arg)
116 read(arg,*,iostat=stat) snapshotnext
117 !if(mype==0)print *,'-snapshotnext has argument=',arg,' passed to snapshotnext=',snapshotnext
118 case("-resume")
120 !if(mype==0)print *,'resume specified: resume_previous_run=T'
121 case("-convert")
122 convert=.true.
123 if(mype==0)print *,'convert specified: convert=T'
124 case("--help","-help")
125 help=.true.
126 EXIT
127 case default
128 unknown_arg=.true.
129 help=.true.
130 EXIT
131 end select
132 i = i+1
133 END DO
134
135 if (unknown_arg) then
136 print*,"======================================="
137 print*,"Error: Command line argument ' ",trim(arg)," ' not recognized"
138 print*,"======================================="
139 help=.true.
140 end if
141
142 ! Show the usage if the help flag was given, or no par file was specified
143 if (help) then
144 if (mype == 0) then
145 print *, 'Usage example:'
146 print *, 'mpirun -np 4 ./amrvac -i file.par [file2.par ...]'
147 print *, ' (later .par files override earlier ones)'
148 print *, ''
149 print *, 'Optional arguments:'
150 print *, '-convert Convert snapshot files'
151 print *, '-if file0001.dat Use this snapshot to restart from'
152 print *, ' (you can modify e.g. output names)'
153 print *, '-resume Automatically resume previous run'
154 print *, ' (useful for long runs on HPC systems)'
155 print *, '-snapshotnext N Manual index for next snapshot'
156 print *, '-slicenext N Manual index for next slice output'
157 print *, '-collapsenext N Manual index for next collapsed output'
158 print *, ''
159 end if
160 call mpi_finalize(ierrmpi)
161 stop
162 end if
163
164 ! Split the input files, in case multiple were given
165 call get_fields_string(all_par_files, " ,'"""//char(9), max_files, &
166 tmp_files, n_par_files)
167
168 allocate(par_files(n_par_files))
169 par_files = tmp_files(1:n_par_files)
170
171 end subroutine read_arguments
172
173 !> Read in the user-supplied parameter-file
174 subroutine read_par_files()
178 use mod_limiter
179 use mod_slice
180 use mod_geometry
181 use mod_source
183
184 double precision, dimension(nsavehi) :: tsave_log, tsave_dat, tsave_slice, &
185 tsave_collapsed, tsave_custom
186 double precision :: dtsave_log, dtsave_dat, dtsave_slice, &
187 dtsave_collapsed, dtsave_custom
188 double precision :: tsavestart_log, tsavestart_dat, tsavestart_slice, &
189 tsavestart_collapsed, tsavestart_custom
190 double precision :: sizeuniformpart^D
191 double precision :: im_delta,im_nu,rka54,rka51,rkb54,rka55
192 double precision :: dx_vec(^ND)
193 integer :: ditsave_log, ditsave_dat, ditsave_slice, &
194 ditsave_collapsed, ditsave_custom
195 integer :: windex, ipower
196 integer :: i, j, k, ifile, io_state
197 integer :: iB, isave, iw, level, idim, islice
198 integer :: nx_vec(^ND), block_nx_vec(^ND)
199 integer :: my_unit, iostate
200 integer :: ilev
201 logical :: fileopen, file_exists
202
203 character :: c_ndim
204 character(len=80) :: fmt_string
205 character(len=std_len) :: err_msg, restart_from_file_arg
206 character(len=std_len) :: basename_full, basename_prev, dummy_file
207 character(len=std_len), dimension(:), allocatable :: &
208 typeboundary_min^D, typeboundary_max^D
209 character(len=std_len), allocatable :: limiter(:)
210 character(len=std_len), allocatable :: gradient_limiter(:)
211 character(len=name_len) :: stretch_dim(ndim)
212 !> How to apply dimensional splitting to the source terms, see
213 !> @ref discretization.md
214 character(len=std_len) :: typesourcesplit
215 !> Which flux scheme of spatial discretization to use (per grid level)
216 character(len=std_len), allocatable :: flux_scheme(:)
217 !> Which type of the maximal bound speed of Riemann fan to use
218 character(len=std_len) :: typeboundspeed
219 !> Which time stepper to use
220 character(len=std_len) :: time_stepper
221 !> Which time integrator to use
222 character(len=std_len) :: time_integrator
223 !> type of curl operator
224 character(len=std_len) :: typecurl
225 !> Limiter used for prolongation to refined grids and ghost cells
226 character(len=std_len) :: typeprolonglimit
227 !> How to compute the CFL-limited time step.
228 !> Options are 'maxsum': max(sum(c/dx)); 'summax': sum(max(c/dx)) and
229 !> 'minimum: max(c/dx), where the summations loop over the grid dimensions and
230 !> c is the velocity. The default 'maxsum' is the conventiontal way of
231 !> computing CFL-limited time steps.
232 character(len=std_len) :: typecourant
233
234
235 namelist /filelist/ base_filename,restart_from_file, &
243
244 namelist /savelist/ tsave,itsave,dtsave,ditsave,nslices,slicedir, &
246 tsave_log, tsave_dat, tsave_slice, tsave_collapsed, tsave_custom, &
247 dtsave_log, dtsave_dat, dtsave_slice, dtsave_collapsed, dtsave_custom, &
248 ditsave_log, ditsave_dat, ditsave_slice, ditsave_collapsed, ditsave_custom,&
249 tsavestart_log, tsavestart_dat, tsavestart_slice, tsavestart_collapsed,&
250 tsavestart_custom, tsavestart
251
254
255 namelist /methodlist/ time_stepper,time_integrator, &
256 source_split_usr,typesourcesplit,local_timestep,&
257 dimsplit,typedimsplit,flux_scheme,&
258 limiter,gradient_limiter,cada3_radius,&
259 loglimit,typeboundspeed, h_correction,&
261 typegrad,typediv,typecurl,&
263 flatcd,flatsh,&
268 schmid_rad^d
269
270 namelist /boundlist/ nghostcells,ghost_copy,&
272
273 namelist /meshlist/ refine_max_level,nbufferx^d,refine_threshold,&
275 stretch_dim, stretch_uncentered, &
279 typeprolonglimit, &
281 namelist /paramlist/ courantpar, dtpar, dtdiffpar, &
282 typecourant, slowsteps
283
284 namelist /emissionlist/ filename_euv,wavelength,&
292
293 ! default maximum number of grid blocks in a processor
294 max_blocks=4000
295
296 ! allocate cell size of all levels
297 allocate(dx(ndim,nlevelshi))
298 {allocate(dg^d(nlevelshi))\}
299 {allocate(ng^d(nlevelshi))\}
300
301 ! default block size excluding ghost cells
302 {block_nx^d = 16\}
303
304 ! default resolution of level-1 mesh (full domain)
305 {domain_nx^d = 32\}
306
307 !! default number of ghost-cell layers at each boundary of a block
308 ! this is now done when the variable is defined in mod_global_parameters
309 ! the physics modules might set this variable in their init subroutine called earlier
310 nghostcells = 2
311
312 ! Allocate boundary conditions arrays in new and old style
313 {
314 allocate(typeboundary_min^d(nwfluxbc))
315 allocate(typeboundary_max^d(nwfluxbc))
316 typeboundary_min^d = undefined
317 typeboundary_max^d = undefined
318 }
319
320 allocate(typeboundary(nw_recon+nwaux,2*ndim))
322
323 ! not save physical boundary in dat files by default
324 save_physical_boundary = .false.
325
326 internalboundary = .false.
327
328 ! defaults for specific options
329 typegrad = 'central'
330 typediv = 'central'
331 typecurl = 'central'
332
333 ! defaults for smallest physical values allowed
334 small_temperature = 0.d0
335 small_pressure = 0.d0
336 small_density = 0.d0
337
338 allocate(small_values_fix_iw(nw))
339 small_values_fix_iw(:) = .true.
340
341 ! defaults for convert behavior
342
343 ! store the -if value from argument in command line
344 restart_from_file_arg = restart_from_file
345 nwauxio = 0
346 nocartesian = .false.
347 saveprim = .false.
348 autoconvert = .false.
349 convert_type = 'vtuBCCmpi'
350 slice_type = 'vtuCC'
351 collapse_type = 'vti'
352 allocate(w_write(nw))
353 w_write(1:nw) = .true.
354 allocate(writelevel(nlevelshi))
355 writelevel(1:nlevelshi) = .true.
356 writespshift(1:ndim,1:2) = zero
357 level_io = -1
358 level_io_min = 1
360 ! endianness: littleendian (default) is 1, bigendian otherwise
361 type_endian = 1
362
363 ! normalization of primitive variables: only for output
364 ! note that length_convert_factor is for length
365 ! this scaling is optional, and must be set consistently if used
366 allocate(w_convert_factor(nw))
367 w_convert_factor(:) = 1.0d0
368 time_convert_factor = 1.0d0
370
371 ! AMR related defaults
373 {nbufferx^d = 0\}
374 allocate(refine_threshold(nlevelshi))
375 refine_threshold(1:nlevelshi) = 0.1d0
376 allocate(derefine_ratio(nlevelshi))
377 derefine_ratio(1:nlevelshi) = 1.0d0/8.0d0
378 typeprolonglimit = 'default'
380 allocate(w_refine_weight(nw+1))
381 w_refine_weight = 0.d0
382 allocate(logflag(nw+1))
383 logflag = .false.
384 allocate(amr_wavefilter(nlevelshi))
385 amr_wavefilter(1:nlevelshi) = 1.0d-2
386 tfixgrid = bigdouble
387 itfixgrid = biginteger
388 ditregrid = 1
389
390 ! Grid stretching defaults
391 stretch_uncentered = .true.
392 stretch_dim(1:ndim) = undefined
393 qstretch_baselevel(1:ndim) = bigdouble
395
396 ! IO defaults
397 it_init = 0
398 it_max = biginteger
399 time_init = 0.d0
400 time_max = bigdouble
401 wall_time_max = bigdouble
402 if(local_timestep) then
403 final_dt_reduction=.false.
404 else
405 final_dt_reduction=.true.
406 endif
407 final_dt_exit=.false.
408 dtmin = 1.0d-10
409 nslices = 0
410 collapse = .false.
411 collapselevel = 1
412 time_between_print = 30.0d0 ! Print status every 30 seconds
413
414 do ifile=1,nfile
415 do isave=1,nsavehi
416 tsave(isave,ifile) = bigdouble ! global_time of saves into the output files
417 itsave(isave,ifile) = biginteger ! it of saves into the output files
418 end do
419 dtsave(ifile) = bigdouble ! time between saves
420 ditsave(ifile) = biginteger ! timesteps between saves
421 isavet(ifile) = 1 ! index for saves by global_time
422 isaveit(ifile) = 1 ! index for saves by it
423 tsavestart(ifile) = 0.0d0
424 end do
425
426 tsave_log = bigdouble
427 tsave_dat = bigdouble
428 tsave_slice = bigdouble
429 tsave_collapsed = bigdouble
430 tsave_custom = bigdouble
431
432 dtsave_log = bigdouble
433 dtsave_dat = bigdouble
434 dtsave_slice = bigdouble
435 dtsave_collapsed = bigdouble
436 dtsave_custom = bigdouble
437
438 ditsave_log = biginteger
439 ditsave_dat = biginteger
440 ditsave_slice = biginteger
441 ditsave_collapsed = biginteger
442 ditsave_custom = biginteger
443
444 tsavestart_log = bigdouble
445 tsavestart_dat = bigdouble
446 tsavestart_slice = bigdouble
447 tsavestart_collapsed = bigdouble
448 tsavestart_custom = bigdouble
449
450 typefilelog = 'default'
451
452 ! defaults for input
453 reset_time = .false.
454 reset_it = .false.
455 firstprocess = .false.
456 reset_grid = .false.
457 base_filename = 'data'
458 usr_filename = ''
459
460
461 ! Defaults for discretization methods
462 typeaverage = 'default'
463 tvdlfeps = one
465 nxdiffusehllc = 0
466 flathllc = .false.
467 slowsteps = -1
468 courantpar = 0.8d0
469 typecourant = 'maxsum'
470 dimsplit = .false.
471 typedimsplit = 'default'
472 if(physics_type=='mhd') then
473 cada3_radius = 0.5d0
474 else
475 cada3_radius = 0.1d0
476 end if
477 {schmid_rad^d = 1.d0\}
478 typetvd = 'roe'
479 typeboundspeed = 'Einfeldt'
480 source_split_usr= .false.
481 time_stepper = 'twostep'
482 time_integrator = 'default'
483 ! default PC or explicit midpoint, hence alfa=0.5
484 rk2_alfa = half
485 ! default IMEX-RK22Ln hence lambda = 1 - 1/sqrt(2)
486 imex222_lambda = 1.0d0 - 1.0d0 / dsqrt(2.0d0)
487 ! default SSPRK(3,3) or Gottlieb-Shu 1998 for threestep
488 ! default SSPRK(4,3) or Spireti-Ruuth for fourstep
489 ! default SSPRK(5,4) using Gottlieb coeffs
490 ssprk_order = 3
491 ! default RK3 butcher table: Heun 3rd order
492 rk3_switch = 3
493 ! default IMEX threestep is IMEX_ARK(232)
494 imex_switch = 1
495
496 ! Defaults for synthesing emission
497 los_theta = 0.d0
498 los_phi = 0.d0
499 image_rotate = 0.d0
500 x_origin = 0.d0
501 big_image = .false.
502 location_slit = 0.d0
503 direction_slit = -1
506 whitelight_instrument='LASCO/C2'
507 r_occultor=-1.d0
508 r_opt_thick=1.d0
509 dat_resolution=.false.
510
511 allocate(flux_scheme(nlevelshi),typepred1(nlevelshi),flux_method(nlevelshi))
512 allocate(limiter(nlevelshi),gradient_limiter(nlevelshi))
513 do level=1,nlevelshi
514 flux_scheme(level) = 'tvdlf'
515 typepred1(level) = 0
516 limiter(level) = 'minmod'
517 gradient_limiter(level) = 'minmod'
518 end do
519
520 flatcd = .false.
521 flatsh = .false.
522 typesourcesplit = 'sfs'
523 allocate(loglimit(nw))
524 loglimit(1:nw) = .false.
525
526 allocate(typeentropy(nw))
527
528 do iw=1,nw
529 typeentropy(iw)='nul' ! Entropy fix type
530 end do
531
532 dtdiffpar = 0.5d0
533 dtpar = -1.d0
534
535 ! problem setup defaults
536 iprob = 1
537
538 ! end defaults
539
540 ! Initialize Kronecker delta, and Levi-Civita tensor
541 do i=1,3
542 do j=1,3
543 if(i==j)then
544 kr(i,j)=1
545 else
546 kr(i,j)=0
547 endif
548 do k=1,3
549 if(i==j.or.j==k.or.k==i)then
550 lvc(i,j,k)=0
551 else if(i+1==j.or.i-2==j)then
552 lvc(i,j,k)=1
553 else
554 lvc(i,j,k)=-1
555 endif
556 enddo
557 enddo
558 enddo
559
560 ! These are used to construct file and log names from multiple par files
561 basename_full = ''
562 basename_prev = ''
563
564 do i = 1, size(par_files)
565 if (mype == 0) print *, "Reading " // trim(par_files(i))
566
567 ! Check whether the file exists
568 inquire(file=trim(par_files(i)), exist=file_exists)
569
570 if (.not. file_exists) then
571 write(err_msg, *) "The parameter file " // trim(par_files(i)) // &
572 " does not exist"
573 call mpistop(trim(err_msg))
574 end if
575
576 open(unitpar, file=trim(par_files(i)), status='old')
577
578 ! Try to read in the namelists. They can be absent or in a different
579 ! order, since we rewind before each read.
580 rewind(unitpar)
581 read(unitpar, filelist, end=101)
582
583101 rewind(unitpar)
584 read(unitpar, savelist, end=102)
585
586102 rewind(unitpar)
587 read(unitpar, stoplist, end=103)
588
589103 rewind(unitpar)
590 read(unitpar, methodlist, end=104)
591
592104 rewind(unitpar)
593 read(unitpar, boundlist, end=105)
594
595105 rewind(unitpar)
596 read(unitpar, meshlist, end=106)
597
598106 rewind(unitpar)
599 read(unitpar, paramlist, end=107)
600
601107 rewind(unitpar)
602 read(unitpar, emissionlist, end=108)
603
604108 close(unitpar)
605
606 ! Append the log and file names given in the par files
607 if (base_filename /= basename_prev) &
608 basename_full = trim(basename_full) // trim(base_filename)
609 basename_prev = base_filename
610 end do
611
612 base_filename = basename_full
613
614 ! Check whether output directory is writable
615 if(mype==0) then
616 dummy_file = trim(base_filename)//"DUMMY"
617 open(newunit=my_unit, file=trim(dummy_file), iostat=iostate)
618 if (iostate /= 0) then
619 call mpistop("Can't write to output directory (" // &
620 trim(base_filename) // ")")
621 else
622 close(my_unit, status='delete')
623 end if
624 end if
625
626 if(source_split_usr) any_source_split=.true.
627
628 ! restart filename from command line overwrites the one in par file
629 if(restart_from_file_arg /= undefined) &
630 restart_from_file=restart_from_file_arg
631
632 ! Root process will search snapshot
633 if (mype == 0) then
634 if(restart_from_file == undefined) then
635 ! search file from highest index
636 file_exists=.false.
637 do index_latest_data = 9999, 0, -1
638 if(snapshot_exists(index_latest_data)) then
639 file_exists=.true.
640 exit
641 end if
642 end do
643 if(.not.file_exists) index_latest_data=-1
644 else
645 ! get index of the given data restarted from
646 index_latest_data=get_snapshot_index(trim(restart_from_file))
647 end if
648 end if
649 call mpi_bcast(index_latest_data, 1, mpi_integer, 0, icomm, ierrmpi)
650
651 if (resume_previous_run) then
652 if (index_latest_data == -1) then
653 if(mype==0) write(*,*) "No snapshots found to resume from, start a new run..."
654 else
655 ! Set file name to restart from
656 write(restart_from_file, "(a,i4.4,a)") trim(base_filename),index_latest_data, ".dat"
657 end if
658 end if
659
660 if (restart_from_file == undefined) then
661 snapshotnext = 0
662 slicenext = 0
663 collapsenext = 0
664 if (firstprocess) &
665 call mpistop("Please restart from a snapshot when firstprocess=T")
666 if (convert) &
667 call mpistop('Change convert to .false. for a new run!')
668 end if
669
670 if (small_pressure < 0.d0) call mpistop("small_pressure should be positive.")
671 if (small_density < 0.d0) call mpistop("small_density should be positive.")
672 ! Give priority to non-zero small temperature
673 if (small_temperature>0.d0) small_pressure=small_density*small_temperature
674
675 if(convert) autoconvert=.false.
676
677 where (tsave_log < bigdouble) tsave(:, 1) = tsave_log
678 where (tsave_dat < bigdouble) tsave(:, 2) = tsave_dat
679 where (tsave_slice < bigdouble) tsave(:, 3) = tsave_slice
680 where (tsave_collapsed < bigdouble) tsave(:, 4) = tsave_collapsed
681 where (tsave_custom < bigdouble) tsave(:, 5) = tsave_custom
682
683 if (dtsave_log < bigdouble) dtsave(1) = dtsave_log
684 if (dtsave_dat < bigdouble) dtsave(2) = dtsave_dat
685 if (dtsave_slice < bigdouble) dtsave(3) = dtsave_slice
686 if (dtsave_collapsed < bigdouble) dtsave(4) = dtsave_collapsed
687 if (dtsave_custom < bigdouble) dtsave(5) = dtsave_custom
688
689 if (tsavestart_log < bigdouble) tsavestart(1) = tsavestart_log
690 if (tsavestart_dat < bigdouble) tsavestart(2) = tsavestart_dat
691 if (tsavestart_slice < bigdouble) tsavestart(3) = tsavestart_slice
692 if (tsavestart_collapsed < bigdouble) tsavestart(4) = tsavestart_collapsed
693 if (tsavestart_custom < bigdouble) tsavestart(5) = tsavestart_custom
694
695 if (ditsave_log < bigdouble) ditsave(1) = ditsave_log
696 if (ditsave_dat < bigdouble) ditsave(2) = ditsave_dat
697 if (ditsave_slice < bigdouble) ditsave(3) = ditsave_slice
698 if (ditsave_collapsed < bigdouble) ditsave(4) = ditsave_collapsed
699 if (ditsave_custom < bigdouble) ditsave(5) = ditsave_custom
700 ! convert hours to seconds for ending wall time
701 if (wall_time_max < bigdouble) wall_time_max=wall_time_max*3600.d0
702
703 if (mype == 0) then
704 write(unitterm, *) ''
705 write(unitterm, *) 'Output type | tsavestart | dtsave | ditsave | itsave(1) | tsave(1)'
706 write(fmt_string, *) '(A12," | ",E9.3E2," | ",E9.3E2," | ",I6," | "'//&
707 ',I6, " | ",E9.3E2)'
708 end if
709
710 do ifile = 1, nfile
711 if (mype == 0) write(unitterm, fmt_string) trim(output_names(ifile)), &
712 tsavestart(ifile), dtsave(ifile), ditsave(ifile), itsave(1, ifile), tsave(1, ifile)
713 end do
714
715 if (mype == 0) write(unitterm, *) ''
716
717 do islice=1,nslices
718 if(slicedir(islice) > ndim) &
719 write(uniterr,*)'Warning in read_par_files: ', &
720 'Slice ', islice,' direction',slicedir(islice),'larger than ndim=',ndim
721 if(slicedir(islice) < 1) &
722 write(uniterr,*)'Warning in read_par_files: ', &
723 'Slice ', islice,' direction',slicedir(islice),'too small, should be [',1,ndim,']'
724 end do
725
726 if(it_max==biginteger .and. time_max==bigdouble.and.mype==0) write(uniterr,*) &
727 'Warning in read_par_files: it_max or time_max not given!'
728
729 select case (typecourant)
730 case ('maxsum')
731 type_courant=type_maxsum
732 case ('summax')
733 type_courant=type_summax
734 case ('minimum')
735 type_courant=type_minimum
736 case default
737 write(unitterm,*)'Unknown typecourant=',typecourant
738 call mpistop("Error from read_par_files: no such typecourant!")
739 end select
740
741 do level=1,nlevelshi
742 select case (flux_scheme(level))
743 case ('hll')
744 flux_method(level)=fs_hll
745 case ('hllc')
746 flux_method(level)=fs_hllc
747 case ('hlld')
748 flux_method(level)=fs_hlld
749 case ('hllcd')
750 flux_method(level)=fs_hllcd
751 case ('tvdlf')
752 flux_method(level)=fs_tvdlf
753 case ('tvdmu')
754 flux_method(level)=fs_tvdmu
755 case ('tvd')
756 flux_method(level)=fs_tvd
757 case ('cd')
758 flux_method(level)=fs_cd
759 case ('cd4')
760 flux_method(level)=fs_cd4
761 case ('fd')
762 flux_method(level)=fs_fd
763 case ('source')
764 flux_method(level)=fs_source
765 case ('nul','null')
766 flux_method(level)=fs_nul
767 case default
768 call mpistop("unkown or bad flux scheme")
769 end select
770 if(flux_scheme(level)=='tvd'.and.time_stepper/='onestep') &
771 call mpistop(" tvd is onestep method, reset time_stepper='onestep'")
772 if(flux_scheme(level)=='tvd')then
773 if(mype==0.and.(.not.dimsplit)) write(unitterm,*) &
774 'Warning: setting dimsplit=T for tvd, as used for level=',level
775 dimsplit=.true.
776 endif
777 if(flux_scheme(level)=='hlld'.and.physics_type/='mhd' .and. physics_type/='twofl') &
778 call mpistop("Cannot use hlld flux if not using MHD or 2FL only charges physics!")
779
780 if(flux_scheme(level)=='hllc'.and.physics_type=='mf') &
781 call mpistop("Cannot use hllc flux if using magnetofriction physics!")
782
783 if(flux_scheme(level)=='tvd'.and.physics_type=='mf') &
784 call mpistop("Cannot use tvd flux if using magnetofriction physics!")
785
786 if(flux_scheme(level)=='tvdmu'.and.physics_type=='mf') &
787 call mpistop("Cannot use tvdmu flux if using magnetofriction physics!")
788
789 if (typepred1(level)==0) then
790 select case (flux_scheme(level))
791 case ('cd')
792 typepred1(level)=fs_cd
793 case ('cd4')
794 typepred1(level)=fs_cd4
795 case ('fd')
796 typepred1(level)=fs_fd
797 case ('tvdlf','tvdmu')
798 typepred1(level)=fs_hancock
799 case ('hll')
800 typepred1(level)=fs_hll
801 case ('hllc')
802 typepred1(level)=fs_hllc
803 case ('hllcd')
804 typepred1(level)=fs_hllcd
805 case ('hlld')
806 typepred1(level)=fs_hlld
807 case ('nul','source','tvd')
808 typepred1(level)=fs_nul
809 case default
810 call mpistop("No default predictor for this full step")
811 end select
812 end if
813 end do
814
815 ! finite difference scheme fd need global maximal speed
816 if(any(flux_scheme=='fd')) need_global_cmax=.true.
817 if(any(limiter=='schmid1')) need_global_a2max=.true.
818
819 ! initialize type_curl
820 select case (typecurl)
821 case ("central")
822 type_curl=central
823 case ("Gaussbased")
824 type_curl=gaussbased
825 case ("Stokesbased")
826 type_curl=stokesbased
827 case default
828 write(unitterm,*) "typecurl=",typecurl
829 call mpistop("unkown type of curl operator in read_par_files")
830 end select
831
832 ! initialize types of time stepper and time integrator
833 select case (time_stepper)
834 case ("onestep")
835 t_stepper=onestep
836 nstep=1
837 if (time_integrator=='default') then
838 time_integrator="Forward_Euler"
839 end if
840 select case (time_integrator)
841 case ("Forward_Euler")
842 t_integrator=forward_euler
843 case ("IMEX_Euler")
844 t_integrator=imex_euler
845 case ("IMEX_SP")
846 t_integrator=imex_sp
847 case default
848 write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
849 call mpistop("unkown onestep time_integrator in read_par_files")
850 end select
851 use_imex_scheme=(t_integrator==imex_euler.or.t_integrator==imex_sp)
852 case ("twostep")
853 t_stepper=twostep
854 nstep=2
855 if (time_integrator=='default') then
856 time_integrator="Predictor_Corrector"
857 endif
858 select case (time_integrator)
859 case ("Predictor_Corrector")
860 t_integrator=predictor_corrector
861 case ("RK2_alfa")
862 t_integrator=rk2_alf
863 case ("ssprk2")
864 t_integrator=ssprk2
865 case ("IMEX_Midpoint")
866 t_integrator=imex_midpoint
867 case ("IMEX_Trapezoidal")
868 t_integrator=imex_trapezoidal
869 case ("IMEX_222")
870 t_integrator=imex_222
871 case default
872 write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
873 call mpistop("unkown twostep time_integrator in read_par_files")
874 end select
875 use_imex_scheme=(t_integrator==imex_midpoint.or.t_integrator==imex_trapezoidal&
876 .or.t_integrator==imex_222)
877 if (t_integrator==rk2_alf) then
878 if(rk2_alfa<smalldouble.or.rk2_alfa>one)call mpistop("set rk2_alfa within [0,1]")
879 rk_a21=rk2_alfa
880 rk_b2=1.0d0/(2.0d0*rk2_alfa)
881 rk_b1=1.0d0-rk_b2
882 endif
883 case ("threestep")
884 t_stepper=threestep
885 nstep=3
886 if (time_integrator=='default') then
887 time_integrator='ssprk3'
888 endif
889 select case (time_integrator)
890 case ("ssprk3")
891 t_integrator=ssprk3
892 case ("RK3_BT")
893 t_integrator=rk3_bt
894 case ("IMEX_ARS3")
895 t_integrator=imex_ars3
896 case ("IMEX_232")
897 t_integrator=imex_232
898 case ("IMEX_CB3a")
899 t_integrator=imex_cb3a
900 case default
901 write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
902 call mpistop("unkown threestep time_integrator in read_par_files")
903 end select
904 if(t_integrator==rk3_bt) then
905 select case(rk3_switch)
906 case(1)
907 ! we code up Ralston 3rd order here
908 rk3_a21=1.0d0/2.0d0
909 rk3_a31=0.0d0
910 rk3_a32=3.0d0/4.0d0
911 rk3_b1=2.0d0/9.0d0
912 rk3_b2=1.0d0/3.0d0
913 case(2)
914 ! we code up RK-Wray 3rd order here
915 rk3_a21=8.0d0/15.0d0
916 rk3_a31=1.0d0/4.0d0
917 rk3_a32=5.0d0/12.0d0
918 rk3_b1=1.0d0/4.0d0
919 rk3_b2=0.0d0
920 case(3)
921 ! we code up Heun 3rd order here
922 rk3_a21=1.0d0/3.0d0
923 rk3_a31=0.0d0
924 rk3_a32=2.0d0/3.0d0
925 rk3_b1=1.0d0/4.0d0
926 rk3_b2=0.0d0
927 case(4)
928 ! we code up Nystrom 3rd order here
929 rk3_a21=2.0d0/3.0d0
930 rk3_a31=0.0d0
931 rk3_a32=2.0d0/3.0d0
932 rk3_b1=1.0d0/4.0d0
933 rk3_b2=3.0d0/8.0d0
934 case default
935 call mpistop("Unknown rk3_switch")
936 end select
937 ! the rest is fixed from above
938 rk3_b3=1.0d0-rk3_b1-rk3_b2
939 rk3_c2=rk3_a21
940 rk3_c3=rk3_a31+rk3_a32
941 endif
942 if(t_integrator==ssprk3) then
943 select case(ssprk_order)
944 case(3) ! this is SSPRK(3,3) Gottlieb-Shu
945 rk_beta11=1.0d0
946 rk_beta22=1.0d0/4.0d0
947 rk_beta33=2.0d0/3.0d0
948 rk_alfa21=3.0d0/4.0d0
949 rk_alfa31=1.0d0/3.0d0
950 rk_c2=1.0d0
951 rk_c3=1.0d0/2.0d0
952 case(2) ! this is SSP(3,2)
953 rk_beta11=1.0d0/2.0d0
954 rk_beta22=1.0d0/2.0d0
955 rk_beta33=1.0d0/3.0d0
956 rk_alfa21=0.0d0
957 rk_alfa31=1.0d0/3.0d0
958 rk_c2=1.0d0/2.0d0
959 rk_c3=1.0d0
960 case default
961 call mpistop("Unknown ssprk3_order")
962 end select
963 rk_alfa22=1.0d0-rk_alfa21
964 rk_alfa33=1.0d0-rk_alfa31
965 endif
966 if(t_integrator==imex_ars3) then
967 ars_gamma=(3.0d0+dsqrt(3.0d0))/6.0d0
968 endif
969 if(t_integrator==imex_232) then
970 select case(imex_switch)
971 case(1) ! this is IMEX_ARK(232)
972 im_delta=1.0d0-1.0d0/dsqrt(2.0d0)
973 im_nu=(3.0d0+2.0d0*dsqrt(2.0d0))/6.0d0
974 imex_a21=2.0d0*im_delta
975 imex_a31=1.0d0-im_nu
976 imex_a32=im_nu
977 imex_b1=1.0d0/(2.0d0*dsqrt(2.0d0))
978 imex_b2=1.0d0/(2.0d0*dsqrt(2.0d0))
979 imex_ha21=im_delta
980 imex_ha22=im_delta
981 case(2) ! this is IMEX_SSP(232)
982 ! doi 10.1002/2017MS001065 Rokhzadi et al
983 imex_a21=0.711664700366941d0
984 imex_a31=0.077338168947683d0
985 imex_a32=0.917273367886007d0
986 imex_b1=0.398930808264688d0
987 imex_b2=0.345755244189623d0
988 imex_ha21=0.353842865099275d0
989 imex_ha22=0.353842865099275d0
990 case default
991 call mpistop("Unknown imex_siwtch")
992 end select
993 imex_c2=imex_a21
994 imex_c3=imex_a31+imex_a32
995 imex_b3=1.0d0-imex_b1-imex_b2
996 endif
997 if(t_integrator==imex_cb3a) then
998 imex_c2 = 0.8925502329346865
999 imex_a22 = imex_c2
1000 imex_ha21 = imex_c2
1001 imex_c3 = imex_c2 / (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0)
1002 imex_ha32 = imex_c3
1003 imex_b2 = (3.0d0*imex_c2 - 1.0d0) / (6.0d0*imex_c2**2)
1004 imex_b3 = (6.0d0*imex_c2**2 - 3.0d0*imex_c2 + 1.0d0) / (6.0d0*imex_c2**2)
1005 imex_a33 = (1.0d0/6.0d0 - imex_b2*imex_c2**2 - imex_b3*imex_c2*imex_c3) / (imex_b3*(imex_c3-imex_c2))
1006 imex_a32 = imex_c3 - imex_a33
1007 ! if (mype == 0) then
1008 ! write(*,*) "================================="
1009 ! write(*,*) "Asserting the order conditions..."
1010 ! ! First order convergence: OK
1011 ! ! Second order convergence
1012 ! write(*,*) -1.0d0/2.0d0 + imex_b2*imex_c2 + imex_b3*imex_c3
1013 ! ! Third order convergence
1014 ! write(*,*) -1.0d0/3.0d0 + imex_b2*imex_c2**2 + imex_b3*imex_c3**2
1015 ! write(*,*) -1.0d0/6.0d0 + imex_b3*imex_ha32*imex_c2
1016 ! write(*,*) -1.0d0/6.0d0 + imex_b2*imex_a22*imex_c2 + imex_b3*imex_a32*imex_c2 + imex_b3*imex_a33*imex_c3
1017 ! write(*,*) "================================="
1018 ! end if
1019 end if
1020 use_imex_scheme=(t_integrator==imex_ars3.or.t_integrator==imex_232.or.t_integrator==imex_cb3a)
1021 case ("fourstep")
1022 t_stepper=fourstep
1023 nstep=4
1024 if (time_integrator=='default') then
1025 time_integrator="ssprk4"
1026 endif
1027 select case (time_integrator)
1028 case ("ssprk4")
1029 t_integrator=ssprk4
1030 case ("rk4")
1031 t_integrator=rk4
1032 case default
1033 write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
1034 call mpistop("unkown fourstep time_integrator in read_par_files")
1035 end select
1036 if(t_integrator==ssprk4) then
1037 select case(ssprk_order)
1038 case(3) ! this is SSPRK(4,3) Spireti-Ruuth
1039 rk_beta11=1.0d0/2.0d0
1040 rk_beta22=1.0d0/2.0d0
1041 rk_beta33=1.0d0/6.0d0
1042 rk_beta44=1.0d0/2.0d0
1043 rk_alfa21=0.0d0
1044 rk_alfa31=2.0d0/3.0d0
1045 rk_alfa41=0.0d0
1046 rk_c2=1.0d0/2.0d0
1047 rk_c3=1.0d0
1048 rk_c4=1.0d0/2.0d0
1049 case(2) ! this is SSP(4,2)
1050 rk_beta11=1.0d0/3.0d0
1051 rk_beta22=1.0d0/3.0d0
1052 rk_beta33=1.0d0/3.0d0
1053 rk_beta44=1.0d0/4.0d0
1054 rk_alfa21=0.0d0
1055 rk_alfa31=0.0d0
1056 rk_alfa41=1.0d0/4.0d0
1057 rk_c2=1.0d0/3.0d0
1058 rk_c3=2.0d0/3.0d0
1059 rk_c4=1.0d0
1060 case default
1061 call mpistop("Unknown ssprk_order")
1062 end select
1063 rk_alfa22=1.0d0-rk_alfa21
1064 rk_alfa33=1.0d0-rk_alfa31
1065 rk_alfa44=1.0d0-rk_alfa41
1066 endif
1067 case ("fivestep")
1068 t_stepper=fivestep
1069 nstep=5
1070 if (time_integrator=='default') then
1071 time_integrator="ssprk5"
1072 end if
1073 select case (time_integrator)
1074 case ("ssprk5")
1075 t_integrator=ssprk5
1076 case default
1077 write(unitterm,*) "time_integrator=",time_integrator,"time_stepper=",time_stepper
1078 call mpistop("unkown fivestep time_integrator in read_par_files")
1079 end select
1080 if(t_integrator==ssprk5) then
1081 select case(ssprk_order)
1082 ! we use ssprk_order to intercompare the different coefficient choices
1083 case(3) ! From Gottlieb 2005
1084 rk_beta11=0.391752226571890d0
1085 rk_beta22=0.368410593050371d0
1086 rk_beta33=0.251891774271694d0
1087 rk_beta44=0.544974750228521d0
1088 rk_beta54=0.063692468666290d0
1089 rk_beta55=0.226007483236906d0
1090 rk_alfa21=0.444370493651235d0
1091 rk_alfa31=0.620101851488403d0
1092 rk_alfa41=0.178079954393132d0
1093 rk_alfa53=0.517231671970585d0
1094 rk_alfa54=0.096059710526147d0
1095
1096 rk_alfa22=0.555629506348765d0
1097 rk_alfa33=0.379898148511597d0
1098 rk_alfa44=0.821920045606868d0
1099 rk_alfa55=0.386708617503269d0
1100 rk_alfa22=1.0d0-rk_alfa21
1101 rk_alfa33=1.0d0-rk_alfa31
1102 rk_alfa44=1.0d0-rk_alfa41
1103 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1104 rk_c2=rk_beta11
1105 rk_c3=rk_alfa22*rk_c2+rk_beta22
1106 rk_c4=rk_alfa33*rk_c3+rk_beta33
1107 rk_c5=rk_alfa44*rk_c4+rk_beta44
1108 case(2) ! From Spireti-Ruuth
1109 rk_beta11=0.39175222700392d0
1110 rk_beta22=0.36841059262959d0
1111 rk_beta33=0.25189177424738d0
1112 rk_beta44=0.54497475021237d0
1113 !rk_beta54=0.06369246925946d0
1114 !rk_beta55=0.22600748319395d0
1115 rk_alfa21=0.44437049406734d0
1116 rk_alfa31=0.62010185138540d0
1117 rk_alfa41=0.17807995410773d0
1118 rk_alfa53=0.51723167208978d0
1119 !rk_alfa54=0.09605971145044d0
1120
1121 rk_alfa22=1.0d0-rk_alfa21
1122 rk_alfa33=1.0d0-rk_alfa31
1123 rk_alfa44=1.0d0-rk_alfa41
1124
1125 rka51=0.00683325884039d0
1126 rka54=0.12759831133288d0
1127 rkb54=0.08460416338212d0
1128 rk_beta54=rkb54-rk_beta44*rka51/rk_alfa41
1129 rk_alfa54=rka54-rk_alfa44*rka51/rk_alfa41
1130
1131 rk_alfa55=1.0d0-rk_alfa53-rk_alfa54
1132 rk_c2=rk_beta11
1133 rk_c3=rk_alfa22*rk_c2+rk_beta22
1134 rk_c4=rk_alfa33*rk_c3+rk_beta33
1135 rk_c5=rk_alfa44*rk_c4+rk_beta44
1136 rk_beta55=1.0d0-rk_beta54-rk_alfa53*rk_c3-rk_alfa54*rk_c4-rk_alfa55*rk_c5
1137 case default
1138 call mpistop("Unknown ssprk_order")
1139 end select
1140 ! the following combinations must be unity
1141 !print *,rk_beta55+rk_beta54+rk_alfa53*rk_c3+rk_alfa54*rk_c4+rk_alfa55*rk_c5
1142 !print *,rk_alfa22+rk_alfa21
1143 !print *,rk_alfa33+rk_alfa31
1144 !print *,rk_alfa44+rk_alfa41
1145 !print *,rk_alfa55+rk_alfa53+rk_alfa54
1146 endif
1147 use_imex_scheme=.false.
1148 case default
1149 call mpistop("Unknown time_stepper in read_par_files")
1150 end select
1151
1152 do i = 1, ndim
1153 select case (stretch_dim(i))
1154 case (undefined, 'none')
1155 stretch_type(i) = stretch_none
1156 stretched_dim(i) = .false.
1157 case ('uni','uniform')
1158 stretch_type(i) = stretch_uni
1159 stretched_dim(i) = .true.
1160 case ('symm', 'symmetric')
1161 stretch_type(i) = stretch_symm
1162 stretched_dim(i) = .true.
1163 case default
1164 stretch_type(i) = stretch_none
1165 stretched_dim(i) = .false.
1166 if (mype == 0) print *, 'Got stretch_type = ', stretch_type(i)
1167 call mpistop('Unknown stretch type')
1168 end select
1169 end do
1170
1171 ! Harmonize the parameters for dimensional splitting and source splitting
1172 if(typedimsplit =='default'.and. dimsplit) typedimsplit='xyyx'
1173 if(typedimsplit =='default'.and..not.dimsplit) typedimsplit='unsplit'
1174 dimsplit = typedimsplit /='unsplit'
1175
1176 ! initialize types of split-source addition
1177 select case (typesourcesplit)
1178 case ('sfs')
1179 sourcesplit=sourcesplit_sfs
1180 case ('sf')
1181 sourcesplit=sourcesplit_sf
1182 case ('ssf')
1183 sourcesplit=sourcesplit_ssf
1184 case ('ssfss')
1185 sourcesplit=sourcesplit_ssfss
1186 case default
1187 write(unitterm,*)'No such typesourcesplit=',typesourcesplit
1188 call mpistop("Error: Unknown typesourcesplit!")
1189 end select
1190
1191 if(coordinate==-1) then
1192 coordinate=cartesian
1193 if(mype==0) then
1194 write(*,*) 'Warning: coordinate system is not specified!'
1195 write(*,*) 'call set_coordinate_system in usr_init in mod_usr.t'
1196 write(*,*) 'Now use Cartesian coordinate'
1197 end if
1198 end if
1199
1200 if(coordinate==cartesian) then
1201 slab=.true.
1202 slab_uniform=.true.
1203 if(any(stretched_dim)) then
1204 coordinate=cartesian_stretched
1205 slab_uniform=.false.
1206 end if
1207 else
1208 slab=.false.
1209 slab_uniform=.false.
1210 end if
1211
1212 if(coordinate==spherical) then
1213 if(dimsplit) then
1214 if(mype==0)print *,'Warning: spherical symmetry needs dimsplit=F, resetting'
1215 dimsplit=.false.
1216 end if
1217 end if
1218
1219 if (ndim==1) dimsplit=.false.
1220
1221 ! type limiter of prolongation
1222 select case(typeprolonglimit)
1223 case('unlimit')
1224 ! unlimited
1225 prolong_limiter=1
1226 case('minmod')
1227 prolong_limiter=2
1228 case('woodward')
1229 prolong_limiter=3
1230 case('koren')
1231 prolong_limiter=4
1232 case default
1233 prolong_limiter=0
1234 end select
1235
1236 ! Type limiter is of integer type for performance
1237 allocate(type_limiter(nlevelshi))
1238 allocate(type_gradient_limiter(nlevelshi))
1239
1240 do level=1,nlevelshi
1241 type_limiter(level) = limiter_type(limiter(level))
1242 type_gradient_limiter(level) = limiter_type(gradient_limiter(level))
1243 end do
1244
1245 if (any(limiter(1:nlevelshi)== 'ppm')&
1246 .and.(flatsh.and.physics_type=='rho')) then
1247 call mpistop(" PPM with flatsh=.true. can not be used with physics_type='rho'!")
1248 end if
1249
1250 ! Copy boundary conditions to typeboundary, which is used internally
1251 {
1252 do iw=1,nwfluxbc
1253 select case(typeboundary_min^d(iw))
1254 case("special")
1255 typeboundary(iw,2*^d-1)=bc_special
1256 case("cont")
1257 typeboundary(iw,2*^d-1)=bc_cont
1258 case("symm")
1259 typeboundary(iw,2*^d-1)=bc_symm
1260 case("asymm")
1261 typeboundary(iw,2*^d-1)=bc_asymm
1262 case("periodic")
1263 typeboundary(iw,2*^d-1)=bc_periodic
1264 case("aperiodic")
1265 typeboundary(iw,2*^d-1)=bc_aperiodic
1266 case("noinflow")
1267 typeboundary(iw,2*^d-1)=bc_noinflow
1268 case("pole")
1269 typeboundary(iw,2*^d-1)=12
1270 case("bc_data")
1271 typeboundary(iw,2*^d-1)=bc_data
1272 case("bc_icarus")
1273 typeboundary(iw,2*^d-1)=bc_icarus
1274 case("character")
1275 typeboundary(iw,2*^d-1)=bc_character
1276 case default
1277 write (unitterm,*) "Undefined boundarytype found in read_par_files", &
1278 typeboundary_min^d(iw),"for variable iw=",iw," and side iB=",2*^d-1
1279 end select
1280 end do
1281 do iw=1,nwfluxbc
1282 select case(typeboundary_max^d(iw))
1283 case("special")
1284 typeboundary(iw,2*^d)=bc_special
1285 case("cont")
1286 typeboundary(iw,2*^d)=bc_cont
1287 case("symm")
1288 typeboundary(iw,2*^d)=bc_symm
1289 case("asymm")
1290 typeboundary(iw,2*^d)=bc_asymm
1291 case("periodic")
1292 typeboundary(iw,2*^d)=bc_periodic
1293 case("aperiodic")
1294 typeboundary(iw,2*^d)=bc_aperiodic
1295 case("noinflow")
1296 typeboundary(iw,2*^d)=bc_noinflow
1297 case("pole")
1298 typeboundary(iw,2*^d)=12
1299 case("bc_data")
1300 typeboundary(iw,2*^d)=bc_data
1301 case("bc_icarus")
1302 typeboundary(iw,2*^d-1)=bc_icarus
1303 case("bc_character")
1304 typeboundary(iw,2*^d)=bc_character
1305 case default
1306 write (unitterm,*) "Undefined boundarytype found in read_par_files", &
1307 typeboundary_max^d(iw),"for variable iw=",iw," and side iB=",2*^d
1308 end select
1309 end do
1310 }
1311
1312 ! psi, tracers take the same boundary type as the first variable
1313 if (nwfluxbc<nw_recon) then
1314 do iw=nwfluxbc+1,nw_recon
1315 typeboundary(iw,:) = typeboundary(1, :)
1316 end do
1317 end if
1318 ! auxiliary variables take the same boundary type as the first variable
1319 if (nwaux>0) then
1320 do iw=nw_recon+1, nw_recon+nwaux
1321 typeboundary(iw,:) = typeboundary(1, :)
1322 end do
1323 end if
1324
1325 if (any(typeboundary == 0)) then
1326 call mpistop("Not all boundary conditions have been defined")
1327 end if
1328
1329 do idim=1,ndim
1330 periodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_periodic))
1331 aperiodb(idim)=(any(typeboundary(:,2*idim-1:2*idim)==bc_aperiodic))
1332 if (periodb(idim).or.aperiodb(idim)) then
1333 do iw=1,nw_recon
1334 if (typeboundary(iw,2*idim-1) .ne. typeboundary(iw,2*idim)) &
1335 call mpistop("Wrong counterpart in periodic boundary")
1336
1337 if (typeboundary(iw,2*idim-1) /= bc_periodic .and. &
1338 typeboundary(iw,2*idim-1) /= bc_aperiodic) then
1339 call mpistop("Each dimension should either have all "//&
1340 "or no variables periodic, some can be aperiodic")
1341 end if
1342 end do
1343 end if
1344 end do
1345 {^nooned
1346 do idim=1,ndim
1347 if(any(typeboundary(:,2*idim-1)==12)) then
1348 if(any(typeboundary(:,2*idim-1)/=12)) typeboundary(:,2*idim-1)=12
1349 if(phys_energy) then
1350 windex=2
1351 else
1352 windex=1
1353 end if
1354 typeboundary(:,2*idim-1)=bc_symm
1355 if(physics_type/='rho') then
1356 select case(coordinate)
1357 case(cylindrical)
1358 typeboundary(phi_+1,2*idim-1)=bc_asymm
1359 if(physics_type=='mhd') typeboundary(ndir+windex+phi_,2*idim-1)=bc_asymm
1360 case(spherical)
1361 typeboundary(3:ndir+1,2*idim-1)=bc_asymm
1362 if(physics_type=='mhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim-1)=bc_asymm
1363 case default
1364 call mpistop('Pole is in cylindrical, polar, spherical coordinates!')
1365 end select
1366 end if
1367 end if
1368 if(any(typeboundary(:,2*idim)==12)) then
1369 if(any(typeboundary(:,2*idim)/=12)) typeboundary(:,2*idim)=12
1370 if(phys_energy) then
1371 windex=2
1372 else
1373 windex=1
1374 end if
1375 typeboundary(:,2*idim)=bc_symm
1376 if(physics_type/='rho') then
1377 select case(coordinate)
1378 case(cylindrical)
1379 typeboundary(phi_+1,2*idim)=bc_asymm
1380 if(physics_type=='mhd') typeboundary(ndir+windex+phi_,2*idim)=bc_asymm
1381 case(spherical)
1382 typeboundary(3:ndir+1,2*idim)=bc_asymm
1383 if(physics_type=='mhd') typeboundary(ndir+windex+2:ndir+windex+ndir,2*idim)=bc_asymm
1384 case default
1385 call mpistop('Pole is in cylindrical, polar, spherical coordinates!')
1386 end select
1387 end if
1388 end if
1389 end do
1390 }
1391
1392 if(.not.phys_energy) then
1393 flatcd=.false.
1394 flatsh=.false.
1395 end if
1396
1397 if(any(limiter(1:nlevelshi)=='mp5')) then
1398 nghostcells=max(nghostcells,3)
1399 end if
1400
1401 if(any(limiter(1:nlevelshi)=='weno5')) then
1402 nghostcells=max(nghostcells,3)
1403 end if
1404
1405 if(any(limiter(1:nlevelshi)=='weno5nm')) then
1406 nghostcells=max(nghostcells,3)
1407 end if
1408
1409 if(any(limiter(1:nlevelshi)=='wenoz5')) then
1410 nghostcells=max(nghostcells,3)
1411 end if
1412
1413 if(any(limiter(1:nlevelshi)=='wenoz5nm')) then
1414 nghostcells=max(nghostcells,3)
1415 end if
1416
1417 if(any(limiter(1:nlevelshi)=='wenozp5')) then
1418 nghostcells=max(nghostcells,3)
1419 end if
1420
1421 if(any(limiter(1:nlevelshi)=='wenozp5nm')) then
1422 nghostcells=max(nghostcells,3)
1423 end if
1424
1425 if(any(limiter(1:nlevelshi)=='teno5ad')) then
1426 nghostcells=max(nghostcells,3)
1427 end if
1428
1429 if(any(limiter(1:nlevelshi)=='weno5cu6')) then
1430 nghostcells=max(nghostcells,3)
1431 end if
1432
1433 if(any(limiter(1:nlevelshi)=='ppm')) then
1434 if(flatsh .or. flatcd) then
1435 nghostcells=max(nghostcells,4)
1436 else
1437 nghostcells=max(nghostcells,3)
1438 end if
1439 end if
1440
1441 if(any(limiter(1:nlevelshi)=='weno7')) then
1442 nghostcells=max(nghostcells,4)
1443 end if
1444
1445 if(any(limiter(1:nlevelshi)=='mpweno7')) then
1446 nghostcells=max(nghostcells,4)
1447 end if
1448
1449 ! If a wider stencil is used, extend the number of ghost cells
1450 nghostcells = nghostcells + phys_wider_stencil
1451
1452 ! prolongation in AMR for constrained transport MHD needs even number ghosts
1453 if(stagger_grid .and. refine_max_level>1 .and. mod(nghostcells,2)/=0) then
1454 nghostcells=nghostcells+1
1455 end if
1456
1457 select case (coordinate)
1458 {^nooned
1459 case (spherical)
1460 xprob^lim^de=xprob^lim^de*two*dpi;
1461 \}
1462 case (cylindrical)
1463 {
1464 if (^d==phi_) then
1465 xprob^lim^d=xprob^lim^d*two*dpi;
1466 end if
1467 \}
1468 end select
1469
1470 ! full block size including ghostcells
1471 {ixghi^d = block_nx^d + 2*nghostcells\}
1472 {ixgshi^d = ixghi^d\}
1473
1474 nx_vec = [{domain_nx^d|, }]
1475 block_nx_vec = [{block_nx^d|, }]
1476
1477 if (any(nx_vec < 4) .or. any(mod(nx_vec, 2) == 1)) &
1478 call mpistop('Grid size (domain_nx^D) has to be even and >= 4')
1479
1480 if (any(block_nx_vec < 4) .or. any(mod(block_nx_vec, 2) == 1)) &
1481 call mpistop('Block size (block_nx^D) has to be even and >= 4')
1482
1483 { if(mod(domain_nx^d,block_nx^d)/=0) &
1484 call mpistop('Grid (domain_nx^D) and block (block_nx^D) must be consistent') \}
1485
1486 if(refine_max_level>nlevelshi.or.refine_max_level<1)then
1487 write(unitterm,*)'Error: refine_max_level',refine_max_level,'>nlevelshi ',nlevelshi
1488 call mpistop("Reset nlevelshi and recompile!")
1489 endif
1490
1491 if (any(stretched_dim)) then
1492 allocate(qstretch(0:nlevelshi,1:ndim),dxfirst(0:nlevelshi,1:ndim),&
1493 dxfirst_1mq(0:nlevelshi,1:ndim),dxmid(0:nlevelshi,1:ndim))
1494 allocate(nstretchedblocks(1:nlevelshi,1:ndim))
1495 qstretch(0:nlevelshi,1:ndim)=0.0d0
1496 dxfirst(0:nlevelshi,1:ndim)=0.0d0
1497 nstretchedblocks(1:nlevelshi,1:ndim)=0
1498 {if (stretch_type(^d) == stretch_uni) then
1499 ! first some sanity checks
1500 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) then
1501 if(mype==0) then
1502 write(*,*) 'stretched grid needs finite qstretch_baselevel>1'
1503 write(*,*) 'will try default value for qstretch_baselevel in dimension', ^d
1504 endif
1505 if(xprobmin^d>smalldouble)then
1506 qstretch_baselevel(^d)=(xprobmax^d/xprobmin^d)**(1.d0/dble(domain_nx^d))
1507 else
1508 call mpistop("can not set qstretch_baselevel automatically")
1509 endif
1510 endif
1511 if(mod(block_nx^d,2)==1) &
1512 call mpistop("stretched grid needs even block size block_nxD")
1513 if(mod(domain_nx^d/block_nx^d,2)/=0) &
1514 call mpistop("number level 1 blocks in D must be even")
1515 qstretch(1,^d)=qstretch_baselevel(^d)
1516 dxfirst(1,^d)=(xprobmax^d-xprobmin^d) &
1517 *(1.0d0-qstretch(1,^d))/(1.0d0-qstretch(1,^d)**domain_nx^d)
1518 qstretch(0,^d)=qstretch(1,^d)**2
1519 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1520 if(refine_max_level>1)then
1521 do ilev=2,refine_max_level
1522 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1523 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1524 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1525 enddo
1526 endif
1527 endif \}
1528 if(mype==0) then
1529 {if(stretch_type(^d) == stretch_uni) then
1530 write(*,*) 'Stretched dimension ', ^d
1531 write(*,*) 'Using stretched grid with qs=',qstretch(0:refine_max_level,^d)
1532 write(*,*) ' and first cell sizes=',dxfirst(0:refine_max_level,^d)
1533 endif\}
1534 end if
1535 {if(stretch_type(^d) == stretch_symm) then
1536 if(mype==0) then
1537 write(*,*) 'will apply symmetric stretch in dimension', ^d
1538 endif
1539 if(mod(block_nx^d,2)==1) &
1540 call mpistop("stretched grid needs even block size block_nxD")
1541 ! checks on the input variable nstretchedblocks_baselevel
1542 if(nstretchedblocks_baselevel(^d)==0) &
1543 call mpistop("need finite even number of stretched blocks at baselevel")
1544 if(mod(nstretchedblocks_baselevel(^d),2)==1) &
1545 call mpistop("need even number of stretched blocks at baselevel")
1546 if(qstretch_baselevel(^d)<1.0d0.or.qstretch_baselevel(^d)==bigdouble) &
1547 call mpistop('stretched grid needs finite qstretch_baselevel>1')
1548 ! compute stretched part to ensure uniform center
1549 ipower=(nstretchedblocks_baselevel(^d)/2)*block_nx^d
1550 if(nstretchedblocks_baselevel(^d)==domain_nx^d/block_nx^d)then
1551 xstretch^d=0.5d0*(xprobmax^d-xprobmin^d)
1552 else
1553 xstretch^d=(xprobmax^d-xprobmin^d) &
1554 /(2.0d0+dble(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d) &
1555 *(1.0d0-qstretch_baselevel(^d))/(1.0d0-qstretch_baselevel(^d)**ipower))
1556 endif
1557 if(xstretch^d>(xprobmax^d-xprobmin^d)*0.5d0) &
1558 call mpistop(" stretched grid part should not exceed full domain")
1559 dxfirst(1,^d)=xstretch^d*(1.0d0-qstretch_baselevel(^d)) &
1560 /(1.0d0-qstretch_baselevel(^d)**ipower)
1561 nstretchedblocks(1,^d)=nstretchedblocks_baselevel(^d)
1562 qstretch(1,^d)=qstretch_baselevel(^d)
1563 qstretch(0,^d)=qstretch(1,^d)**2
1564 dxfirst(0,^d)=dxfirst(1,^d)*(1.0d0+qstretch(1,^d))
1565 dxmid(1,^d)=dxfirst(1,^d)
1566 dxmid(0,^d)=dxfirst(1,^d)*2.0d0
1567 if(refine_max_level>1)then
1568 do ilev=2,refine_max_level
1569 nstretchedblocks(ilev,^d)=2*nstretchedblocks(ilev-1,^d)
1570 qstretch(ilev,^d)=dsqrt(qstretch(ilev-1,^d))
1571 dxfirst(ilev,^d)=dxfirst(ilev-1,^d) &
1572 /(1.0d0+dsqrt(qstretch(ilev-1,^d)))
1573 dxmid(ilev,^d)=dxmid(ilev-1,^d)/2.0d0
1574 enddo
1575 endif
1576 ! sanity check on total domain size:
1577 sizeuniformpart^d=dxfirst(1,^d) &
1578 *(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
1579 if(mype==0) then
1580 print *,'uniform part of size=',sizeuniformpart^d
1581 print *,'setting of domain is then=',2*xstretch^d+sizeuniformpart^d
1582 print *,'versus=',xprobmax^d-xprobmin^d
1583 endif
1584 if(dabs(xprobmax^d-xprobmin^d-2*xstretch^d-sizeuniformpart^d)>smalldouble) then
1585 call mpistop('mismatch in domain size!')
1586 endif
1587 endif \}
1588 dxfirst_1mq(0:refine_max_level,1:ndim)=dxfirst(0:refine_max_level,1:ndim) &
1589 /(1.0d0-qstretch(0:refine_max_level,1:ndim))
1590 end if
1591
1592 dx_vec = [{xprobmax^d-xprobmin^d|, }] / nx_vec
1593
1594 if (mype==0) then
1595 write(c_ndim, '(I1)') ^nd
1596 write(unitterm, '(A30,' // c_ndim // '(I0," "))') &
1597 ' Domain size (cells): ', nx_vec
1598 write(unitterm, '(A30,' // c_ndim // '(E9.3," "))') &
1599 ' Level one dx: ', dx_vec
1600 end if
1601
1602 if (any(dx_vec < smalldouble)) &
1603 call mpistop("Incorrect domain size (too small grid spacing)")
1604
1605 dx(:, 1) = dx_vec
1606
1607 if(sum(w_refine_weight(:))==0) w_refine_weight(1) = 1.d0
1608 if(dabs(sum(w_refine_weight(:))-1.d0)>smalldouble) then
1609 write(unitterm,*) "Sum of all elements in w_refine_weight be 1.d0"
1610 call mpistop("Reset w_refine_weight so the sum is 1.d0")
1611 end if
1612
1613 select case (typeboundspeed)
1614 case('Einfeldt')
1615 boundspeed=1
1616 case('cmaxmean')
1617 boundspeed=2
1618 case('cmaxleftright')
1619 boundspeed=3
1620 case default
1621 call mpistop("set typeboundspeed='Einfieldt' or 'cmaxmean' or 'cmaxleftright'")
1622 end select
1623
1624 if (mype==0) write(unitterm, '(A30)', advance='no') 'Refine estimation: '
1625
1626 select case (refine_criterion)
1627 case (0)
1628 if (mype==0) write(unitterm, '(A)') "user defined"
1629 case (1)
1630 if (mype==0) write(unitterm, '(A)') "relative error"
1631 case (2)
1632 if (mype==0) write(unitterm, '(A)') "Lohner's original scheme"
1633 case (3)
1634 if (mype==0) write(unitterm, '(A)') "Lohner's scheme"
1635 case default
1636 call mpistop("Unknown error estimator, change refine_criterion")
1637 end select
1638
1639 if (tfixgrid<bigdouble/2.0d0) then
1640 if(mype==0)print*,'Warning, at time=',tfixgrid,'the grid will be fixed'
1641 end if
1642 if (itfixgrid<biginteger/2) then
1643 if(mype==0)print*,'Warning, at iteration=',itfixgrid,'the grid will be fixed'
1644 end if
1645 if (ditregrid>1) then
1646 if(mype==0)print*,'Note, Grid is reconstructed once every',ditregrid,'iterations'
1647 end if
1648
1649
1650 do islice=1,nslices
1651 select case(slicedir(islice))
1652 {case(^d)
1653 if(slicecoord(islice)<xprobmin^d.or.slicecoord(islice)>xprobmax^d) &
1654 write(uniterr,*)'Warning in read_par_files: ', &
1655 'Slice ', islice, ' coordinate',slicecoord(islice),'out of bounds for dimension ',slicedir(islice)
1656 \}
1657 end select
1658 end do
1659
1660 if (mype==0) then
1661 write(unitterm, '(A30,A,A)') 'restart_from_file: ', ' ', trim(restart_from_file)
1662 write(unitterm, '(A30,L1)') 'converting: ', convert
1663 write(unitterm, '(A)') ''
1664 endif
1665
1666 deallocate(flux_scheme)
1667
1668 end subroutine read_par_files
1669
1670 !> Routine to find entries in a string
1671 subroutine get_fields_string(line, delims, n_max, fields, n_found, fully_read)
1672 !> The line from which we want to read
1673 character(len=*), intent(in) :: line
1674 !> A string with delimiters. For example delims = " ,'"""//char(9)
1675 character(len=*), intent(in) :: delims
1676 !> Maximum number of entries to read in
1677 integer, intent(in) :: n_max
1678 !> Number of entries found
1679 integer, intent(inout) :: n_found
1680 !> Fields in the strings
1681 character(len=*), intent(inout) :: fields(n_max)
1682 logical, intent(out), optional :: fully_read
1683
1684 integer :: ixs_start(n_max)
1685 integer :: ixs_end(n_max)
1686 integer :: ix, ix_prev
1687
1688 ix_prev = 0
1689 n_found = 0
1690
1691 do while (n_found < n_max)
1692 ! Find the starting point of the next entry (a non-delimiter value)
1693 ix = verify(line(ix_prev+1:), delims)
1694 if (ix == 0) exit
1695
1696 n_found = n_found + 1
1697 ixs_start(n_found) = ix_prev + ix ! This is the absolute position in 'line'
1698
1699 ! Get the end point of the current entry (next delimiter index minus one)
1700 ix = scan(line(ixs_start(n_found)+1:), delims) - 1
1701
1702 if (ix == -1) then ! If there is no last delimiter,
1703 ixs_end(n_found) = len(line) ! the end of the line is the endpoint
1704 else
1705 ixs_end(n_found) = ixs_start(n_found) + ix
1706 end if
1707
1708 fields(n_found) = line(ixs_start(n_found):ixs_end(n_found))
1709 ix_prev = ixs_end(n_found) ! We continue to search from here
1710 end do
1711
1712 if (present(fully_read)) then
1713 ix = verify(line(ix_prev+1:), delims)
1714 fully_read = (ix == 0) ! Are there only delimiters?
1715 end if
1716
1717 end subroutine get_fields_string
1718
1719 subroutine saveamrfile(ifile)
1720
1723 use mod_particles, only: write_particles_snapshot
1724 use mod_slice, only: write_slice
1725 use mod_collapse, only: write_collapsed
1727 integer:: ifile
1728
1729 select case (ifile)
1730 case (fileout_)
1731 ! Write .dat snapshot
1732 call write_snapshot()
1733
1734 ! Generate formatted output (e.g., VTK)
1736
1737 if(use_particles) call write_particles_snapshot()
1738
1740 case (fileslice_)
1741 call write_slice
1742 case (filecollapse_)
1743 call write_collapsed
1744 case (filelog_)
1745 select case (typefilelog)
1746 case ('default')
1747 call printlog_default
1748 case ('regression_test')
1750 case ('special')
1751 if (.not. associated(usr_print_log)) then
1752 call mpistop("usr_print_log not defined")
1753 else
1754 call usr_print_log()
1755 end if
1756 case default
1757 call mpistop("Error in SaveFile: Unknown typefilelog")
1758 end select
1759 case (fileanalysis_)
1760 if (associated(usr_write_analysis)) then
1761 call usr_write_analysis()
1762 end if
1763 case default
1764 write(*,*) 'No save method is defined for ifile=',ifile
1765 call mpistop("")
1766 end select
1767
1768 ! opedit: Flush stdout and stderr from time to time.
1769 flush(unit=unitterm)
1770
1771 end subroutine saveamrfile
1772
1773
1774 ! Check if a snapshot exists
1775 logical function snapshot_exists(ix)
1777 integer, intent(in) :: ix !< Index of snapshot
1778 character(len=std_len) :: filename
1779
1780 write(filename, "(a,i4.4,a)") trim(base_filename), ix, ".dat"
1781 inquire(file=trim(filename), exist=snapshot_exists)
1782 end function snapshot_exists
1783
1784 integer function get_snapshot_index(filename)
1785 character(len=*), intent(in) :: filename
1786 integer :: i
1787
1788 ! Try to parse index in restart_from_file string (e.g. basename0000.dat)
1789 i = len_trim(filename) - 7
1790 read(filename(i:i+3), '(I4)') get_snapshot_index
1791 end function get_snapshot_index
1792
1793
1794
1795 !> Write header for a snapshot
1796 !>
1797 !> If you edit the header, don't forget to update: snapshot_write_header(),
1798 !> snapshot_read_header(), doc/fileformat.md, tools/python/dat_reader.py
1799 subroutine snapshot_write_header(fh, offset_tree, offset_block)
1800 use mod_forest
1801 use mod_physics
1803 use mod_slice, only: slicenext
1805 integer, intent(in) :: fh !< File handle
1806 integer(kind=MPI_OFFSET_KIND), intent(in) :: offset_tree !< Offset of tree info
1807 integer(kind=MPI_OFFSET_KIND), intent(in) :: offset_block !< Offset of block data
1808 call snapshot_write_header1(fh, offset_tree, offset_block, cons_wnames, nw)
1809 end subroutine snapshot_write_header
1810
1811 !> Read header for a snapshot
1812 !>
1813 !> If you edit the header, don't forget to update: snapshot_write_header(),
1814 !> snapshot_read_header(), doc/fileformat.md, tools/python/dat_reader.py
1815 subroutine snapshot_read_header(fh, offset_tree, offset_block)
1816 use mod_forest
1818 use mod_physics, only: physics_type
1819 use mod_slice, only: slicenext
1820 integer, intent(in) :: fh !< File handle
1821 integer(MPI_OFFSET_KIND), intent(out) :: offset_tree !< Offset of tree info
1822 integer(MPI_OFFSET_KIND), intent(out) :: offset_block !< Offset of block data
1823
1824 double precision :: rbuf(ndim)
1825 double precision, allocatable :: params(:)
1826 integer :: i, version
1827 integer :: ibuf(ndim), iw
1828 integer :: er, n_par, tmp_int
1829 integer, dimension(MPI_STATUS_SIZE) :: st
1830 logical :: periodic(ndim)
1831 character(len=name_len), allocatable :: var_names(:), param_names(:)
1832 character(len=name_len) :: phys_name, geom_name
1833
1834 ! Version number
1835 call mpi_file_read(fh, version, 1, mpi_integer, st, er)
1836 if (all(compatible_versions /= version)) then
1837 call mpistop("Incompatible file version (maybe old format?)")
1838 end if
1839
1840 ! offset_tree
1841 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1842 offset_tree = ibuf(1)
1843
1844 ! offset_block
1845 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1846 offset_block = ibuf(1)
1847
1848 ! nw
1849 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1850 nw_found=ibuf(1)
1851 if (nw /= ibuf(1)) then
1852 write(*,*) "nw=",nw," and nw found in restart file=",ibuf(1)
1853 write(*,*) "Please be aware of changes in w at restart."
1854 !call mpistop("currently, changing nw at restart is not allowed")
1855 end if
1856
1857 ! ndir
1858 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1859 if (ibuf(1) /= ndir) then
1860 write(*,*) "ndir in restart file = ",ibuf(1)
1861 write(*,*) "ndir = ",ndir
1862 call mpistop("reset ndir to ndir in restart file")
1863 end if
1864
1865 ! ndim
1866 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1867 if (ibuf(1) /= ndim) then
1868 write(*,*) "ndim in restart file = ",ibuf(1)
1869 write(*,*) "ndim = ",ndim
1870 call mpistop("reset ndim to ndim in restart file")
1871 end if
1872
1873 ! levmax
1874 call mpi_file_read(fh, ibuf(1), 1, mpi_integer, st, er)
1875 if (ibuf(1) > refine_max_level) then
1876 write(*,*) "number of levels in restart file = ",ibuf(1)
1877 write(*,*) "refine_max_level = ",refine_max_level
1878 call mpistop("refine_max_level < num. levels in restart file")
1879 end if
1880
1881 ! nleafs
1882 call mpi_file_read(fh, nleafs, 1, mpi_integer, st, er)
1883
1884 ! nparents
1885 call mpi_file_read(fh, nparents, 1, mpi_integer, st, er)
1886
1887 ! it
1888 call mpi_file_read(fh, it, 1, mpi_integer, st, er)
1889
1890 ! global time
1891 call mpi_file_read(fh, global_time, 1, mpi_double_precision, st, er)
1892
1893 ! xprobmin^D
1894 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1895 if (maxval(abs(rbuf(1:ndim) - [ xprobmin^d ])) > 0) then
1896 write(*,*) "Error: xprobmin differs from restart data: ", rbuf(1:ndim)
1897 call mpistop("change xprobmin^D in par file")
1898 end if
1899
1900 ! xprobmax^D
1901 call mpi_file_read(fh,rbuf(1:ndim),ndim,mpi_double_precision,st,er)
1902 if (maxval(abs(rbuf(1:ndim) - [ xprobmax^d ])) > 0) then
1903 write(*,*) "Error: xprobmax differs from restart data: ", rbuf(1:ndim)
1904 call mpistop("change xprobmax^D in par file")
1905 end if
1906
1907 ! domain_nx^D
1908 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
1909 if (any(ibuf(1:ndim) /= [ domain_nx^d ])) then
1910 write(*,*) "Error: mesh size differs from restart data: ", ibuf(1:ndim)
1911 call mpistop("change domain_nx^D in par file")
1912 end if
1913
1914 ! block_nx^D
1915 call mpi_file_read(fh,ibuf(1:ndim), ndim, mpi_integer,st,er)
1916 if (any(ibuf(1:ndim) /= [ block_nx^d ])) then
1917 write(*,*) "Error: block size differs from restart data:", ibuf(1:ndim)
1918 call mpistop("change block_nx^D in par file")
1919 end if
1920
1921 ! From version 5, read more info about the grid
1922 if (version > 4) then
1923 call mpi_file_read(fh, periodic, ndim, mpi_logical, st, er)
1924 if ({periodic(^d) .and. .not.periodb(^d) .or. .not.periodic(^d) .and. periodb(^d)| .or. }) &
1925 call mpistop("change in periodicity in par file")
1926
1927 call mpi_file_read(fh, geom_name, name_len, mpi_character, st, er)
1928
1929 if (geom_name /= geometry_name(1:name_len)) then
1930 write(*,*) "type of coordinates in data is: ", geom_name
1931 call mpistop("select the correct coordinates in mod_usr.t file")
1932 end if
1933
1934 call mpi_file_read(fh, stagger_mark_dat, 1, mpi_logical, st, er)
1935 if (stagger_grid .and. .not. stagger_mark_dat .or. .not.stagger_grid.and.stagger_mark_dat) then
1936 write(*,*) "Warning: stagger grid flag differs from restart data:", stagger_mark_dat
1937 !call mpistop("change parameter to use stagger grid")
1938 end if
1939 end if
1940
1941 ! From version 4 onwards, the later parts of the header must be present
1942 if (version > 3) then
1943 ! w_names (not used here)
1944 allocate(var_names(nw_found))
1945 do iw = 1, nw_found
1946 call mpi_file_read(fh, var_names(iw), name_len, mpi_character, st, er)
1947 end do
1948
1949 ! Physics related information
1950 call mpi_file_read(fh, phys_name, name_len, mpi_character, st, er)
1951
1952 if (phys_name /= physics_type) then
1953! call mpistop("Cannot restart with a different physics type")
1954 end if
1955
1956 call mpi_file_read(fh, n_par, 1, mpi_integer, st, er)
1957 allocate(params(n_par))
1958 allocate(param_names(n_par))
1959 call mpi_file_read(fh, params, n_par, mpi_double_precision, st, er)
1960 call mpi_file_read(fh, param_names, name_len * n_par, mpi_character, st, er)
1961
1962 ! Read snapshotnext etc. for restarting
1963 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1964
1965 ! Only set snapshotnext if the user hasn't specified it
1966 if (snapshotnext == -1) snapshotnext = tmp_int
1967
1968 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1969 if (slicenext == -1) slicenext = tmp_int
1970
1971 call mpi_file_read(fh, tmp_int, 1, mpi_integer, st, er)
1972 if (collapsenext == -1) collapsenext = tmp_int
1973 else
1974 ! Guess snapshotnext from file name if not set
1975 if (snapshotnext == -1) &
1977 ! Set slicenext and collapsenext if not set
1978 if (slicenext == -1) slicenext = 0
1979 if (collapsenext == -1) collapsenext = 0
1980 end if
1981
1982 ! Still used in convert
1984
1985 end subroutine snapshot_read_header
1986
1988 use mod_forest
1990 use mod_physics
1993
1994 double precision, allocatable :: w_buffer(:)
1995 integer :: file_handle, igrid, Morton_no, iwrite
1996 integer :: ipe, ix_buffer(2*ndim+1), n_values
1997 integer :: ixO^L, n_ghost(2*ndim)
1998 integer :: ixOs^L,n_values_stagger
1999 integer :: iorecvstatus(MPI_STATUS_SIZE)
2000 integer :: ioastatus(MPI_STATUS_SIZE)
2001 integer :: igrecvstatus(MPI_STATUS_SIZE)
2002 integer :: istatus(MPI_STATUS_SIZE)
2003 integer(kind=MPI_OFFSET_KIND) :: offset_tree_info
2004 integer(kind=MPI_OFFSET_KIND) :: offset_block_data
2005 integer(kind=MPI_OFFSET_KIND) :: offset_offsets
2006 integer, allocatable :: block_ig(:, :)
2007 integer, allocatable :: block_lvl(:)
2008 integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)
2009 type(tree_node), pointer :: pnode
2010
2011 call mpi_barrier(icomm, ierrmpi)
2012
2013 ! Allocate send/receive buffer
2014 n_values = count_ix(ixg^ll) * nw
2015 if(stagger_grid) then
2016 n_values = n_values + count_ix(ixgs^ll) * nws
2017 end if
2018 allocate(w_buffer(n_values))
2019
2020 ! Allocate arrays with information about grid blocks
2021 allocate(block_ig(ndim, nleafs))
2022 allocate(block_lvl(nleafs))
2023 allocate(block_offset(nleafs+1))
2024
2025 ! master processor
2026 if (mype==0) then
2027 call create_output_file(file_handle, snapshotnext, ".dat")
2028
2029 ! Don't know offsets yet, we will write header again later
2030 offset_tree_info = -1
2031 offset_block_data = -1
2032 call snapshot_write_header(file_handle, offset_tree_info, &
2033 offset_block_data)
2034
2035 call mpi_file_get_position(file_handle, offset_tree_info, ierrmpi)
2036
2037 call write_forest(file_handle)
2038
2039 ! Collect information about the spatial index (ig^D) and refinement level
2040 ! of leaves
2041 do morton_no = morton_start(0), morton_stop(npe-1)
2042 igrid = sfc(1, morton_no)
2043 ipe = sfc(2, morton_no)
2044 pnode => igrid_to_node(igrid, ipe)%node
2045
2046 block_ig(:, morton_no) = [ pnode%ig^d ]
2047 block_lvl(morton_no) = pnode%level
2048 block_offset(morton_no) = 0 ! Will be determined later
2049 end do
2050
2051 call mpi_file_write(file_handle, block_lvl, size(block_lvl), &
2052 mpi_integer, istatus, ierrmpi)
2053
2054 call mpi_file_write(file_handle, block_ig, size(block_ig), &
2055 mpi_integer, istatus, ierrmpi)
2056
2057 ! Block offsets are currently unknown, but will be overwritten later
2058 call mpi_file_get_position(file_handle, offset_offsets, ierrmpi)
2059 call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
2060 mpi_offset, istatus, ierrmpi)
2061
2062 call mpi_file_get_position(file_handle, offset_block_data, ierrmpi)
2063
2064 ! Check whether data was written as expected
2065 if (offset_block_data - offset_tree_info /= &
2066 (nleafs + nparents) * size_logical + &
2067 nleafs * ((1+ndim) * size_int + 2 * size_int)) then
2068 if (mype == 0) then
2069 print *, "Warning: MPI_OFFSET type /= 8 bytes"
2070 print *, "This *could* cause problems when reading .dat files"
2071 end if
2072 end if
2073
2074 block_offset(1) = offset_block_data
2075 iwrite = 0
2076 end if
2077
2078 do morton_no=morton_start(mype), morton_stop(mype)
2079 igrid = sfc_to_igrid(morton_no)
2080 itag = morton_no
2081
2082 call block_shape_io(igrid, n_ghost, ixo^l, n_values)
2083 if(stagger_grid) then
2084 w_buffer(1:n_values) = pack(ps(igrid)%w(ixo^s, 1:nw), .true.)
2085 {ixosmin^d = ixomin^d -1\}
2086 {ixosmax^d = ixomax^d \}
2087 n_values_stagger= count_ix(ixos^l)*nws
2088 w_buffer(n_values+1:n_values+n_values_stagger) = pack(ps(igrid)%ws(ixos^s, 1:nws), .true.)
2089 n_values=n_values+n_values_stagger
2090 else
2091 w_buffer(1:n_values) = pack(ps(igrid)%w(ixo^s, 1:nw), .true.)
2092 end if
2093 ix_buffer(1) = n_values
2094 ix_buffer(2:) = n_ghost
2095
2096 if (mype /= 0) then
2097 call mpi_send(ix_buffer, 2*ndim+1, &
2098 mpi_integer, 0, itag, icomm, ierrmpi)
2099 call mpi_send(w_buffer, n_values, &
2100 mpi_double_precision, 0, itag, icomm, ierrmpi)
2101 else
2102 iwrite = iwrite+1
2103 call mpi_file_write(file_handle, ix_buffer(2:), &
2104 2*ndim, mpi_integer, istatus, ierrmpi)
2105 call mpi_file_write(file_handle, w_buffer, &
2106 n_values, mpi_double_precision, istatus, ierrmpi)
2107
2108 ! Set offset of next block
2109 block_offset(iwrite+1) = block_offset(iwrite) + &
2110 int(n_values, mpi_offset_kind) * size_double + &
2111 2 * ndim * size_int
2112 end if
2113 end do
2114
2115 ! Write data communicated from other processors
2116 if (mype == 0) then
2117 do ipe = 1, npe-1
2118 do morton_no=morton_start(ipe), morton_stop(ipe)
2119 iwrite=iwrite+1
2120 itag=morton_no
2121
2122 call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, ipe, itag, icomm,&
2123 igrecvstatus, ierrmpi)
2124 n_values = ix_buffer(1)
2125
2126 call mpi_recv(w_buffer, n_values, mpi_double_precision,&
2127 ipe, itag, icomm, iorecvstatus, ierrmpi)
2128
2129 call mpi_file_write(file_handle, ix_buffer(2:), &
2130 2*ndim, mpi_integer, istatus, ierrmpi)
2131 call mpi_file_write(file_handle, w_buffer, &
2132 n_values, mpi_double_precision, istatus, ierrmpi)
2133
2134 ! Set offset of next block
2135 block_offset(iwrite+1) = block_offset(iwrite) + &
2136 int(n_values, mpi_offset_kind) * size_double + &
2137 2 * ndim * size_int
2138 end do
2139 end do
2140
2141 ! Write block offsets (now we know them)
2142 call mpi_file_seek(file_handle, offset_offsets, mpi_seek_set, ierrmpi)
2143 call mpi_file_write(file_handle, block_offset(1:nleafs), nleafs, &
2144 mpi_offset, istatus, ierrmpi)
2145
2146 ! Write header again, now with correct offsets
2147 call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
2148 call snapshot_write_header(file_handle, offset_tree_info, &
2149 offset_block_data)
2150
2151 call mpi_file_close(file_handle, ierrmpi)
2152 end if
2153
2154 call mpi_barrier(icomm, ierrmpi)
2155 end subroutine write_snapshot
2156
2157 !> Routine to read in snapshots (.dat files). When it cannot recognize the
2158 !> file version, it will automatically try the 'old' reader.
2159 subroutine read_snapshot
2162 use mod_forest
2164 use mod_slice, only: slicenext
2167
2168 double precision :: ws(ixGs^T,1:ndim)
2169 double precision, allocatable :: w_buffer(:)
2170 double precision, dimension(:^D&,:), allocatable :: w
2171 integer :: ix_buffer(2*ndim+1), n_values, n_values_stagger
2172 integer :: ixO^L, ixOs^L
2173 integer :: file_handle, amode, igrid, Morton_no, iread
2174 integer :: istatus(MPI_STATUS_SIZE)
2175 integer :: iorecvstatus(MPI_STATUS_SIZE)
2176 integer :: ipe,inrecv,nrecv, file_version
2177 integer(MPI_OFFSET_KIND) :: offset_tree_info
2178 integer(MPI_OFFSET_KIND) :: offset_block_data
2179 logical :: fexist
2180
2181 if (mype==0) then
2182 inquire(file=trim(restart_from_file), exist=fexist)
2183 if (.not.fexist) call mpistop(trim(restart_from_file)//" not found!")
2184
2185 call mpi_file_open(mpi_comm_self,restart_from_file,mpi_mode_rdonly, &
2186 mpi_info_null,file_handle,ierrmpi)
2187 call mpi_file_read(file_handle, file_version, 1, mpi_integer, &
2188 istatus, ierrmpi)
2189 end if
2190
2191 call mpi_bcast(file_version,1,mpi_integer,0,icomm,ierrmpi)
2192
2193 if (all(compatible_versions /= file_version)) then
2194 if (mype == 0) print *, "Unknown version, trying old snapshot reader..."
2195 call mpi_file_close(file_handle,ierrmpi)
2196 call read_snapshot_old()
2197
2198 ! Guess snapshotnext from file name if not set
2199 if (snapshotnext == -1) &
2201 ! Set slicenext and collapsenext if not set
2202 if (slicenext == -1) slicenext = 0
2203 if (collapsenext == -1) collapsenext = 0
2204
2205 ! Still used in convert
2207
2208 return ! Leave this routine
2209 else if (mype == 0) then
2210 call mpi_file_seek(file_handle, 0_mpi_offset_kind, mpi_seek_set, ierrmpi)
2211 call snapshot_read_header(file_handle, offset_tree_info, &
2212 offset_block_data)
2213 end if
2214
2215 ! Share information about restart file
2216 call mpi_bcast(nw_found,1,mpi_integer,0,icomm,ierrmpi)
2217 call mpi_bcast(nleafs,1,mpi_integer,0,icomm,ierrmpi)
2218 call mpi_bcast(nparents,1,mpi_integer,0,icomm,ierrmpi)
2219 call mpi_bcast(it,1,mpi_integer,0,icomm,ierrmpi)
2220 call mpi_bcast(global_time,1,mpi_double_precision,0,icomm,ierrmpi)
2221
2222 call mpi_bcast(snapshotnext,1,mpi_integer,0,icomm,ierrmpi)
2223 call mpi_bcast(slicenext,1,mpi_integer,0,icomm,ierrmpi)
2224 call mpi_bcast(collapsenext,1,mpi_integer,0,icomm,ierrmpi)
2225 call mpi_bcast(stagger_mark_dat,1,mpi_logical,0,icomm,ierrmpi)
2226
2227 ! Allocate send/receive buffer
2228 n_values = count_ix(ixg^ll) * nw_found
2229 if(stagger_mark_dat) then
2230 n_values = n_values + count_ix(ixgs^ll) * nws
2231 end if
2232 allocate(w_buffer(n_values))
2233 allocate(w(ixg^t,1:nw_found))
2234
2236
2237 if (mype == 0) then
2238 call mpi_file_seek(file_handle, offset_tree_info, &
2239 mpi_seek_set, ierrmpi)
2240 end if
2241
2242 call read_forest(file_handle)
2243
2244 do morton_no=morton_start(mype),morton_stop(mype)
2245 igrid=sfc_to_igrid(morton_no)
2246 call alloc_node(igrid)
2247 end do
2248
2249 if (mype==0) then
2250 call mpi_file_seek(file_handle, offset_block_data, mpi_seek_set, ierrmpi)
2251
2252 iread = 0
2253 do ipe = 0, npe-1
2254 do morton_no=morton_start(ipe),morton_stop(ipe)
2255 iread=iread+1
2256 itag=morton_no
2257
2258 call mpi_file_read(file_handle,ix_buffer(1:2*ndim), 2*ndim, &
2259 mpi_integer, istatus,ierrmpi)
2260
2261 ! Construct ixO^L array from number of ghost cells
2262 {ixomin^d = ixmlo^d - ix_buffer(^d)\}
2263 {ixomax^d = ixmhi^d + ix_buffer(ndim+^d)\}
2264 n_values = count_ix(ixo^l) * nw_found
2265 if(stagger_mark_dat) then
2266 {ixosmin^d = ixomin^d - 1\}
2267 {ixosmax^d = ixomax^d\}
2268 n_values_stagger = n_values
2269 n_values = n_values + count_ix(ixos^l) * nws
2270 end if
2271
2272 call mpi_file_read(file_handle, w_buffer, n_values, &
2273 mpi_double_precision, istatus, ierrmpi)
2274
2275 if (mype == ipe) then ! Root task
2276 igrid=sfc_to_igrid(morton_no)
2277 block=>ps(igrid)
2278 if(stagger_mark_dat) then
2279 w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values_stagger), &
2280 shape(w(ixo^s, 1:nw_found)))
2281 if(stagger_grid) &
2282 ps(igrid)%ws(ixos^s,1:nws)=reshape(w_buffer(n_values_stagger+1:n_values), &
2283 shape(ws(ixos^s, 1:nws)))
2284 else
2285 w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values), &
2286 shape(w(ixo^s, 1:nw_found)))
2287 end if
2288 if (nw_found<nw) then
2289 if (associated(usr_transform_w)) then
2290 call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
2291 else
2292 ps(igrid)%w(ixo^s,1:nw_found)=w(ixo^s,1:nw_found)
2293 end if
2294 else if (nw_found>nw) then
2295 if (associated(usr_transform_w)) then
2296 call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
2297 else
2298 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2299 end if
2300 else
2301 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2302 end if
2303 else
2304 call mpi_send([ ixo^l, n_values ], 2*ndim+1, &
2305 mpi_integer, ipe, itag, icomm, ierrmpi)
2306 call mpi_send(w_buffer, n_values, &
2307 mpi_double_precision, ipe, itag, icomm, ierrmpi)
2308 end if
2309 end do
2310 end do
2311
2312 call mpi_file_close(file_handle,ierrmpi)
2313
2314 else ! mype > 0
2315
2316 do morton_no=morton_start(mype),morton_stop(mype)
2317 igrid=sfc_to_igrid(morton_no)
2318 block=>ps(igrid)
2319 itag=morton_no
2320
2321 call mpi_recv(ix_buffer, 2*ndim+1, mpi_integer, 0, itag, icomm,&
2322 iorecvstatus, ierrmpi)
2323 {ixomin^d = ix_buffer(^d)\}
2324 {ixomax^d = ix_buffer(ndim+^d)\}
2325 n_values = ix_buffer(2*ndim+1)
2326
2327 call mpi_recv(w_buffer, n_values, mpi_double_precision,&
2328 0, itag, icomm, iorecvstatus, ierrmpi)
2329
2330 if(stagger_mark_dat) then
2331 n_values_stagger = count_ix(ixo^l) * nw_found
2332 {ixosmin^d = ixomin^d - 1\}
2333 {ixosmax^d = ixomax^d\}
2334 w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values_stagger), &
2335 shape(w(ixo^s, 1:nw_found)))
2336 if(stagger_grid) &
2337 ps(igrid)%ws(ixos^s,1:nws)=reshape(w_buffer(n_values_stagger+1:n_values), &
2338 shape(ws(ixos^s, 1:nws)))
2339 else
2340 w(ixo^s, 1:nw_found) = reshape(w_buffer(1:n_values), &
2341 shape(w(ixo^s, 1:nw_found)))
2342 end if
2343 if (nw_found<nw) then
2344 if (associated(usr_transform_w)) then
2345 call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
2346 else
2347 ps(igrid)%w(ixo^s,1:nw_found)=w(ixo^s,1:nw_found)
2348 end if
2349 else if (nw_found>nw) then
2350 if (associated(usr_transform_w)) then
2351 call usr_transform_w(ixg^ll,ixm^ll,nw_found,w,ps(igrid)%x,ps(igrid)%w)
2352 else
2353 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2354 end if
2355 else
2356 ps(igrid)%w(ixo^s,1:nw)=w(ixo^s,1:nw)
2357 end if
2358 end do
2359 end if
2360
2361 call mpi_barrier(icomm,ierrmpi)
2362
2363 end subroutine read_snapshot
2364
2366 use mod_forest
2370
2371 double precision :: wio(ixG^T,1:nw)
2372 double precision :: eqpar_dummy(100)
2373 integer :: fh, igrid, Morton_no, iread
2374 integer :: levmaxini, ndimini, ndirini
2375 integer :: nwini, neqparini, nxini^D
2376 integer(kind=MPI_OFFSET_KIND) :: offset
2377 integer :: istatus(MPI_STATUS_SIZE)
2378 integer, allocatable :: iorecvstatus(:,:)
2379 integer :: ipe,inrecv,nrecv
2380 integer :: sendini(7+^ND)
2381 logical :: fexist
2382 character(len=80) :: filename
2383
2384 if (mype==0) then
2385 call mpi_file_open(mpi_comm_self,trim(restart_from_file), &
2386 mpi_mode_rdonly,mpi_info_null,fh,ierrmpi)
2387
2388 offset=-int(7*size_int+size_double,kind=mpi_offset_kind)
2389 call mpi_file_seek(fh,offset,mpi_seek_end,ierrmpi)
2390
2391 call mpi_file_read(fh,nleafs,1,mpi_integer,istatus,ierrmpi)
2393 call mpi_file_read(fh,levmaxini,1,mpi_integer,istatus,ierrmpi)
2394 call mpi_file_read(fh,ndimini,1,mpi_integer,istatus,ierrmpi)
2395 call mpi_file_read(fh,ndirini,1,mpi_integer,istatus,ierrmpi)
2396 call mpi_file_read(fh,nwini,1,mpi_integer,istatus,ierrmpi)
2397 call mpi_file_read(fh,neqparini,1,mpi_integer,istatus,ierrmpi)
2398 call mpi_file_read(fh,it,1,mpi_integer,istatus,ierrmpi)
2399 call mpi_file_read(fh,global_time,1,mpi_double_precision,istatus,ierrmpi)
2400
2401 ! check if settings are suitable for restart
2402 if (levmaxini>refine_max_level) then
2403 write(*,*) "number of levels in restart file = ",levmaxini
2404 write(*,*) "refine_max_level = ",refine_max_level
2405 call mpistop("refine_max_level < number of levels in restart file")
2406 end if
2407 if (ndimini/=ndim) then
2408 write(*,*) "ndim in restart file = ",ndimini
2409 write(*,*) "ndim = ",ndim
2410 call mpistop("reset ndim to ndim in restart file")
2411 end if
2412 if (ndirini/=ndir) then
2413 write(*,*) "ndir in restart file = ",ndirini
2414 write(*,*) "ndir = ",ndir
2415 call mpistop("reset ndir to ndir in restart file")
2416 end if
2417 if (nw/=nwini) then
2418 write(*,*) "nw=",nw," and nw in restart file=",nwini
2419 call mpistop("currently, changing nw at restart is not allowed")
2420 end if
2421
2422 offset=offset-int(ndimini*size_int+neqparini*size_double,kind=mpi_offset_kind)
2423 call mpi_file_seek(fh,offset,mpi_seek_end,ierrmpi)
2424
2425 {call mpi_file_read(fh,nxini^d,1,mpi_integer,istatus,ierrmpi)\}
2426 if (ixghi^d/=nxini^d+2*nghostcells|.or.) then
2427 write(*,*) "Error: reset resolution to ",nxini^d+2*nghostcells
2428 call mpistop("change with setamrvac")
2429 end if
2430
2431 call mpi_file_read(fh,eqpar_dummy,neqparini, &
2432 mpi_double_precision,istatus,ierrmpi)
2433 end if
2434
2435 ! broadcast the global parameters first
2436 if (npe>1) then
2437 if (mype==0) then
2438 sendini=(/nleafs,levmaxini,ndimini,ndirini,nwini,neqparini,it ,^d&nxini^d /)
2439 end if
2440 call mpi_bcast(sendini,7+^nd,mpi_integer,0,icomm,ierrmpi)
2441 nleafs=sendini(1);levmaxini=sendini(2);ndimini=sendini(3);
2442 ndirini=sendini(4);nwini=sendini(5);
2443 neqparini=sendini(6);it=sendini(7);
2444 nxini^d=sendini(7+^d);
2446 call mpi_bcast(global_time,1,mpi_double_precision,0,icomm,ierrmpi)
2447 end if
2448
2449 if (mype == 0) then
2450 offset = int(size_block_io,kind=mpi_offset_kind) * &
2451 int(nleafs,kind=mpi_offset_kind)
2452 call mpi_file_seek(fh,offset,mpi_seek_set,ierrmpi)
2453 end if
2454
2455 call read_forest(fh)
2456
2457 do morton_no=morton_start(mype),morton_stop(mype)
2458 igrid=sfc_to_igrid(morton_no)
2459 call alloc_node(igrid)
2460 end do
2461
2462 if (mype==0)then
2463 iread=0
2464
2465 do morton_no=morton_start(0),morton_stop(0)
2466 igrid=sfc_to_igrid(morton_no)
2467 iread=iread+1
2468 offset=int(size_block_io,kind=mpi_offset_kind) &
2469 *int(morton_no-1,kind=mpi_offset_kind)
2470 call mpi_file_read_at(fh,offset,ps(igrid)%w,1,type_block_io, &
2471 istatus,ierrmpi)
2472 end do
2473 if (npe>1) then
2474 do ipe=1,npe-1
2475 do morton_no=morton_start(ipe),morton_stop(ipe)
2476 iread=iread+1
2477 itag=morton_no
2478 offset=int(size_block_io,kind=mpi_offset_kind)&
2479 *int(morton_no-1,kind=mpi_offset_kind)
2480 call mpi_file_read_at(fh,offset,wio,1,type_block_io,&
2481 istatus,ierrmpi)
2482 call mpi_send(wio,1,type_block_io,ipe,itag,icomm,ierrmpi)
2483 end do
2484 end do
2485 end if
2486 call mpi_file_close(fh,ierrmpi)
2487 else
2488 nrecv=(morton_stop(mype)-morton_start(mype)+1)
2489 allocate(iorecvstatus(mpi_status_size,nrecv))
2490 inrecv=0
2491 do morton_no=morton_start(mype),morton_stop(mype)
2492 igrid=sfc_to_igrid(morton_no)
2493 itag=morton_no
2494 inrecv=inrecv+1
2495 call mpi_recv(ps(igrid)%w,1,type_block_io,0,itag,icomm,&
2496 iorecvstatus(:,inrecv),ierrmpi)
2497 end do
2498 deallocate(iorecvstatus)
2499 end if
2500
2501 call mpi_barrier(icomm,ierrmpi)
2502
2503 end subroutine read_snapshot_old
2504
2505 !> Write volume-averaged values and other information to the log file
2507
2508 use mod_timing
2511
2512 double precision :: dtTimeLast, now, cellupdatesPerSecond
2513 double precision :: activeBlocksPerCore, wctPerCodeTime, timeToFinish
2514 double precision :: wmean(1:nw), total_volume
2515 double precision :: volume_coverage(refine_max_level)
2516 integer :: i, iw, level
2517 integer :: nx^D, nc, ncells, dit
2518 integer :: amode, istatus(MPI_STATUS_SIZE)
2519 integer, parameter :: my_unit = 20
2520 logical, save :: opened = .false.
2521 logical :: fileopen
2522 character(len=40) :: fmt_string
2523 character(len=80) :: filename
2524 character(len=2048) :: line
2525
2526 ! Compute the volume-average of w**1 = w
2527 call get_volume_average(1, wmean, total_volume)
2528
2529 ! Compute the volume coverage
2530 call get_volume_coverage(volume_coverage)
2531
2532 if (mype == 0) then
2533
2534 ! To compute cell updates per second, we do the following:
2535 nx^d=ixmhi^d-ixmlo^d+1;
2536 nc={nx^d*}
2537 ncells = nc * nleafs_active
2538
2539 ! assumes the number of active leafs haven't changed since last compute.
2540 now = mpi_wtime()
2541 dit = it - ittimelast
2542 dttimelast = now - timelast
2543 ittimelast = it
2544 timelast = now
2545 cellupdatespersecond = dble(ncells) * dble(nstep) * &
2546 dble(dit) / (dttimelast * dble(npe))
2547
2548 ! blocks per core:
2549 activeblockspercore = dble(nleafs_active) / dble(npe)
2550
2551 ! Wall clock time per code time unit in seconds:
2552 wctpercodetime = dttimelast / max(dit * dt, epsilon(1.0d0))
2553
2554 ! Wall clock time to finish in hours:
2555 timetofinish = (time_max - global_time) * wctpercodetime / 3600.0d0
2556
2557 ! On first entry, open the file and generate the header
2558 if (.not. opened) then
2559
2560 filename = trim(base_filename) // ".log"
2561
2562 ! Delete the log when not doing a restart run
2563 if (restart_from_file == undefined) then
2564 open(unit=my_unit,file=trim(filename),status='replace')
2565 close(my_unit, status='delete')
2566 end if
2567
2568 amode = ior(mpi_mode_create,mpi_mode_wronly)
2569 amode = ior(amode,mpi_mode_append)
2570
2571 call mpi_file_open(mpi_comm_self, filename, amode, &
2572 mpi_info_null, log_fh, ierrmpi)
2573
2574 opened = .true.
2575
2576 ! Start of file headern
2577 line = "it global_time dt"
2578 do level=1,nw
2579 i = len_trim(line) + 2
2580 write(line(i:),"(a,a)") trim(cons_wnames(level)), " "
2581 end do
2582
2583 ! Volume coverage per level
2584 do level = 1, refine_max_level
2585 i = len_trim(line) + 2
2586 write(line(i:), "(a,i0)") "c", level
2587 end do
2588
2589 ! Cell counts per level
2590 do level=1,refine_max_level
2591 i = len_trim(line) + 2
2592 write(line(i:), "(a,i0)") "n", level
2593 end do
2594
2595 ! Rest of file header
2596 line = trim(line) // " | Xload Xmemory 'Cell_Updates /second/core'"
2597 line = trim(line) // " 'Active_Blocks/Core' 'Wct Per Code Time [s]'"
2598 line = trim(line) // " 'TimeToFinish [hrs]'"
2599
2600 ! Only write header if not restarting
2601 if (restart_from_file == undefined .or. reset_time) then
2602 call mpi_file_write(log_fh, trim(line) // new_line('a'), &
2603 len_trim(line)+1, mpi_character, istatus, ierrmpi)
2604 end if
2605 end if
2606
2607 ! Construct the line to be added to the log
2608
2609 fmt_string = '(' // fmt_i // ',2' // fmt_r // ')'
2610 write(line, fmt_string) it, global_time, dt
2611 i = len_trim(line) + 2
2612
2613 write(fmt_string, '(a,i0,a)') '(', nw, fmt_r // ')'
2614 write(line(i:), fmt_string) wmean(1:nw)
2615 i = len_trim(line) + 2
2616
2617 write(fmt_string, '(a,i0,a)') '(', refine_max_level, fmt_r // ')'
2618 write(line(i:), fmt_string) volume_coverage(1:refine_max_level)
2619 i = len_trim(line) + 2
2620
2621 write(fmt_string, '(a,i0,a)') '(', refine_max_level, fmt_i // ')'
2622 write(line(i:), fmt_string) nleafs_level(1:refine_max_level)
2623 i = len_trim(line) + 2
2624
2625 fmt_string = '(a,6' // fmt_r2 // ')'
2626 write(line(i:), fmt_string) '| ', xload, xmemory, cellupdatespersecond, &
2627 activeblockspercore, wctpercodetime, timetofinish
2628
2629 call mpi_file_write(log_fh, trim(line) // new_line('a') , &
2630 len_trim(line)+1, mpi_character, istatus, ierrmpi)
2631 end if
2632
2633 end subroutine printlog_default
2634
2635 !> Print a log that can be used to check whether the code still produces the
2636 !> same output (regression test)
2639
2640 double precision :: modes(nw, 2), volume
2641 integer, parameter :: n_modes = 2
2642 integer :: power
2643 integer :: amode, istatus(MPI_STATUS_SIZE)
2644 logical, save :: file_open = .false.
2645 character(len=40) :: fmt_string
2646 character(len=2048) :: line
2647 character(len=80) :: filename
2648
2649 do power = 1, n_modes
2650 call get_volume_average(power, modes(:, power), volume)
2651 end do
2652
2653 if (mype == 0) then
2654 if (.not. file_open) then
2655 filename = trim(base_filename) // ".log"
2656 amode = ior(mpi_mode_create,mpi_mode_wronly)
2657 amode = ior(amode,mpi_mode_append)
2658
2659 call mpi_file_open(mpi_comm_self, filename, amode, &
2660 mpi_info_null, log_fh, ierrmpi)
2661 file_open = .true.
2662
2663 line= "# time mean(w) mean(w**2)"
2664 call mpi_file_write(log_fh, trim(line) // new_line('a'), &
2665 len_trim(line)+1, mpi_character, istatus, ierrmpi)
2666 end if
2667
2668 write(fmt_string, "(a,i0,a)") "(", nw * n_modes + 1, fmt_r // ")"
2669 write(line, fmt_string) global_time, modes
2670 call mpi_file_write(log_fh, trim(line) // new_line('a') , &
2671 len_trim(line)+1, mpi_character, istatus, ierrmpi)
2672 end if
2673 end subroutine printlog_regression_test
2674
2675 !> Compute mean(w**power) over the leaves of the grid. The first mode
2676 !> (power=1) corresponds to the mean, the second to the mean squared values
2677 !> and so on.
2678 subroutine get_volume_average(power, mode, volume)
2680
2681 integer, intent(in) :: power !< Which mode to compute
2682 double precision, intent(out) :: mode(nw) !< The computed mode
2683 double precision, intent(out) :: volume !< The total grid volume
2684
2685 double precision :: wsum(nw+1)
2686 double precision :: dsum_recv(1:nw+1)
2687 integer :: iigrid, igrid, iw
2688
2689 wsum(:) = 0
2690
2691 ! Loop over all the grids
2692 do iigrid = 1, igridstail
2693 igrid = igrids(iigrid)
2694
2695 ! Store total volume in last element
2696 wsum(nw+1) = wsum(nw+1) + sum(ps(igrid)%dvolume(ixm^t))
2697
2698 ! Compute the modes of the cell-centered variables, weighted by volume
2699 do iw = 1, nw
2700 wsum(iw) = wsum(iw) + &
2701 sum(ps(igrid)%dvolume(ixm^t)*ps(igrid)%w(ixm^t,iw)**power)
2702 end do
2703 end do
2704
2705 ! Make the information available on all tasks
2706 call mpi_allreduce(wsum, dsum_recv, nw+1, mpi_double_precision, &
2707 mpi_sum, icomm, ierrmpi)
2708
2709 ! Set the volume and the average
2710 volume = dsum_recv(nw+1)
2711 mode = dsum_recv(1:nw) / volume
2712
2713 end subroutine get_volume_average
2714
2715 !> Compute how much of the domain is covered by each grid level. This routine
2716 !> does not take a non-Cartesian geometry into account.
2717 subroutine get_volume_coverage(vol_cov)
2719
2720 double precision, intent(out) :: vol_cov(1:refine_max_level)
2721 double precision :: dsum_recv(1:refine_max_level)
2722 integer :: iigrid, igrid, iw, level
2723
2724 ! First determine the total 'flat' volume in each level
2725 vol_cov(1:refine_max_level)=zero
2726
2727 do iigrid = 1, igridstail
2728 igrid = igrids(iigrid);
2729 level = node(plevel_,igrid)
2730 vol_cov(level) = vol_cov(level)+ &
2731 {(rnode(rpxmax^d_,igrid)-rnode(rpxmin^d_,igrid))|*}
2732 end do
2733
2734 ! Make the information available on all tasks
2735 call mpi_allreduce(vol_cov, dsum_recv, refine_max_level, mpi_double_precision, &
2736 mpi_sum, icomm, ierrmpi)
2737
2738 ! Normalize
2739 vol_cov = dsum_recv / sum(dsum_recv)
2740 end subroutine get_volume_coverage
2741
2742 !> Compute the volume average of func(w) over the leaves of the grid.
2743 subroutine get_volume_average_func(func, f_avg, volume)
2745
2746 interface
2747 pure function func(w_vec, w_size) result(val)
2748 integer, intent(in) :: w_size
2749 double precision, intent(in) :: w_vec(w_size)
2750 double precision :: val
2751 end function func
2752 end interface
2753 double precision, intent(out) :: f_avg !< The volume average of func
2754 double precision, intent(out) :: volume !< The total grid volume
2755 double precision :: wsum(2)
2756 double precision :: dsum_recv(2)
2757 integer :: iigrid, igrid, i^D
2758
2759 wsum(:) = 0
2760
2761 ! Loop over all the grids
2762 do iigrid = 1, igridstail
2763 igrid = igrids(iigrid)
2764
2765 ! Store total volume in last element
2766 wsum(2) = wsum(2) + sum(ps(igrid)%dvolume(ixm^t))
2767
2768 ! Compute the modes of the cell-centered variables, weighted by volume
2769 {do i^d = ixmlo^d, ixmhi^d\}
2770 wsum(1) = wsum(1) + ps(igrid)%dvolume(i^d) * &
2771 func(ps(igrid)%w(i^d, :), nw)
2772 {end do\}
2773 end do
2774
2775 ! Make the information available on all tasks
2776 call mpi_allreduce(wsum, dsum_recv, 2, mpi_double_precision, &
2777 mpi_sum, icomm, ierrmpi)
2778
2779 ! Set the volume and the average
2780 volume = dsum_recv(2)
2781 f_avg = dsum_recv(1) / volume
2782
2783 end subroutine get_volume_average_func
2784
2785 !> Compute global maxima of iw variables over the leaves of the grid.
2786 subroutine get_global_maxima(wmax)
2788
2789 double precision, intent(out) :: wmax(nw) !< The global maxima
2790
2791 double precision :: wmax_mype(nw),wmax_recv(nw)
2792 integer :: iigrid, igrid, iw
2793
2794 wmax_mype(1:nw) = -bigdouble
2795
2796 ! Loop over all the grids
2797 do iigrid = 1, igridstail
2798 igrid = igrids(iigrid)
2799 do iw = 1, nw
2800 wmax_mype(iw)=max(wmax_mype(iw),maxval(ps(igrid)%w(ixm^t,iw)))
2801 end do
2802 end do
2803
2804 ! Make the information available on all tasks
2805 call mpi_allreduce(wmax_mype, wmax_recv, nw, mpi_double_precision, &
2806 mpi_max, icomm, ierrmpi)
2807
2808 wmax(1:nw)=wmax_recv(1:nw)
2809
2810 end subroutine get_global_maxima
2811
2812 !> Compute global minima of iw variables over the leaves of the grid.
2813 subroutine get_global_minima(wmin)
2815
2816 double precision, intent(out) :: wmin(nw) !< The global maxima
2817
2818 double precision :: wmin_mype(nw),wmin_recv(nw)
2819 integer :: iigrid, igrid, iw
2820
2821 wmin_mype(1:nw) = bigdouble
2822
2823 ! Loop over all the grids
2824 do iigrid = 1, igridstail
2825 igrid = igrids(iigrid)
2826 do iw = 1, nw
2827 wmin_mype(iw)=min(wmin_mype(iw),minval(ps(igrid)%w(ixm^t,iw)))
2828 end do
2829 end do
2830
2831 ! Make the information available on all tasks
2832 call mpi_allreduce(wmin_mype, wmin_recv, nw, mpi_double_precision, &
2833 mpi_min, icomm, ierrmpi)
2834
2835 wmin(1:nw)=wmin_recv(1:nw)
2836
2837 end subroutine get_global_minima
2838
2839end module mod_input_output
subroutine, public alloc_node(igrid)
allocate arrays on igrid node
Collapses D-dimensional output to D-1 view by line-of-sight integration.
Definition mod_collapse.t:4
subroutine write_collapsed
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine generate_plotfile
Module with basic grid data structures.
Definition mod_forest.t:2
integer, dimension(:), allocatable, save nleafs_level
How many leaves are present per refinement level.
Definition mod_forest.t:81
integer, dimension(:), allocatable, save sfc_to_igrid
Go from a Morton number to an igrid index (for a single processor)
Definition mod_forest.t:53
integer nleafs_active
Definition mod_forest.t:78
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition mod_forest.t:62
integer, save nleafs
Number of leaf block.
Definition mod_forest.t:76
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition mod_forest.t:65
integer, dimension(:,:), allocatable, save sfc
Array to go from a Morton number to an igrid and processor index. Sfc(1:3, MN) contains [igrid,...
Definition mod_forest.t:43
integer, save nparents
Number of parent blocks.
Definition mod_forest.t:73
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition mod_forest.t:32
subroutine, public write_forest(file_handle)
subroutine, public read_forest(file_handle)
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
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
double precision r_opt_thick
for spherical coordinate, region below it (unit=Rsun) is treated as not transparent
character(len=std_len) filename_sxr
Base file name for synthetic SXR emission output.
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.
integer ixghi
Upper index of grid block arrays.
character(len=std_len) geometry_name
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
logical activate_unit_arcsec
use arcsec as length unit of images/spectra
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) filename_whitelight
Base file name for synthetic white light.
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.
double precision, dimension(:), allocatable amr_wavefilter
refinement: lohner estimate wavefilter setting
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.
logical coarsenprimitive
coarsen primitive variables in level-jump ghost cells
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.
character(len=std_len) whitelight_instrument
white light observation instrument
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.
logical, dimension(ndim) collapse
If collapse(DIM) is true, generate output integrated over DIM.
double precision, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
double precision dt
global time step
integer refine_criterion
select types of refine criterion
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
Conversion factor for length unit.
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 of a physical block without 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
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
logical dat_resolution
resolution of the images
double precision r_occultor
the white light emission below it (unit=Rsun) is not visible
integer, dimension(ndim) nstretchedblocks_baselevel
(even) number of (symmetrically) stretched blocks at AMR level 1, per dimension
integer npe
The number of MPI tasks.
double precision, dimension(^nd) qstretch_baselevel
stretch factor between cells at AMR level 1, per dimension
integer nwauxio
Number of auxiliary variables that are only included in output.
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
if true write the w variable in output
logical prolongprimitive
prolongate primitive variables in level-jump ghost cells
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 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.
logical fix_small_values
fix small values with average or replace methods
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.
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.
double precision instrument_resolution_factor
times for enhancing spatial resolution for EUV image/spectra
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
Error tolerance ratio for derefinement decision.
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 (for dat resolution only)
double precision, dimension(1:3) x_origin
where the is the origin (X=0,Y=0) of image
double precision, dimension(^nd, 2) writespshift
domain percentage cut off shifted from each boundary when converting data
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.
subroutine, public create_output_file(fh, ix, extension, suffix)
Standard method for creating a new output file.
subroutine, public block_shape_io(igrid, n_ghost, ixol, n_values)
Determine the shape of a block for output (whether to include ghost cells, and on which sides)
subroutine, public snapshot_write_header1(fh, offset_tree, offset_block, dataset_names, nw_vars)
Write header for a snapshot, generalize cons_wnames and nw.
character(len=name_len) function, dimension(1:nwc), public get_names_from_string(aux_variable_names, nwc)
pure integer function, public count_ix(ixol)
Compute number of elements in index range.
Module for reading input and writing output.
subroutine saveamrfile(ifile)
logical function snapshot_exists(ix)
integer function get_snapshot_index(filename)
subroutine read_par_files()
Read in the user-supplied parameter-file.
character(len= *), parameter fmt_r
subroutine get_volume_average_func(func, f_avg, volume)
Compute the volume average of func(w) over the leaves of the grid.
subroutine get_global_minima(wmin)
Compute global minima of iw variables over the leaves of the grid.
subroutine get_global_maxima(wmax)
Compute global maxima of iw variables over the leaves of the grid.
subroutine printlog_regression_test()
Print a log that can be used to check whether the code still produces the same output (regression tes...
subroutine read_arguments()
Read the command line arguments passed to amrvac.
subroutine get_volume_average(power, mode, volume)
Compute mean(w**power) over the leaves of the grid. The first mode (power=1) corresponds to the mean,...
subroutine read_snapshot
Routine to read in snapshots (.dat files). When it cannot recognize the file version,...
integer nw_found
number of w found in dat files
character(len= *), parameter fmt_r2
subroutine get_fields_string(line, delims, n_max, fields, n_found, fully_read)
Routine to find entries in a string.
subroutine get_volume_coverage(vol_cov)
Compute how much of the domain is covered by each grid level. This routine does not take a non-Cartes...
character(len= *), parameter fmt_i
subroutine printlog_default
Write volume-averaged values and other information to the log file.
subroutine, public snapshot_write_header(fh, offset_tree, offset_block)
Write header for a snapshot.
subroutine snapshot_read_header(fh, offset_tree, offset_block)
Read header for a snapshot.
integer, dimension(3), parameter compatible_versions
List of compatible versions.
subroutine read_snapshot_old()
Module with slope/flux limiters.
Definition mod_limiter.t:2
double precision cada3_radius
radius of the asymptotic region [0.001, 10], larger means more accurate in smooth region but more ove...
Definition mod_limiter.t:13
double precision schmid_rad
Definition mod_limiter.t:14
Module containing all the particle routines.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
Definition mod_physics.t:17
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition mod_physics.t:50
logical phys_energy
Solve energy equation or not.
Definition mod_physics.t:35
Writes D-1 slice, can do so in various formats, depending on slice_type.
Definition mod_slice.t:2
integer slicenext
the file number of slices
Definition mod_slice.t:14
subroutine write_slice
Definition mod_slice.t:31
double precision, dimension(1000) slicecoord
Slice coordinates, see Slice output.
Definition mod_slice.t:8
character(len=std_len) slice_type
choose data type of slice: vtu, vtuCC, dat, or csv
Definition mod_slice.t:26
integer, dimension(nslicemax) slicedir
The slice direction for each slice.
Definition mod_slice.t:20
integer nslices
Number of slices to output.
Definition mod_slice.t:17
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)
Definition mod_source.t:2
integer ittimelast
Definition mod_timing.t:13
double precision timelast
Definition mod_timing.t:12
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.
Definition mod_forest.t:11