!!******************************************************************************
!!
!! 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-2020 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/>.
!!
!!******************************************************************************
!!
!! module: EVOLUTION
!!
!! This module provides an interface for temporal integration with
!! the stability handling. New integration methods can be added by
!! implementing more evolve_* subroutines.
!!
!!******************************************************************************
!
module evolution
#ifdef PROFILE
! import external subroutines
!
use timers, only : set_timer, start_timer, stop_timer
#endif /* PROFILE */
! module variables are not implicit by default
!
implicit none
#ifdef PROFILE
! timer indices
!
integer , save :: imi, ima, imt, imu, imf, iui, imv
#endif /* PROFILE */
! interfaces for procedure pointers
!
abstract interface
subroutine evolve_iface()
end subroutine
end interface
! pointer to the temporal integration subroutine
!
procedure(evolve_iface), pointer, save :: evolve => null()
! evolution parameters
!
character(len=255), save :: name_int = ""
integer , save :: stages = 2
real(kind=8) , save :: cfl = 5.0d-01
! coefficient controlling the decay of scalar potential ѱ
!
real(kind=8) , save :: alpha = 2.0d+00
real(kind=8) , save :: decay = 1.0d+00
! time variables
!
integer , save :: step = 0
real(kind=8) , save :: time = 0.0d+00
real(kind=8) , save :: dt = 1.0d+00
real(kind=8) , save :: dtn = 1.0d+00
! the increment array
!
real(kind=8), dimension(:,:,:,:), allocatable :: du
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_evolution, finalize_evolution, print_evolution
public :: advance, new_time_step
! declare public variables
!
public :: cfl, step, time, dt, dtn
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!!
!!*** PUBLIC SUBROUTINES *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_EVOLUTION:
! -------------------------------
!
! Subroutine initializes module EVOLUTION by setting its parameters.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
! status - an integer flag for error return value;
!
!===============================================================================
!
subroutine initialize_evolution(verbose, status)
! include external procedures
!
use coordinates , only : nn => bcells
use equations , only : nv
use iso_fortran_env, only : error_unit
use parameters , only : get_parameter
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
logical, intent(in) :: verbose
integer, intent(out) :: status
! local variables
!
character(len=255) :: integration = "rk2"
integer :: n
! local parameters
!
character(len=*), parameter :: loc = 'EVOLUTION::initialize_evolution()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! set timer descriptions
!
call set_timer('evolution:: initialization' , imi)
call set_timer('evolution:: solution advance', ima)
call set_timer('evolution:: new time step' , imt)
call set_timer('evolution:: solution update' , imu)
call set_timer('evolution:: flux update' , imf)
call set_timer('evolution:: increment update', iui)
call set_timer('evolution:: variable update' , imv)
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! get the integration method and the value of the CFL coefficient
!
call get_parameter("time_advance", integration)
call get_parameter("stages" , stages )
call get_parameter("cfl" , cfl )
call get_parameter("alpha" , alpha )
! select the integration method, check the correctness of the integration
! parameters and adjust the CFL coefficient if necessary
!
select case(trim(integration))
case ("euler", "EULER", "rk1", "RK1")
name_int = "1st order Euler"
evolve => evolve_euler
case ("rk2", "RK2")
name_int = "2nd order Runge-Kutta"
evolve => evolve_rk2
case ("ssprk(m,2)", "SSPRK(m,2)")
stages = max(2, min(9, stages))
write(name_int, "('2nd order SSPRK(',i0,',2)')") stages
evolve => evolve_ssprk2
cfl = (stages - 1) * cfl
case ("rk3", "RK3")
name_int = "3rd order Runge-Kutta"
evolve => evolve_rk3
case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)")
name_int = "3rd order SSPRK(4,3)"
evolve => evolve_ssprk34
cfl = 2.0d+00 * cfl
case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)")
name_int = "3rd order SSPRK(5,3)"
evolve => evolve_ssprk35
cfl = 2.65062919143939d+00 * cfl
case ("rk3.m", "ssprk(m,3)", "SSPRK(m,3)")
n = 2
do while(n**2 <= stages)
n = n + 1
end do
n = n - 1
stages = max(4, n**2)
write(name_int, "('Optimal 3rd order SSPRK(',i0,',3)')") stages
evolve => evolve_ssprk3_m
cfl = (n - 1) * n * cfl
case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)")
name_int = "Optimal 4th order SSPRK(10,4)"
evolve => evolve_ssprk4_10
cfl = 6.0d+00 * cfl
case default
if (verbose) then
write(*,*)
write(*,"(1x,a)") "ERROR!"
write(*,"(1x,a)") "The selected time advance method is not " // &
"implemented: " // trim(integration)
write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," // &
" 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," // &
" 'ssprk(m,3)', 'ssprk(10,4)'."
end if
status = 1
end select
! proceed if everything is fine
!
if (status == 0) then
! calculate the decay factor for magnetic field divergence scalar source term
!
decay = exp(- alpha * cfl)
! allocate space for the increment array
!
#if NDIMS == 3
allocate(du(nv,nn,nn,nn), stat = status)
#else /* NDIMS == 3 */
allocate(du(nv,nn,nn, 1), stat = status)
#endif /* NDIMS == 3 */
if (status == 0) then
du(:,:,:,:) = 0.0d+00
else
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot allocate memory for the increment array!"
end if
end if ! status
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine initialize_evolution
!
!===============================================================================
!
! subroutine FINALIZE_EVOLUTION:
! -----------------------------
!
! Subroutine releases memory used by the module.
!
! Arguments:
!
! status - an integer flag for error return value;
!
!===============================================================================
!
subroutine finalize_evolution(status)
! include external procedures
!
use iso_fortran_env, only : error_unit
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
integer, intent(out) :: status
! local parameters
!
character(len=*), parameter :: loc = 'EVOLUTION::finalize_evolution()'
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for module initialization/finalization
!
call start_timer(imi)
#endif /* PROFILE */
! reset the status flag
!
status = 0
! nullify pointer to integration subroutine
!
nullify(evolve)
! deallocate the increment array
!
if (allocated(du)) then
deallocate(du, stat = status)
if (status > 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc) &
, "Cannot deallocate memory for the increment array!"
end if
end if
#ifdef PROFILE
! stop accounting time for module initialization/finalization
!
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine finalize_evolution
!
!===============================================================================
!
! subroutine PRINT_EVOLUTION:
! --------------------------
!
! Subroutine prints module parameters.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
!
!===============================================================================
!
subroutine print_evolution(verbose)
! import external procedures and variables
!
use equations, only : magnetized
use helpers , only : print_section, print_parameter
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
logical, intent(in) :: verbose
!
!-------------------------------------------------------------------------------
!
if (verbose) then
call print_section(verbose, "Evolution")
call print_parameter(verbose, "time advance" , name_int)
call print_parameter(verbose, "CFL coefficient" , cfl )
if (magnetized) then
call print_parameter(verbose, "GLM alpha coefficient", alpha )
end if
end if
!-------------------------------------------------------------------------------
!
end subroutine print_evolution
!
!===============================================================================
!
! subroutine ADVANCE:
! ------------------
!
! Subroutine advances the solution by one time step using the selected time
! integration method.
!
! Arguments:
!
! dtnext - the next time step;
! status - the subroutine call status: 0 for success, otherwise failure;
!
!===============================================================================
!
subroutine advance(dtnext, status)
! references
!
use blocks , only : set_blocks_update
use coordinates, only : toplev
use forcing , only : update_forcing, forcing_enabled
use mesh , only : update_mesh
! local variables are not implicit by default
!
implicit none
! input variables
!
real(kind=8), intent(in) :: dtnext
integer , intent(out) :: status
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for solution advance
!
call start_timer(ima)
#endif /* PROFILE */
! find new time step
!
call new_time_step(dtnext)
! advance the solution using the selected method
!
call evolve()
! add forcing contribution
!
if (forcing_enabled) call update_forcing(time, dt)
! check if we need to perform the refinement step
!
if (toplev > 1) then
! set all meta blocks to not be updated
!
call set_blocks_update(.false.)
! check refinement and refine
!
call update_mesh(status)
if (status /= 0) go to 100
! update primitive variables
!
call update_variables(time + dt, 0.0d+00)
! set all meta blocks to be updated
!
call set_blocks_update(.true.)
end if ! toplev > 1
#ifdef DEBUG
! check variables for NaNs
!
call check_variables()
#endif /* DEBUG */
! error entry point
!
100 continue
#ifdef PROFILE
! stop accounting time for solution advance
!
call stop_timer(ima)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine advance
!
!===============================================================================
!
! subroutine NEW_TIME_STEP:
! ------------------------
!
! Subroutine estimates the new time step from the maximum speed in the system
! and source term constraints.
!
! Arguments:
!
! dtnext - next time step;
!
!===============================================================================
!
subroutine new_time_step(dtnext)
! include external procedures
!
use equations , only : maxspeed, cmax, cmax2
#ifdef MPI
use mpitools , only : reduce_maximum_real, reduce_maximum_integer
#endif /* MPI */
! include external variables
!
use blocks , only : block_data, list_data
use coordinates , only : adx, ady, adz
use coordinates , only : toplev
use sources , only : viscosity, resistivity
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: dtnext
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
integer :: iret, mlev
real(kind=8) :: cm, dx_min
! local parameters
!
real(kind=8), parameter :: eps = tiny(cmax)
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for new time step estimation
!
call start_timer(imt)
#endif /* PROFILE */
! reset the maximum speed, and the highest level
!
cmax = eps
mlev = 1
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks in order to find the maximum speed among them
! and the highest level which is required to obtain the minimum spacial step
!
do while (associated(pdata))
! update the maximum level
!
mlev = max(mlev, pdata%meta%level)
! get the maximum speed for the current block
!
cm = maxspeed(pdata%q(:,:,:,:))
! update the maximum speed
!
cmax = max(cmax, cm)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
#ifdef MPI
! reduce maximum speed and level over all processors
!
call reduce_maximum_real (cmax, iret)
call reduce_maximum_integer(mlev, iret)
#endif /* MPI */
! calculate the squared cmax
!
cmax2 = cmax * cmax
! find the smallest spacial step
!
#if NDIMS == 2
dx_min = min(adx(mlev), ady(mlev))
#endif /* NDIMS == 2 */
#if NDIMS == 3
dx_min = min(adx(mlev), ady(mlev), adz(mlev))
#endif /* NDIMS == 3 */
! calculate the new time step
!
dtn = cfl * dx_min / max(cmax &
+ 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)
! round the time
!
if (dtnext > 0.0d+00) dt = min(dtn, dtnext)
#ifdef PROFILE
! stop accounting time for new time step estimation
!
call stop_timer(imt)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine new_time_step
!
!===============================================================================
!!
!!*** PRIVATE SUBROUTINES ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine EVOLVE_EULER:
! -----------------------
!
! Subroutine advances the solution by one time step using the 1st order
! Euler integration method.
!
! References:
!
! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
! "Numerical Recipes in Fortran",
! Cambridge University Press, Cambridge, 1992
!
!===============================================================================
!
subroutine evolve_euler()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
! prepare times
!
tm = time + dt
dtm = dt
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the solution
!
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + dt * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u0
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_euler
!
!===============================================================================
!
! subroutine EVOLVE_RK2:
! ---------------------
!
! Subroutine advances the solution by one time step using the 2nd order
! Runge-Kutta time integration method.
!
! References:
!
! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
! "Numerical Recipes in Fortran",
! Cambridge University Press, Cambridge, 1992
!
!===============================================================================
!
subroutine evolve_rk2()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
!= 1st step: U(1) = U(n) + dt * F[U(n)]
!
! prepare times
!
tm = time + dt
dtm = dt
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + dt * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u1
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)]
!
! prepare times
!
tm = time + dt
dtm = 0.5d+00 * dt
! update fluxes from the intermediate stage
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = 0.5d+00 * (pdata%u0(:,:,:,:) &
+ pdata%u1(:,:,:,:) + dt * du(:,:,:,:))
! update the conservative variable pointer
!
pdata%u => pdata%u0
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_rk2
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK2:
! ------------------------
!
! Subroutine advances the solution by one time step using the 2nd order
! m-stage Strong Stability Preserving Runge-Kutta time integration method.
! Up to 9 stages are allowed, due to stability problems with more stages.
!
! References:
!
! [1] Gottlieb, S. and Gottlieb, L.-A., J.
! "Strong Stability Preserving Properties of Runge-Kutta Time
! Discretization Methods for Linear Constant Coefficient Operators",
! Journal of Scientific Computing,
! 2003, vol. 18, no. 1, pp. 83-109
!
!===============================================================================
!
subroutine evolve_ssprk2()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
integer :: n
real(kind=8) :: ds
real(kind=8) :: tm, dtm
! local saved variables
!
logical , save :: first = .true.
real(kind=8), save :: ft, fl, fr
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
! prepare things which don't change later
!
if (first) then
! calculate integration coefficients
!
ft = 1.0d+00 / (stages - 1)
fl = 1.0d+00 / stages
fr = 1.0d+00 - fl
! update first flag
!
first = .false.
end if
! calculate the fractional time step
!
ds = ft * dt
!= 1st step: U(0) = U(n)
!
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! copy conservative array u0 to u1
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u1
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
!= 2nd step: U(i) = [1 + dt/(m-1) L] U(i-1), for i = 1, ..., m-1
!
! integrate the intermediate steps
!
do n = 1, stages - 1
! prepare times
!
tm = time + n * ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
end do ! n = 1, stages - 1
!= the final step: U(n+1) = 1/m U(0) + (m-1)/m [1 + dt/(m-1) L] U(m-1)
!
! prepare times
!
tm = time + dt
dtm = fl * dt
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = fl * pdata%u0(:,:,:,:) &
+ fr * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:))
! update the conservative variable pointer
!
pdata%u => pdata%u0
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pointer to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk2
!
!===============================================================================
!
! subroutine EVOLVE_RK3:
! ---------------------
!
! Subroutine advances the solution by one time step using the 3rd order
! Runge-Kutta time integration method.
!
! References:
!
! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.,
! "Numerical Recipes in Fortran",
! Cambridge University Press, Cambridge, 1992
!
!
!===============================================================================
!
subroutine evolve_rk3()
! include external procedures
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
real(kind=8) :: ds
real(kind=8) :: tm, dtm
! local integration parameters
!
real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00
real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
!= 1st substep: U(1) = U(n) + dt F[U(n)]
!
! prepare the fractional time step
!
ds = dt
! prepare times
!
tm = time + ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the first intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u1
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)]
!
! prepare the fractional time step
!
ds = f22 * dt
! prepare times
!
tm = time + 0.5d+00 * dt
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the second intermediate solution
!
pdata%u1(:,:,:,:) = f21 * pdata%u0(:,:,:,:) + f22 * pdata%u1(:,:,:,:) &
+ ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)]
!
! prepare the fractional time step
!
ds = f32 * dt
! prepare times
!
tm = time + dt
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = f31 * pdata%u0(:,:,:,:) + f32 * pdata%u1(:,:,:,:) &
+ ds * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u0
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_rk3
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK34:
! -------------------------
!
! Subroutine advances the solution by one time step using the 3rd order
! 4-stage Strong Stability Preserving Runge-Kutta time integration method.
!
! References:
!
! [1] Ruuth, S. J.,
! "Global Optimization of Explicit Strong-Stability-Preserving
! Runge-Kutta methods",
! Mathematics of Computation,
! 2006, vol. 75, no. 253, pp. 183-207
!
!===============================================================================
!
subroutine evolve_ssprk34()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
real(kind=8) :: ds
real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
!
! calculate the fractional time step
!
ds = dt / 2.0d+00
! prepare times
!
tm = time + ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + dtm * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u1
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 2nd step: U(2) = U(2) + 1/2 dt F[U(1)]
!
! prepare times
!
tm = time + dt
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + dtm * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)])
!
! prepare times
!
tm = time + ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = (2.0d+00 * pdata%u0(:,:,:,:) + pdata%u1(:,:,:,:) &
+ dtm * du(:,:,:,:)) / 3.0d+00
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= the final step: U(n+1) = U(3) + 1/2 dt F[U(3)]
!
! prepare times
!
tm = time + dt
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u0
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk34
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK35:
! -------------------------
!
! Subroutine advances the solution by one time step using the 3rd order
! 5-stage Strong Stability Preserving Runge-Kutta time integration method.
!
! References:
!
! [1] Ruuth, S. J.,
! "Global Optimization of Explicit Strong-Stability-Preserving
! Runge-Kutta methods",
! Mathematics of Computation,
! 2006, vol. 75, no. 253, pp. 183-207
!
!===============================================================================
!
subroutine evolve_ssprk35()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
real(kind=8) :: ds
real(kind=8) :: tm, dtm
! local integration parameters
!
real(kind=8), parameter :: b1 = 3.77268915331368d-01
real(kind=8), parameter :: b3 = 2.42995220537396d-01
real(kind=8), parameter :: b4 = 2.38458932846290d-01
real(kind=8), parameter :: b5 = 2.87632146308408d-01
real(kind=8), parameter :: a31 = 3.55909775063327d-01
real(kind=8), parameter :: a33 = 6.44090224936673d-01
real(kind=8), parameter :: a41 = 3.67933791638137d-01
real(kind=8), parameter :: a44 = 6.32066208361863d-01
real(kind=8), parameter :: a53 = 2.37593836598569d-01
real(kind=8), parameter :: a55 = 7.62406163401431d-01
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
!= 1st step: U(1) = U(n) + b1 dt F[U(n)]
!
! calculate the fractional time step
!
ds = b1 * dt
! prepare times
!
tm = time + ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u1
! assign pdata to the next block
!
pdata => pdata%next
end do
! update primitive variables
!
call update_variables(tm, dtm)
!= 2nd step: U(2) = U(1) + b1 dt F[U(1)]
!
! prepare times
!
tm = time + 2.0d+00 * ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)]
!
! calculate the fractional time step
!
ds = b3 * dt
! prepare times
!
tm = time + (2.0d+00 * a33 * b1 + b3) * dt
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u1(:,:,:,:) = a31 * pdata%u0(:,:,:,:) + a33 * pdata%u1(:,:,:,:) &
+ ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)]
!
! calculate the fractional time step
!
ds = b4 * dt
! prepare times
!
tm = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u0(:,:,:,:) = a41 * pdata%u0(:,:,:,:) + a44 * pdata%u1(:,:,:,:) &
+ ds * du(:,:,:,:)
! update the conservative variable pointer
!
pdata%u => pdata%u0
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)]
!
! calculate the fractional time step
!
ds = b5 * dt
! prepare times
!
tm = time + dt
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = a53 * pdata%u1(:,:,:,:) + a55 * pdata%u0(:,:,:,:) &
+ ds * du(:,:,:,:)
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk35
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK3_M:
! --------------------------
!
! Subroutine advances the solution by one time step using the 3rd order
! m-stage Strong Stability Preserving Runge-Kutta time integration method.
!
! References:
!
! [1] Gottlieb, S., Ketcheson, D., Shu C. - W.,
! "Strong stability preserving Runge-Kutta and multistep
! time discretizations",
! World Scientific Publishing, 2011
!
!===============================================================================
!
subroutine evolve_ssprk3_m()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local saved variables
!
logical , save :: first = .true.
integer , save :: n, n1, n2, n3, n4
real(kind=8), save :: r, f1, f2
! local variables
!
integer :: i
real(kind=8) :: ds
real(kind=8) :: tm, dtm
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
! prepare coefficients
!
if (first) then
n = int(sqrt(1.0d+00 * stages))
n1 = (n - 1) * (n - 2) / 2
n2 = n1 + 1
n3 = n * (n + 1) / 2 - 1
n4 = n * (n + 1) / 2 + 1
r = (1.0d+00 * (2 * n - 1))
f1 = (1.0d+00 * n ) / r
f2 = (1.0d+00 * (n - 1)) / r
r = 1.0d+00 * (stages - n)
! update first flag
!
first = .false.
end if
!= 1st step: U(1) = [1 + dt/r) L] U(1), for i = 1, ..., n1
!
! calculate the fractional time step
!
ds = dt / r
! integrate the intermediate steps
!
do i = 1, n1
! prepare times
!
tm = time + i * ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
end do ! n = 1, n1
!= 2nd step: U(2) = U(1)
!
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! update the final solution
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
!= 3rd step: U(1) = [1 + dt/r) L] U(1), for i = n2, ..., n3
!
! integrate the intermediate steps
!
do i = n2, n3
! prepare times
!
tm = time + i * ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
end do ! i = n2, n3
!= 4th step: U(1) = n * U(2) + (n - 1) * [1 + dt/r L] U(1)) / (2 * n - 1)
!
! prepare times
!
tm = time + dt
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = f1 * pdata%u1(:,:,:,:) + f2 * (pdata%u0(:,:,:,:) &
+ ds * du(:,:,:,:))
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
!= the final step: U(1) = [1 + dt/r) L] U(1), for i = n4, ..., m, U(n+1) = U(1)
!
! integrate the intermediate steps
!
do i = n4, stages
! prepare times
!
tm = time + (i - n) * ds
dtm = ds
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
end do ! i = n4, stages
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk3_m
!
!===============================================================================
!
! subroutine EVOLVE_SSPRK4_10:
! ---------------------------
!
! Subroutine advances the solution by one time step using the 4rd order
! 10-stage Strong Stability Preserving Runge-Kutta time integration method.
!
! References:
!
! [1] Gottlieb, S., Ketcheson, D., Shu C. - W.,
! "Strong stability preserving Runge-Kutta and multistep
! time discretizations",
! World Scientific Publishing, 2011
!
!===============================================================================
!
subroutine evolve_ssprk4_10()
! include external variables
!
use blocks , only : block_data, list_data
use equations, only : ibp
use sources , only : update_sources
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local variables
!
integer :: n
real(kind=8) :: tm, dtm
! local vectors
!
real(kind=8), dimension(9) :: ds
! local parameters
!
real(kind=8), dimension(9), parameter :: c = &
(/ 1.0d+00, 2.0d+00, 3.0d+00, 4.0d+00, 2.0d+00, &
3.0d+00, 4.0d+00, 5.0d+00, 6.0d+00 /) / 6.0d+00
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for one step update
!
call start_timer(imu)
#endif /* PROFILE */
! prepare time step fractions
!
ds(:) = c(:) * dt
!= 1st step: U(2) = U(1)
!
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! update the final solution
!
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
!= 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5
!
! integrate the intermediate steps
!
do n = 1, 5
! prepare times
!
tm = time + ds(n)
dtm = ds(1)
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + dtm * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
end do ! n = 1, 5
!= 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1)
!
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! update the final solution
!
pdata%u1(:,:,:,:) = (pdata%u1(:,:,:,:) &
+ 9.0d+00 * pdata%u0(:,:,:,:)) / 2.5d+01
pdata%u0(:,:,:,:) = 1.5d+01 * pdata%u1(:,:,:,:) &
- 5.0d+00 * pdata%u0(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
!= 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9
!
! integrate the intermediate steps
!
do n = 6, 9
! prepare times
!
tm = time + ds(n)
dtm = ds(1)
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the intermediate solution
!
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + dtm * du(:,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
end do ! n = 6, 9
!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)]
!
! prepare times
!
tm = time + dt
dtm = dt / 1.0d+01
! update fluxes
!
call update_fluxes()
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! calculate the variable increment
!
call update_increment(pdata)
! add the source terms
!
call update_sources(pdata, tm, dtm, du(:,:,:,:))
! update the final solution
!
pdata%u0(:,:,:,:) = pdata%u1(:,:,:,:) + 6.0d-01 * pdata%u0(:,:,:,:) &
+ dtm * du(:,:,:,:)
! update ψ with its source term
!
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! update primitive variables
!
call update_variables(tm, dtm)
#ifdef PROFILE
! stop accounting time for one step update
!
call stop_timer(imu)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_ssprk4_10
!
!===============================================================================
!
! subroutine UPDATE_FLUXES:
! ------------------------
!
! Subroutine iterates over all data blocks and calculates the numerical
! fluxes for each block. After the fluxes are updated, they are corrected
! for blocks which have neighbours at higher refinement level.
!
!
!===============================================================================
!
subroutine update_fluxes()
! include external procedures
!
use boundaries , only : boundary_fluxes
use schemes , only : update_flux
! include external variables
!
use blocks , only : block_data, list_data
use coordinates , only : adx, ady, adz
! local variables are not implicit by default
!
implicit none
! local pointers
!
type(block_data), pointer :: pdata
! local vectors
!
real(kind=8), dimension(NDIMS) :: dx
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for fluxe update
!
call start_timer(imf)
#endif /* PROFILE */
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! obtain dx, dy, and dz for the current block
!
dx(1) = adx(pdata%meta%level)
dx(2) = ady(pdata%meta%level)
#if NDIMS == 3
dx(3) = adz(pdata%meta%level)
#endif /* NDIMS == 3 */
! update fluxes for the current block
!
call update_flux(dx(1:NDIMS), pdata%q(:,:,:,:), pdata%f(1:NDIMS,:,:,:,:))
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
! correct the numerical fluxes of the blocks which have neighbours at higher
! levels
!
call boundary_fluxes()
#ifdef PROFILE
! stop accounting time for flux update
!
call stop_timer(imf)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_fluxes
!
!===============================================================================
!
! subroutine UPDATE_INCREMENT:
! ---------------------------
!
! Subroutine calculates the conservative variable increment from
! directional fluxes.
!
! Arguments:
!
! pdata - the point to data block storing the directional fluxes;
!
!===============================================================================
!
subroutine update_increment(pdata)
! include external variables
!
use blocks , only : block_data
use coordinates, only : nbl, neu
use coordinates, only : adxi, adyi, adzi
use equations , only : nf, ns
use equations , only : idn, isl, isu
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
type(block_data), pointer, intent(inout) :: pdata
! local variables
!
integer :: i , j , k = 1, p
integer :: im1, jm1, km1
integer :: ip1, jp1, kp1
real(kind=8) :: dxi, dyi, dzi
real(kind=8) :: df
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for the increment update
!
call start_timer(iui)
#endif /* PROFILE */
! prepare the coordinate intervals
!
dxi = adxi(pdata%meta%level)
dyi = adyi(pdata%meta%level)
#if NDIMS == 3
dzi = adzi(pdata%meta%level)
#endif /* NDIMS == 3 */
! calculate the variable update from the directional fluxes
!
#if NDIMS == 3
do k = nbl, neu
km1 = k - 1
#endif /* NDIMS == 3 */
do j = nbl, neu
jm1 = j - 1
do i = nbl, neu
im1 = i - 1
#if NDIMS == 3
du(1:nf,i,j,k) = - dxi * (pdata%f(1,:,i,j,k) - pdata%f(1,:,im1,j,k)) &
- dyi * (pdata%f(2,:,i,j,k) - pdata%f(2,:,i,jm1,k)) &
- dzi * (pdata%f(3,:,i,j,k) - pdata%f(3,:,i,j,km1))
#else /* NDIMS == 3 */
du(1:nf,i,j,k) = - dxi * (pdata%f(1,:,i,j,k) - pdata%f(1,:,im1,j,k)) &
- dyi * (pdata%f(2,:,i,j,k) - pdata%f(2,:,i,jm1,k))
#endif /* NDIMS == 3 */
end do ! i = nbl, neu
end do ! j = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
! update passive scalars
!
if (ns > 0) then
! reset passive scalar increments
!
du(isl:isu,:,:,:) = 0.0d+00
#if NDIMS == 3
do k = nbl - 1, neu
kp1 = k + 1
#endif /* NDIMS == 3 */
do j = nbl - 1, neu
jp1 = j + 1
do i = nbl - 1, neu
ip1 = i + 1
! X-face
!
if (pdata%f(1,idn,i,j,k) >= 0.0d+00) then
do p = isl, isu
df = dxi * pdata%f(1,idn,i,j,k) * pdata%q(p,i ,j,k)
du(p,i ,j,k) = du(p,i ,j,k) - df
du(p,ip1,j,k) = du(p,ip1,j,k) + df
end do
else
do p = isl, isu
df = dxi * pdata%f(1,idn,i,j,k) * pdata%q(p,ip1,j,k)
du(p,i ,j,k) = du(p,i ,j,k) - df
du(p,ip1,j,k) = du(p,ip1,j,k) + df
end do
end if
! Y-face
!
if (pdata%f(2,idn,i,j,k) >= 0.0d+00) then
do p = isl, isu
df = dyi * pdata%f(2,idn,i,j,k) * pdata%q(p,i,j ,k)
du(p,i,j ,k) = du(p,i,j ,k) - df
du(p,i,jp1,k) = du(p,i,jp1,k) + df
end do
else
do p = isl, isu
df = dyi * pdata%f(2,idn,i,j,k) * pdata%q(p,i,jp1,k)
du(p,i,j ,k) = du(p,i,j ,k) - df
du(p,i,jp1,k) = du(p,i,jp1,k) + df
end do
end if
#if NDIMS == 3
! Z-face
!
if (pdata%f(3,idn,i,j,k) >= 0.0d+00) then
do p = isl, isu
df = dzi * pdata%f(3,idn,i,j,k) * pdata%q(p,i,j,k )
du(p,i,j,k ) = du(p,i,j,k ) - df
du(p,i,j,kp1) = du(p,i,j,kp1) + df
end do
else
do p = isl, isu
df = dzi * pdata%f(3,idn,i,j,k) * pdata%q(p,i,j,kp1)
du(p,i,j,k ) = du(p,i,j,k ) - df
du(p,i,j,kp1) = du(p,i,j,kp1) + df
end do
end if
#endif /* NDIMS == 3 */
end do ! i = nbl, neu
end do ! j = nbl, neu
#if NDIMS == 3
end do ! k = nbl, neu
#endif /* NDIMS == 3 */
end if
#ifdef PROFILE
! stop accounting time for the increment update
!
call stop_timer(iui)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_increment
!
!===============================================================================
!
! subroutine UPDATE_VARIABLES:
! ---------------------------
!
! Subroutine iterates over all data blocks and converts the conservative
! variables to their primitive representation.
!
! Arguments:
!
! tm - time at the moment of update;
! dtm - time step since the last update;
!
!===============================================================================
!
subroutine update_variables(tm, dtm)
! include external procedures
!
use blocks , only : set_neighbors_update
use boundaries , only : boundary_variables
use equations , only : update_primitive_variables
use equations , only : fix_unphysical_cells, correct_unphysical_states
use shapes , only : update_shapes
! include external variables
!
use blocks , only : block_meta, list_meta
use blocks , only : block_data, list_data
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
real(kind=8), intent(in) :: tm, dtm
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
! start accounting time for variable update
!
call start_timer(imv)
#endif /* PROFILE */
! update primitive variables in the changed blocks
!
pdata => list_data
do while (associated(pdata))
pmeta => pdata%meta
if (pmeta%update) call update_primitive_variables(pdata%u, pdata%q)
pdata => pdata%next
end do
! update boundaries
!
call boundary_variables(tm, dtm)
! correct unphysical states if detected
!
if (fix_unphysical_cells) then
! if an unphysical cell appeared in a block while updating its primitive
! variables it could be propagated to its neighbors through boundary update;
! mark all neighbors of such a block to be verified and corrected for
! unphysical cells too
!
pmeta => list_meta
do while (associated(pmeta))
if (pmeta%update) call set_neighbors_update(pmeta)
pmeta => pmeta%next
end do
! verify and correct, if necessary, unphysical cells in recently updated blocks
!
pdata => list_data
do while (associated(pdata))
pmeta => pdata%meta
if (pmeta%update) &
call correct_unphysical_states(step, pmeta%id, pdata%q, pdata%u)
pdata => pdata%next
end do
end if
! apply shapes in blocks which need it
!
pdata => list_data
do while (associated(pdata))
pmeta => pdata%meta
if (pmeta%update) call update_shapes(pdata, tm, dtm)
pdata => pdata%next
end do
#ifdef PROFILE
! stop accounting time for variable update
!
call stop_timer(imv)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine update_variables
#ifdef DEBUG
!
!===============================================================================
!
! subroutine CHECK_VARIABLES:
! --------------------------
!
! Subroutine iterates over all data blocks and converts the conservative
! variables to their primitive representation.
!
!
!===============================================================================
!
subroutine check_variables()
! include external procedures
!
use coordinates , only : nn => bcells
use equations , only : nv, pvars, cvars
use ieee_arithmetic, only : ieee_is_nan
! include external variables
!
use blocks , only : block_meta, list_meta
use blocks , only : block_data, list_data
! local variables are not implicit by default
!
implicit none
! local variables
!
integer :: i, j, k = 1, p
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
! associate the pointer with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! associate pmeta with the corresponding meta block
!
pmeta => pdata%meta
! check if there are NaNs in primitive variables
!
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nn
do i = 1, nn
do p = 1, nv
if (ieee_is_nan(pdata%u(p,i,j,k))) then
print *, 'U NaN:', cvars(p), pdata%meta%id, i, j, k
end if
if (ieee_is_nan(pdata%q(p,i,j,k))) then
print *, 'Q NaN:', pvars(p), pdata%meta%id, i, j, k
end if
end do ! p = 1, nv
end do ! i = 1, nn
end do ! j = 1, nn
#if NDIMS == 3
end do ! k = 1, nn
#endif /* NDIMS == 3 */
! assign pointer to the next block
!
pdata => pdata%next
end do
!-------------------------------------------------------------------------------
!
end subroutine check_variables
#endif /* DEBUG */
!===============================================================================
!
end module evolution