!!******************************************************************************
!!
!! module: EVOLUTION - handling the time evolution of the block structure
!!
!! Copyright (C) 2008-2011 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!!******************************************************************************
!!
!! This file is part of the AMUN code.
!!
!! 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 2
!! 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, write to the Free Software Foundation,
!! Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
!!
!!******************************************************************************
!!
!
module evolution
implicit none
integer, save :: n
real , save :: t, dt, dtn, dxmin
contains
!
!===============================================================================
!
! evolve: subroutine sweeps over all leaf blocks and performs one step time
! evolution for each according to the selected integration scheme
!
!===============================================================================
!
subroutine evolve()
use blocks , only : block_data, list_data
use boundaries, only : boundary_variables
#ifdef CONSERVATIVE
use boundaries, only : boundary_correct_fluxes
#endif /* CONSERVATIVE */
#ifdef REFINE
use config , only : maxlev
#endif /* REFINE */
use mesh , only : update_mesh
use timer , only : start_timer, stop_timer
#ifdef FORCE
use forcing , only : fourier_transform, evolve_forcing
use variables , only : idn, imz
#endif /* FORCE */
implicit none
! local variables
!
type(block_data), pointer :: pblock
real :: cm
!
!-------------------------------------------------------------------------------
!
! start the evolution timer
!
call start_timer(2)
#ifdef FORCE
! perform the Fourier transform of the velocity field
!
pblock => list_data
do while (associated(pblock))
if (pblock%meta%leaf) &
call fourier_transform(pblock%meta%level &
, pblock%meta%xmin, pblock%meta%ymin, pblock%meta%zmin &
, pblock%u(idn:imz,:,:,:))
pblock => pblock%next ! assign pointer to the next block
end do
! evolve the forcing source terms by the time interval dt
!
call evolve_forcing(dt)
#endif /* FORCE */
! iterate over all data blocks and perform one step of time evolution
!
pblock => list_data
do while (associated(pblock))
! check if this block is a leaf
!
#ifdef CONSERVATIVE
if (pblock%meta%leaf) &
#ifdef EULER
call flux_euler(pblock)
#endif /* EULER */
#ifdef RK2
call flux_rk2(pblock)
#endif /* RK2 */
#ifdef RK3
call flux_rk3(pblock)
#endif /* RK3 */
#else /* CONSERVATIVE */
if (pblock%meta%leaf) &
#ifdef EULER
call evolve_euler(pblock)
#endif /* EULER */
#ifdef RK2
call evolve_rk2(pblock)
#endif /* RK2 */
#ifdef RK3
call evolve_rk3(pblock)
#endif /* RK3 */
#endif /* CONSERVATIVE */
! assign pointer to the next block
!
pblock => pblock%next
end do
#ifdef CONSERVATIVE
! correct the numerical fluxes between neighboring blocks which are at different
! levels
!
call boundary_correct_fluxes()
! update solution using numerical fluxes stored in data blocks
!
pblock => list_data
do while (associated(pblock))
! check if this block is a leaf and update its conserved variables using
! corrected numerical fluxes
!
if (pblock%meta%leaf) call update_solution(pblock)
! assign pointer to the next block
!
pblock => pblock%next
end do
#endif /* CONSERVATIVE */
! update boundaries
!
call boundary_variables()
#ifdef REFINE
! chec if we need to perform the refinement step
!
if (maxlev .gt. 1) then
! check refinement and refine
!
call update_mesh()
! update boundaries
!
call boundary_variables()
end if ! maxlev > 1
#endif /* REFINE */
! find new time step
!
call find_new_timestep()
! stop the evolution timer
!
call stop_timer(2)
!
!-------------------------------------------------------------------------------
!
end subroutine evolve
!
!===============================================================================
!
! find_new_timestep: subroutine updates the maximum speed among the leafs and
! calculates new time step
!
!===============================================================================
!
subroutine find_new_timestep()
use blocks , only : block_meta, block_data, list_meta, list_data
use config , only : maxlev
#ifdef MPI
use mpitools, only : mallreducemaxr
#endif /* MPI */
use coords , only : adx, ady, adz
use scheme , only : maxspeed, cmax
use timer , only : start_timer, stop_timer
#ifdef VISCOSITY
use config , only : visc
#endif /* VISCOSITY */
#if defined MHD && defined RESISTIVITY
use config , only : ueta
#endif /* MHD & RESISTIVITY */
implicit none
! local variables
!
real :: cm
integer(kind=4) :: lev
! local pointers
!
type(block_meta), pointer :: pmeta
type(block_data), pointer :: pdata
!
!-------------------------------------------------------------------------------
!
! start the evolution timer
!
call start_timer(6)
! reset the maximum speed, and highest level
!
cmax = 1.0d-16
lev = 1
! if maxlev > 1, find the highest level
!
if (maxlev .gt. 1) then
! iterate over all meta blocks and find the highest level with leafs
!
pmeta => list_meta
do while (associated(pmeta))
! check if the metablock is a leaf, if so obtaind the highest level
!
if (pmeta%leaf) lev = max(lev, pmeta%level)
! associate the pointer with the next block
!
pmeta => pmeta%next
end do ! meta blocks
! find the smallest spacial step
!
#if NDIMS == 2
dxmin = min(adx(lev), ady(lev))
#endif /* NDIMS == 2 */
#if NDIMS == 3
dxmin = min(adx(lev), ady(lev), adz(lev))
#endif /* NDIMS == 3 */
end if ! maxlev > 1
! iterate over all data blocks in order to find the maximum speed among them
!
pdata => list_data
do while (associated(pdata))
! check if this block is a leaf
!
if (pdata%meta%leaf) then
! find the maximum level occupied by blocks (can be smaller than maxlev)
!
! obtain the maximum speed for the current block
!
cm = maxspeed(pdata%u(:,:,:,:))
! compare global and local maximum speeds
!
cmax = max(cmax, cm)
end if ! leaf
! assiociate the pointer with the next block
!
pdata => pdata%next
end do
#ifdef MPI
! reduce the maximum speed over all processes
!
call mallreducemaxr(cmax)
#endif /* MPI */
! calcilate new time step
!
dtn = dxmin / max(cmax, 1.0d-16)
#ifdef VISCOSITY
dtn = min(dtn, 0.5d0 * dxmin * dxmin / max(1.0d-16, visc))
#endif /* VISCOSITY */
#if defined MHD && defined RESISTIVITY
dtn = min(dtn, 0.5d0 * dxmin * dxmin / max(1.0d-16, ueta))
#endif /* MHD & RESISTIVITY */
! stop the evolution timer
!
call stop_timer(6)
!
!-------------------------------------------------------------------------------
!
end subroutine find_new_timestep
#ifdef CONSERVATIVE
!
!===============================================================================
!
! update_solution: subroutine performs an one step update of the conserved
! variables for the given data block using the integrated
! numerical fluxes stored in the same data block
!
!===============================================================================
!
subroutine update_solution(pblock)
use blocks , only : block_data
use config , only : im, jm, km
use coords , only : adxi, adyi, adzi
#if defined MHD && defined GLM
use config , only : alpha_p
use scheme , only : cmax
use variables, only : iph
#endif /* MHD & GLM */
#ifdef FORCE
use forcing , only : real_forcing
use variables, only : inx, iny, inz
use variables, only : idn, imx, imy, imz
#ifdef ADI
use variables, only : ien
#endif /* ADI */
#endif /* FORCE */
#ifdef SHAPE
use problem , only : update_shapes
#endif /* SHAPE */
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
#if defined MHD && defined GLM
real :: decay
#endif /* MHD & GLM */
#ifdef FORCE
! local arrays
!
real, dimension( 3,im,jm,km) :: f
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
#endif /* FORCE */
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
#if NDIMS == 3
dzi = adzi(pblock%meta%level)
#endif /* NDIMS == 3 */
! perform update of conserved variables of the given block using its fluxes
!
call advance_solution(pblock%u(:,:,:,:), pblock%f(:,:,:,:,:), dxi, dyi, dzi)
#if defined MHD && defined GLM
! evolve Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dxmin)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* MHD & GLM */
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) &
+ pblock%u(idn,:,:,:) * f(inx,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) &
+ pblock%u(idn,:,:,:) * f(iny,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) &
+ pblock%u(idn,:,:,:) * f(inz,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
#endif /* ADI */
#endif /* FORCE */
#ifdef SHAPE
! update solid shapes
!
call update_shapes(pblock, t)
#endif /* SHAPE */
!-------------------------------------------------------------------------------
!
end subroutine update_solution
!
!===============================================================================
!
! advance_solution: subroutine performs an one step update of the conserved
! variables array using the numerical fluxes passed as an
! argument
!
!===============================================================================
!
subroutine advance_solution(u, f, dxi, dyi, dzi)
use config , only : im, jm, km
use variables, only : nqt, nfl
use variables, only : inx, iny, inz
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
use scheme , only : cmax
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
implicit none
! input arguments
!
real, dimension( nqt,im,jm,km), intent(inout) :: u
real, dimension(NDIMS,nqt,im,jm,km), intent(in) :: f
real :: dxi, dyi, dzi
! local variables
!
integer :: i, j, k, im1, jm1, km1
real :: dhx, dhy, dhz
#if defined MHD && defined GLM
real :: ch2
#endif /* MHD & GLM */
! local arrays
!
real, dimension(nqt,im,jm,km) :: du
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dhx = dt * dxi
dhy = dt * dyi
#if NDIMS == 3
dhz = dt * dzi
#endif /* NDIMS == 3 */
! reset the increment array du
!
du(:,:,:,:) = 0.0d0
! perform update along the X direction
!
do i = 1, im
im1 = max(1, i - 1)
du(:,i,:,:) = du(:,i,:,:) + dhx * (f(inx,:,im1,:,:) - f(inx,:,i,:,:))
end do
! perform update along the Y direction
!
do j = 1, jm
jm1 = max(1, j - 1)
du(:,:,j,:) = du(:,:,j,:) + dhy * (f(iny,:,:,jm1,:) - f(iny,:,:,j,:))
end do
#if NDIMS == 3
! perform update along the Z direction
!
do k = 1, km
km1 = max(1, k - 1)
du(:,:,:,k) = du(:,:,:,k) + dhz * (f(inz,:,:,:,km1) - f(inz,:,:,:,k))
end do
#endif /* NDIMS == 3 */
! update the solution for the fluid variables
!
u( 1:nfl,:,:,:) = u( 1:nfl,:,:,:) + du( 1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
u(ibx:ibz,:,:,:) = u(ibx:ibz,:,:,:) + du(ibx:ibz,:,:,:)
#ifdef GLM
! calculate c_h^2
!
ch2 = cmax * cmax
! update the solution for the scalar potential Psi
!
u(iph,:,:,:) = u(iph,:,:,:) + ch2 * du(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
!-------------------------------------------------------------------------------
!
end subroutine advance_solution
#ifdef EULER
!
!===============================================================================
!
! flux_euler: subroutine performs the first order integration of the numerical
! flux
!
!===============================================================================
!
subroutine flux_euler(pblock)
use blocks , only : block_data
use coords , only : adxi, adyi, adzi
use scheme , only : update_flux
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! 1st step of integration
!
call update_flux(pblock%u(:,:,:,:), pblock%f(:,:,:,:,:), dxi, dyi, dzi)
!-------------------------------------------------------------------------------
!
end subroutine flux_euler
#endif /* EULER */
#ifdef RK2
!
!===============================================================================
!
! flux_rk2: subroutine performs integration of the numerical flux using
! the second order Runge-Kutta method
!
!===============================================================================
!
subroutine flux_rk2(pblock)
use blocks , only : block_data
use config , only : im, jm, km
use coords , only : adxi, adyi, adzi
use scheme , only : update_flux
use variables, only : nqt
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
! local arrays
!
real, dimension( nqt,im,jm,km) :: u
real, dimension(NDIMS,nqt,im,jm,km) :: f0, f1
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! copy the initial state to the local array u
!
u(:,:,:,:) = pblock%u(:,:,:,:)
! calculate fluxes at the moment t
!
call update_flux(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! advance the solution to (t + dt) using computed fluxes in this substep
!
call advance_solution(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! calculate fluxes at the moment (t + dt)
!
call update_flux(u(:,:,:,:), f1(:,:,:,:,:), dxi, dyi, dzi)
! calculate the time averaged flux
!
pblock%f(:,:,:,:,:) = 0.5d0 * (f0(:,:,:,:,:) + f1(:,:,:,:,:))
!-------------------------------------------------------------------------------
!
end subroutine flux_rk2
#endif /* RK2 */
#ifdef RK3
!
!===============================================================================
!
! flux_rk3: subroutine performs integration of the numerical flux using
! the third order Runge-Kutta method
!
!===============================================================================
!
subroutine flux_rk3(pblock)
use blocks , only : block_data
use config , only : im, jm, km
use coords , only : adxi, adyi, adzi
use scheme , only : update_flux
use variables, only : nqt
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
! local arrays
!
real, dimension( nqt,im,jm,km) :: u
real, dimension(NDIMS,nqt,im,jm,km) :: f0, f1, f2
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! copy the initial state to the local array u
!
u(:,:,:,:) = pblock%u(:,:,:,:)
! calculate fluxes at the moment t
!
call update_flux(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! advance the solution to (t + dt) using computed fluxes in this substep
!
call advance_solution(u(:,:,:,:), f0(:,:,:,:,:), dxi, dyi, dzi)
! calculate fluxes at the moment (t + dt)
!
call update_flux(u(:,:,:,:), f1(:,:,:,:,:), dxi, dyi, dzi)
! copy the initial state to the local array u
!
u(:,:,:,:) = pblock%u(:,:,:,:)
! average fluxes from t and t + dt and prepare for half step update
!
f2(:,:,:,:,:) = 0.25d0 * (f0(:,:,:,:,:) + f1(:,:,:,:,:))
! advance the solution to (t + dt/2) using computed flux
!
call advance_solution(u(:,:,:,:), f2(:,:,:,:,:), dxi, dyi, dzi)
! calculate fluxes at the moment (t + dt/2)
!
call update_flux(u(:,:,:,:), f2(:,:,:,:,:), dxi, dyi, dzi)
! calculate the time averaged flux using Gauss formula
!
pblock%f(:,:,:,:,:) = (f0(:,:,:,:,:) + f1(:,:,:,:,:) &
+ 4.0d0 * f2(:,:,:,:,:)) / 6.0d0
!-------------------------------------------------------------------------------
!
end subroutine flux_rk3
#endif /* RK3 */
#else /* CONSERVATIVE */
#ifdef EULER
!
!===============================================================================
!
! evolve_euler: subroutine evolves the current block using Euler integration
!
!===============================================================================
!
subroutine evolve_euler(pblock)
use blocks , only : block_data
use config , only : im, jm, km
#ifdef FORCE
use forcing , only : real_forcing
#endif /* FORCE */
use coords , only : adxi, adyi, adzi
#ifdef SHAPE
use problem , only : update_shapes
#endif /* SHAPE */
use scheme , only : update, cmax
use variables, only : nqt, nfl
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
use config , only : alpha_p
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
use variables, only : idn, imx, imy, imz
#ifdef ADI
use variables, only : ien
#endif /* ADI */
#endif /* FORCE */
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi, ch2
#if defined MHD && defined GLM
real :: decay
#endif /* MHD & GLM */
! local arrays
!
real, dimension(nqt,im,jm,km) :: du
#ifdef FORCE
real, dimension( 3,im,jm,km) :: f
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
#endif /* FORCE */
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
! 1st step of integration
!
call update(pblock%u(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
! update the solution for the fluid variables
!
pblock%u(1:nfl,:,:,:) = pblock%u(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
pblock%u(ibx:ibz,:,:,:) = pblock%u(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)
#ifdef GLM
! calculate c_h^2
!
ch2 = cmax * cmax
! update the solution for the scalar potential Psi
!
pblock%u(iph,:,:,:) = pblock%u(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)
! evolve Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dxmin)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) + pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) + pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) + pblock%u(idn,:,:,:) * f(3,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
#endif /* ADI */
#endif /* FORCE */
#ifdef SHAPE
! restrict update in a defined shape
!
call update_shapes(pblock, t)
#endif /* SHAPE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_euler
#endif /* EULER */
#ifdef RK2
!
!===============================================================================
!
! evolve_rk2: subroutine evolves the current block using the 2nd order
! Runge-Kutta method
!
!===============================================================================
!
subroutine evolve_rk2(pblock)
use blocks , only : block_data
use config , only : im, jm, km
#ifdef FORCE
use forcing , only : real_forcing
#endif /* FORCE */
use coords , only : adxi, adyi, adzi
#ifdef SHAPE
use problem , only : update_shapes
#endif /* SHAPE */
use scheme , only : update, cmax
use variables, only : nqt, nfl
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
use config , only : alpha_p
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
use variables, only : idn, imx, imy, imz
#ifdef ADI
use variables, only : ien
#endif /* ADI */
#endif /* FORCE */
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi, ch2
#if defined MHD && defined GLM
real :: decay
#endif /* MHD & GLM */
! local arrays
!
real, dimension(nqt,im,jm,km) :: u1, du
#ifdef FORCE
real, dimension( 3,im,jm,km) :: f
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
#endif /* FORCE */
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
!! 1st step of integration
!!
call update(pblock%u(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
! update the solution for the fluid variables
!
u1(1:nfl,:,:,:) = pblock%u(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
u1(ibx:ibz,:,:,:) = pblock%u(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)
#ifdef GLM
! calculate c_h^2
!
ch2 = cmax * cmax
! update the solution for the scalar potential Psi
!
u1(iph,:,:,:) = pblock%u(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
! 2nd step of integration
!
call update(u1(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
! update the solution for the fluid variables
!
pblock%u(1:nfl,:,:,:) = 0.5d0 * (pblock%u(1:nfl,:,:,:) &
+ u1(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:))
#ifdef MHD
! update the solution for the magnetic variables
!
pblock%u(ibx:ibz,:,:,:) = 0.5d0 * (pblock%u(ibx:ibz,:,:,:) &
+ u1(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:))
#ifdef GLM
! update the solution for the scalar potential Psi
!
pblock%u(iph,:,:,:) = 0.5d0 * (pblock%u(iph,:,:,:) &
+ u1(iph,:,:,:) + ch2 * dt * du(iph,:,:,:))
! evolve Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dxmin)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) + pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) + pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) + pblock%u(idn,:,:,:) * f(3,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
#endif /* ADI */
#endif /* FORCE */
#ifdef SHAPE
! restrict update in a defined shape
!
call update_shapes(pblock, t)
#endif /* SHAPE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_rk2
#endif /* RK2 */
#ifdef RK3
!
!===============================================================================
!
! evolve_rk3: subroutine evolves the current block using the 3rd order
! Runge-Kutta method
!
!===============================================================================
!
subroutine evolve_rk3(pblock)
use blocks , only : block_data
use config , only : im, jm, km
#ifdef FORCE
use forcing , only : real_forcing
#endif /* FORCE */
use coords , only : adxi, adyi, adzi
#ifdef SHAPE
use problem , only : update_shapes
#endif /* SHAPE */
use scheme , only : update, cmax
use variables, only : nqt, nfl
#ifdef MHD
use variables, only : ibx, ibz
#ifdef GLM
use config , only : alpha_p
use variables, only : iph
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
use variables, only : idn, imx, imy, imz
#ifdef ADI
use variables, only : ien
#endif /* ADI */
#endif /* FORCE */
implicit none
! input arguments
!
type(block_data), pointer, intent(inout) :: pblock
! local variables
!
real :: dxi, dyi, dzi
#if defined MHD && defined GLM
real :: decay, ch2
#endif /* MHD & GLM */
! local arrays
!
real, dimension(nqt,im,jm,km) :: u1, du
#ifdef FORCE
real, dimension( 3,im,jm,km) :: f
#ifdef ADI
real, dimension( im,jm,km) :: ek, dek
#endif /* ADI */
#endif /* FORCE */
! parameters
!
real, parameter :: f4 = 1.0d0 / 4.0d0, f3 = 1.0d0 / 3.0d0
!
!-------------------------------------------------------------------------------
!
! prepare dxi, dyi, and dzi
!
dxi = adxi(pblock%meta%level)
dyi = adyi(pblock%meta%level)
dzi = adzi(pblock%meta%level)
!! 1st step of integration
!!
call update(pblock%u(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
! update the solution for the fluid variables
!
u1(1:nfl,:,:,:) = pblock%u(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)
#ifdef MHD
! update the solution for the magnetic variables
!
u1(ibx:ibz,:,:,:) = pblock%u(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)
#ifdef GLM
! calculate c_h^2
!
ch2 = cmax * cmax
! update the solution for the scalar potential Psi
!
u1(iph,:,:,:) = pblock%u(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
!! 2nd step of integration
!!
call update(u1(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
! update the solution for the fluid variables
!
u1(1:nfl,:,:,:) = f4 * (3.0d0 * pblock%u(1:nfl,:,:,:) &
+ u1(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:))
#ifdef MHD
! update the solution for the magnetic variables
!
u1(ibx:ibz,:,:,:) = f4 * (3.0d0 * pblock%u(ibx:ibz,:,:,:) &
+ u1(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:))
#ifdef GLM
! update the solution for the scalar potential Psi
!
u1(iph,:,:,:) = f4 * (3.0d0 * pblock%u(iph,:,:,:) &
+ u1(iph,:,:,:) + ch2 * dt * du(iph,:,:,:))
#endif /* GLM */
#endif /* MHD */
!! 3rd step of integration
!!
call update(u1(:,:,:,:), du(:,:,:,:), dxi, dyi, dzi)
! update the solution for the fluid variables
!
pblock%u(1:nfl,:,:,:) = f3 * (pblock%u(1:nfl,:,:,:) &
+ 2.0d0 * (u1(1:nfl,:,:,:) + dt * du(1:nfl,:,:,:)))
#ifdef MHD
! update the solution for the magnetic variables
!
pblock%u(ibx:ibz,:,:,:) = f3 * (pblock%u(ibx:ibz,:,:,:) &
+ 2.0d0 * (u1(ibx:ibz,:,:,:) + dt * du(ibx:ibz,:,:,:)))
#ifdef GLM
! update the solution for the scalar potential Psi
!
pblock%u(iph,:,:,:) = f3 * (pblock%u(iph,:,:,:) &
+ 2.0d0 * (u1(iph,:,:,:) + ch2 * dt * du(iph,:,:,:)))
! evolve analytically Psi due to the source term
!
decay = exp(- alpha_p * cmax * dt / dxmin)
pblock%u(iph,:,:,:) = decay * pblock%u(iph,:,:,:)
#endif /* GLM */
#endif /* MHD */
#ifdef FORCE
! obtain the forcing terms in real space
!
call real_forcing(pblock%meta%level, pblock%meta%xmin, pblock%meta%ymin &
, pblock%meta%zmin, f(:,:,:,:))
#ifdef ADI
! calculate kinetic energy before adding the forcing term
!
ek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:)
#endif /* ADI */
! update momenta due to the forcing terms
!
pblock%u(imx,:,:,:) = pblock%u(imx,:,:,:) + pblock%u(idn,:,:,:) * f(1,:,:,:)
pblock%u(imy,:,:,:) = pblock%u(imy,:,:,:) + pblock%u(idn,:,:,:) * f(2,:,:,:)
pblock%u(imz,:,:,:) = pblock%u(imz,:,:,:) + pblock%u(idn,:,:,:) * f(3,:,:,:)
#ifdef ADI
! calculate kinetic energy after adding the forcing term
!
dek(:,:,:) = 0.5d0 * (pblock%u(imx,:,:,:)**2 + pblock%u(imy,:,:,:)**2 &
+ pblock%u(imz,:,:,:)**2) / pblock%u(idn,:,:,:) - ek(:,:,:)
! update total energy with the injected one
!
pblock%u(ien,:,:,:) = pblock%u(ien,:,:,:) + dek(:,:,:)
#endif /* ADI */
#endif /* FORCE */
#ifdef SHAPE
! restrict update in a defined shape
!
call update_shapes(pblock, t)
#endif /* SHAPE */
!-------------------------------------------------------------------------------
!
end subroutine evolve_rk3
#endif /* RK3 */
#endif /* CONSERVATIVE */
!
!===============================================================================
!
end module