!!****************************************************************************** !! !! 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 !! !! 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 . !! !!****************************************************************************** !! !! 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 coordinates , only : print_coordinates use coordinates , only : im, jm, km 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 : 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, 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 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 ! 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 ! 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 */ ! 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 400 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 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("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 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 call initialize_blocks((/ nv, nv, im, jm, km /), master, iret) 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_equations(master) call print_coordinates(master) if (master) then write (*,*) write (*,"(1x,a)" ) "Snapshots:" write (*,"(4x,a22,1x,'=',1x,a)") "precise snapshot times", trim(prec_snap) end if 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() 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) ! 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() 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) ! 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 ! 300 continue call finalize_parameters() ! finalize module MPITOOLS ! 400 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 */ !=============================================================================== !