gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/amun.F90
Grzegorz Kowal dd3da22460 MESH: Move mesh statistics to be handled in INTEGRALS. 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2021-11-18 12:21:06 -03:00

558 lines
17 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-2021 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
use blocks , only : get_nleafs, build_leaf_list
use boundaries , only : boundary_variables
use evolution , only : advance, initialize_time_step, new_time_step
use evolution , only : step, time, dt, errtol
use helpers , only : print_welcome, print_section, print_parameter
use integrals , only : store_integrals
use io , only : initialize_io, finalize_io
use io , only : read_restart_snapshot, write_restart_snapshot
use io , only : write_snapshot, update_dtp
use iso_fortran_env, only : error_unit
use mesh , only : generate_mesh
use mpitools , only : initialize_mpitools, finalize_mpitools
#ifdef MPI
use mpitools , only : reduce_sum
#endif /* MPI */
use mpitools , only : master, nprocs, nproc, check_status
use parameters , only : read_parameters, finalize_parameters
use problem , only : initialize_problem, finalize_problem
use problem , only : print_problem_info, execute_problem
use problem , only : resumed, nrun, name, rngtype
use problem , only : tmax, trun, nwork, nmax
use random , only : initialize_random, finalize_random
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 : user_time_statistics
use workspace , only : initialize_workspace, finalize_workspace
implicit none
! the initialization and execution status flags
!
logical :: initialization_succeeded
logical :: verbose = .true.
logical :: proceed = .true.
integer :: quit = 0
integer :: status = 0
! timer indices, iteration number and other time variables
!
integer :: iin, iev, itm
integer :: i, ed, eh, em, es, ec
integer :: nsteps = 1
real(kind=8) :: tbeg = 0.0d+00, thrs
real(kind=8) :: tm_curr, tm_exec, tm_conv, tm_last = 0.0d+00
! the format string
!
character(len=80) :: sfmt
! vectors for improved estimation of the remaining simulation time
!
real(kind=8), dimension(8) :: tprv, tm_prev
! an array pointer for timings
!
real(kind=8), dimension(:), allocatable :: tm
#ifndef __GFORTRAN__
! signal is a subroutine only for the GNU compiler, otherwise it is a function
!
integer(kind=4) :: signal
#endif /* __GFORTRAN__ */
#ifdef SIGNALS
#ifdef __GFORTRAN__
intrinsic signal
#endif /* __GFORTRAN__ */
integer, parameter :: SIGERR = -1
integer, parameter :: SIGINT = 2, SIGABRT = 6, SIGTERM = 15
#endif /* SIGNALS */
character(len=*), parameter :: loc = 'AMUN::amun()'
!-------------------------------------------------------------------------------
!
#ifdef SIGNALS
! signals handling
!
#ifdef __GFORTRAN__
if (signal(SIGINT , signal_handler) == SIGERR) status = 1
if (signal(SIGABRT, signal_handler) == SIGERR) status = 2
if (signal(SIGTERM, signal_handler) == SIGERR) status = 3
#else /* __GFORTRAN__ */
if (signal(SIGINT , signal_handler, -1) == SIGERR) status = 1
if (signal(SIGABRT, signal_handler, -1) == SIGERR) status = 2
if (signal(SIGTERM, signal_handler, -1) == SIGERR) status = 3
#endif /* __GFORTRAN__ */
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not initialize the signal handling!"
call exit(status)
end if
#endif /* SIGNALS */
! timers
!
call initialize_timers()
call set_timer('INITIALIZATION', iin)
call set_timer('EVOLUTION' , iev)
call set_timer('TERMINATION' , itm)
! parallelization
!
call initialize_mpitools(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not initialize module MPITOOLS!"
end if
if (check_status(status /= 0)) call exit(status)
verbose = master
call start_timer(iin)
call print_welcome(verbose)
#ifdef MPI
call print_section(verbose, "Parallelization")
call print_parameter(verbose, "MPI processes", nprocs)
#endif /* MPI */
initialization_succeeded = .false.
! read parameters
!
call read_parameters(verbose, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not read parameters file!"
end if
if (check_status(status /= 0)) go to 6000
! initialize a few basic modules: IO, PROBLEM, RNG, WORKSPACE
!
call initialize_io(verbose, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not initialize module IO!"
end if
if (check_status(status /= 0)) go to 5000
call initialize_problem(verbose, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not initialize module PROBLEM!"
end if
if (check_status(status /= 0)) go to 4000
call initialize_random(rngtype, 1, 0, nproc, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not initialize module RANDOM!"
end if
if (check_status(status /= 0)) go to 3000
call initialize_workspace(nwork, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem initializing WORKSPACE module!"
end if
if (check_status(status /= 0)) go to 2000
! show some info before the simulation starts
!
call print_problem_info(verbose)
! initiate or resume the simulation
!
if (resumed) then
call read_restart_snapshot(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem reading the restart snapshot!"
end if
if (check_status(status /= 0)) go to 1000
call build_leaf_list(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem building the list of leafs!"
end if
if (check_status(status /= 0)) go to 1000
call boundary_variables(time, 0.0d+00)
else
call generate_mesh(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem generating the initial mesh!"
end if
if (check_status(status /= 0)) go to 1000
call boundary_variables(0.0d+00, 0.0d+00)
call initialize_time_step()
call store_integrals()
call user_time_statistics(time)
call write_snapshot(name)
end if
call stop_timer(iin)
initialization_succeeded = .true.
! time evolution starts here
!
call start_timer(iev)
! initiate the progress info
!
if (verbose) then
tbeg = time
tprv(:) = time
tm_prev(:) = get_timer_total()
ed = 9999
eh = 23
em = 59
es = 59
write(*,*)
write(*,"(1x,a)" ) "Evolving the system:"
write(*,"(4x,'step',5x,'time',11x,'timestep',6x,'err/tol',4x," // &
"'blocks',7x,'ETA')")
#ifdef __INTEL_COMPILER
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1,$)") &
step, time, dt, errtol, get_nleafs(), ed, eh, em, es, char(13)
#else /* __INTEL_COMPILER */
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1)",advance="no") &
step, time, dt, errtol, get_nleafs(), ed, eh, em, es, char(13)
#endif /* __INTEL_COMPILER */
end if
proceed = (nsteps <= nmax) .and. (time < tmax) .and. proceed
! the integration loop
!
do while(proceed)
call advance(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Advancing to the next step failed!"
end if
if (check_status(status /= 0)) go to 1000
! increase the iteration number and time
!
time = time + dt
step = step + 1
nsteps = nsteps + 1
! update the time step for the precise snapshots
!
call update_dtp()
! estimate the next time step
!
call new_time_step()
! store the time statistics and the run snapshot
!
call store_integrals()
call user_time_statistics(time)
call write_snapshot(name)
! estimate the elapsed time to trigger the restart snapshot storage
!
tm_curr = get_timer_total()
thrs = tm_curr / 3.6d+03
call write_restart_snapshot(thrs, name, nrun, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not store the restart snapshot!"
end if
if (check_status(status /= 0)) go to 1000
! check conditions for the integration interruption
!
proceed = (nsteps <= nmax) .and. (time < tmax) .and. &
(thrs <= trun) .and. .not. check_status(quit /= 0)
if (verbose) then
! update the vectors for the remaining time estimation
!
tm_prev(1) = tm_curr
tprv(1) = time
tm_prev = cshift(tm_prev, 1)
tprv = cshift(tprv , 1)
if (time >= tmax .or. (tm_curr - tm_last) >= 1.0d+00) then
ec = int((tm_curr - tm_prev(1)) * (tmax - time) &
/ max(1.0d-03, time - tprv(1)), kind = 4)
es = max(0, min(863999999, ec))
ed = es / 86400
es = es - 86400 * ed
eh = es / 3600
es = es - 3600 * eh
em = es / 60
es = es - 60 * em
! print the progress info
!
#ifdef __INTEL_COMPILER
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1,$)") &
step, time, dt, errtol, get_nleafs(), ed, eh, em, es, char(13)
#else /* __INTEL_COMPILER */
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1)",advance="no")&
step, time, dt, errtol, get_nleafs(), ed, eh, em, es, char(13)
#endif /* __INTEL_COMPILER */
tm_last = tm_curr
end if
end if
end do ! main loop
if (verbose) write(*,*)
call stop_timer(iev)
! store the final restart snapshot
!
call write_restart_snapshot(1.0d+16, name, nrun, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not store the restart snapshot!"
end if
1000 continue
call start_timer(itm)
! finalize all modules
!
call finalize_workspace(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem finalizing WORKSPACE module!"
end if
2000 continue
call finalize_random(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not finalize module RANDOM!"
end if
3000 continue
call finalize_problem(verbose, status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not finalize module PROBLEM!"
end if
4000 continue
call finalize_io(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem finalizing IO module!"
end if
5000 continue
call finalize_parameters()
call stop_timer(itm)
! print the execution time summary (only if the initialization was successful)
!
if (initialization_succeeded) then
allocate(tm(ntimers), stat=status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not allocate memory for timings!"
go to 6000
end if
tm(1) = get_timer_total()
do i = 2, ntimers
tm(i) = get_timer(i)
end do
#ifdef MPI
call reduce_sum(tm(1:ntimers))
#endif /* MPI */
if (verbose) then
write(* ,'(a)') ''
write(sfmt,"(a)") "(2x,a32,1x,':',3x,1f16.3,' secs = ',f6.2,' %')"
write(*,'(1x,a)') 'EXECUTION TIMINGS'
tm_conv = 1.0d+02 / tm(1)
do i = 2, ntimers
if (timer_enabled(i)) then
if (get_count(i) > 0) then
write(*,sfmt) timer_description(i), tm(i), tm_conv * tm(i)
end if
end if
end do
write(sfmt,"(a)") "(1x,a14,20x,':',3x,1f16.3,' secs = ',f6.2,' %')"
write(*,sfmt) 'TOTAL CPU TIME', tm(1) , 1.0d+02
write(*,sfmt) 'TIME PER STEP ', tm(1) / nsteps, 1.0d+02 / nsteps
#ifdef MPI
write(*,sfmt) 'TIME PER CPU ', tm(1) / nprocs, 1.0d+02 / nprocs
#endif /* MPI */
tm_exec = get_timer_total()
es = int(tm_exec)
ec = int(1000 * (tm_exec - es))
ed = es / 86400
es = es - 86400 * ed
eh = es / 3600
es = es - 3600 * eh
em = es / 60
es = es - 60 * em
write(sfmt,"(a)") "(1x,a14,20x,':',3x,i14,'d'" // &
",i3.2,'h',i3.2,'m',i3.2,'.',i3.3,'s')"
write(*,sfmt) 'EXECUTION TIME', ed, eh, em, es, ec
end if
deallocate(tm, stat=status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Could not deallocate memory for timings!"
end if
end if
6000 continue
call finalize_mpitools(status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Problem finalizing MPITOOLS module!"
end if
call finalize_timers()
#ifdef SIGNALS
contains
#ifdef __INTEL_COMPILER
!
!===============================================================================
!
! function SIGNAL_HANDLER:
! -----------------------
!
! Function sets the variable quit after receiving a signal.
!
! Arguments:
!
! sig_num - the number of the signal to be handled;
!
!===============================================================================
!
integer(kind=4) function signal_handler(sig_num)
use iso_fortran_env, only : error_unit
implicit none
integer(kind=4), intent(in) :: sig_num
!-------------------------------------------------------------------------------
!
quit = sig_num
write(error_unit,*)
write(error_unit,*) "Received signal:", iterm
write(error_unit,*) "Closing program..."
signal_handler = 1
!-------------------------------------------------------------------------------
!
end function
#else /* __INTEL_COMPILER */
!
!===============================================================================
!
! subroutine SIGNAL_HANDLER:
! -------------------------
!
! Subroutine sets the variable quit after receiving a signal.
!
!===============================================================================
!
subroutine signal_handler()
use iso_fortran_env, only : error_unit
implicit none
!-------------------------------------------------------------------------------
!
quit = 15
write(error_unit,*)
write(error_unit,*) "Received signal: 2, 9, or 15"
write(error_unit,*) "Closing program..."
return
!-------------------------------------------------------------------------------
!
end subroutine
#endif /* __INTEL_COMPILER */
#endif /* SIGNALS */
!===============================================================================
!
end program
Reference in New Issue View Git Blame Copy Permalink