gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/driver.F90
Grzegorz Kowal 2a6062cbd5 DRIVER, COORDINATES, IO: Pass more fixed parameters to initialize_coordinates(). 
Parameters such block dimensions, the base level dimensions, or the
domain extrema cannot be modified when the job is restarted, so read
them from the restart snapshots.

Also add new subroutine to IO read integer vector attribute from the
HDF5 snapshots.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2019-01-29 10:38:16 -02:00

950 lines
24 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 : boundary_variables
use coordinates , only : initialize_coordinates, finalize_coordinates
use equations , only : initialize_equations, finalize_equations
use evolution , only : initialize_evolution, finalize_evolution
use evolution , only : advance, new_time_step
use evolution , only : step, time, dt
use gravity , only : initialize_gravity, finalize_gravity
use integrals , only : initialize_integrals, finalize_integrals
use integrals , only : store_integrals
use interpolations , only : initialize_interpolations, finalize_interpolations
use io , only : initialize_io, finalize_io
use io , only : restart_from_snapshot, read_snapshot_parameter
use io , only : read_restart_snapshot, write_restart_snapshot
use io , only : write_snapshot, next_tout
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 schemes , only : initialize_schemes, finalize_schemes
use shapes , only : initialize_shapes, finalize_shapes
use sources , only : initialize_sources, finalize_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
! flat to identify if the problem is run from scratch or restarted
!
logical :: job_restart = .false.
integer :: nres = 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
! flag to adjust time precisely to the snapshots
!
logical , save :: precise_snapshots = .false.
character(len=255) :: prec_snap = "off"
! the termination and status flags
!
integer :: iterm, iret
! timer indices
!
integer :: iin, iev, itm
#ifdef PROFILE
integer :: ipr, ipi
#endif /* PROFILE */
! local snapshot file counters
!
integer :: nrun = 1
! 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)
! start time accounting for the initialization
!
call start_timer(iin)
#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 stop_timer(iin)
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) then
call stop_timer(iin)
go to 500
end if
! print the welcome message
!
if (master) then
write (*,"(1x,78('-'))")
write (*,"(1x,18('='),17x,a,17x,19('='))") 'A M U N'
write (*,"(1x,16('='),4x,a,4x,16('='))") &
'Copyright (C) 2008-2019 Grzegorz Kowal'
write (*,"(1x,18('='),9x,a,9x,19('='))") &
'under GNU GPLv3 license'
write (*,"(1x,78('-'))")
#ifdef MPI
! print the parallelization type and the number of parallel processes
!
write (*,*)
write (*,"(1x,a)" ) "Parallelization:"
write (*,"(4x,a,1x,i6 )" ) "MPI processes =", nprocs
#endif /* MPI */
end if
! 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 400
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 400
end if
#endif /* MPI */
! check if the job is restarted
!
call get_parameter("restart_number", nres)
job_restart = nres > 0
! if the run is from a restarted job, read the fixed parameters from
! the restart snapshot, otherwise, read them from the parameter file
!
if (job_restart) 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("nghosts", nghosts, 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("nghosts" , nghosts)
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 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 the precise snapshot flag
!
call get_parameter("precise_snapshots", prec_snap)
! set the precise snapshot flag
!
if (prec_snap == "on" ) precise_snapshots = .true.
if (prec_snap == "ON" ) precise_snapshots = .true.
if (prec_snap == "true") precise_snapshots = .true.
if (prec_snap == "TRUE") precise_snapshots = .true.
if (prec_snap == "yes" ) precise_snapshots = .true.
if (prec_snap == "YES" ) precise_snapshots = .true.
! get integral calculation interval
!
call get_parameter("ndat" , ndat)
! initialize the random number generator (passes the number of OpenMP threads
! and the current thread number)
!
call initialize_random(1, 0)
! initialize module USER_PROBLEM
!
call initialize_user_problem(problem, master, iret)
if (iret > 0) go to 340
! initialize module PROBLEMS
!
call initialize_problems(problem, master, iret)
if (iret > 0) go to 320
! initialize module EQUATIONS
!
call initialize_equations(eqsys, eos, master, iret)
if (iret > 0) go to 300
! initialize module SOURCES
!
call initialize_sources(master, iret)
if (iret > 0) go to 280
! initialize module GRAVITY
!
call initialize_gravity(master, iret)
if (iret > 0) go to 260
! initialize module COORDINATES
!
call initialize_coordinates(ncells, nghosts, toplev, bdims, xmin, xmax, &
ymin, ymax, zmin, zmax, master, iret)
if (iret > 0) go to 240
! initialize refinement module and print info
!
if (master) then
write (*,*)
write (*,"(1x,a)" ) "Refinement:"
end if
! jump to the end if the refinement could not be initialized
!
if (iret > 0) go to 50
! initialize module REFINEMENT
!
call initialize_refinement(master, iret)
! initialize methods modules and print info
!
if (master) then
write (*,*)
write (*,"(1x,a)" ) "Methods:"
end if
! initialize evolution
!
call initialize_evolution(master, iret)
! jump to the end if the schemes could not be initialized
!
if (iret > 0) go to 40
! initialize module SCHEMES
!
call initialize_schemes(master, iret)
! jump to the end if the schemes could not be initialized
!
if (iret > 0) go to 30
! initialize module INTERPOLATIONS
!
call initialize_interpolations(master, iret)
! initialize module OPERATORS
!
call initialize_operators(master, iret)
! initialize block module
!
call initialize_blocks(master, iret)
! initialize boundaries module and print info
!
if (master) then
write (*,*)
write (*,"(1x,a)" ) "Boundaries:"
end if
! initialize boundaries
!
call initialize_boundaries(master, iret)
! initialize module SHAPES
!
call initialize_shapes(master, iret)
! initialize boundaries module and print info
!
if (master) then
write (*,*)
write (*,"(1x,a)" ) "Snapshots:"
write (*,"(4x,a22,1x,'=',1x,a)") "precise snapshot times", trim(prec_snap)
end if
! initialize module IO
!
call initialize_io(master, nrun, iret)
! check if we initiate new problem or restart previous job
!
if (restart_from_snapshot()) then
! increase the run number
!
nrun = nrun + 1
! initialize the mesh module
!
call initialize_mesh(nrun, master, iret)
! quit if mesh couldn't be initialized
!
if (iret > 0) go to 10
! 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()
else
! initialize the mesh module
!
call initialize_mesh(nrun, master, iret)
! quit if mesh couldn't be initialized
!
if (iret > 0) go to 10
! 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
! initialize the integrals module
!
call initialize_integrals(master, nrun, iret)
! check if the module was successfully initialized
!
if (iret > 0) go to 10
! 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)
! 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 integrals module
!
call finalize_integrals()
! finalize I/O module
!
call finalize_io(iret)
! finalize module SHAPES
!
call finalize_shapes(iret)
! finalize the mesh module
!
call finalize_mesh(iret)
! finalize boundary module
!
call finalize_boundaries(iret)
! deallocate block structure
!
call finalize_blocks(iret)
! finalize the random number generator
!
call finalize_random()
! finalize module OPERATORS
!
call finalize_operators(iret)
! finalize module INTERPOLATIONS
!
call finalize_interpolations(iret)
! finalize module SCHEMES
!
call finalize_schemes(iret)
! jump point
!
30 continue
! finalize module EVOLUTION
!
call finalize_evolution(iret)
! jump point
!
40 continue
! finalize module REFINEMENT
!
call finalize_refinement(iret)
! jump point
!
50 continue
! finalize module COORDINATES
!
240 continue
call finalize_coordinates(iret)
! finalize module GRAVITY
!
260 continue
call finalize_gravity(iret)
! finalize module SOURCES
!
280 continue
call finalize_sources(iret)
! finalize module EQUATIONS
!
300 continue
call finalize_equations(iret)
! finalize module PROBLEMS
!
320 continue
call finalize_problems(iret)
! finalize module USER_PROBLEMS
!
340 continue
call finalize_user_problem(iret)
! 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
! finalize modules PARAMETERS
!
400 continue
call finalize_parameters()
! finalize module MPITOOLS
!
500 continue
call finalize_mpitools()
! 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