Merge branch 'master' into reconnection
This commit is contained in:
commit
8c6ec027d1
@ -41,7 +41,8 @@ module evolution
|
||||
|
||||
! evolution parameters
|
||||
!
|
||||
real(kind=8), save :: cfl = 0.5d+00
|
||||
real(kind=8), save :: cfl = 5.0d-01
|
||||
integer , save :: stages = 2
|
||||
|
||||
! coefficient controlling the decay of scalar potential ѱ
|
||||
!
|
||||
@ -96,7 +97,8 @@ module evolution
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use parameters, only : get_parameter_string, get_parameter_real
|
||||
use parameters, only : get_parameter_string, get_parameter_real &
|
||||
, get_parameter_integer
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
@ -116,9 +118,10 @@ module evolution
|
||||
!
|
||||
! get the integration method and the value of the CFL coefficient
|
||||
!
|
||||
call get_parameter_string("time_advance", integration)
|
||||
call get_parameter_real ("cfl" , cfl )
|
||||
call get_parameter_real ("alpha" , alpha )
|
||||
call get_parameter_string ("time_advance", integration)
|
||||
call get_parameter_integer("stages" , stages )
|
||||
call get_parameter_real ("cfl" , cfl )
|
||||
call get_parameter_real ("alpha" , alpha )
|
||||
|
||||
! select the integration method, check the correctness of the integration
|
||||
! parameters and adjust the CFL coefficient if necessary
|
||||
@ -135,17 +138,30 @@ module evolution
|
||||
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(',i1,',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")
|
||||
case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)")
|
||||
|
||||
name_int = "3rd order 4-stage Runge-Kutta"
|
||||
evolve => evolve_rk34
|
||||
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 default
|
||||
|
||||
if (verbose) then
|
||||
@ -626,6 +642,208 @@ module evolution
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! 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 procedures
|
||||
!
|
||||
use boundaries , only : boundary_variables
|
||||
use sources , only : update_sources
|
||||
|
||||
! include external variables
|
||||
!
|
||||
use blocks , only : block_data, list_data
|
||||
use coordinates , only : im, jm, km
|
||||
use equations , only : nv, ibp
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: n
|
||||
real(kind=8) :: ds
|
||||
|
||||
! local saved variables
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
real(kind=8), save :: ft, fl, fr
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(nv,im,jm,km) :: du
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! 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(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! 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
|
||||
|
||||
! 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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
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)
|
||||
!
|
||||
! 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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = fl * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ fr * (pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
|
||||
! assign pointer to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk2
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine EVOLVE_RK3:
|
||||
! ---------------------
|
||||
!
|
||||
@ -820,11 +1038,11 @@ module evolution
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine EVOLVE_RK34:
|
||||
! ----------------------
|
||||
! subroutine EVOLVE_SSPRK34:
|
||||
! -------------------------
|
||||
!
|
||||
! Subroutine advances the solution by one time step using the 3rd order
|
||||
! 4-stage Runge-Kutta time integration method.
|
||||
! 4-stage Strong Stability Preserving Runge-Kutta time integration method.
|
||||
!
|
||||
! References:
|
||||
!
|
||||
@ -832,12 +1050,11 @@ module evolution
|
||||
! "Global Optimization of Explicit Strong-Stability-Preserving
|
||||
! Runge-Kutta methods",
|
||||
! Mathematics of Computation,
|
||||
! 2006, volume 75, number 253, pp. 183-207
|
||||
!
|
||||
! 2006, vol. 75, no. 253, pp. 183-207
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine evolve_rk34()
|
||||
subroutine evolve_ssprk34()
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
@ -856,7 +1073,7 @@ module evolution
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pblock
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local variables
|
||||
!
|
||||
@ -873,42 +1090,46 @@ module evolution
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
! = 1st step of the integration: U(1) = U(n) + 1/2 dt F[U(n)] =
|
||||
!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
|
||||
!
|
||||
! calculate fractional time step
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = b1 * dt
|
||||
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u1
|
||||
pdata%u => pdata%u1
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -918,35 +1139,38 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
! = 2nd step of the integration: U(2) = U(1) + 1/2 dt F[U(1)] =
|
||||
! (we can copy U(2) directly to U(1), since U(1) is not used anymore)
|
||||
!= 2nd step: U(2) = U(1) + 1/2 dt F[U(1)]
|
||||
!
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u1(1:nv,:,:,:) = pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! assign pointer to the next block
|
||||
! update the intermediate solution
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
end do
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -956,40 +1180,43 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
! = 3rd step of the integration: U(3) = 2/3 U(n) + 1/3 U(2) + 1/6 dt F[U(2)] =
|
||||
! (we can copy U(3) directly to U(1), since U(1) is not used anymore)
|
||||
!= 3rd step: U(3) = 2/3 U(n) + 1/3 U(2) + 1/6 dt F[U(2)]
|
||||
!
|
||||
! calculate fractional time step
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = b3 * dt
|
||||
|
||||
! update fluxes for the first step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u1(1:nv,:,:,:) = a31 * pblock%u0(1:nv,:,:,:) &
|
||||
+ a33 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! assign pointer to the next block
|
||||
! update the intermediate solution
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = a31 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ a33 * pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
end do
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -999,46 +1226,51 @@ module evolution
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
! = update the final solution: U(n+1) = U(3) + 1/2 dt F[U(3)]
|
||||
!= the final step: U(n+1) = U(3) + 1/2 dt F[U(3)]
|
||||
!
|
||||
! calculate fractional time step
|
||||
!
|
||||
ds = b1 * dt
|
||||
|
||||
! update fluxes for the second step of the RK2 integration
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! update the solution using numerical fluxes stored in the data blocks
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pblock => list_data
|
||||
do while (associated(pblock))
|
||||
pdata => list_data
|
||||
|
||||
! calculate variable increment for the current block
|
||||
! iterate over all data blocks
|
||||
!
|
||||
call update_increment(pblock, du(:,:,:,:))
|
||||
do while (associated(pdata))
|
||||
|
||||
! add source terms
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_sources(pblock, du(:,:,:,:))
|
||||
call update_increment(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the solution for the fluid variables
|
||||
! add the source terms
|
||||
!
|
||||
pblock%u0(1:nv,:,:,:) = pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:)
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pblock%u => pblock%u0
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ by its source term
|
||||
!
|
||||
if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:)
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
|
||||
! assign pointer to the next block
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pblock => pblock%next
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
@ -1050,7 +1282,310 @@ module evolution
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_rk34
|
||||
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 procedures
|
||||
!
|
||||
use boundaries , only : boundary_variables
|
||||
use sources , only : update_sources
|
||||
|
||||
! include external variables
|
||||
!
|
||||
use blocks , only : block_data, list_data
|
||||
use coordinates , only : im, jm, km
|
||||
use equations , only : nv, ibp
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local variables
|
||||
!
|
||||
real(kind=8) :: ds
|
||||
|
||||
! local arrays
|
||||
!
|
||||
real(kind=8), dimension(nv,im,jm,km) :: du
|
||||
|
||||
! 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.44090224936674d-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
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
!= 1st step: U(1) = U(n) + b1 dt F[U(n)]
|
||||
!
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = b1 * 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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! 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()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
!= 2nd step: U(2) = U(1) + b1 dt F[U(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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)]
|
||||
!
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = b3 * 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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(1:nv,1:im,1:jm,1:km) = a31 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ a33 * pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)]
|
||||
!
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = b4 * 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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = a41 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ a44 * pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! 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()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
!= 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
|
||||
|
||||
! 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, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, du(1:nv,1:im,1:jm,1:km))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(1:nv,1:im,1:jm,1:km) = a53 * pdata%u1(1:nv,1:im,1:jm,1:km) &
|
||||
+ a55 * pdata%u0(1:nv,1:im,1:jm,1:km) &
|
||||
+ ds * du(1:nv,1:im,1:jm,1:km)
|
||||
|
||||
! update ψ by its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = &
|
||||
decay * pdata%u(ibp,1:im,1:jm,1:km)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables()
|
||||
|
||||
! update boundaries
|
||||
!
|
||||
call boundary_variables()
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk35
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user