!!****************************************************************************** !! !! This file is part of the AMUN source code, a program to perform !! Newtonian or relativistic magnetohydrodynamical simulations on uniform or !! adaptive mesh. !! !! Copyright (C) 2008-2021 Grzegorz Kowal !! !! This program is free software: you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation, either version 3 of the License, or !! (at your option) any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program. If not, see . !! !!****************************************************************************** !! !! module: EVOLUTION !! !! This module provides an interface for temporal integration with !! the stability handling. New integration methods can be added by !! implementing more evolve_* subroutines. !! !!****************************************************************************** ! module evolution #ifdef PROFILE ! import external subroutines ! use timers, only : set_timer, start_timer, stop_timer #endif /* PROFILE */ ! module variables are not implicit by default ! implicit none #ifdef PROFILE ! timer indices ! integer , save :: imi, ima, imt, imu, imf, iui, imv #endif /* PROFILE */ ! interfaces for procedure pointers ! abstract interface subroutine evolve_iface() end subroutine end interface ! pointer to the temporal integration subroutine ! procedure(evolve_iface), pointer, save :: evolve => null() ! evolution parameters ! character(len=255), save :: name_int = "" integer , save :: stages = 2 integer , save :: registers = 1 real(kind=8) , save :: cfl = 5.0d-01 ! coefficients controlling the decay of scalar potential ѱ ! real(kind=8) , save :: glm_alpha = 5.0d-01 real(kind=8) , save :: decay = 7.788007830714049d-01 ! time variables ! integer , save :: step = 0 real(kind=8) , save :: time = 0.0d+00 real(kind=8) , save :: dt = 1.0d+00 real(kind=8) , save :: dtn = 1.0d+00 real(kind=8) , save :: dte = 0.0d+00 ! the absolute and relative tolerances, limiting factors, the maximum error, ! the maximum number of passes for the adaptive step, ! the total number of integration iterations, the number of rejected iterations ! real(kind=8) , save :: atol = 1.0d-04 real(kind=8) , save :: rtol = 1.0d-04 real(kind=8) , save :: fac = 9.0d-01 real(kind=8) , save :: facmin = 1.0d-01 real(kind=8) , save :: facmax = 5.0d+00 real(kind=8) , save :: maxerr = 0.0d+00 integer , save :: mrej = 5 integer , save :: niterations = 0 integer , save :: nrejections = 0 ! errors in three recent steps ! real(kind=8), dimension(3), save :: errs ! by default everything is private ! private ! declare public subroutines ! public :: initialize_evolution, finalize_evolution, print_evolution public :: advance, new_time_step ! declare public variables ! public :: step, time, dt, dtn, dte, cfl, glm_alpha, registers public :: atol, rtol, mrej, niterations, nrejections, errs, maxerr !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine INITIALIZE_EVOLUTION: ! ------------------------------- ! ! Subroutine initializes module EVOLUTION by setting its parameters. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine initialize_evolution(verbose, status) ! include external procedures ! use parameters, only : get_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose integer, intent(out) :: status ! local variables ! character(len=255) :: integration = "rk2" integer :: n ! local parameters ! character(len=*), parameter :: loc = 'EVOLUTION::initialize_evolution()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! set timer descriptions ! call set_timer('evolution:: initialization' , imi) call set_timer('evolution:: solution advance', ima) call set_timer('evolution:: new time step' , imt) call set_timer('evolution:: solution update' , imu) call set_timer('evolution:: flux update' , imf) call set_timer('evolution:: increment update', iui) call set_timer('evolution:: variable update' , imv) ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! reset the status flag ! status = 0 ! get the integration method and the value of the CFL coefficient ! call get_parameter("time_advance", integration) call get_parameter("stages" , stages ) call get_parameter("cfl" , cfl ) call get_parameter("glm_alpha" , glm_alpha ) ! get the absolute and relative tolerances, the maximum number of passes for ! the adaptive step, and limiting factors ! call get_parameter("absolute_tolerance", atol) call get_parameter("relative_tolerance", rtol) call get_parameter("maximum_rejections", mrej) call get_parameter("factor" , fac ) call get_parameter("minimum_factor" , facmin) call get_parameter("maximum_factor" , facmax) ! select the integration method, check the correctness of the integration ! parameters and adjust the CFL coefficient if necessary ! select case(trim(integration)) case ("euler", "EULER", "rk1", "RK1") evolve => evolve_euler registers = 1 name_int = "1st order Euler" case ("rk2", "RK2") evolve => evolve_rk2 registers = 2 name_int = "2nd order Runge-Kutta" case ("ssprk(m,2)", "SSPRK(m,2)") evolve => evolve_ssprk2_m registers = 3 stages = max(2, min(9, stages)) cfl = (stages - 1) * cfl write(name_int, "('2nd order SSPRK(',i0,',2)')") stages case ("rk3", "RK3") evolve => evolve_rk3 registers = 2 name_int = "3rd order Runge-Kutta" case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)") evolve => evolve_ssprk34 registers = 2 cfl = 2.0d+00 * cfl name_int = "3rd order SSPRK(4,3)" case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)") evolve => evolve_ssprk35 registers = 2 cfl = 2.65062919143939d+00 * cfl name_int = "3rd order SSPRK(5,3)" case ("rk3.m", "ssprk(m,3)", "SSPRK(m,3)") evolve => evolve_ssprk3_m registers = 3 n = 2 do while(n**2 <= stages) n = n + 1 end do n = n - 1 stages = max(4, n**2) cfl = (n - 1) * n * cfl write(name_int, "('Optimal 3rd order SSPRK(',i0,',3)')") stages case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") evolve => evolve_ssprk4_10 registers = 2 cfl = 6.0d+00 * cfl name_int = "Optimal 4th order SSPRK(10,4)" case default if (verbose) then write(*,*) write(*,"(1x,a)") "ERROR!" write(*,"(1x,a)") "The selected time advance method is not " // & "implemented: " // trim(integration) write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," // & " 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," // & " 'ssprk(m,3)', 'ssprk(10,4)'." end if status = 1 end select ! proceed if everything is fine ! if (status == 0) then ! calculate the decay factor for magnetic field divergence scalar source term ! decay = exp(- glm_alpha * cfl) end if ! status ! reset recent error history ! errs = 1.0d+00 #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine initialize_evolution ! !=============================================================================== ! ! subroutine FINALIZE_EVOLUTION: ! ----------------------------- ! ! Subroutine releases memory used by the module. ! ! Arguments: ! ! status - an integer flag for error return value; ! !=============================================================================== ! subroutine finalize_evolution(verbose, status) ! include external procedures ! use helpers , only : print_section, print_parameter use iso_fortran_env, only : error_unit ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose integer, intent(out) :: status ! local parameters ! character(len=*), parameter :: loc = 'EVOLUTION::finalize_evolution()' ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for module initialization/finalization ! call start_timer(imi) #endif /* PROFILE */ ! reset the status flag ! status = 0 ! nullify pointer to integration subroutine ! nullify(evolve) ! print integration summary ! if (niterations > 0) then call print_section(verbose, "Integration summary") call print_parameter(verbose, "Number of iterations", niterations) call print_parameter(verbose, "Number of rejections", nrejections) end if #ifdef PROFILE ! stop accounting time for module initialization/finalization ! call stop_timer(imi) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine finalize_evolution ! !=============================================================================== ! ! subroutine PRINT_EVOLUTION: ! -------------------------- ! ! Subroutine prints module parameters. ! ! Arguments: ! ! verbose - a logical flag turning the information printing; ! !=============================================================================== ! subroutine print_evolution(verbose) ! import external procedures and variables ! use equations, only : magnetized use helpers , only : print_section, print_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose ! !------------------------------------------------------------------------------- ! if (verbose) then call print_section(verbose, "Evolution") call print_parameter(verbose, "time advance", name_int) call print_parameter(verbose, "no. of registers", registers) call print_parameter(verbose, "CFL coefficient", cfl) if (magnetized) then call print_parameter(verbose, "GLM alpha coefficient", glm_alpha) end if call print_parameter(verbose, "absolute tolerance", atol) call print_parameter(verbose, "relative tolerance", rtol) call print_parameter(verbose, "maximum rejections", mrej) call print_parameter(verbose, "factor" , fac) call print_parameter(verbose, "minimum factor" , facmin) call print_parameter(verbose, "maximum factor" , facmax) end if !------------------------------------------------------------------------------- ! end subroutine print_evolution ! !=============================================================================== ! ! subroutine ADVANCE: ! ------------------ ! ! Subroutine advances the solution by one time step using the selected time ! integration method. ! ! Arguments: ! ! dtnext - the next time step; ! status - the subroutine call status: 0 for success, otherwise failure; ! !=============================================================================== ! subroutine advance(dtnext, status) ! references ! use blocks , only : set_blocks_update use coordinates, only : toplev use forcing , only : update_forcing, forcing_enabled use mesh , only : update_mesh ! local variables are not implicit by default ! implicit none ! input variables ! real(kind=8), intent(in) :: dtnext integer , intent(out) :: status ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for solution advance ! call start_timer(ima) #endif /* PROFILE */ ! find new time step ! call new_time_step(dtnext) ! advance the solution using the selected method ! call evolve() ! add forcing contribution ! if (forcing_enabled) call update_forcing(dt) ! check if we need to perform the refinement step ! if (toplev > 1) then ! set all meta blocks to not be updated ! call set_blocks_update(.false.) ! check refinement and refine ! call update_mesh(status) if (status /= 0) go to 100 ! update primitive variables ! call update_variables(time + dt, 0.0d+00) ! set all meta blocks to be updated ! call set_blocks_update(.true.) end if ! toplev > 1 #ifdef DEBUG ! check variables for NaNs ! call check_variables() #endif /* DEBUG */ ! error entry point ! 100 continue #ifdef PROFILE ! stop accounting time for solution advance ! call stop_timer(ima) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine advance ! !=============================================================================== ! ! subroutine NEW_TIME_STEP: ! ------------------------ ! ! Subroutine estimates the new time step from the maximum speed in the system ! and source term constraints. ! ! Arguments: ! ! dtnext - next time step; ! !=============================================================================== ! subroutine new_time_step(dtnext) ! include external procedures ! use equations , only : maxspeed, cmax, cmax2 #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ ! include external variables ! use blocks , only : block_data, list_data use coordinates , only : adx, ady #if NDIMS == 3 use coordinates , only : adz #endif /* NDIMS == 3 */ use sources , only : viscosity, resistivity ! local variables are not implicit by default ! implicit none ! subroutine arguments ! real(kind=8), intent(in) :: dtnext ! local pointers ! type(block_data), pointer :: pdata ! local variables ! integer :: mlev real(kind=8) :: cm, dx_min ! local parameters ! real(kind=8), parameter :: eps = tiny(cmax) ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for new time step estimation ! call start_timer(imt) #endif /* PROFILE */ ! reset the maximum speed, and the highest level ! cmax = eps mlev = 1 ! assign pdata with the first block on the data block list ! pdata => list_data ! iterate over all data blocks in order to find the maximum speed among them ! and the highest level which is required to obtain the minimum spacial step ! do while (associated(pdata)) ! update the maximum level ! mlev = max(mlev, pdata%meta%level) ! get the maximum speed for the current block ! cm = maxspeed(pdata%q(:,:,:,:)) ! update the maximum speed ! cmax = max(cmax, cm) ! assign pdata to the next block ! pdata => pdata%next end do ! over data blocks #ifdef MPI ! reduce maximum speed and level over all processors ! call reduce_maximum(cmax) call reduce_maximum(mlev) #endif /* MPI */ ! calculate the squared cmax ! cmax2 = cmax * cmax ! find the smallest spacial step ! #if NDIMS == 2 dx_min = min(adx(mlev), ady(mlev)) #endif /* NDIMS == 2 */ #if NDIMS == 3 dx_min = min(adx(mlev), ady(mlev), adz(mlev)) #endif /* NDIMS == 3 */ ! calculate the new time step ! dtn = cfl * dx_min / max(cmax & + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps) ! round the time ! if (dte <= 0.0d+00) dte = dtn if (dtnext > 0.0d+00) dt = min(dtn, dte, dtnext) #ifdef PROFILE ! stop accounting time for new time step estimation ! call stop_timer(imt) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine new_time_step ! !=============================================================================== !! !!*** PRIVATE SUBROUTINES **************************************************** !! !=============================================================================== ! !=============================================================================== ! ! subroutine EVOLVE_EULER: ! ----------------------- ! ! Subroutine advances the solution by one time step using the 1st order ! Euler integration method. ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine evolve_euler() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp use sources , only : update_sources implicit none type(block_data), pointer :: pdata real(kind=8) :: tm, dtm ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ tm = time + dt dtm = dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_euler ! !=============================================================================== ! ! subroutine EVOLVE_RK2: ! --------------------- ! ! Subroutine advances the solution by one time step using the 2nd order ! Runge-Kutta time integration method. ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! !=============================================================================== ! subroutine evolve_rk2() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp use sources , only : update_sources implicit none type(block_data), pointer :: pdata real(kind=8) :: tm, dtm ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ != 1st step: U(1) = U(n) + dt * F[U(n)] ! tm = time + dt dtm = dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm) != 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)] ! tm = time + dt dtm = 0.5d+00 * dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = 0.5d+00 * (pdata%uu(:,:,:,:,1) & + pdata%uu(:,:,:,:,2) & + dt * pdata%du(:,:,:,:)) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_rk2 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK2_M: ! -------------------------- ! ! 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_m() use blocks , only : block_data, list_data, get_dblocks use boundaries , only : boundary_fluxes use coordinates, only : nc => ncells, nb, ne use equations , only : errors, nf, ibp #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ use sources , only : update_sources implicit none type(block_data), pointer :: pdata logical :: test integer :: n, l, nrej real(kind=8) :: tm, dtm, ds, umax, utol, emax real(kind=8) :: fc, fcmn, fcmx logical , save :: first = .true. real(kind=8), save :: ft, fl, fr, gt, gl real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ if (first) then ft = 1.0d+00 / (stages - 1) fl = 1.0d+00 / stages fr = 1.0d+00 - fl gt = fl - ft gl = fl first = .false. end if fc = fac fcmn = facmin fcmx = facmax l = get_dblocks() #if NDIMS == 3 allocate(lerr(l,nf,nc,nc,nc)) #else /* NDIMS == 3 */ allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ test = .true. nrej = 0 do while(test) lerr(:,:,:,:,:) = 0.0d+00 ds = ft * dt != 1st step: U(1) = U(n), U(2) = U(n) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 ! do n = 1, stages - 1 tm = time + n * ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() l = 1 pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do call update_variables(tm, dtm) end do ! n = 1, stages - 1 != the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] ! tm = time + dt dtm = fl * dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() l = 1 pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = fr * pdata%uu(:,:,:,:,2) & + fl * (pdata%uu(:,:,:,:,3) & + dt * pdata%du(:,:,:,:)) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do call update_variables(tm, dtm) ! find umax ! umax = 0.0d+00 pdata => list_data do while (associated(pdata)) #if NDIMS == 3 umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))), & maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,2)))) #else /* NDIMS == 3 */ umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,1))), & maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,2)))) #endif /* NDIMS == 3 */ pdata => pdata%next end do ! get the maximum error of each variable ! do l = 1, nf errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) end do #ifdef MPI call reduce_maximum(umax) call reduce_maximum(errors) #endif /* MPI */ ! calculate tolerance and time step ! maxerr = maxval(errors) utol = atol + rtol * umax emax = maxerr / utol if (emax <= 1.0d+00 .or. nrej >= mrej) then test = .false. errs(3) = errs(2) errs(2) = errs(1) errs(1) = emax maxerr = maxerr / umax dte = dt * min(fcmx, max(fcmn, & fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3)) fcmx = facmax else errs(1) = emax dte = dt * min(fcmx, max(fcmn, & fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3)) dt = dte fcmx = fac nrej = nrej + 1 ! rejection count in the current step nrejections = nrejections + 1 end if niterations = niterations + 1 end do if (allocated(lerr)) deallocate(lerr) != final step: U(n+1) = U(1) ! tm = time + dt dtm = ft * dt pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) ! update ψ with its source term ! if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk2_m ! !=============================================================================== ! ! subroutine EVOLVE_RK3: ! --------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! Runge-Kutta time integration method. ! ! References: ! ! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., ! "Numerical Recipes in Fortran", ! Cambridge University Press, Cambridge, 1992 ! ! !=============================================================================== ! subroutine evolve_rk3() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp use sources , only : update_sources implicit none type(block_data), pointer :: pdata real(kind=8) :: ds real(kind=8) :: tm, dtm real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00 real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ != 1st substep: U(1) = U(n) + dt F[U(n)] ! ds = dt tm = time + ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm) != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)] ! ds = f22 * dt tm = time + 0.5d+00 * dt dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) & + f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) != 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)] ! ds = f32 * dt tm = time + dt dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = f31 * pdata%uu(:,:,:,:,1) & + f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_rk3 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK34: ! ------------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! 4-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Ruuth, S. J., ! "Global Optimization of Explicit Strong-Stability-Preserving ! Runge-Kutta methods", ! Mathematics of Computation, ! 2006, vol. 75, no. 253, pp. 183-207 ! !=============================================================================== ! subroutine evolve_ssprk34() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp use sources , only : update_sources implicit none type(block_data), pointer :: pdata real(kind=8) :: ds real(kind=8) :: tm, dtm ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ != 1st step: U(1) = U(n) + 1/2 dt F[U(n)] ! ds = dt / 2.0d+00 tm = time + ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm) != 2nd step: U(2) = U(2) + 1/2 dt F[U(1)] ! tm = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dtm * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) != 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)]) ! tm = time + ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) & + pdata%uu(:,:,:,:,2) & + dtm * pdata%du(:,:,:,:)) / 3.0d+00 pdata => pdata%next end do call update_variables(tm, dtm) != the final step: U(n+1) = U(3) + 1/2 dt F[U(3)] ! tm = time + dt pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk34 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK35: ! ------------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! 5-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Ruuth, S. J., ! "Global Optimization of Explicit Strong-Stability-Preserving ! Runge-Kutta methods", ! Mathematics of Computation, ! 2006, vol. 75, no. 253, pp. 183-207 ! !=============================================================================== ! subroutine evolve_ssprk35() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp use sources , only : update_sources implicit none type(block_data), pointer :: pdata real(kind=8) :: ds real(kind=8) :: tm, dtm real(kind=8), parameter :: b1 = 3.77268915331368d-01 real(kind=8), parameter :: b3 = 2.42995220537396d-01 real(kind=8), parameter :: b4 = 2.38458932846290d-01 real(kind=8), parameter :: b5 = 2.87632146308408d-01 real(kind=8), parameter :: a31 = 3.55909775063327d-01 real(kind=8), parameter :: a33 = 6.44090224936673d-01 real(kind=8), parameter :: a41 = 3.67933791638137d-01 real(kind=8), parameter :: a44 = 6.32066208361863d-01 real(kind=8), parameter :: a53 = 2.37593836598569d-01 real(kind=8), parameter :: a55 = 7.62406163401431d-01 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ != 1st step: U(1) = U(n) + b1 dt F[U(n)] ! ds = b1 * dt tm = time + ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm) != 2nd step: U(2) = U(1) + b1 dt F[U(1)] ! tm = time + 2.0d+00 * ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) != 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] ! ds = b3 * dt tm = time + (2.0d+00 * a33 * b1 + b3) * dt dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) & + a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) != 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] ! ds = b4 * dt tm = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = a41 * pdata%uu(:,:,:,:,1) & + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do call update_variables(tm, dtm) != the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] ! ds = b5 * dt tm = time + dt dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = a53 * pdata%uu(:,:,:,:,2) & + a55 * pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) pdata => pdata%next end do if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk35 ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK3_M: ! -------------------------- ! ! Subroutine advances the solution by one time step using the 3rd order ! m-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Gottlieb, S., Ketcheson, D., Shu C. - W., ! "Strong stability preserving Runge-Kutta and multistep ! time discretizations", ! World Scientific Publishing, 2011 ! !=============================================================================== ! subroutine evolve_ssprk3_m() use blocks , only : block_data, list_data, get_dblocks use boundaries , only : boundary_fluxes use coordinates, only : nc => ncells, nb, ne use equations , only : errors, nf, ibp #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ use sources , only : update_sources implicit none type(block_data), pointer :: pdata logical , save :: first = .true. integer , save :: n, n1, n2, n3, n4 real(kind=8), save :: r, f1, f2, g1, g2 logical :: test integer :: i, l, nrej real(kind=8) :: tm, dtm, ds, umax, utol, emax real(kind=8) :: fc, fcmn, fcmx real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr real(kind=8), parameter :: k1 = -5.8d-01 / 3.0d+00, k2 = 7.0d-02, & k3 = -1.0d-01 / 3.0d+00 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ if (first) then n = int(sqrt(1.0d+00 * stages)) n1 = (n - 1) * (n - 2) / 2 n2 = n1 + 1 n3 = n * (n + 1) / 2 - 1 n4 = n * (n + 1) / 2 + 1 r = (1.0d+00 * (2 * n - 1)) f1 = (1.0d+00 * n ) / r f2 = (1.0d+00 * (n - 1)) / r r = 1.0d+00 * (stages - n) g1 = 1.0d+00 / ( stages - n) - 1.0d+00 / stages g2 = 1.0d+00 / (2 * stages - n) - 1.0d+00 / stages first = .false. end if fc = fac fcmn = facmin fcmx = facmax l = get_dblocks() #if NDIMS == 3 allocate(lerr(l,nf,nc,nc,nc)) #else /* NDIMS == 3 */ allocate(lerr(l,nf,nc,nc, 1)) #endif /* NDIMS == 3 */ test = .true. nrej = 0 do while(test) lerr(:,:,:,:,:) = 0.0d+00 ds = dt / r != 1st step: U(1) = U(n) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2 ! do i = 1, n1 tm = time + i * ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() l = 1 pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do call update_variables(tm, dtm) end do ! n = 1, n1 != 3rd step: U(2) = U(1) ! ! iterate over all data blocks ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,2) pdata => pdata%next end do != 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1 ! do i = n2, n3 tm = time + i * ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() l = 1 pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do call update_variables(tm, dtm) end do ! i = n2, n3 != 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1) ! tm = time + dt dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() l = 1 pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,3) & + f2 * (pdata%uu(:,:,:,:,2) & + ds * pdata%du(:,:,:,:)) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do call update_variables(tm, dtm) != 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m ! do i = n4, stages tm = time + (i - n) * ds dtm = ds pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() l = 1 pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) #if NDIMS == 3 lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ l = l + 1 pdata => pdata%next end do call update_variables(tm, dtm) end do ! i = n4, stages ! find umax ! umax = 0.0d+00 pdata => list_data do while (associated(pdata)) #if NDIMS == 3 umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))), & maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,2)))) #else /* NDIMS == 3 */ umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,1))), & maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,2)))) #endif /* NDIMS == 3 */ pdata => pdata%next end do ! get the maximum error of each variable ! do l = 1, nf errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) end do #ifdef MPI call reduce_maximum(umax) call reduce_maximum(errors) #endif /* MPI */ ! calculate tolerance and time step ! maxerr = maxval(errors) utol = atol + rtol * umax emax = maxerr / utol if (emax <= 1.0d+00 .or. nrej >= mrej) then test = .false. errs(3) = errs(2) errs(2) = errs(1) errs(1) = emax maxerr = maxerr / umax dte = dt * min(fcmx, max(fcmn, & fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3)) fcmx = facmax else errs(1) = emax dte = dt * min(fcmx, max(fcmn, & fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3)) dt = dte fcmx = fac nrej = nrej + 1 ! rejection count in the current step nrejections = nrejections + 1 end if niterations = niterations + 1 end do if (allocated(lerr)) deallocate(lerr) != final step: U(n+1) = U(1) ! tm = time + dt dtm = dt / r pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,1) pdata => pdata%next end do call update_variables(tm, dtm) if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk3_m ! !=============================================================================== ! ! subroutine EVOLVE_SSPRK4_10: ! --------------------------- ! ! Subroutine advances the solution by one time step using the 4rd order ! 10-stage Strong Stability Preserving Runge-Kutta time integration method. ! ! References: ! ! [1] Gottlieb, S., Ketcheson, D., Shu C. - W., ! "Strong stability preserving Runge-Kutta and multistep ! time discretizations", ! World Scientific Publishing, 2011 ! !=============================================================================== ! subroutine evolve_ssprk4_10() use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes use equations , only : ibp use sources , only : update_sources implicit none type(block_data), pointer :: pdata integer :: n real(kind=8) :: tm, dtm real(kind=8), dimension(9) :: ds real(kind=8), dimension(9), parameter :: c = & (/ 1.0d+00, 2.0d+00, 3.0d+00, 4.0d+00, 2.0d+00, & 3.0d+00, 4.0d+00, 5.0d+00, 6.0d+00 /) / 6.0d+00 ! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ ds(:) = c(:) * dt != 1st step: U(2) = U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata => pdata%next end do != 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5 ! do n = 1, 5 tm = time + ds(n) dtm = ds(1) pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) end do ! n = 1, 5 != 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1) ! pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,2) = (pdata%uu(:,:,:,:,2) & + 9.0d+00 * pdata%uu(:,:,:,:,1)) / 2.5d+01 pdata%uu(:,:,:,:,1) = 1.5d+01 * pdata%uu(:,:,:,:,2) & - 5.0d+00 * pdata%uu(:,:,:,:,1) pdata => pdata%next end do != 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9 ! ! integrate the intermediate steps ! do n = 6, 9 tm = time + ds(n) dtm = ds(1) pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) end do ! n = 6, 9 != the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] ! tm = time + dt dtm = dt / 1.0d+01 pdata => list_data do while (associated(pdata)) call update_increment(pdata) call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) pdata => pdata%next end do call boundary_fluxes() pdata => list_data do while (associated(pdata)) pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) & + 6.0d-01 * pdata%uu(:,:,:,:,1) & + dtm * pdata%du(:,:,:,:) pdata => pdata%next end do if (ibp > 0) then pdata => list_data do while (associated(pdata)) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) pdata => pdata%next end do end if call update_variables(tm, dtm) #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk4_10 ! !=============================================================================== ! ! subroutine UPDATE_INCREMENT: ! --------------------------- ! ! Subroutine calculates the conservative variable increment from ! directional fluxes. ! ! Arguments: ! ! pdata - the point to data block storing the directional fluxes; ! !=============================================================================== ! subroutine update_increment(pdata) ! include external variables ! use blocks , only : block_data use coordinates, only : nbl, ne, neu, nn => bcells use coordinates, only : adx, ady, adxi, adyi #if NDIMS == 3 use coordinates, only : adz, adzi #endif /* NDIMS == 3 */ use equations , only : nf, ns use equations , only : idn, isl, isu use schemes , only : update_flux ! local variables are not implicit by default ! implicit none ! subroutine arguments ! type(block_data), pointer, intent(inout) :: pdata ! local variables ! integer :: i , j , k = 1, p integer :: im1, ip1 integer :: jm1, jp1 #if NDIMS == 3 integer :: km1, kp1 #endif /* NDIMS == 3 */ real(kind=8) :: df real(kind=8), dimension(NDIMS) :: dh, dhi #if NDIMS == 3 real(kind=8), dimension(nf,nn,nn,nn,3) :: f #else /* NDIMS == 3 */ real(kind=8), dimension(nf,nn,nn, 1,2) :: f #endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for the increment update ! call start_timer(iui) #endif /* PROFILE */ ! prepare the coordinate intervals ! dh(1) = adx(pdata%meta%level) dh(2) = ady(pdata%meta%level) #if NDIMS == 3 dh(3) = adz(pdata%meta%level) #endif /* NDIMS == 3 */ dhi(1) = adxi(pdata%meta%level) dhi(2) = adyi(pdata%meta%level) #if NDIMS == 3 dhi(3) = adzi(pdata%meta%level) #endif /* NDIMS == 3 */ ! calculate the interface fluxes ! call update_flux(dh(1:NDIMS), pdata%q(:,:,:,:), f(:,:,:,:,1:NDIMS)) ! store block interface fluxes ! pdata%fx(:,1,:,:) = f(:,nbl,:,:,1) pdata%fx(:,2,:,:) = f(:,ne ,:,:,1) pdata%fy(:,:,1,:) = f(:,:,nbl,:,2) pdata%fy(:,:,2,:) = f(:,:,ne ,:,2) #if NDIMS == 3 pdata%fz(:,:,:,1) = f(:,:,:,nbl,3) pdata%fz(:,:,:,2) = f(:,:,:,ne ,3) #endif /* NDIMS == 3 */ ! calculate the variable update from the directional fluxes ! #if NDIMS == 3 do k = nbl, neu km1 = k - 1 #endif /* NDIMS == 3 */ do j = nbl, neu jm1 = j - 1 do i = nbl, neu im1 = i - 1 #if NDIMS == 3 pdata%du(1:nf,i,j,k) = - dhi(1) * (f(:,i,j,k,1) - f(:,im1,j,k,1)) & - dhi(2) * (f(:,i,j,k,2) - f(:,i,jm1,k,2)) & - dhi(3) * (f(:,i,j,k,3) - f(:,i,j,km1,3)) #else /* NDIMS == 3 */ pdata%du(1:nf,i,j,k) = - dhi(1) * (f(:,i,j,k,1) - f(:,im1,j,k,1)) & - dhi(2) * (f(:,i,j,k,2) - f(:,i,jm1,k,2)) #endif /* NDIMS == 3 */ end do ! i = nbl, neu end do ! j = nbl, neu #if NDIMS == 3 end do ! k = nbl, neu #endif /* NDIMS == 3 */ ! update passive scalars ! if (ns > 0) then ! reset passive scalar increments ! pdata%du(isl:isu,:,:,:) = 0.0d+00 #if NDIMS == 3 do k = nbl - 1, neu kp1 = k + 1 #endif /* NDIMS == 3 */ do j = nbl - 1, neu jp1 = j + 1 do i = nbl - 1, neu ip1 = i + 1 ! X-face ! if (f(idn,i,j,k,1) >= 0.0d+00) then do p = isl, isu df = dhi(1) * f(idn,i,j,k,1) * pdata%q(p,i ,j,k) pdata%du(p,i ,j,k) = pdata%du(p,i ,j,k) - df pdata%du(p,ip1,j,k) = pdata%du(p,ip1,j,k) + df end do else do p = isl, isu df = dhi(1) * f(idn,i,j,k,1) * pdata%q(p,ip1,j,k) pdata%du(p,i ,j,k) = pdata%du(p,i ,j,k) - df pdata%du(p,ip1,j,k) = pdata%du(p,ip1,j,k) + df end do end if ! Y-face ! if (f(idn,i,j,k,2) >= 0.0d+00) then do p = isl, isu df = dhi(2) * f(idn,i,j,k,2) * pdata%q(p,i,j ,k) pdata%du(p,i,j ,k) = pdata%du(p,i,j ,k) - df pdata%du(p,i,jp1,k) = pdata%du(p,i,jp1,k) + df end do else do p = isl, isu df = dhi(2) * f(idn,i,j,k,2) * pdata%q(p,i,jp1,k) pdata%du(p,i,j ,k) = pdata%du(p,i,j ,k) - df pdata%du(p,i,jp1,k) = pdata%du(p,i,jp1,k) + df end do end if #if NDIMS == 3 ! Z-face ! if (f(idn,i,j,k,3) >= 0.0d+00) then do p = isl, isu df = dhi(3) * f(idn,i,j,k,3) * pdata%q(p,i,j,k ) pdata%du(p,i,j,k ) = pdata%du(p,i,j,k ) - df pdata%du(p,i,j,kp1) = pdata%du(p,i,j,kp1) + df end do else do p = isl, isu df = dhi(3) * f(idn,i,j,k,3) * pdata%q(p,i,j,kp1) pdata%du(p,i,j,k ) = pdata%du(p,i,j,k ) - df pdata%du(p,i,j,kp1) = pdata%du(p,i,j,kp1) + df end do end if #endif /* NDIMS == 3 */ end do ! i = nbl, neu end do ! j = nbl, neu #if NDIMS == 3 end do ! k = nbl, neu #endif /* NDIMS == 3 */ end if #ifdef PROFILE ! stop accounting time for the increment update ! call stop_timer(iui) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine update_increment ! !=============================================================================== ! ! subroutine UPDATE_VARIABLES: ! --------------------------- ! ! Subroutine iterates over all data blocks and converts the conservative ! variables to their primitive representation. ! ! Arguments: ! ! tm - time at the moment of update; ! dtm - time step since the last update; ! !=============================================================================== ! subroutine update_variables(tm, dtm) ! include external procedures ! use blocks , only : set_neighbors_update use boundaries , only : boundary_variables use equations , only : update_primitive_variables use equations , only : fix_unphysical_cells, correct_unphysical_states use shapes , only : update_shapes ! include external variables ! use blocks , only : block_meta, list_meta use blocks , only : block_data, list_data ! local variables are not implicit by default ! implicit none ! subroutine arguments ! real(kind=8), intent(in) :: tm, dtm ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! !------------------------------------------------------------------------------- ! #ifdef PROFILE ! start accounting time for variable update ! call start_timer(imv) #endif /* PROFILE */ ! update primitive variables in the changed blocks ! pdata => list_data do while (associated(pdata)) pmeta => pdata%meta if (pmeta%update) call update_primitive_variables(pdata%u, pdata%q) pdata => pdata%next end do ! update boundaries ! call boundary_variables(tm, dtm) ! apply shapes in blocks which need it ! pdata => list_data do while (associated(pdata)) pmeta => pdata%meta if (pmeta%update) call update_shapes(pdata, tm, dtm) pdata => pdata%next end do ! correct unphysical states if detected ! if (fix_unphysical_cells) then ! if an unphysical cell appeared in a block while updating its primitive ! variables it could be propagated to its neighbors through boundary update; ! mark all neighbors of such a block to be verified and corrected for ! unphysical cells too ! pmeta => list_meta do while (associated(pmeta)) if (pmeta%update) call set_neighbors_update(pmeta) pmeta => pmeta%next end do ! verify and correct, if necessary, unphysical cells in recently updated blocks ! pdata => list_data do while (associated(pdata)) pmeta => pdata%meta if (pmeta%update) & call correct_unphysical_states(step, pmeta%id, pdata%q, pdata%u) pdata => pdata%next end do end if #ifdef PROFILE ! stop accounting time for variable update ! call stop_timer(imv) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! end subroutine update_variables #ifdef DEBUG ! !=============================================================================== ! ! subroutine CHECK_VARIABLES: ! -------------------------- ! ! Subroutine iterates over all data blocks and converts the conservative ! variables to their primitive representation. ! ! !=============================================================================== ! subroutine check_variables() ! include external procedures ! use coordinates , only : nn => bcells use equations , only : nv, pvars, cvars use ieee_arithmetic, only : ieee_is_nan ! include external variables ! use blocks , only : block_meta, list_meta use blocks , only : block_data, list_data ! local variables are not implicit by default ! implicit none ! local variables ! integer :: i, j, k = 1, p ! local pointers ! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata ! !------------------------------------------------------------------------------- ! ! associate the pointer with the first block on the data block list ! pdata => list_data ! iterate over all data blocks ! do while (associated(pdata)) ! associate pmeta with the corresponding meta block ! pmeta => pdata%meta ! check if there are NaNs in primitive variables ! #if NDIMS == 3 do k = 1, nn #endif /* NDIMS == 3 */ do j = 1, nn do i = 1, nn do p = 1, nv if (ieee_is_nan(pdata%u(p,i,j,k))) then print *, 'U NaN:', cvars(p), pdata%meta%id, i, j, k end if if (ieee_is_nan(pdata%q(p,i,j,k))) then print *, 'Q NaN:', pvars(p), pdata%meta%id, i, j, k end if end do ! p = 1, nv end do ! i = 1, nn end do ! j = 1, nn #if NDIMS == 3 end do ! k = 1, nn #endif /* NDIMS == 3 */ ! assign pointer to the next block ! pdata => pdata%next end do !------------------------------------------------------------------------------- ! end subroutine check_variables #endif /* DEBUG */ !=============================================================================== ! end module evolution