dnsk Subroutine

subroutine dnsk(t, y, ydot, neq, res, psol, wt, rpar, ipar, savr, delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, s, confac, tolnew, muldel, maxit, ires, iersl, iernew)

Uses

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

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

for the unknown . The method used is a modified 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.

integer, intent(in) :: neq

Problem size.

procedure(res_t) :: res

User-defined residuals routine.

procedure(psol_t) :: psol

User-defined preconditioner routine.

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

Weights for error control.

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

User real workspace.

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

User integer workspace.

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(in) :: sqrtn

Square root of neq.

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

Reciprocal of square root of neq.

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

Tolerance for linear system solver.

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

Tolerance for convergence of the Newton iteration.

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

Convergence factor for the Newton iteration. s=rate/(1 - rate), where rate is the estimated rate of convergence of the Newton iteration. The closer rate is to 0, the faster the Newton iteration is converging. the closer rate is to 1, the slower the Newton iteration is converging.

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

Residual scale factor to improve convergence.

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

Tolerance on the norm of the Newton correction in the alternative Newton convergence test.

integer, intent(in) :: muldel

Flag indicating whether or not to multiply delta by confac.

integer, intent(in) :: maxit

Maximum number of Newton iterations allowed.

integer, intent(out) :: ires

Error flag from res. If ires < -1, then iernew = -1.

integer, intent(out) :: iersl

Error flag from the linear system solver. See dslvk for details. If iersl < 0, then iernew = -1. If iersl == 1, then iernew = 1.

integer, intent(out) :: iernew

Error flag for the Newton iteration. 0: the Newton iteration converged. 1: recoverable error inside Newton iteration. -1: unrecoverable error inside Newton iteration.


Calls

proc~~dnsk~~CallsGraph proc~dnsk dnsk ddwnrm ddwnrm proc~dnsk->ddwnrm dslvk dslvk proc~dnsk->dslvk

Variables

Type Visibility Attributes Name Initial
integer, public :: m
real(kind=rk), public :: deltanorm
real(kind=rk), public :: oldnorm
real(kind=rk), public :: rate
real(kind=rk), public :: rhok
logical, public :: converged

Source Code

subroutine dnsk( &
   t, y, ydot, neq, res, psol, wt, rpar, ipar, &
   savr, delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, &
   s, confac, tolnew, muldel, maxit, ires, iersl, iernew)
!! This routines 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.

   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(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(inout) :: wt(neq) ! @todo: confusing notation: wt ?= ewt ?= whgt
      !! Weights for error control.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   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(in) :: sqrtn
      !! Square root of `neq`.
   real(rk), intent(in) :: rsqrtn
      !! Reciprocal of square root of `neq`.
   real(rk), intent(in) :: epslin
      !! Tolerance for linear system solver.
   real(rk), intent(in) :: epscon
      !! Tolerance for convergence of the Newton iteration.
   real(rk), intent(out) :: s
      !! Convergence factor for the Newton iteration. `s=rate/(1 - rate)`, where
      !! `rate` is the estimated rate of convergence of the Newton iteration.
      !! The closer `rate` is to 0, the faster the Newton iteration is converging.
      !! the closer `rate` is to 1, the slower the Newton iteration is converging.
   real(rk), intent(in) :: confac
      !! Residual scale factor to improve convergence.
   real(rk), intent(in) :: tolnew
      !! Tolerance on the norm of the Newton correction in the alternative Newton
      !! convergence test.
   integer, intent(in) :: muldel ! @todo: convert to logical
      !! Flag indicating whether or not to multiply `delta` by `confac`.
   integer, intent(in) :: maxit
      !! Maximum number of Newton iterations allowed.
   integer, intent(out) :: ires
      !! Error flag from `res`. If `ires < -1`, then `iernew = -1`.
   integer, intent(out) :: iersl
      !! Error flag from the linear system solver. See [[dslvk]] for details.
      !! If `iersl < 0`, then `iernew = -1`.
      !! If `iersl == 1`, then `iernew = 1`.
   integer, intent(out) :: iernew
      !! Error flag for the Newton iteration.
      !!  `0`: the Newton iteration converged.
      !!  `1`: recoverable error inside Newton iteration.
      !! `-1`: unrecoverable error inside Newton iteration.

   integer :: m
   real(rk) :: deltanorm, oldnorm, rate, rhok
   real(rk), external :: ddwnrm ! @todo: remove this once inside module
   logical :: converged

   ! Initialize Newton counter M and accumulation vector E.
   m = 0
   e = zero

   ! Corrector loop.
   converged = .false.
   do
      iwm(iwloc_nni) = iwm(iwloc_nni) + 1

      ! If necessary, multiply residual by convergence factor.
      if (muldel == 1) then
         delta = delta*confac
      end if

      ! Save residual in SAVR.
      savr = delta

      ! Compute a new iterate. Store the correction in DELTA.
      call dslvk(neq, y, t, ydot, savr, delta, wt, rwm, iwm, &
                 res, ires, psol, iersl, cj, epslin, sqrtn, rsqrtn, rhok, &
                 rpar, ipar)
      if ((ires /= 0) .or. (iersl /= 0)) exit

      ! Update Y, E, and YDOT.
      y = y - delta
      e = e - delta
      ydot = ydot - cj*delta

      ! Test for convergence of the iteration.
      deltanorm = ddwnrm(neq, delta, wt, rpar, ipar)
      if (m == 0) then
         oldnorm = deltanorm
         if (deltanorm <= tolnew) then
            converged = .true.
            exit
         end if
      else
         rate = (deltanorm/oldnorm)**(one/m)
         if (rate > 0.9_rk) exit
         s = rate/(one - rate)
      end if

      if (s*deltanorm <= epscon) then
         converged = .true.
         exit
      end if

      ! The corrector has not yet converged. Update M and test whether
      ! the maximum number of iterations have been tried.
      m = m + 1
      if (m >= maxit) exit

      ! Evaluate the residual, and go back to do another iteration.
      ires = 0
      call res(t, y, ydot, cj, delta, ires, rpar, ipar)
      iwm(iwloc_nre) = iwm(iwloc_nre) + 1
      if (ires < 0) exit

   end do

   if (converged) then
      iernew = 0
   else
      if ((ires <= -2) .or. (iersl < 0)) then
         iernew = -1
      else
         iernew = 1
      end if
   end if

end subroutine dnsk