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)

Uses

  • proc~~dnedk~~UsesGraph proc~dnedk dnedk module~daskr daskr proc~dnedk->module~daskr module~daskr_kinds daskr_kinds proc~dnedk->module~daskr_kinds module~daskr->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

This routine solves a nonlinear system of algebraic equations of the form:

for the unknown . The method used is a matrix-free Newton scheme.

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: t

Independent variable.

real(kind=rk), intent(inout) :: y(neq)

Solution vector.

real(kind=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(kind=rk), intent(in) :: h

Step size.

real(kind=rk), intent(inout) :: wt(neq)

Weights for error control.

integer, intent(in) :: jstart

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(kind=rk), intent(inout) :: rpar(*)

User real workspace.

integer, intent(inout) :: ipar(*)

User integer workspace.

real(kind=rk), intent(in) :: phi(neq,kp1)

Array of divided differences used by the Newton iteration.

real(kind=rk), intent(in) :: gama(kp1)

Array used to predict y and ydot.

real(kind=rk), intent(out) :: savr(neq)

Saved residual vector.

real(kind=rk), intent(inout) :: delta(neq)

Real workspace.

real(kind=rk), intent(out) :: e(neq)

Error accumulation vector.

real(kind=rk), intent(inout) :: rwm(*)

Real workspace for the linear system solver.

integer, intent(inout) :: iwm(*)

Integer workspace for the linear system solver.

real(kind=rk), intent(in) :: cj

Scalar used in forming the system Jacobian.

real(kind=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(kind=rk), intent(in) :: cjlast

Previous value of cj.

real(kind=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(kind=rk), intent(in) :: uround

Unit roundoff (not used).

real(kind=rk), intent(in) :: epli

Convergence test constant for the iterative solution of linear systems. See daskr for details.

real(kind=rk), intent(in) :: sqrtn

Square root of neq.

real(kind=rk), intent(in) :: rsqrtn

Reciprocal of square root of neq.

real(kind=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.


Calls

proc~~dnedk~~CallsGraph proc~dnedk dnedk ddwnrm ddwnrm proc~dnedk->ddwnrm dnsk dnsk proc~dnedk->dnsk

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: muldel = 0
integer, public, parameter :: maxit = 4
real(kind=rk), public, parameter :: xrate = 0.25_rk
integer, public :: iernew
integer, public :: ierpj
integer, public :: iersl
integer, public :: iertyp
integer, public :: ires
integer, public :: j
integer, public :: liwp
integer, public :: lrwp
real(kind=rk), public :: deltanorm
real(kind=rk), public :: epslin
real(kind=rk), public :: temp1
real(kind=rk), public :: temp2
real(kind=rk), public :: tolnew

Source Code

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