Rewrite module EVOLUTION. Move update_interval() to SCHEME.
Improve subroutine comments and remove unused code.
This commit is contained in:
parent
2413de249f
commit
bbed8c41e5
@ -23,7 +23,10 @@
|
|||||||
!!
|
!!
|
||||||
!! module: EVOLUTION
|
!! module: EVOLUTION
|
||||||
!!
|
!!
|
||||||
!! This module performs the time integration using different methods.
|
!! This module constains subroutines to perform the time integration using
|
||||||
|
!! different methods preserving the conservation and estimates the new time
|
||||||
|
!! step from the stability condition.
|
||||||
|
!!
|
||||||
!!
|
!!
|
||||||
!!******************************************************************************
|
!!******************************************************************************
|
||||||
!
|
!
|
||||||
@ -35,17 +38,14 @@ module evolution
|
|||||||
|
|
||||||
! evolution parameters
|
! evolution parameters
|
||||||
!
|
!
|
||||||
real , save :: cfl = 0.5d0
|
real , save :: cfl = 0.5d+0
|
||||||
|
|
||||||
#if defined MHD && defined GLM
|
! time variables
|
||||||
! coefficient controlling the decay of scalar potential Psi
|
|
||||||
!
|
!
|
||||||
real , save :: alpha_p = 0.5d0
|
integer, save :: n = 0
|
||||||
real , save :: decay = 1.0d0
|
real , save :: t = 0.0d+0
|
||||||
#endif /* MHD & GLM */
|
real , save :: dt = 1.0d+0
|
||||||
|
real , save :: dtn = 1.0d+0
|
||||||
integer, save :: n
|
|
||||||
real , save :: t, dt, dtn, dxmin
|
|
||||||
|
|
||||||
! by default everything is private
|
! by default everything is private
|
||||||
!
|
!
|
||||||
@ -57,7 +57,7 @@ module evolution
|
|||||||
|
|
||||||
! declare public variables
|
! declare public variables
|
||||||
!
|
!
|
||||||
public :: n, cfl, t, dt, dtn, dxmin
|
public :: cfl, n, t, dt, dtn
|
||||||
|
|
||||||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||||||
!
|
!
|
||||||
@ -72,16 +72,16 @@ module evolution
|
|||||||
! subroutine INITIALIZE_EVOLUTION:
|
! subroutine INITIALIZE_EVOLUTION:
|
||||||
! -------------------------------
|
! -------------------------------
|
||||||
!
|
!
|
||||||
! Subroutine initializes the EVOLUTION module by setting its parameters.
|
! Subroutine initializes module EVOLUTION by setting its parameters.
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
subroutine initialize_evolution()
|
subroutine initialize_evolution()
|
||||||
|
|
||||||
! include external procedures and variables
|
! include external procedures
|
||||||
!
|
!
|
||||||
use parameters, only : get_parameter_real
|
use parameters , only : get_parameter_real
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
@ -93,19 +93,9 @@ module evolution
|
|||||||
!
|
!
|
||||||
call get_parameter_real("cfl", cfl)
|
call get_parameter_real("cfl", cfl)
|
||||||
|
|
||||||
#if defined MHD && defined GLM
|
|
||||||
! get the value of the stability coefficient
|
|
||||||
!
|
|
||||||
call get_parameter_real("alpha_p", alpha_p)
|
|
||||||
|
|
||||||
! calculate decay factor for magnetic field divergence scalar
|
|
||||||
!
|
|
||||||
decay = exp(- alpha_p * cfl)
|
|
||||||
#endif /* MHD & GLM */
|
|
||||||
|
|
||||||
! calculate the initial time step
|
! calculate the initial time step
|
||||||
!
|
!
|
||||||
call find_new_timestep()
|
call new_time_step()
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
@ -116,22 +106,23 @@ module evolution
|
|||||||
! subroutine ADVANCE:
|
! subroutine ADVANCE:
|
||||||
! ------------------
|
! ------------------
|
||||||
!
|
!
|
||||||
! Subroutine advances the solution by one time step using the selected
|
! Subroutine advances the solution by one time step using the selected time
|
||||||
! time integration method.
|
! integration method.
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
subroutine advance()
|
subroutine advance()
|
||||||
|
|
||||||
! include external procedures and variables
|
! include external procedures
|
||||||
|
!
|
||||||
|
use boundaries , only : boundary_variables
|
||||||
|
use mesh , only : update_mesh
|
||||||
|
|
||||||
|
! include external variables
|
||||||
!
|
!
|
||||||
use blocks , only : block_data, list_data
|
use blocks , only : block_data, list_data
|
||||||
use boundaries , only : boundary_variables
|
use coordinates , only : toplev
|
||||||
#ifdef REFINE
|
|
||||||
use coordinates, only : toplev
|
|
||||||
#endif /* REFINE */
|
|
||||||
use mesh , only : update_mesh
|
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
@ -148,7 +139,6 @@ module evolution
|
|||||||
call advance_rk2()
|
call advance_rk2()
|
||||||
#endif /* RK2 */
|
#endif /* RK2 */
|
||||||
|
|
||||||
#ifdef REFINE
|
|
||||||
! chec if we need to perform the refinement step
|
! chec if we need to perform the refinement step
|
||||||
!
|
!
|
||||||
if (toplev > 1) then
|
if (toplev > 1) then
|
||||||
@ -162,7 +152,6 @@ module evolution
|
|||||||
call boundary_variables()
|
call boundary_variables()
|
||||||
|
|
||||||
end if ! toplev > 1
|
end if ! toplev > 1
|
||||||
#endif /* REFINE */
|
|
||||||
|
|
||||||
! update primitive variables
|
! update primitive variables
|
||||||
!
|
!
|
||||||
@ -170,7 +159,7 @@ module evolution
|
|||||||
|
|
||||||
! find new time step
|
! find new time step
|
||||||
!
|
!
|
||||||
call find_new_timestep()
|
call new_time_step()
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
@ -193,12 +182,16 @@ module evolution
|
|||||||
!
|
!
|
||||||
subroutine advance_rk2()
|
subroutine advance_rk2()
|
||||||
|
|
||||||
! include external procedures and variables
|
! include external procedures
|
||||||
|
!
|
||||||
|
use boundaries , only : boundary_variables
|
||||||
|
use scheme , only : update_interval
|
||||||
|
|
||||||
|
! include external variables
|
||||||
!
|
!
|
||||||
use blocks , only : block_data, list_data
|
use blocks , only : block_data, list_data
|
||||||
use boundaries , only : boundary_variables
|
use coordinates , only : adx, ady, adz
|
||||||
use coordinates, only : im, jm, km
|
use coordinates , only : im, jm, km
|
||||||
use coordinates, only : adx, ady, adz
|
|
||||||
use variables , only : nfl
|
use variables , only : nfl
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
@ -313,12 +306,15 @@ module evolution
|
|||||||
!
|
!
|
||||||
subroutine update_fluxes()
|
subroutine update_fluxes()
|
||||||
|
|
||||||
! include external procedures and variables
|
! include external procedures
|
||||||
|
!
|
||||||
|
use boundaries , only : boundary_correct_fluxes
|
||||||
|
use scheme , only : update_flux
|
||||||
|
|
||||||
|
! include external variables
|
||||||
!
|
!
|
||||||
use blocks , only : block_data, list_data
|
use blocks , only : block_data, list_data
|
||||||
use boundaries , only : boundary_correct_fluxes
|
use coordinates , only : adx, ady, adz
|
||||||
use coordinates, only : adx, ady, adz
|
|
||||||
use scheme , only : update_flux
|
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
@ -383,10 +379,13 @@ module evolution
|
|||||||
!
|
!
|
||||||
subroutine update_variables()
|
subroutine update_variables()
|
||||||
|
|
||||||
! include external procedures and variables
|
! include external procedures
|
||||||
|
!
|
||||||
|
use equations , only : update_primitive_variables
|
||||||
|
|
||||||
|
! include external variables
|
||||||
!
|
!
|
||||||
use blocks , only : block_data, list_data
|
use blocks , only : block_data, list_data
|
||||||
use equations , only : update_primitive_variables
|
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
@ -419,162 +418,77 @@ module evolution
|
|||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
! subroutine UPDATE_INTERVAL:
|
! subroutine NEW_TIME_STEP:
|
||||||
! --------------------------
|
! ------------------------
|
||||||
!
|
!
|
||||||
! Subroutine calculate the conservative variable interval from fluxes.
|
! Subroutine estimates the new time step from the maximum speed in the system
|
||||||
|
! and source term constraints.
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
subroutine update_interval(dh, f, du)
|
subroutine new_time_step()
|
||||||
|
|
||||||
! include external procedures and variables
|
! include external procedures
|
||||||
!
|
!
|
||||||
use coordinates, only : im, jm, km
|
use equations , only : maxspeed
|
||||||
use variables , only : nfl
|
#ifdef MPI
|
||||||
|
use mpitools , only : reduce_maximum_real
|
||||||
|
#endif /* MPI */
|
||||||
|
|
||||||
|
! include external variables
|
||||||
|
!
|
||||||
|
use blocks , only : block_data, list_data
|
||||||
|
use coordinates , only : adx, ady, adz
|
||||||
|
use coordinates , only : toplev
|
||||||
|
use variables , only : cmax
|
||||||
|
|
||||||
! local variables are not implicit by default
|
! local variables are not implicit by default
|
||||||
!
|
!
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
! subroutine arguments
|
! local pointers
|
||||||
!
|
!
|
||||||
real, dimension(3) , intent(in) :: dh
|
type(block_data), pointer :: pblock
|
||||||
real, dimension(NDIMS,nfl,im,jm,km), intent(in) :: f
|
|
||||||
real, dimension( nfl,im,jm,km), intent(inout) :: du
|
|
||||||
|
|
||||||
! local variables
|
|
||||||
!
|
|
||||||
integer :: i, j, k
|
|
||||||
!
|
|
||||||
!-------------------------------------------------------------------------------
|
|
||||||
!
|
|
||||||
! reset the increment array du
|
|
||||||
!
|
|
||||||
du(:,:,:,:) = 0.0d0
|
|
||||||
|
|
||||||
! perform update along the X direction
|
|
||||||
!
|
|
||||||
do i = 2, im
|
|
||||||
du(:,i,:,:) = du(:,i,:,:) - dh(1) * (f(1,:,i,:,:) - f(1,:,i-1,:,:))
|
|
||||||
end do
|
|
||||||
du(:,1,:,:) = du(:,1,:,:) - dh(1) * f(1,:,1,:,:)
|
|
||||||
|
|
||||||
! perform update along the Y direction
|
|
||||||
!
|
|
||||||
do j = 2, jm
|
|
||||||
du(:,:,j,:) = du(:,:,j,:) - dh(2) * (f(2,:,:,j,:) - f(2,:,:,j-1,:))
|
|
||||||
end do
|
|
||||||
du(:,:,1,:) = du(:,:,1,:) - dh(2) * f(2,:,:,1,:)
|
|
||||||
|
|
||||||
#if NDIMS == 3
|
|
||||||
! perform update along the Z direction
|
|
||||||
!
|
|
||||||
do k = 2, km
|
|
||||||
du(:,:,:,k) = du(:,:,:,k) - dh(3) * (f(3,:,:,:,k) - f(3,:,:,:,k-1))
|
|
||||||
end do
|
|
||||||
du(:,:,:,1) = du(:,:,:,1) - dh(3) * f(3,:,:,:,1)
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
|
||||||
!
|
|
||||||
end subroutine update_interval
|
|
||||||
!
|
|
||||||
!===============================================================================
|
|
||||||
!
|
|
||||||
! 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 coordinates, only : toplev
|
|
||||||
use coordinates, only : adx, ady, adz
|
|
||||||
use equations , only : maxspeed
|
|
||||||
#ifdef MPI
|
|
||||||
use mpitools , only : reduce_maximum_real
|
|
||||||
#endif /* MPI */
|
|
||||||
use variables , only : cmax
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! local variables
|
! local variables
|
||||||
!
|
!
|
||||||
integer :: iret
|
integer :: iret
|
||||||
real :: cm
|
|
||||||
integer(kind=4) :: lev
|
integer(kind=4) :: lev
|
||||||
|
real :: cm, dx_min
|
||||||
|
|
||||||
! local pointers
|
! local parameters
|
||||||
!
|
!
|
||||||
type(block_meta), pointer :: pmeta
|
real, parameter :: eps = tiny(cmax)
|
||||||
type(block_data), pointer :: pdata
|
|
||||||
!
|
!
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
! reset the maximum speed, and highest level
|
! reset the maximum speed, and the highest level
|
||||||
!
|
!
|
||||||
cmax = 1.0d-16
|
cmax = eps
|
||||||
lev = 1
|
lev = 1
|
||||||
|
|
||||||
! if toplev > 1, find the highest level
|
|
||||||
!
|
|
||||||
if (toplev .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
|
|
||||||
|
|
||||||
end if ! toplev > 1
|
|
||||||
|
|
||||||
! 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 */
|
|
||||||
|
|
||||||
! iterate over all data blocks in order to find the maximum speed among them
|
! iterate over all data blocks in order to find the maximum speed among them
|
||||||
|
! and the highest level which is required to obtain the spatial step
|
||||||
!
|
!
|
||||||
pdata => list_data
|
pblock => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pblock))
|
||||||
|
|
||||||
! check if this block is a leaf
|
|
||||||
!
|
|
||||||
if (pdata%meta%leaf) then
|
|
||||||
|
|
||||||
! find the maximum level occupied by blocks (can be smaller than toplev)
|
! find the maximum level occupied by blocks (can be smaller than toplev)
|
||||||
!
|
!
|
||||||
|
lev = max(lev, pblock%meta%level)
|
||||||
|
|
||||||
! obtain the maximum speed for the current block
|
! obtain the maximum speed for the current block
|
||||||
!
|
!
|
||||||
cm = maxspeed(pdata%q(:,:,:,:))
|
cm = maxspeed(pblock%q(:,:,:,:))
|
||||||
|
|
||||||
! compare global and local maximum speeds
|
! compare global and local maximum speeds
|
||||||
!
|
!
|
||||||
cmax = max(cmax, cm)
|
cmax = max(cmax, cm)
|
||||||
|
|
||||||
end if ! leaf
|
|
||||||
|
|
||||||
! assiociate the pointer with the next block
|
! assiociate the pointer with the next block
|
||||||
!
|
!
|
||||||
pdata => pdata%next
|
pblock => pblock%next
|
||||||
|
|
||||||
end do
|
end do
|
||||||
|
|
||||||
@ -584,13 +498,22 @@ module evolution
|
|||||||
call reduce_maximum_real(cmax, iret)
|
call reduce_maximum_real(cmax, iret)
|
||||||
#endif /* MPI */
|
#endif /* MPI */
|
||||||
|
|
||||||
! calcilate new time step
|
! find the smallest spatial step
|
||||||
!
|
!
|
||||||
dtn = dxmin / max(cmax, 1.0d-16)
|
#if NDIMS == 2
|
||||||
|
dx_min = min(adx(lev), ady(lev))
|
||||||
|
#endif /* NDIMS == 2 */
|
||||||
|
#if NDIMS == 3
|
||||||
|
dx_min = min(adx(lev), ady(lev), adz(lev))
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
! calcilate the new time step
|
||||||
|
!
|
||||||
|
dtn = dx_min / max(cmax, eps)
|
||||||
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
end subroutine find_new_timestep
|
end subroutine new_time_step
|
||||||
|
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
!
|
!
|
||||||
|
@ -213,8 +213,8 @@ driver.o : driver.F90 blocks.o coordinates.o equations.o evolution.o \
|
|||||||
equations.o : equations.F90 coordinates.o parameters.o variables.o
|
equations.o : equations.F90 coordinates.o parameters.o variables.o
|
||||||
error.o : error.F90
|
error.o : error.F90
|
||||||
evolution.o : evolution.F90 blocks.o boundaries.o coordinates.o \
|
evolution.o : evolution.F90 blocks.o boundaries.o coordinates.o \
|
||||||
equations.o forcing.o interpolations.o mesh.o mpitools.o \
|
equations.o mesh.o mpitools.o parameters.o scheme.o \
|
||||||
problems.o scheme.o variables.o
|
variables.o
|
||||||
domains.o : domains.F90 blocks.o boundaries.o coordinates.o parameters.o
|
domains.o : domains.F90 blocks.o boundaries.o coordinates.o parameters.o
|
||||||
forcing.o : forcing.F90 constants.o coordinates.o mpitools.o \
|
forcing.o : forcing.F90 constants.o coordinates.o mpitools.o \
|
||||||
parameters.o random.o variables.o
|
parameters.o random.o variables.o
|
||||||
|
@ -310,6 +310,70 @@ module scheme
|
|||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
end subroutine update_flux
|
end subroutine update_flux
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
! subroutine UPDATE_INTERVAL:
|
||||||
|
! --------------------------
|
||||||
|
!
|
||||||
|
! Subroutine calculate the conservative variable interval from fluxes.
|
||||||
|
!
|
||||||
|
!
|
||||||
|
!===============================================================================
|
||||||
|
!
|
||||||
|
subroutine update_interval(dh, f, du)
|
||||||
|
|
||||||
|
! include external procedures and variables
|
||||||
|
!
|
||||||
|
use coordinates , only : im, jm, km
|
||||||
|
use variables , only : nfl
|
||||||
|
|
||||||
|
! local variables are not implicit by default
|
||||||
|
!
|
||||||
|
implicit none
|
||||||
|
|
||||||
|
! subroutine arguments
|
||||||
|
!
|
||||||
|
real, dimension(3) , intent(in) :: dh
|
||||||
|
real, dimension(NDIMS,nfl,im,jm,km), intent(in) :: f
|
||||||
|
real, dimension( nfl,im,jm,km), intent(inout) :: du
|
||||||
|
|
||||||
|
! local variables
|
||||||
|
!
|
||||||
|
integer :: i, j, k
|
||||||
|
!
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
! reset the increment array du
|
||||||
|
!
|
||||||
|
du(:,:,:,:) = 0.0d0
|
||||||
|
|
||||||
|
! perform update along the X direction
|
||||||
|
!
|
||||||
|
do i = 2, im
|
||||||
|
du(:,i,:,:) = du(:,i,:,:) - dh(1) * (f(1,:,i,:,:) - f(1,:,i-1,:,:))
|
||||||
|
end do
|
||||||
|
du(:,1,:,:) = du(:,1,:,:) - dh(1) * f(1,:,1,:,:)
|
||||||
|
|
||||||
|
! perform update along the Y direction
|
||||||
|
!
|
||||||
|
do j = 2, jm
|
||||||
|
du(:,:,j,:) = du(:,:,j,:) - dh(2) * (f(2,:,:,j,:) - f(2,:,:,j-1,:))
|
||||||
|
end do
|
||||||
|
du(:,:,1,:) = du(:,:,1,:) - dh(2) * f(2,:,:,1,:)
|
||||||
|
|
||||||
|
#if NDIMS == 3
|
||||||
|
! perform update along the Z direction
|
||||||
|
!
|
||||||
|
do k = 2, km
|
||||||
|
du(:,:,:,k) = du(:,:,:,k) - dh(3) * (f(3,:,:,:,k) - f(3,:,:,:,k-1))
|
||||||
|
end do
|
||||||
|
du(:,:,:,1) = du(:,:,:,1) - dh(3) * f(3,:,:,:,1)
|
||||||
|
#endif /* NDIMS == 3 */
|
||||||
|
|
||||||
|
!-------------------------------------------------------------------------------
|
||||||
|
!
|
||||||
|
end subroutine update_interval
|
||||||
#else /* CONSERVATIVE */
|
#else /* CONSERVATIVE */
|
||||||
!
|
!
|
||||||
!===============================================================================
|
!===============================================================================
|
||||||
|
Loading…
x
Reference in New Issue
Block a user