dnedk Subroutine

subroutine dnedk( &
   t, y, ydot, neq, res, jack, psol, &
   h, wt, jstart, idid, rpar, ipar, phi, gama, savr, delta, e, &
   rwm, iwm, cj, cjold, cjlast, s, uround, epli, sqrtn, rsqrtn, &
   epscon, jcalc, jflag, kp1, nonneg, ntype, iernls)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown \(y\). The method used is a matrix-free Newton scheme.

   use daskr_kinds, only: rk, zero, one
   use daskr
   implicit none

   real(rk), intent(in) :: t
      !! Independent variable.
   real(rk), intent(inout) :: y(neq)
      !! Solution vector.
   real(rk), intent(inout) :: ydot(neq)
      !! Derivative of solution vector after successful step.
   integer, intent(in) :: neq
      !! Problem size.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jack_t) :: jack
      !! User-defined Jacobian routine.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(in) :: h
      !! Step size.
   real(rk), intent(inout) :: wt(neq) ! @todo: confusing notation: wt ?= ewt ?= whgt
      !! Weights for error control.
   integer, intent(in) :: jstart ! @todo: ? convert to logical
      !! Flag indicating whether this is the first call to this routine.
      !! If `jstart = 0`, then this is the first call, otherwise it is not.
   integer, intent(out) :: idid
      !! Completion flag.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(in) :: phi(neq, kp1)
      !! Array of divided differences used by the Newton iteration.
   real(rk), intent(in) :: gama(kp1)
      !! Array used to predict `y` and `ydot`.
   real(rk), intent(out) :: savr(neq)
      !! Saved residual vector.
   real(rk), intent(inout) :: delta(neq)
      !! Real workspace.
   real(rk), intent(out) :: e(neq)
      !! Error accumulation vector.
   real(rk), intent(inout) :: rwm(*)
      !! Real workspace for the linear system solver.
   integer, intent(inout) :: iwm(*)
      !! Integer workspace for the linear system solver.
   real(rk), intent(in) :: cj
      !! Scalar used in forming the system Jacobian.
   real(rk), intent(out) :: cjold
      !! Value of `cj` as of the last call to [[ditmd]]. Accounts for changes in `cj`
      !! needed to decide whether to call [[ditmd]].
   real(rk), intent(in) :: cjlast
      !! Previous value of `cj`.
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration. See [[dnsk]] for details. On the
      !! first Newton iteration with an updated preconditioner, `s = 100`. The value of
      !! `s` is preserved from call to call so that the rate estimate from a previous
      !! step can be applied to the current step.
   real(rk), intent(in) :: uround ! @todo: ? remove
      !! Unit roundoff (not used).
   real(rk), intent(in) :: epli ! @note: not the same as 'epslin'!
      !! Convergence test constant for the iterative solution of linear systems. See
      !! [[daskr]] for details.
   real(rk), intent(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   integer, intent(out) :: jcalc
      !! Flag indicating whether the Jacobian matrix needs to be updated.
      !! `-1`: call [[ditmd]] to update the Jacobian matrix.
      !!  `0`: Jacobian matrix is up-to-date.
      !!  `1`: Jacobian matrix is outdated, but [[ditmd]] will not be called unless `jcalc` is set to `-1`.
   integer, intent(in) :: jflag
      !! Flag indicating whether a `jac` routine is supplied by the user.
   integer, intent(in) :: kp1
      !! The current order + 1. Updated across calls.
   integer, intent(in) :: nonneg
      !! Flag indicating whether nonnegativity constraints are imposed on the solution.
      !! `0`: no nonnegativity constraints.
      !! `1`: nonnegativity constraints are imposed.
   integer, intent(in) :: ntype
      !! Type of the nonlinear solver.
      !! `1`: modified Newton with iterative linear solver.
      !! `2`: modified Newton with user-supplied linear solver.
   integer, intent(out) :: iernls
      !! Error flag for the nonlinear solver.
      !!  `0`: nonlinear solver converged.
      !!  `1`: recoverable error inside non-linear solver.
      !! `-1`: unrecoverable error inside non-linear solver.

   integer, parameter :: muldel = 0, maxit = 4
   real(rk), parameter :: xrate = 0.25_rk

   integer :: iernew, ierpj, iersl, iertyp, ires, j, liwp, lrwp
   real(rk) :: deltanorm, epslin, temp1, temp2, tolnew
   real(rk), external :: ddwnrm ! @todo: remove this once inside module

   ! Verify that this is the correct subroutine.
   iertyp = 0
   if (ntype /= 1) then
      iertyp = 1
      goto 380
   end if

   ! If this is the first step, perform initializations.
   if (jstart == 0) then
      cjold = cj
      jcalc = -1
      s = 1e2_rk
   end if

   ! Perform all other initializations.
   iernls = 0
   lrwp = iwm(iwloc_lrwp)
   liwp = iwm(iwloc_liwp)

   ! Decide whether to update the preconditioner.
   if (jflag /= 0) then
      temp1 = (one - xrate)/(one + xrate)
      temp2 = 1/temp1
      if ((cj/cjold < temp1) .or. (cj/cjold > temp2)) jcalc = -1
      if (cj /= cjlast) s = 1e2_rk
   else
      jcalc = 0
   end if

   ! Looping point for updating preconditioner with current stepsize.
300 continue

   ! Initialize all error flags to zero.
   ierpj = 0
   ires = 0
   iersl = 0
   iernew = 0

   ! Predict the solution and derivative and compute the tolerance for the Newton iteration.
   y = phi(:, 1)
   ydot = zero
   do j = 2, kp1
      y = y + phi(:, j)
      ydot = ydot + gama(j)*phi(:, j)
   end do

330 continue
   epslin = epli*epscon
   tolnew = epslin

   ! Call RES to initialize DELTA.
   call res(t, y, ydot, cj, delta, ires, rpar, ipar)
   iwm(iwloc_nre) = iwm(iwloc_nre) + 1
   if (ires < 0) goto 380

   ! If indicated, update the preconditioner.
   ! Set JCALC to 0 as an indicator that this has been done.
   if (jcalc == -1) then
      jcalc = 0
      call jack(res, ires, neq, t, y, ydot, wt, delta, e, h, cj, rwm(lrwp), iwm(liwp), ierpj, rpar, ipar)
      iwm(iwloc_nje) = iwm(iwloc_nje) + 1
      cjold = cj
      s = 1e2_rk
      if (ires < 0) goto 380
      if (ierpj /= 0) goto 380
   end if

   ! Call the nonlinear Newton solver.
   call dnsk(t, y, ydot, neq, res, psol, wt, rpar, ipar, savr, &
             delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, &
             s, temp1, tolnew, muldel, maxit, ires, iersl, iernew)

   if (iernew > 0 .and. jcalc /= 0) then
      ! The Newton iteration had a recoverable failure with an old
      ! preconditioner. Retry the step with a new preconditioner.
      jcalc = -1
      goto 300
   end if

   if (iernew /= 0) goto 380

   ! The Newton iteration has converged. If nonnegativity of solution is required, set
   ! the solution nonnegative, if the perturbation to do it is small enough. If the
   ! change is too large, then consider the corrector iteration to have failed.
   if (nonneg == 0) goto 390
   delta = min(y, zero)
   deltanorm = ddwnrm(neq, delta, wt, rpar, ipar)
   if (deltanorm > epscon) goto 380
   e = e - delta
   goto 390

   ! Exits from nonlinear solver.
   ! No convergence with current preconditioner.
   ! Compute IERNLS and IDID accordingly.
380 continue
   if ((ires <= -2) .or. (iersl < 0) .or. (iertyp /= 0)) then
      iernls = -1
      if (ires <= -2) idid = -11
      if (iersl < 0) idid = -13
      if (iertyp /= 0) idid = -15
   else
      iernls = 1
      if (ires == -1) idid = -10
      if (ierpj /= 0) idid = -5
      if (iersl > 0) idid = -14
   end if

390 continue
   jcalc = 1

end subroutine dnedk