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.
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. |
||
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 |
Step size. |
||
real(kind=rk), | intent(inout) | :: | wt(neq) |
Weights for error control. |
||
integer, | intent(in) | :: | jstart |
Flag indicating whether this is the first call to this routine.
If |
||
integer, | intent(out) | :: | idid |
Completion flag. |
||
real(kind=rk), | intent(inout) | :: | rpar(*) |
User real workspace. |
||
integer, | intent(inout) | :: | ipar(*) |
User integer workspace. |
||
real(kind=rk), | intent(in) | :: | phi(neq,kp1) |
Array of divided differences used by the Newton iteration. |
||
real(kind=rk), | intent(in) | :: | gama(kp1) |
Array used to predict |
||
real(kind=rk), | intent(in) | :: | dum2(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(out) | :: | cjold |
Value of |
||
real(kind=rk), | intent(in) | :: | cjlast |
Previous value of |
||
real(kind=rk), | intent(out) | :: | s |
Convergence factor for the Newton iteration. See dnsd and dnsk for
details. On the first Newton iteration with an updated preconditioner, |
||
real(kind=rk), | intent(in) | :: | uround |
Unit roundoff. |
||
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. |
||
integer, | intent(out) | :: | jcalc |
Flag indicating whether the Jacobian matrix needs to be updated.
|
||
integer, | intent(in) | :: | dum6 |
Dummy argument. |
||
integer, | intent(in) | :: | kp1 |
The current order + 1. Updated across calls. |
||
integer, | intent(in) | :: | nonneg |
Flag indicating whether nonnegativity constraints are imposed on the solution.
|
||
integer, | intent(in) | :: | ntype |
Type of the nonlinear solver.
|
||
integer, | intent(out) | :: | iernls |
Error flag for the nonlinear solver.
|
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public, | parameter | :: | muldel | = | 1 | |
integer, | public, | parameter | :: | maxit | = | 4 | |
real(kind=rk), | public, | parameter | :: | xrate | = | 0.25_rk | |
integer, | public | :: | idum | ||||
integer, | public | :: | ierj | ||||
integer, | public | :: | iernew | ||||
integer, | public | :: | iertyp | ||||
integer, | public | :: | ires | ||||
integer, | public | :: | j | ||||
real(kind=rk), | public | :: | delnorm | ||||
real(kind=rk), | public | :: | pnorm | ||||
real(kind=rk), | public | :: | temp1 | ||||
real(kind=rk), | public | :: | temp2 | ||||
real(kind=rk), | public | :: | tolnew |
subroutine dnedd( & t, y, ydot, neq, res, jac, dum1, & h, wt, jstart, idid, rpar, ipar, phi, gama, dum2, delta, e, & rwm, iwm, cj, cjold, cjlast, s, uround, dum3, dum4, dum5, & epscon, jcalc, dum6, kp1, nonneg, ntype, iernls) !! 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(jac_t) :: jac !! User-defined Jacobian routine. procedure(psol_t) :: dum1 !! Dummy argument. real(rk), intent(in) :: h !! Step size. real(rk), intent(inout) :: wt(neq) ! @todo: confusing notation: wt ?= ewt ?= whgt !! Weights for error control. integer, intent(in) :: jstart ! @todo: ? convert to logical !! Flag indicating whether this is the first call to this routine. !! If `jstart = 0`, then this is the first call, otherwise it is not. integer, intent(out) :: idid !! Completion flag. real(rk), intent(inout) :: rpar(*) !! User real workspace. integer, intent(inout) :: ipar(*) !! User integer workspace. real(rk), intent(in) :: phi(neq, kp1) !! Array of divided differences used by the Newton iteration. real(rk), intent(in) :: gama(kp1) !! Array used to predict `y` and `ydot`. real(rk), intent(in) :: dum2(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(out) :: cjold !! Value of `cj` as of the last call to [[dmatd]]. Accounts for changes in `cj` !! needed to decide whether to call [[dmatd]]. real(rk), intent(in) :: cjlast !! Previous value of `cj`. real(rk), intent(out) :: s !! Convergence factor for the Newton iteration. See [[dnsd]] and [[dnsk]] for !! details. On the first Newton iteration with an updated preconditioner, `s = 100`. !! The value of `s` is preserved from call to call so that the rate estimate from !! a previous step can be applied to the current step. real(rk), intent(in) :: uround ! @todo: ? remove !! Unit roundoff. 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. integer, intent(out) :: jcalc !! Flag indicating whether the Jacobian matrix needs to be updated. !! `-1`: call [[dmatd]] to update the Jacobian matrix. !! `0`: Jacobian matrix is up-to-date. !! `1`: Jacobian matrix is outdated, but [[dmatd]] will not be called unless `jcalc` is set to `-1`. integer, intent(in) :: dum6 !! Dummy argument. integer, intent(in) :: kp1 !! The current order + 1. Updated across calls. integer, intent(in) :: nonneg !! Flag indicating whether nonnegativity constraints are imposed on the solution. !! `0`: no nonnegativity constraints. !! `1`: nonnegativity constraints are imposed. integer, intent(in) :: ntype !! Type of the nonlinear solver. !! `0`: modified Newton with direct solver. !! `1`: modified Newton with iterative linear solver. !! `2`: modified Newton with user-supplied linear solver. integer, intent(out) :: iernls !! Error flag for the nonlinear solver. !! `0`: nonlinear solver converged. !! `1`: recoverable error inside non-linear solver. !! `-1`: unrecoverable error inside non-linear solver. integer, parameter :: muldel = 1, maxit = 4 real(rk), parameter :: xrate = 0.25_rk integer :: idum, ierj, iernew, iertyp, ires, j real(rk) :: delnorm, pnorm, temp1, temp2, tolnew real(rk), external :: ddwnrm ! @todo: remove this once inside module ! Verify that this is the correct subroutine. iertyp = 0 if (ntype .ne. 0) then iertyp = 1 goto 380 end if ! If this is the first step, perform initializations. if (jstart .eq. 0) then cjold = cj jcalc = -1 end if ! Perform all other initializations. iernls = 0 ! Decide whether new Jacobian is needed. temp1 = (one - xrate)/(one + xrate) temp2 = 1/temp1 if ((cj/cjold .lt. temp1) .or. (cj/cjold .gt. temp2)) jcalc = -1 if (cj .ne. cjlast) s = 1e2_rk ! Entry point for updating the Jacobian with current stepsize. 300 continue ! Initialize all error flags to zero. ierj = 0 ires = 0 iernew = 0 ! Predict the solution and derivative and compute the tolerance for the Newton iteration. y = phi(:, 1) ydot = zero do j = 2, kp1 y = y + phi(:, j) ydot = ydot + gama(j)*phi(:, j) end do pnorm = ddwnrm(neq, y, wt, rpar, ipar) tolnew = 100*uround*pnorm ! Call RES to initialize DELTA. call res(t, y, ydot, cj, delta, ires, rpar, ipar) iwm(iwloc_nre) = iwm(iwloc_nre) + 1 if (ires .lt. 0) goto 380 ! If indicated, reevaluate the iteration matrix ! J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0). ! Set JCALC to 0 as an indicator that this has been done. if (jcalc .eq. -1) then iwm(iwloc_nje) = iwm(iwloc_nje) + 1 jcalc = 0 call dmatd(neq, t, y, ydot, delta, cj, h, ierj, wt, e, rwm, iwm, res, ires, uround, jac, rpar, ipar) cjold = cj s = 1e2_rk if (ires .lt. 0) goto 380 if (ierj .ne. 0) goto 380 end if ! Call the nonlinear Newton solver. temp1 = 2/(one + cj/cjold) call dnsd(t, y, ydot, neq, res, dum1, wt, rpar, ipar, dum2, & delta, e, rwm, iwm, cj, dum4, dum5, dum3, epscon, s, temp1, & tolnew, muldel, maxit, ires, idum, iernew) if (iernew .gt. 0 .and. jcalc .ne. 0) then ! The Newton iteration had a recoverable failure with an old ! iteration matrix. Retry the step with a new iteration matrix. jcalc = -1 goto 300 end if if (iernew .ne. 0) goto 380 ! The Newton iteration has converged. If nonnegativity of ! solution is required, set the solution nonnegative, if the ! perturbation to do it is small enough. If the change is too ! large, then consider the corrector iteration to have failed. if (nonneg .eq. 0) goto 390 delta = min(y, zero) delnorm = ddwnrm(neq, delta, wt, rpar, ipar) if (delnorm .gt. epscon) goto 380 e = e - delta ! Exits from nonlinear solver. ! No convergence with current iteration matrix, or singular iteration matrix. ! Compute IERNLS and IDID accordingly. 380 continue if (ires .le. -2 .or. iertyp .ne. 0) then iernls = -1 if (ires .le. -2) idid = -11 if (iertyp .ne. 0) idid = -15 else iernls = 1 if (ires .lt. 0) idid = -10 if (ierj .ne. 0) idid = -8 end if 390 continue jcalc = 1 end subroutine dnedd