dnedd Subroutine

subroutine dnedd(t, y, ydot, neq, res, jac, dum1, h, wt, jstart, idid, rpar, ipar, phi, gama, dum2, delta, e, rwm, iwm, cj, cjold, cjlast, s, uround, dum3, dum4, dum5, epscon, jcalc, dum6, kp1, nonneg, ntype, iernls)

Uses

  • proc~~dnedd~~UsesGraph proc~dnedd dnedd module~daskr daskr proc~dnedd->module~daskr module~daskr_kinds daskr_kinds proc~dnedd->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 modified Newton scheme.

All arguments with "dum" in their names are dummy arguments which are not used in this routine.

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.

integer, intent(in) :: neq

Problem size.

procedure(res_t) :: res

User-defined residuals routine.

procedure(jac_t) :: jac

User-defined Jacobian routine.

procedure(psol_t) :: dum1

Dummy argument.

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(in) :: dum2(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 dmatd. Accounts for changes in cj needed to decide whether to call dmatd.

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

Previous value of cj.

real(kind=rk), intent(out) :: s

Convergence factor for the Newton iteration. See dnsd and 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.

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

Dummy argument.

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

Dummy argument.

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

Dummy argument.

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 dmatd to update the Jacobian matrix. 0: Jacobian matrix is up-to-date. 1: Jacobian matrix is outdated, but dmatd will not be called unless jcalc is set to -1.

integer, intent(in) :: dum6

Dummy argument.

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. 0: modified Newton with direct 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~~dnedd~~CallsGraph proc~dnedd dnedd ddwnrm ddwnrm proc~dnedd->ddwnrm dmatd dmatd proc~dnedd->dmatd dnsd dnsd proc~dnedd->dnsd

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: muldel = 1
integer, public, parameter :: maxit = 4
real(kind=rk), public, parameter :: xrate = 0.25_rk
integer, public :: idum
integer, public :: ierj
integer, public :: iernew
integer, public :: iertyp
integer, public :: ires
integer, public :: j
real(kind=rk), public :: delnorm
real(kind=rk), public :: pnorm
real(kind=rk), public :: temp1
real(kind=rk), public :: temp2
real(kind=rk), public :: tolnew

Source Code

subroutine dnedd( &
   t, y, ydot, neq, res, jac, dum1, &
   h, wt, jstart, idid, rpar, ipar, phi, gama, dum2, delta, e, &
   rwm, iwm, cj, cjold, cjlast, s, uround, dum3, dum4, dum5, &
   epscon, jcalc, dum6, 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 modified Newton scheme.
!!
!! All arguments with "dum" in their names are dummy arguments which are not used in
!! this routine.

   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.
   integer, intent(in) :: neq
      !! Problem size.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jac_t) :: jac
      !! User-defined Jacobian routine.
   procedure(psol_t) :: dum1
      !! Dummy argument.
   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(in) :: dum2(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 [[dmatd]]. Accounts for changes in `cj`
      !! needed to decide whether to call [[dmatd]].
   real(rk), intent(in) :: cjlast
      !! Previous value of `cj`.
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration. See [[dnsd]] and [[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.
   real(rk), intent(in) :: dum3
      !! Dummy argument.
   real(rk), intent(in) :: dum4
      !! Dummy argument.
   real(rk), intent(in) :: dum5
      !! Dummy argument.
   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 [[dmatd]] to update the Jacobian matrix.
      !!  `0`: Jacobian matrix is up-to-date.
      !!  `1`: Jacobian matrix is outdated, but [[dmatd]] will not be called unless `jcalc` is set to `-1`.
   integer, intent(in) :: dum6
      !! Dummy argument.
   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.
      !! `0`: modified Newton with direct 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 = 1, maxit = 4
   real(rk), parameter :: xrate = 0.25_rk
   integer :: idum, ierj, iernew, iertyp, ires, j
   real(rk) :: delnorm, pnorm, temp1, temp2, tolnew
   real(rk), external :: ddwnrm ! @todo: remove this once inside module

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

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

   ! Perform all other initializations.
   iernls = 0

   ! Decide whether new Jacobian is needed.
   temp1 = (one - xrate)/(one + xrate)
   temp2 = 1/temp1
   if ((cj/cjold .lt. temp1) .or. (cj/cjold .gt. temp2)) jcalc = -1
   if (cj .ne. cjlast) s = 1e2_rk

   ! Entry point for updating the Jacobian with current stepsize.
300 continue

   ! Initialize all error flags to zero.
   ierj = 0
   ires = 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

   pnorm = ddwnrm(neq, y, wt, rpar, ipar)
   tolnew = 100*uround*pnorm

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

   ! If indicated, reevaluate the iteration matrix
   ! J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
   ! Set JCALC to 0 as an indicator that this has been done.
   if (jcalc .eq. -1) then
      iwm(iwloc_nje) = iwm(iwloc_nje) + 1
      jcalc = 0
      call dmatd(neq, t, y, ydot, delta, cj, h, ierj, wt, e, rwm, iwm, res, ires, uround, jac, rpar, ipar)
      cjold = cj
      s = 1e2_rk
      if (ires .lt. 0) goto 380
      if (ierj .ne. 0) goto 380
   end if

   ! Call the nonlinear Newton solver.
   temp1 = 2/(one + cj/cjold)
   call dnsd(t, y, ydot, neq, res, dum1, wt, rpar, ipar, dum2, &
             delta, e, rwm, iwm, cj, dum4, dum5, dum3, epscon, s, temp1, &
             tolnew, muldel, maxit, ires, idum, iernew)

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

   if (iernew .ne. 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 .eq. 0) goto 390
   delta = min(y, zero)
   delnorm = ddwnrm(neq, delta, wt, rpar, ipar)
   if (delnorm .gt. epscon) goto 380
   e = e - delta

   ! Exits from nonlinear solver.
   ! No convergence with current iteration matrix, or singular iteration matrix.
   ! Compute IERNLS and IDID accordingly.
380 continue
   if (ires .le. -2 .or. iertyp .ne. 0) then
      iernls = -1
      if (ires .le. -2) idid = -11
      if (iertyp .ne. 0) idid = -15
   else
      iernls = 1
      if (ires .lt. 0) idid = -10
      if (ierj .ne. 0) idid = -8
   end if

390 continue
   jcalc = 1

end subroutine dnedd