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