gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/driver.F90
Grzegorz Kowal 601adc8037 IO, DRIVER: Allow for restart with different number of ghost zones. 
Now we can restart from restart snapshots with different number of ghost
zones. Also, remove the old restart snapshot format support from
read_datablocks_h5().

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2019-02-11 11:02:36 -02:00

811 lines
22 KiB
Fortran
Raw Blame History

!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2008-2019 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!******************************************************************************
!!
!! program: AMUN
!!
!! AMUN is a code to perform numerical simulations in a fluid approximation on
!! adaptive mesh for Newtonian and relativistic environments with or without
!! magnetic field.
!!
!!
!!******************************************************************************
!
program amun
! include external subroutines used in this module
!
use iso_fortran_env, only : error_unit
use blocks , only : initialize_blocks, finalize_blocks, get_nleafs
use blocks , only : build_leaf_list
use boundaries , only : initialize_boundaries, finalize_boundaries
use boundaries , only : print_boundaries, boundary_variables
use coordinates , only : initialize_coordinates, finalize_coordinates
use coordinates , only : print_coordinates
use coordinates , only : nn => bcells
use domains , only : initialize_domains, finalize_domains
use equations , only : initialize_equations, finalize_equations
use equations , only : print_equations
use equations , only : nv
use evolution , only : initialize_evolution, finalize_evolution
use evolution , only : print_evolution
use evolution , only : advance, new_time_step
use evolution , only : step, time, dt
use gravity , only : initialize_gravity, finalize_gravity
use helpers , only : print_welcome, print_section, print_parameter
use integrals , only : initialize_integrals, finalize_integrals
use integrals , only : store_integrals
use interpolations , only : initialize_interpolations, finalize_interpolations
use interpolations , only : print_interpolations
use io , only : initialize_io, finalize_io, print_io
use io , only : restart_snapshot_number, restart_from_snapshot
use io , only : read_snapshot_parameter
use io , only : read_restart_snapshot, write_restart_snapshot
use io , only : write_snapshot, next_tout, precise_snapshots
use mesh , only : initialize_mesh, finalize_mesh
use mesh , only : generate_mesh, store_mesh_stats
use mpitools , only : initialize_mpitools, finalize_mpitools
#ifdef MPI
use mpitools , only : bcast_integer_variable
use mpitools , only : reduce_maximum_integer, reduce_sum_real_array
#endif /* MPI */
use mpitools , only : master, nprocs, nproc
use operators , only : initialize_operators, finalize_operators
use parameters , only : read_parameters, finalize_parameters
#ifdef MPI
use parameters , only : redistribute_parameters
#endif /* MPI */
use parameters , only : get_parameter
use problems , only : initialize_problems, finalize_problems
use random , only : initialize_random, finalize_random
use refinement , only : initialize_refinement, finalize_refinement
use refinement , only : print_refinement
use schemes , only : initialize_schemes, finalize_schemes
use schemes , only : print_schemes
use shapes , only : initialize_shapes, finalize_shapes, print_shapes
use sources , only : initialize_sources, finalize_sources
use sources , only : print_sources
use timers , only : initialize_timers, finalize_timers
use timers , only : start_timer, stop_timer, set_timer, get_timer
use timers , only : get_timer_total, timer_enabled, timer_description
use timers , only : get_count, ntimers
use user_problem , only : initialize_user_problem, finalize_user_problem
! module variables are not implicit by default
!
implicit none
! the success, termination and status flags
!
logical :: initialization_succeeded
integer :: iterm, iret
! the number of restarted runs
!
integer :: nrun = 0
! default parameters
!
character(len=64) :: problem = "none"
character(len=32) :: eqsys = "hydrodynamic"
character(len=32) :: eos = "adiabatic"
integer :: ncells = 8
integer :: nghosts = 2
integer :: toplev = 1
integer, dimension(3) :: bdims = 1
integer :: nmax = huge(1), ndat = 1
real(kind=8) :: xmin = 0.0d+00, xmax = 1.0d+00
real(kind=8) :: ymin = 0.0d+00, ymax = 1.0d+00
real(kind=8) :: zmin = 0.0d+00, zmax = 1.0d+00
real(kind=8) :: tmax = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01
real(kind=8) :: dtnext = 0.0d+00
! timer indices
!
integer :: iin, iev, itm
#ifdef PROFILE
integer :: ipr, ipi
#endif /* PROFILE */
! iteration and time variables
!
integer :: i, ed, eh, em, es, ec
integer :: nsteps = 1
character(len=80) :: fmt, tmp
real(kind=8) :: tbeg, thrs
real(kind=8) :: tm_curr, tm_exec, tm_conv, tm_last = 0.0d+00
#ifdef INTEL
! the type of the function SIGNAL should be defined for Intel compiler
!
integer(kind=4) :: signal
#endif /* INTEL */
#ifdef SIGNALS
! references to functions handling signals
!
#ifdef GNU
intrinsic signal
#endif /* GNU */
external terminate
! signal definitions
!
integer, parameter :: SIGERR = -1
integer, parameter :: SIGINT = 2, SIGABRT = 6, SIGTERM = 15
#endif /* SIGNALS */
! an array to store execution times
!
real(kind=8), dimension(ntimers) :: tm
! common block
!
common /termination/ iterm
!
!-------------------------------------------------------------------------------
!
! initialize the termination and return flags
!
iterm = 0
iret = 0
! initialize module TIMERS
!
call initialize_timers()
! set timer descriptions
!
call set_timer('INITIALIZATION', iin)
call set_timer('EVOLUTION' , iev)
call set_timer('TERMINATION' , itm)
#ifdef SIGNALS
! assign function terminate() to handle signals
!
#ifdef GNU
if (signal(SIGINT , terminate) == SIGERR) iret = 1
if (signal(SIGABRT, terminate) == SIGERR) iret = 2
if (signal(SIGTERM, terminate) == SIGERR) iret = 3
#endif /* GNU */
#ifdef INTEL
if (signal(SIGINT , terminate, -1) == SIGERR) iret = 1
if (signal(SIGABRT, terminate, -1) == SIGERR) iret = 2
if (signal(SIGTERM, terminate, -1) == SIGERR) iret = 3
#endif /* INTEL */
! in the case of problems with signal handler assignment, quit the program
!
if (iret > 0) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem initializing the signal handler!"
call finalize_timers()
call exit(1)
end if
#endif /* SIGNALS */
! initialize module MPITOOLS
!
call initialize_mpitools(iret)
! quit, in the case of problems with the MPI initialization
!
if (iret > 0) go to 400
! set the initiallization success flag
!
initialization_succeeded = .false.
! start time accounting for the initialization
!
call start_timer(iin)
! print welcome messages
!
call print_welcome(master)
call print_section(master, "Parallelization")
#ifdef MPI
call print_parameter(master, "MPI processes", nprocs)
#endif /* MPI */
! initialize and read parameters from the parameter file
!
if (master) call read_parameters(iterm)
#ifdef MPI
! broadcast the termination flag
!
call bcast_integer_variable(iterm, iret)
#endif /* MPI */
! check if the parameters were read successfully
!
if (iterm > 0 .or. iret > 0) then
if (master) then
write(error_unit,"('[AMUN::program]: ', a)") "Problem reading parameters!"
end if
go to 300
end if
#ifdef MPI
! redistribute parameters among all processors
!
call redistribute_parameters(iterm)
! reduce the termination flag over all processors to check if everything is fine
!
call reduce_maximum_integer(iterm, iret)
! check if the parameters were broadcasted successfully
!
if (iterm > 0 .or. iret > 0) then
if (master) then
write(error_unit,"('[AMUN::program]: ', a)") &
"Problem broadcasting parameters!"
end if
go to 300
end if
#endif /* MPI */
! initialize IO to handle restart snapshots if necessary
!
call initialize_io(master, iret)
if (iret > 0) go to 200
! get the run number
!
nrun = max(1, restart_snapshot_number() + 1)
! if the run is from a restarted job, read the fixed parameters from
! the restart snapshot, otherwise, read them from the parameter file
!
if (restart_from_snapshot()) then
call read_snapshot_parameter("problem", problem, iret)
call read_snapshot_parameter("eqsys" , eqsys , iret)
call read_snapshot_parameter("eos" , eos , iret)
call read_snapshot_parameter("ncells" , ncells , iret)
call read_snapshot_parameter("maxlev" , toplev , iret)
call read_snapshot_parameter("rdims" , bdims , iret)
call read_snapshot_parameter("xmin" , xmin , iret)
call read_snapshot_parameter("xmax" , xmax , iret)
call read_snapshot_parameter("ymin" , ymin , iret)
call read_snapshot_parameter("ymax" , ymax , iret)
#if NDIMS == 3
call read_snapshot_parameter("zmin" , zmin , iret)
call read_snapshot_parameter("zmax" , zmax , iret)
#endif /* NDIMS == 3 */
else
call get_parameter("problem" , problem)
call get_parameter("equation_system" , eqsys )
call get_parameter("equation_of_state", eos )
call get_parameter("ncells" , ncells )
call get_parameter("maxlev" , toplev )
call get_parameter("xblocks" , bdims(1))
call get_parameter("yblocks" , bdims(2))
#if NDIMS == 3
call get_parameter("zblocks" , bdims(3))
#endif /* NDIMS == 3 */
call get_parameter("xmin" , xmin )
call get_parameter("xmax" , xmax )
call get_parameter("ymin" , ymin )
call get_parameter("ymax" , ymax )
#if NDIMS == 3
call get_parameter("zmin" , zmin )
call get_parameter("zmax" , zmax )
#endif /* NDIMS == 3 */
end if
! get the number of ghost zones
!
call get_parameter("nghosts", nghosts)
! get the execution termination parameters
!
call get_parameter("nmax", nmax)
call get_parameter("tmax", tmax)
call get_parameter("trun", trun)
call get_parameter("tsav", tsav)
! correct the run time by the save time
!
trun = trun - tsav / 6.0d+01
! initialize dtnext
!
dtnext = 2.0d+00 * tmax
! get integral calculation interval
!
call get_parameter("ndat", ndat)
! initialize the remaining modules
!
call initialize_random(1, 0)
if (iret > 0) go to 190
call initialize_equations(eqsys, eos, master, iret)
if (iret > 0) go to 180
call initialize_coordinates(ncells, nghosts, toplev, bdims, xmin, xmax, &
ymin, ymax, zmin, zmax, master, iret)
if (iret > 0) go to 170
#if NDIMS == 3
call initialize_blocks((/ nv, nv, nn, nn, nn /), master, iret)
#else /* NDIMS == 3 */
call initialize_blocks((/ nv, nv, nn, nn, 1 /), master, iret)
#endif /* NDIMS == 3 */
if (iret > 0) go to 160
call initialize_operators(master, iret)
if (iret > 0) go to 150
call initialize_sources(master, iret)
if (iret > 0) go to 140
call initialize_user_problem(problem, master, iret)
if (iret > 0) go to 130
call initialize_problems(problem, master, iret)
if (iret > 0) go to 120
call initialize_domains(problem, master, iret)
if (iret > 0) go to 110
call initialize_boundaries(master, iret)
if (iret > 0) go to 100
call initialize_refinement(master, iret)
if (iret > 0) go to 90
call initialize_mesh(nrun, master, iret)
if (iret > 0) go to 80
call initialize_shapes(master, iret)
if (iret > 0) go to 70
call initialize_gravity(master, iret)
if (iret > 0) go to 60
call initialize_interpolations(master, iret)
if (iret > 0) go to 50
call initialize_schemes(master, iret)
if (iret > 0) go to 40
call initialize_evolution(master, iret)
if (iret > 0) go to 30
call initialize_integrals(master, nrun, iret)
if (iret > 0) go to 20
! print module information
!
call print_section(master, "Problem")
call print_parameter(master, "problem name", trim(problem))
call print_equations(master)
call print_sources(master)
call print_coordinates(master)
call print_boundaries(master)
call print_shapes(master)
call print_refinement(master)
call print_evolution(master)
call print_schemes(master)
call print_interpolations(master)
call print_io(master)
! check if we initiate new problem or restart previous job
!
if (restart_from_snapshot()) then
! reconstruct the meta and data block structures from a given restart file
!
call read_restart_snapshot(iterm)
#ifdef MPI
! reduce termination flag over all processors
!
call reduce_maximum_integer(iterm, iret)
#endif /* MPI */
! quit if there was a problem with reading restart snapshots
!
if (iterm > 0) go to 10
! update the list of leafs
!
call build_leaf_list()
! update boundaries
!
call boundary_variables(time, dtnext)
else
! generate the initial mesh, refine that mesh to the desired level according to
! the initialized problem
!
call generate_mesh()
! update boundaries
!
call boundary_variables(0.0d+00, dtnext)
! calculate new timestep
!
call new_time_step(dtnext)
end if
! store mesh statistics
!
call store_mesh_stats(step, time)
#ifdef MPI
! reduce termination flag over all processors
!
call reduce_maximum_integer(iterm, iret)
#endif /* MPI */
! check if the problem was successfully initialized or restarted
!
if (iterm > 0) go to 10
! store integrals and data to a file
!
if (.not. restart_from_snapshot()) then
call store_integrals()
call write_snapshot()
#ifdef MPI
! reduce termination flag over all processors
!
call reduce_maximum_integer(iterm, iret)
#endif /* MPI */
end if
! if the initial data were not successfully writen, exit the program
!
if (iterm > 0) go to 10
! reset the termination flag
!
iterm = 0
iret = 0
! stop time accounting for the initialization
!
call stop_timer(iin)
! set the initiallization success flag
!
initialization_succeeded = .true.
! start time accounting for the evolution
!
call start_timer(iev)
! get current time in seconds
!
if (master) &
tbeg = time
! print progress info on master processor
!
if (master) then
! initialize estimated remaining time of calculations
!
ed = 9999
eh = 23
em = 59
es = 59
! print progress info
!
write(*,*)
write(*,"(1x,a)" ) "Evolving the system:"
write(*,'(4x,a4,5x,a4,11x,a2,12x,a6,7x,a3)') 'step', 'time', 'dt' &
, 'blocks', 'ETA'
#ifdef INTEL
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
',1i2.2,"s",15x,a1,$)') &
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#else /* INTEL */
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
',1i2.2,"s",15x,a1)',advance="no") &
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#endif /* INTEL */
end if
! main loop
!
do while((nsteps <= nmax) .and. (time < tmax) .and. (iterm == 0))
! get the next snapshot time
!
if (precise_snapshots) dtnext = next_tout() - time
! performe one step evolution
!
call advance(dtnext)
! advance the iteration number and time
!
time = time + dt
step = step + 1
nsteps = nsteps + 1
! get current time in seconds
!
tm_curr = get_timer_total()
! compute elapsed time
!
thrs = tm_curr / 3.6d+03
! store mesh statistics
!
call store_mesh_stats(step, time)
! store integrals
!
call store_integrals()
! write down the restart snapshot
!
call write_restart_snapshot(thrs, nrun, iret)
! store data
!
call write_snapshot()
! check if the time exceeds execution time limit
!
if (thrs > trun) iterm = 100
! print progress info to console, but not too often
!
if (master) then
if (time >= tmax .or. (tm_curr - tm_last) >= 1.0d+00) then
! calculate days, hours, seconds
!
ec = int(tm_curr * (tmax - time) / max(1.0d-08, time - tbeg), kind = 4)
es = max(0, int(mod(ec, 60)))
em = int(mod(ec / 60, 60))
eh = int(ec / 3600)
ed = int(eh / 24)
eh = int(mod(eh, 24))
ed = min(9999,ed)
#ifdef INTEL
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
',1i2.2,"s",15x,a1,$)') &
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#else /* INTEL */
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
',1i2.2,"s",15x,a1)',advance="no") &
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
#endif /* INTEL */
! update the timestamp
!
tm_last = tm_curr
end if
end if
! prepare iterm
!
iterm = max(iterm, iret)
#ifdef MPI
! reduce termination flag over all processors
!
call reduce_maximum_integer(iterm, iret)
#endif /* MPI */
end do
! add one empty line
!
if (master) write(*,*)
! stop time accounting for the evolution
!
call stop_timer(iev)
! write down the restart snapshot
!
call write_restart_snapshot(1.0d+16, nrun, iret)
! a label to go to if there are any problems, but since all modules have been
! initialized, we have to finalize them first
!
10 continue
! start time accounting for the termination
!
call start_timer(itm)
! finalize modules
!
20 continue
call finalize_integrals(iret)
30 continue
call finalize_evolution(iret)
40 continue
call finalize_schemes(iret)
50 continue
call finalize_interpolations(iret)
60 continue
call finalize_gravity(iret)
70 continue
call finalize_shapes(iret)
80 continue
call finalize_mesh(iret)
90 continue
call finalize_refinement(iret)
100 continue
call finalize_boundaries(iret)
110 continue
call finalize_domains(iret)
120 continue
call finalize_problems(iret)
130 continue
call finalize_user_problem(iret)
140 continue
call finalize_sources(iret)
150 continue
call finalize_operators(iret)
160 continue
call finalize_blocks(iret)
170 continue
call finalize_coordinates(iret)
180 continue
call finalize_equations(iret)
190 continue
call finalize_random()
200 continue
call finalize_io(iret)
300 continue
call finalize_parameters()
! print execution times only if the initialization was fine
!
if (initialization_succeeded) then
! stop time accounting for the termination
!
call stop_timer(itm)
! get total time
!
tm(1) = get_timer_total()
! get subtasks timers
!
do i = 2, ntimers
tm(i) = get_timer(i)
end do
#ifdef MPI
! sum up timers from all processes
!
call reduce_sum_real_array(ntimers, tm(:), iret)
#endif /* MPI */
! print timings only on master processor
!
if (master) then
! calculate the conversion factor
!
tm_conv = 1.0d+02 / tm(1)
! print one empty line
!
write (*,'(a)') ''
! print the execution times
!
write (tmp,"(a)") "(2x,a32,1x,':',3x,1f16.3,' secs = ',f6.2,' %')"
write (*,'(1x,a)') 'EXECUTION TIMINGS'
do i = 2, ntimers
if (timer_enabled(i)) then
if (get_count(i) > 0) then
write (*,tmp) timer_description(i), tm(i), tm_conv * tm(i)
end if ! timer counter > 0
end if ! enabled
end do
! print the execution times
!
write (tmp,"(a)") "(1x,a14,20x,':',3x,1f16.3,' secs = ',f6.2,' %')"
write (*,tmp) 'TOTAL CPU TIME', tm(1) , 1.0d+02
write (*,tmp) 'TIME PER STEP ', tm(1) / nsteps, 1.0d+02 / nsteps
#ifdef MPI
write (*,tmp) 'TIME PER CPU ', tm(1) / nprocs, 1.0d+02 / nprocs
#endif /* MPI */
! get the execution time
!
tm_exec = get_timer_total()
! convert the execution time to days, hours, minutes, and seconds and print it
!
tm(1) = tm_exec / 8.64d+04
tm(2) = mod(tm_exec / 3.6d+03, 2.4d+01)
tm(3) = mod(tm_exec / 6.0d+01, 6.0d+01)
tm(4) = mod(tm_exec, 6.0d+01)
tm(5) = nint((tm_exec - floor(tm_exec)) * 1.0d+03)
write (tmp,"(a)") "(1x,a14,20x,':',3x,i14,'d'" // &
",i3.2,'h',i3.2,'m',i3.2,'.',i3.3,'s')"
write (*,tmp) 'EXECUTION TIME', int(tm(1:5))
end if
end if
! finalize module MPITOOLS
!
400 continue
call finalize_mpitools(iret)
! finalize module TIMERS
!
call finalize_timers()
!-------------------------------------------------------------------------------
!
end program
#ifdef SIGNALS
!
!===============================================================================
!
! function TERMINATE:
! ------------------
!
! Function sets variable iterm after receiving a signal.
!
!
!===============================================================================
!
integer(kind=4) function terminate(sig_num)
implicit none
! input arguments
!
integer(kind=4), intent(in) :: sig_num
integer :: iterm
! common block
!
common /termination/ iterm
!-------------------------------------------------------------------------------
!
#ifdef INTEL
iterm = sig_num
#else /* INTEL */
iterm = 15
#endif /* INTEL */
terminate = 1
!-------------------------------------------------------------------------------
!
end
#endif /* SIGNALS */
!===============================================================================
!
Reference in New Issue View Git Blame Copy Permalink