quasilinear_solve Subroutine

  subroutine quasilinear_solve(L, f, jac_prod, solution, order, resid_norm, &
                               flag, info, tol, precond, differentiate,     &
                               iter_max, gmres_iter_max, krylov_dim)
    !* Author: Chris MacMackin
    !  Date: March 2017
    !
    ! This is an iterative method to solve nonlinear systems of
    ! ODEs. Consider $$ L^{(n)}\vec{u}(x) = \vec{f}(\vec{u}(x),
    ! \vec{u}^{(1)}(x), \ldots, \vec{u}^{n-1}(x), x), \qquad \vec{u}
    ! \in \mathbb{R}^{m}. $$ For a domain from \( (0,b) \), boundary
    ! conditions are specified by $$ g_k(\vec{u}(0) = 0,
    ! \vec{u}^{(1)}(0), \ldots, \vec{u}^{n-1}(0)), \qquad k = 1,
    ! \ldots, l $$ $$ g_k(\vec{u}(b), \vec{u}^{(1)}(b), \ldots,
    ! \vec{u}^{n-1}(b)) = 0, \qquad k = l+1, \ldots, mn. $$ Here,
    ! \(L^{(n)}\) is an \(n\)th order ordinary differential operator,
    ! and \(\vec{f}, g_1, g_2, \ldots, g_{nm}\) are nonlinear
    ! functions of \(\vec{u}(x)\) and its first \(n-1\) derivatives.
    ! 
    ! With quasilinearisation method, we iterate according to $$
    ! L^{(n)}\vec{u}_{r+1}(x) = \vec{f}(\vec{u}_r(x),
    ! \vec{u}^{(1)}_r(x), \ldots, \vec{u}^{n-1}_r(x), x) + \\
    ! \qquad\sum_{s=0}^{n-1}\vec{f}_{\vec{u}^{(s)}}(\vec{u}_r(x),
    ! \vec{u}^{(1)}_r(x), \ldots, \vec{u}^{n-1}_r(x),
    ! x)(\vec{u}^{(s)}_{r+1}(x) - \vec{u}^{(s)}_r(x)). $$ The boundary
    ! conditions are iterated according to $$ \sum_{s=0}^{n-1}
    ! g_{k,\vec{u}^{(s)}}(\vec{u}_r(0), \vec{u}^{(1)}_r(0), \ldots,
    ! \vec{u}^{n-1}_r(0), x)\cdot(\vec{u}^{(s)}_{r+1}(0) -
    ! \vec{u}^{(s)}_r(0)) = 0, k=1,\ldots,l$$ and $$ \sum_{s=0}^{n-1}
    ! g_{k,\vec{u}^{(s)}}(\vec{u}_r(b), \vec{u}^{(1)}_r(b), \ldots,
    ! \vec{u}^{n-1}_r(b), x)\cdot(\vec{u}^{(s)}_{r+1}(b) -
    ! \vec{u}^{(s)}_r(b)) = 0, k=l+1,\ldots,mn.$$ Here,
    ! \(\vec{f}_{\vec{u}^{(s)}}\) is the derivative of \(\vec{f}\)
    ! with respect to \(\vec{u}^{(s)}\) and is a tensor, while
    ! \(g_{k,\vec{u}^{(s)}}\) is the derivative of \(g_k\) with
    ! respect to \(\vec{u}^{(s)}\) and is a vector.
    ! 
    ! The user must provide functions for \(L\) and
    ! \(\vec{f}\). Currently this implementation only handles linear
    ! boundary conditions, which remain the same between
    ! iterations. The boundary conditions should be handled in the
    ! functions passed as arguments. In most situations, a
    ! preconditioner will also be needed. A simple one is to
    ! approximate \(L^{-1}\). The derivatives of \(\vec{f}\) and
    ! \(g_k\) are estimated using a finite-difference.
    !
    !####Output parameters
    !
    ! On output, the components of the `info` argument are as follows:
    !
    !     info(1)   = nLe   (number of evaluations of `L`)
    !     info(2)   = nfe   (number of evaluations of `f`)
    !     info(3)   = nrpre (number of preconditioner evaluations)
    !     info(4)   = nli   (number of linear iterations)
    !     info(5)   = nni   (number of nonlinear iterations)
    !
    procedure(L_intr)                            :: L
      !! A function providing the linear, left-hand-side of the ODE
      !! being solved.
    procedure(f_intr)                            :: f
      !! A function providing the nonlinear, right-hand-side of the
      !! ODE being solved.
    procedure(jac_intr)                          :: jac_prod
      !! A function providing the product of the Jacobian of the
      !! nonlinear, right-hand-side of the ODE being solved and
      !! another vector.
    real(r8), dimension(:), intent(inout)        :: solution
      !! On input, an estimate of the solution to the ODE. On output,
      !! the actual solution.
    integer, intent(in)                          :: order
      !! The order of the derivative taken by `L`
    real(r8), intent(out)                        :: resid_norm
      !! Norm of the residual of the final solution.
    integer, intent(out)                         :: flag
      !! Status flag indicating whether the iterations ended succesfully.
      !!
      !!< 0:
      !!:    `|flag|` is [[gmres_solve]] return code
      !!
      !!0
      !!:    Normal termination: acceptable solution found
      !!
      !!1
      !!:    Convergence, but residual greater than required tolerance
      !!
      !!2
      !!:    Solution did not converge within `iter_max` iterations
      !!
      !!3
      !!:    Solution began to diverge
      !!
      !!4
      !!:    No `diff` procedure provided when `order > 1`
      !!
    integer, dimension(5), intent(out), optional :: info
      !! Array containing various outputs; see above
    real(r8), intent(in), optional               :: tol
      !! The required reduction in the solution residual. Default is
      !! `size(solution) * 1e-8`.
    procedure(pre_intr), optional                :: precond
      !! A right-preconditioner which may be used to improve
      !! convergence of the solution.
    procedure(diff_intr), optional               :: differentiate
      !! A procedure which will evaluate the `n`th derivative of the
      !! state vector, when `n` is less than `order`.
    integer, intent(in), optional                :: iter_max
      !! Maximum allowable number of quasilinearised
      !! iterations. Default is 15.
    integer, intent(in), optional                :: gmres_iter_max
      !! Maximum allowable number of GMRES iterations. Default is
      !! 1000.
    integer, intent(in), optional                :: krylov_dim
      !! Maximum Krylov subspace dimension; default 10. Larger values
      !! will allow for faster convergence (and in some cases be the
      !! difference between whether or not convergence is possible),
      !! but require more memory.

    integer :: npoints, itmax, gitmax, kdim
    real(r8) :: eta, gmres_eta

    integer :: i, stagnant_iters, gmres_flag
    integer :: nlhs, nrpre, nli, tnlhs, tnrpre, tnli
    real(r8) :: init_resid, old_resid, gmres_norm
    real(r8), dimension(size(solution),order) :: u_prev
    real(r8), dimension(size(solution))       :: f_prev, rhs

    tnlhs  = 0
    tnrpre = 0
    tnli   = 0

    if (.not. present(differentiate) .and. order > 1) then
      flag = 3
      call logger%error('quasilinear_solve','Did not provide routine '//&
           'to get derivative values.')
      return
    end if

    npoints = size(solution)
    if (present(tol)) then
      eta = tol
    else
      eta = 1.e-8_r8 * npoints
    end if
    if (present(iter_max)) then
      itmax = iter_max
    else
      itmax = 15
    end if
    if (present(gmres_iter_max)) then
      gitmax = gmres_iter_max
    else
      gitmax = 1000
    end if
    if (present(krylov_dim)) then
      kdim = krylov_dim
    else
      kdim = 10
    end if

    i = 0
    stagnant_iters = 0

    u_prev = get_derivs(solution)
    f_prev = f(u_prev)
    resid_norm = dnrm2(npoints, L(solution) - f_prev, 1)
    init_resid = resid_norm
    old_resid = resid_norm * 1e3_r8

    do while(resid_norm > eta)
      call logger%trivia('quasilinear_solve','Nonlinear iteration '//str(i)//&
           ', with '//str(tnli)//'linear iterations, and residual '//str(resid_norm))
      i = i + 1
      if (abs(old_resid - resid_norm)/resid_norm < 1e-2_r8) then
        stagnant_iters = stagnant_iters + 1
      else
        stagnant_iters = 0
      end if
      if (stagnant_iters > 3) then
        flag = 1
        return
      end if
      if (i > itmax) then
        flag = 2
        return
      end if
      if (20*old_resid < resid_norm) then
        flag = 3
        return
      end if

      rhs = f_prev - jac_prod(u_prev, u_prev)
      gmres_eta = max(min(eta*10._r8**min(i+2,6),1e-4_r8),1e-10_r8)
      gmres_eta = gmres_eta * 10._r8**(-2*stagnant_iters)

      call gmres_solve(solution, lin_op, rhs, gmres_norm, gmres_flag, &
                       nlhs, nrpre, nli, gmres_eta, preconditioner,   &
                       resid_update=0, iter_max=gitmax, krylov_dim=kdim)

      tnlhs  = tnlhs  + nlhs
      tnrpre = tnrpre + nrpre
      tnli   = tnli   + nli
      if(gmres_flag > 0) call logger%warning('quasilinear_solve', &
           'GMRES solver returned with flag '//str(gmres_flag))

      u_prev = get_derivs(solution)
      f_prev = f(u_prev)
      old_resid = resid_norm
      resid_norm = dnrm2(npoints, L(solution) - f_prev, 1)


      if (gmres_flag /= 0 .and. (resid_norm > 20*old_resid .or. any(isnan(solution)))) then
        if (present(info)) then
          info(1) = i + tnlhs + tnrpre
          info(2) = 2*i + tnlhs
          info(3) = tnrpre
          info(4) = tnli
          info(5) = i
        end if
        if (init_resid < resid_norm) then
          flag = 3
        else
          flag = -gmres_flag
        end if
        return
      end if
    end do

    if (present(info)) then
      info(1) = 1 + i + tnlhs + tnrpre
      info(2) = 1 + 2*i + tnlhs
      info(3) = tnrpre
      info(4) = tnli
      info(5) = i
    end if
    flag = 0
    call logger%trivia('quasilinear_solve','Nonlinear iteration '//str(i)//&
         ', with '//str(tnli)//' linear iterations, and residual '//str(resid_norm))

  contains

    function get_derivs(v)
      !! Calculates the necessary number of derivatives and assembles
      !! them in 2-D array.
      real(r8), dimension(:), intent(in) :: v
      real(r8), dimension(size(v), order) :: get_derivs
      integer :: j
      get_derivs(:,1) = v
      do j = 2, order
        get_derivs(:,j) = differentiate(v, i-1)
      end do
    end function get_derivs

    function lin_op(v, xcur, rhs, rpar, ipar, success)
      !! The linear operator for the quasilinearised system.
      real(r8), dimension(:), intent(in)    :: v
        !! The vector to be operated upon
      real(r8), dimension(:), intent(in)    :: xcur
        !! Array containing the current estimate of the independent
        !! variables in the linear system. This may not be needed, but
        !! is provided just in case.
      real(r8), dimension(:), intent(in)    :: rhs
        !! Array containing the right hand side of the linear
        !! system. This may not be needed, but is provided just in
        !! case.
      real(r8), dimension(*), intent(inout) :: rpar
        !! Parameter/work array 
      integer, dimension(*), intent(inout)  :: ipar
        !! Parameter/work array
      logical, intent(out)                  :: success
        !! Indicates whether operation was completed succesfully
      real(r8), dimension(size(xcur))       :: lin_op
        !! Result of the operation
      real(r8), dimension(size(xcur),order) :: v_derivs
      v_derivs = get_derivs(v)
      lin_op = L(v) - jac_prod(u_prev, v_derivs)
      success = .true.
    end function lin_op

    function preconditioner(v, xcur, rhs, rpar, ipar, success)
      !! The preconditioner for the quasilinearised system.
      real(r8), dimension(:), intent(in)    :: v
        !! The vector to be preconditioned
      real(r8), dimension(:), intent(in)    :: xcur
        !! Array containing the current estimate of the independent
        !! variables in the linear system. This may not be needed, but
        !! is provided just in case.
      real(r8), dimension(:), intent(in)    :: rhs
        !! Array containing the right hand side of the linear
        !! system. This may not be needed, but is provided just in
        !! case.
      real(r8), dimension(*), intent(inout) :: rpar
        !! Parameter/work array 
      integer, dimension(*), intent(inout)  :: ipar
        !! Parameter/work array
      logical, intent(out)                  :: success
        !! Indicates whether operation was completed succesfully
      real(r8), dimension(size(xcur))       :: preconditioner
        !! Result of the operation
      preconditioner = precond(v, reshape(xcur, [size(xcur), 1]), L, f, f_prev, rhs)
      success = .true.
    end function preconditioner

  end subroutine quasilinear_solve