hrweno_tvdode.f90

module hrweno_tvdode
!!   This module contains two total variation diminishing (TVD) high-order schemes for solving
!! initial value problems. It is very important to use TVD schemes for time integration. Even
!! with a TVD spacial discretization, if the time discretization is done by a non-TVD method,
!! the result may be oscillatory.
!!   Source: ICASE 97-65 by Shu, 1997.
   use hrweno_kinds, only: rk
   use stdlib_optval, only: optval
   implicit none
   private

   public :: rktvd, mstvd

   type, abstract :: tvdode
   !! Abstract class for TVD ODE solvers.
      procedure(integrand), pointer, nopass, private :: fu => null()
         !! subroutine with the derivative \( u'(t,u) \)
      integer :: neq
         !! number of equations
      integer :: order
         !! order of the method
      integer :: fevals = 0
         !! number of function evaluations
      integer :: istate = 0
         !! flag indicating the state of the integration:
         !! 1 first call for a problem,
         !! 2 subsequent call for a problem.
      character(:), allocatable:: msg
         !! error message
      real(rk), allocatable, private :: ui(:)
      real(rk), allocatable, private :: udot(:)
   contains
      procedure, pass(self) :: error_msg
   end type tvdode

   type, extends(tvdode) :: rktvd
   !! Runge-Kutta TVD ODE solver class.
   contains
      procedure, pass(self) :: integrate => rktvd_integrate
   end type rktvd

   type, extends(tvdode) :: mstvd
   !! Multi-step TVD ODE solver class.
      real(rk), allocatable, private :: uold(:, :)
      real(rk), allocatable, private :: udotold(:, :)
   contains
      procedure, pass(self) :: integrate => mstvd_integrate
   end type mstvd

   abstract interface
      subroutine integrand(t, u, udot)
      !!  Integrand for `tvdode` class
         import :: rk
         real(rk), intent(in) :: t, u(:)
         real(rk), intent(out) :: udot(:)
      end subroutine
   end interface

   interface rktvd
      module procedure :: rktvd_init
   end interface rktvd

   interface mstvd
      module procedure :: mstvd_init
   end interface mstvd

contains

   type(rktvd) function rktvd_init(fu, neq, order) result(self)
   !! Initialize `rktvd` object.
      procedure(integrand) :: fu
         !! subroutine with the derivative \( u'(t,u) \)
      integer, intent(in) :: neq
         !! number of equations
      integer, intent(in) :: order
         !! order of the method (1, 2 or 3)

      self%fu => fu

      if (neq > 0) then
         self%neq = neq
      else
         call self%error_msg("Invalid input 'neq'. Valid range: neq >= 1.")
      end if

      if ((order >= 1) .and. (order <= 3)) then
         self%order = order
      else
         call self%error_msg("Invalid input 'order' in 'rktvd'. Valid range: 1 <= k <= 3.")
      end if

      allocate (self%ui(self%neq), self%udot(self%neq))
      self%istate = 1

   end function rktvd_init

   subroutine rktvd_integrate(self, u, t, tout, dt, itask)
   !!   This subroutine implements the optimal 1st, 2nd and 3rd order TVD RK methods described
   !! in ICASE 97-65 (Shu, 1997).
   !!   The routine was built to work similarly to LSODE.
   !!
   !! @note
   !!   There are also 4th and 5th order methods, but they have lower CFL coefficients and
   !! are more difficult to implement. See Equation 4.15, page 44.
   !!
   !! @todo
   !! - Adjust dt in final step to avoid overshoting tout by some fraction of dt.
      class(rktvd), intent(inout) :: self
         !! object
      real(rk), intent(inout) :: u(:)
         !! vector(neq) with the variables to integrate \( u(t) \)
      real(rk), intent(inout) :: t
         !! time; on return it will be the current value of \(t\) (close to tout)
      real(rk), intent(in) :: tout
         !! time where next output is desired
      real(rk), intent(in) :: dt
         !! time step
      integer, intent(in), optional :: itask
         !! flag indicating the task to be performed:
         !! 1 normal integration until tout;
         !! 2 single `dt` step.

      integer :: itask_

      ! Check input conditions
      if (self%istate < 1) return
      if (is_done(t, tout, dt)) return
      itask_ = optval(itask, 1)

      ! Algorthm selection
      associate (ui => self%ui, udot => self%udot)
         select case (self%order)

            ! ------------------------------- 1st-order RK (Euler) ----------------------------
            ! Equation (4.10), page 43.
         case (1)
            do
               call self%fu(t, u, udot)
               u = u + dt*udot
               t = t + dt
               self%fevals = self%fevals + 1
               if (is_done(t, tout, dt) .or. itask_ == 2) exit
            end do

            ! --------------------------------- 2nd-order RK ----------------------------------
            ! Equation (4.10), page 43.
         case (2)
            do
               call self%fu(t, u, udot)
               ui = u + dt*udot
               call self%fu(t + dt, ui, udot)
               u = (u + ui + dt*udot)/2
               t = t + dt
               self%fevals = self%fevals + 2
               if (is_done(t, tout, dt) .or. itask_ == 2) exit
            end do

            ! --------------------------------- 3rd-order RK -------------------------------------
            ! Equation (4.11), page 43.
         case (3)
            do
               call self%fu(t, u, udot)
               ui = u + dt*udot
               call self%fu(t + dt, ui, udot)
               ui = (3*u + ui + dt*udot)/4
               call self%fu(t + dt/2, ui, udot)
               u = (u + 2*ui + 2*dt*udot)/3
               t = t + dt
               self%fevals = self%fevals + 3
               if (is_done(t, tout, dt) .or. itask_ == 2) exit
            end do

         end select
      end associate

      if (self%istate == 1) self%istate = 2

   end subroutine rktvd_integrate

   type(mstvd) function mstvd_init(fu, neq) result(self)
   !! Initialize `mstvd` object.
      procedure(integrand) :: fu
         !! subroutine with the derivative \( u'(t,u) \)
      integer, intent(in) :: neq
         !! number of equations

      self%fu => fu

      if (neq > 0) then
         self%neq = neq
      else
         call self%error_msg("Invalid input 'neq'. Valid range: neq >= 1.")
      end if

      self%order = 3

      allocate (self%ui(self%neq), self%udot(self%neq), &
                self%uold(self%neq, self%order + 1), self%udotold(self%neq, self%order + 1))
      self%istate = 1

   end function mstvd_init

   subroutine mstvd_integrate(self, u, t, tout, dt)
   !!   This subroutine implements a 5-step, 3rd order TVD multi-step method described
   !! in ICASE 97-65 (Shu, 1997). In theory, this method should have an efficiency 1.5 times
   !! higher than the RK method of the same order. However, in practice they appear to be
   !! almost identical.
   !!   The routine was built to work similarly to LSODE.
   !!
   !! @note
   !!   There is a 2nd order multi-step method, but the corresponding CFL value is half that
   !! of the 2nd order RK method. Thus, there is no reason to implement it.
      class(mstvd), intent(inout) :: self
         !! object
      real(rk), intent(inout) :: u(:)
         !! vector(neq) with the variables to integrate \( u(t) \)
      real(rk), intent(inout) :: t
         !! time; on return it will be the current value of \(t\) (close to tout)
      real(rk), intent(in) :: tout
         !! time where next output is desired
      real(rk), intent(in) :: dt
         !! time step

      type(rktvd) :: ode_start
      integer :: i

      ! Check input conditions
      if (self%istate < 1) return
      if (is_done(t, tout, dt)) return

      ! The first starting values must be computed with a single-step method: we chose
      ! the RK method of the same order.
      ! The factor 2 in 't+2*dt' is not important, it just needs to be larger than 1.0
      ! so that one full 'dt' step can be computed.
      associate (ui => self%ui, udot => self%udot, uold => self%uold, udotold => self%udotold)
      if (self%istate == 1) then

         ode_start = rktvd(self%fu, self%neq, self%order)

         do i = (self%order + 1), 1, -1
            uold(:, i) = u
            call self%fu(t, u, udotold(:, i))
            call ode_start%integrate(u, t, t + 2*dt, dt, itask=2)
         end do

         self%fevals = ode_start%fevals
         self%istate = 2

      end if

      ! Equation (4.26), page 48.
      do

         if (is_done(t, tout, dt)) exit

         call self%fu(t, u, udot)
         ui = (25*u + 50*dt*udot + 7*uold(:, 4) + 10*dt*udotold(:, 4))/32
         t = t + dt
         self%fevals = self%fevals + 1

         ! Shift u and udot values one step into the past
         udotold = eoshift(udotold, shift=-1, dim=2)
         uold = eoshift(uold, shift=-1, dim=2)
         udotold(:, 1) = udot
         uold(:, 1) = u
         u = ui

      end do
      end associate

   end subroutine mstvd_integrate

   pure logical function is_done(t, tout, dt)
    !! Aux function to check if the integration is finished.
      real(rk), intent(in) :: t
        !! current time
      real(rk), intent(in) :: tout
        !! time where next output is desired
      real(rk), intent(in) :: dt
        !! time step

      is_done = (t - tout)*sign(1._rk, dt) > 0._rk

   end function is_done

   pure subroutine error_msg(self, msg)
   !! Error method.
      class(tvdode), intent(inout) :: self
         !! object
      character(*), intent(in) :: msg
         !! message

      self%msg = msg
      self%istate = -1
      error stop self%msg

   end subroutine

end module hrweno_tvdode