dnsd Subroutine

subroutine dnsd(t, y, ydot, neq, res, dum1, wt, rpar, ipar, dum2, delta, e, rwm, iwm, cj, dum3, dum4, dum5, epscon, s, confac, tolnew, muldel, maxit, ires, dum6, iernew)

Uses

  • proc~~dnsd~~UsesGraph proc~dnsd dnsd module~daskr daskr proc~dnsd->module~daskr module~daskr_kinds daskr_kinds proc~dnsd->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(psol_t) :: dum1

Dummy argument.

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

Weights for error criterion.

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

User real workspace.

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

User integer workspace.

real(kind=rk), intent(in) :: dum2(neq)

Dummy argument.

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

Real work array.

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) :: 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.

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

Convergence factor for the Newton iteration. s=rate/(1 - rate), where rate is the estimated rate of convergence of the Newton iteration.

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. If ires < -1, then iernew = -1.

integer, intent(in) :: dum6

Dummy argument.

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~~dnsd~~CallsGraph proc~dnsd dnsd ddwnrm ddwnrm proc~dnsd->ddwnrm dslvd dslvd proc~dnsd->dslvd

Variables

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

Source Code

subroutine dnsd( &
   t, y, ydot, neq, res, dum1, wt, rpar, ipar, &
   dum2, delta, e, rwm, iwm, cj, dum3, dum4, dum5, epscon, &
   s, confac, tolnew, muldel, maxit, ires, dum6, iernew)
!! 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(psol_t) :: dum1
      !! Dummy argument.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error criterion.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(in) :: dum2(neq)
      !! Dummy argument.
   real(rk), intent(inout) :: delta(neq)
      !! Real work array.
   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) :: 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.
   real(rk), intent(inout) :: s
      !! Convergence factor for the Newton iteration. `s=rate/(1 - rate)`, where
      !! `rate` is the estimated rate of convergence of the Newton iteration.
   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`.
      !! If `ires < -1`, then `iernew = -1`.
   integer, intent(in) :: dum6
      !! Dummy argument.
   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
   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 .eq. 1) then
         delta = delta*confac
      end if

      ! Compute a new iterate (back-substitution).
      ! Store the correction in DELTA.
      call dslvd(neq, delta, rwm, iwm)

      ! 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 .eq. 0) then
         oldnorm = deltanorm
         if (deltanorm .le. tolnew) then
            converged = .true.
            exit
         end if
      else
         rate = (deltanorm/oldnorm)**(one/m)
         if (rate .gt. 0.9_rk) exit
         s = rate/(1 - rate)
      end if
      if (s*deltanorm .le. 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 .ge. 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 .lt. 0) exit

   end do

   if (converged) then
      iernew = 0
   else
      if (ires .le. -2) then
         iernew = -1
      else
         iernew = 1
      end if
   end if

end subroutine dnsd