diff --git a/src/evolution.F90 b/src/evolution.F90 index f75e9bb..f07d0ed 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -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 ! !=============================================================================== !