This routine solves a nonlinear system of algebraic equations of the form:
for the unknown . The method used is a matrix-free Newton scheme.
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 after successful step. |
||
integer, | intent(in) | :: | neq |
Problem size. |
||
procedure(res_t) | :: | res |
User-defined residuals routine. |
|||
procedure(jack_t) | :: | jack |
User-defined Jacobian routine. |
|||
procedure(psol_t) | :: | psol |
User-defined preconditioner routine. |
|||
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(out) | :: | savr(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 dnsk for details. On the
first Newton iteration with an updated preconditioner, |
||
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 |
||
real(kind=rk), | intent(in) | :: | rsqrtn |
Reciprocal of square root of |
||
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) | :: | jflag |
Flag indicating whether a |
||
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 | = | 0 | |
integer, | public, | parameter | :: | maxit | = | 4 | |
real(kind=rk), | public, | parameter | :: | xrate | = | 0.25_rk | |
integer, | public | :: | iernew | ||||
integer, | public | :: | ierpj | ||||
integer, | public | :: | iersl | ||||
integer, | public | :: | iertyp | ||||
integer, | public | :: | ires | ||||
integer, | public | :: | j | ||||
integer, | public | :: | liwp | ||||
integer, | public | :: | lrwp | ||||
real(kind=rk), | public | :: | deltanorm | ||||
real(kind=rk), | public | :: | epslin | ||||
real(kind=rk), | public | :: | temp1 | ||||
real(kind=rk), | public | :: | temp2 | ||||
real(kind=rk), | public | :: | tolnew |
subroutine dnedk( & t, y, ydot, neq, res, jack, psol, & h, wt, jstart, idid, rpar, ipar, phi, gama, savr, delta, e, & rwm, iwm, cj, cjold, cjlast, s, uround, epli, sqrtn, rsqrtn, & epscon, jcalc, jflag, 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 matrix-free Newton scheme. 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 after successful step. integer, intent(in) :: neq !! Problem size. procedure(res_t) :: res !! User-defined residuals routine. procedure(jack_t) :: jack !! User-defined Jacobian routine. procedure(psol_t) :: psol !! User-defined preconditioner routine. 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(out) :: savr(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 [[ditmd]]. Accounts for changes in `cj` !! needed to decide whether to call [[ditmd]]. real(rk), intent(in) :: cjlast !! Previous value of `cj`. real(rk), intent(out) :: s !! Convergence factor for the Newton iteration. See [[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 (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. integer, intent(out) :: jcalc !! Flag indicating whether the Jacobian matrix needs to be updated. !! `-1`: call [[ditmd]] to update the Jacobian matrix. !! `0`: Jacobian matrix is up-to-date. !! `1`: Jacobian matrix is outdated, but [[ditmd]] will not be called unless `jcalc` is set to `-1`. integer, intent(in) :: jflag !! Flag indicating whether a `jac` routine is supplied by the user. 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. !! `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 = 0, maxit = 4 real(rk), parameter :: xrate = 0.25_rk integer :: iernew, ierpj, iersl, iertyp, ires, j, liwp, lrwp real(rk) :: deltanorm, epslin, temp1, temp2, tolnew real(rk), external :: ddwnrm ! @todo: remove this once inside module ! Verify that this is the correct subroutine. iertyp = 0 if (ntype /= 1) then iertyp = 1 goto 380 end if ! If this is the first step, perform initializations. if (jstart == 0) then cjold = cj jcalc = -1 s = 1e2_rk end if ! Perform all other initializations. iernls = 0 lrwp = iwm(iwloc_lrwp) liwp = iwm(iwloc_liwp) ! Decide whether to update the preconditioner. if (jflag /= 0) then temp1 = (one - xrate)/(one + xrate) temp2 = 1/temp1 if ((cj/cjold < temp1) .or. (cj/cjold > temp2)) jcalc = -1 if (cj /= cjlast) s = 1e2_rk else jcalc = 0 end if ! Looping point for updating preconditioner with current stepsize. 300 continue ! Initialize all error flags to zero. ierpj = 0 ires = 0 iersl = 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 330 continue epslin = epli*epscon tolnew = epslin ! Call RES to initialize DELTA. call res(t, y, ydot, cj, delta, ires, rpar, ipar) iwm(iwloc_nre) = iwm(iwloc_nre) + 1 if (ires < 0) goto 380 ! If indicated, update the preconditioner. ! Set JCALC to 0 as an indicator that this has been done. if (jcalc == -1) then jcalc = 0 call jack(res, ires, neq, t, y, ydot, wt, delta, e, h, cj, rwm(lrwp), iwm(liwp), ierpj, rpar, ipar) iwm(iwloc_nje) = iwm(iwloc_nje) + 1 cjold = cj s = 1e2_rk if (ires < 0) goto 380 if (ierpj /= 0) goto 380 end if ! Call the nonlinear Newton solver. call dnsk(t, y, ydot, neq, res, psol, wt, rpar, ipar, savr, & delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, & s, temp1, tolnew, muldel, maxit, ires, iersl, iernew) if (iernew > 0 .and. jcalc /= 0) then ! The Newton iteration had a recoverable failure with an old ! preconditioner. Retry the step with a new preconditioner. jcalc = -1 goto 300 end if if (iernew /= 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 == 0) goto 390 delta = min(y, zero) deltanorm = ddwnrm(neq, delta, wt, rpar, ipar) if (deltanorm > epscon) goto 380 e = e - delta goto 390 ! Exits from nonlinear solver. ! No convergence with current preconditioner. ! Compute IERNLS and IDID accordingly. 380 continue if ((ires <= -2) .or. (iersl < 0) .or. (iertyp /= 0)) then iernls = -1 if (ires <= -2) idid = -11 if (iersl < 0) idid = -13 if (iertyp /= 0) idid = -15 else iernls = 1 if (ires == -1) idid = -10 if (ierpj /= 0) idid = -5 if (iersl > 0) idid = -14 end if 390 continue jcalc = 1 end subroutine dnedk