ddasid Subroutine

subroutine ddasid(t, y, ydot, neq, icopt, idy, res, jac, dum1, h, tscale, wt, dum2, rpar, ipar, dum3, delta, r, y0, ydot0, dum4, rwm, iwm, cj, uround, dum5, dum6, dum7, epscon, ratemax, stptol, dum8, icnflg, icnstr, iernls)

Uses

  • proc~~ddasid~~UsesGraph proc~ddasid ddasid module~daskr daskr proc~ddasid->module~daskr module~daskr_kinds daskr_kinds proc~ddasid->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. 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.

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(jac_t) :: jac

User-defined Jacobian routine.

procedure(psol_t) :: dum1

Dummy argument.

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

Dummy argument.

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

User real workspace.

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

User integer workspace.

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

Dummy argument.

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) :: dum4(neq)

Dummy argument.

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.

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

Dummy argument.

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

Dummy argument.

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

Dummy argument.

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

Dummy argument.

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~~ddasid~~CallsGraph proc~ddasid ddasid dmatd dmatd proc~ddasid->dmatd dnsid dnsid proc~ddasid->dnsid

Variables

Type Visibility Attributes Name Initial
integer, public :: ierj
integer, public :: iernew
integer, public :: ires
integer, public :: mxnit
integer, public :: mxnj
integer, public :: nj
logical, public :: failed

Source Code

subroutine ddasid( &
   t, y, ydot, neq, icopt, idy, res, jac, dum1, h, tscale, &
   wt, dum2, rpar, ipar, dum3, delta, r, y0, ydot0, dum4, rwm, iwm, cj, uround, &
   dum5, dum6, dum7, epscon, ratemax, stptol, dum8, 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. 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
   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(jac_t) :: jac
      !! User-defined Jacobian routine.
   procedure(psol_t) :: dum1
      !! Dummy argument.
   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) :: dum2
      !! Dummy argument.
   real(rk), intent(inout) :: rpar(*)
      !! User real workspace.
   integer, intent(inout) :: ipar(*)
      !! User integer workspace.
   real(rk), intent(out) :: dum3(neq)
      !! Dummy argument.
   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) :: dum4(neq)
      !! Dummy argument.
   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 ! @todo: used?
      !! Unit roundoff.
   real(rk), intent(in) :: dum5
      !! Dummy argument.
   real(rk), intent(in) :: dum6
      !! Dummy argument.
   real(rk), intent(in) :: dum7
      !! Dummy argument.
   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) :: dum8
      !! Dummy argument.
   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 :: ierj, iernew, ires, mxnit, mxnj, nj
   logical :: failed

   ! Perform initializations.
   mxnit = iwm(iwloc_mxnit)
   mxnj = iwm(iwloc_mxnj)
   iernls = 0
   nj = 0
   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) then
      failed = .true.
   end if

   ! Looping point for updating the Jacobian.
   if (.not. failed) then
      do

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

         ! Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYDOT.
         call dmatd(neq, t, y, ydot, delta, cj, h, ierj, wt, r, rwm, iwm, res, ires, &
                    uround, jac, rpar, ipar)
         nj = nj + 1
         iwm(iwloc_nje) = iwm(iwloc_nje) + 1

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

         ! Call the nonlinear Newton solver.
         call dnsid(t, y, ydot, neq, icopt, idy, res, wt, rpar, ipar, delta, r, &
                    y0, ydot0, rwm, iwm, cj, tscale, epscon, ratemax, mxnit, stptol, &
                    icnflg, icnstr, iernew)

         if ((iernew .eq. 1) .and. (nj .lt. mxnj)) then
            ! MXNIT iterations were done, the convergence rate is < 1,
            ! and the number of Jacobian evaluations is less than MXNJ.
            ! Call RES, reevaluate the Jacobian, and try again.
            ires = 0
            call res(t, y, ydot, cj, delta, ires, rpar, ipar)
            iwm(iwloc_nre) = iwm(iwloc_nre) + 1
            if (ires .lt. 0) then
               failed = .true.
               exit
            end if
            cycle
         else
            exit ! Either success or convergence failure
         end if

      end do
   end if

   ! Unsuccessful exits from nonlinear solver.
   ! Compute IERNLS accordingly.
   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 ddasid