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.
Type | Intent | Optional | 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 |
||
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 ( |
||
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.
|
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 |
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