ddasik Subroutine

subroutine ddasik(t, y, ydot, neq, icopt, idy, res, jack, psol, h, tscale, wt, jskip, rpar, ipar, savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, uround, epli, sqrtn, rsqrtn, epscon, ratemax, stptol, jflg, icnflg, icnstr, iernls)

Uses

  • proc~~ddasik~~UsesGraph proc~ddasik ddasik module~daskr daskr proc~ddasik->module~daskr module~daskr_kinds daskr_kinds proc~ddasik->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 parts of and in the initial conditions. An initial value for and initial guess for are input. The method used is a Newton scheme with Krylov iteration and a linesearch algorithm.

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.

integer, intent(in) :: icopt

Initial condition flag.

integer, intent(in) :: idy(neq)

Array indicating which variables are differential and which are algebraic.

procedure(res_t) :: res

User-defined residuals routine.

procedure(jack_t) :: jack

User-supplied routine to update the preconditioner.

procedure(psol_t) :: psol

User-defined preconditioner routine.

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

Scaling factor for this initial condition calculation.

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

Scale factor in t; used for stopping tests if nonzero.

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

Weights for error criterion.

integer, intent(inout) :: jskip

Flag to signal if initial jack call is to be skipped. 1: skip the call, 0: do not skip call.

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

User real workspace.

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

User integer workspace.

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

Real work array.

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

Real work array.

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

Real work array.

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

Real work array.

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

Real work array.

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

Real work array.

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

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

Maximum convergence rate for which Newton iteration is considered converging.

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

Tolerance used in calculating the minimum lambda (rl) value allowed.

integer, intent(in) :: jflg

Flag indicating whether a jac routine is supplied by the user.

integer, intent(in) :: icnflg

Constraint flag. If nonzero, then constraint violations in the proposed new approximate solution will be checked for, and the maximum step length will be adjusted accordingly.

integer, intent(out) :: icnstr(neq)

Flags for checking constraints.

integer, intent(out) :: iernls

Error flag for nonlinear solver. 0: nonlinear solver converged. 1, 2: recoverable error inside nonlinear solver. 1: retry with current (y,y'). 2: retry with original (y,y'). -1: unrecoverable error in nonlinear solver.


Calls

proc~~ddasik~~CallsGraph proc~ddasik ddasik dnsik dnsik proc~ddasik->dnsik

Variables

Type Visibility Attributes Name Initial
integer, public :: iernew
integer, public :: ierpj
integer, public :: ires
integer, public :: liwp
integer, public :: lrwp
integer, public :: maxnit
integer, public :: maxnj
integer, public :: nj
real(kind=rk), public :: epslin
logical, public :: failed

Source Code

subroutine ddasik( &
   t, y, ydot, neq, icopt, idy, res, jack, psol, h, tscale, &
   wt, jskip, rpar, ipar, savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, uround, &
   epli, sqrtn, rsqrtn, epscon, ratemax, stptol, jflg, icnflg, icnstr, iernls)
!! This routine solves a nonlinear system of algebraic equations of the form:
!!
!!  $$ G(t, y, \dot{y}) = 0 $$
!!
!! for the unknown parts of \(y\) and \(\dot{y}\) in the initial conditions. An initial
!! value for \(y\) and initial guess for \(\dot{y}\) are input. The method used is a
!! Newton scheme with Krylov iteration and a linesearch algorithm.

   use daskr_kinds, only: rk
   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.
   integer, intent(in) :: icopt
      !! Initial condition flag.
   integer, intent(in) :: idy(neq)
      !! Array indicating which variables are differential and which are algebraic.
   procedure(res_t) :: res
      !! User-defined residuals routine.
   procedure(jack_t) :: jack
      !! User-supplied routine to update the preconditioner.
   procedure(psol_t) :: psol
      !! User-defined preconditioner routine.
   real(rk), intent(in) :: h
      !! Scaling factor for this initial condition calculation.
   real(rk), intent(in) :: tscale ! @todo: what is "t"?
      !! Scale factor in `t`; used for stopping tests if nonzero.
   real(rk), intent(inout) :: wt(neq)
      !! Weights for error criterion.
   integer, intent(inout) :: jskip
      !! Flag to signal if initial `jack` call is to be skipped.
      !! `1`: skip the call,
      !! `0`: do not skip call.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: savr(neq)
      !! Real work array.
   real(rk), intent(inout) :: delta(neq)
      !! Real work array.
   real(rk), intent(out) :: r(neq)
      !! Real work array.
   real(rk), intent(out) :: y0(neq)
      !! Real work array.
   real(rk), intent(out) :: ydot0(neq)
      !! Real work array.
   real(rk), intent(inout) :: pwk(neq)
      !! Real work array.
   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) :: uround
      !! 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.
   real(rk), intent(in) :: ratemax
      !! Maximum convergence rate for which Newton iteration is considered converging.
   real(rk), intent(in) :: stptol
      !! Tolerance used in calculating the minimum lambda (`rl`) value allowed.
   integer, intent(in) :: jflg
      !! Flag indicating whether a `jac` routine is supplied by the user.
   integer, intent(in) :: icnflg
      !! Constraint flag. If nonzero, then constraint violations in the proposed new
      !! approximate solution will be checked for, and the maximum step length will be
      !! adjusted accordingly.
   integer, intent(out) :: icnstr(neq)
      !! Flags for checking constraints.
   integer, intent(out) :: iernls
      !! Error flag for nonlinear solver.
      !! `0`: nonlinear solver converged.
      !! `1`, `2`: recoverable error inside nonlinear solver.
      !! `1`: retry with current (y,y').
      !! `2`: retry with original (y,y').
      !! `-1`: unrecoverable error in nonlinear solver.

   integer :: iernew, ierpj, ires, liwp, lrwp, maxnit, maxnj, nj
   real(rk) :: epslin
   logical :: failed

   ! Perform initializations.
   lrwp = iwm(iwloc_lrwp)
   liwp = iwm(iwloc_liwp)
   maxnit = iwm(iwloc_mxnit)
   maxnj = iwm(iwloc_mxnj)
   nj = 0
   epslin = epli*epscon
   failed = .false.

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

   ! Preconditioner update loop
   if (.not. failed) then
      do

         ! Set all error flags to zero.
         iernls = 0
         ierpj = 0
         ires = 0
         iernew = 0

         ! If a Jacobian routine was supplied, call it.
         if ((jflg .eq. 1) .and. (jskip .eq. 0)) then
            nj = nj + 1
            iwm(iwloc_nje) = iwm(iwloc_nje) + 1
            call jack(res, ires, neq, t, y, ydot, wt, delta, r, h, cj, &
                      rwm(lrwp), iwm(liwp), ierpj, rpar, ipar)

            if ((ires .lt. 0) .or. (ierpj .ne. 0)) then
               failed = .true.
               exit
            end if
         end if

         jskip = 0

         ! Call the nonlinear Newton solver.
         call dnsik(t, y, ydot, neq, icopt, idy, res, psol, wt, rpar, ipar, &
                    savr, delta, r, y0, ydot0, pwk, rwm, iwm, cj, tscale, sqrtn, rsqrtn, &
                    epslin, epscon, ratemax, maxnit, stptol, icnflg, icnstr, iernew)

         if ((iernew .eq. 1) .and. (nj .lt. maxnj) .and. (jflg .eq. 1)) then
            ! Up to MXNIT iterations were done, the convergence rate is < 1,
            ! a Jacobian routine is supplied, and the number of JACK calls is less than MXNJ.
            ! Copy the residual SAVR to DELTA, call JACK, and try again.
            delta = savr
            cycle
         else
            exit ! Either success or convergence failure
         end if

      end do
   end if

   if (failed) then
      iernls = 2
      if (ires .le. -2) iernls = -1
      return
   end if

   if (iernew .ne. 0) then
      iernls = min(iernew, 2)
   end if

end subroutine ddasik