This routines solves a nonlinear system of algebraic equations of the form:
for the unknown . The method used is a modified 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. |
||
integer, | intent(in) | :: | neq |
Problem size. |
||
procedure(res_t) | :: | res |
User-defined residuals routine. |
|||
procedure(psol_t) | :: | psol |
User-defined preconditioner routine. |
|||
real(kind=rk), | intent(inout) | :: | wt(neq) |
Weights for error control. |
||
real(kind=rk), | intent(inout) | :: | rpar(*) |
User real workspace. |
||
integer, | intent(inout) | :: | ipar(*) |
User integer workspace. |
||
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(in) | :: | sqrtn |
Square root of |
||
real(kind=rk), | intent(in) | :: | rsqrtn |
Reciprocal of square root of |
||
real(kind=rk), | intent(in) | :: | epslin |
Tolerance for linear system solver. |
||
real(kind=rk), | intent(in) | :: | epscon |
Tolerance for convergence of the Newton iteration. |
||
real(kind=rk), | intent(out) | :: | s |
Convergence factor for the Newton iteration. |
||
real(kind=rk), | intent(in) | :: | confac |
Residual scale factor to improve convergence. |
||
real(kind=rk), | intent(in) | :: | tolnew |
Tolerance on the norm of the Newton correction in the alternative Newton convergence test. |
||
integer, | intent(in) | :: | muldel |
Flag indicating whether or not to multiply |
||
integer, | intent(in) | :: | maxit |
Maximum number of Newton iterations allowed. |
||
integer, | intent(out) | :: | ires |
Error flag from |
||
integer, | intent(out) | :: | iersl |
Error flag from the linear system solver. See dslvk for details.
If |
||
integer, | intent(out) | :: | iernew |
Error flag for the Newton iteration.
|
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public | :: | m | ||||
real(kind=rk), | public | :: | deltanorm | ||||
real(kind=rk), | public | :: | oldnorm | ||||
real(kind=rk), | public | :: | rate | ||||
real(kind=rk), | public | :: | rhok | ||||
logical, | public | :: | converged |
subroutine dnsk( & t, y, ydot, neq, res, psol, wt, rpar, ipar, & savr, delta, e, rwm, iwm, cj, sqrtn, rsqrtn, epslin, epscon, & s, confac, tolnew, muldel, maxit, ires, iersl, iernew) !! This routines 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. 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(psol_t) :: psol !! User-defined preconditioner routine. real(rk), intent(inout) :: wt(neq) ! @todo: confusing notation: wt ?= ewt ?= whgt !! Weights for error control. real(rk), intent(inout) :: rpar(*) !! User real workspace. integer, intent(inout) :: ipar(*) !! User integer workspace. 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(in) :: sqrtn !! Square root of `neq`. real(rk), intent(in) :: rsqrtn !! Reciprocal of square root of `neq`. real(rk), intent(in) :: epslin !! Tolerance for linear system solver. real(rk), intent(in) :: epscon !! Tolerance for convergence of the Newton iteration. real(rk), intent(out) :: s !! Convergence factor for the Newton iteration. `s=rate/(1 - rate)`, where !! `rate` is the estimated rate of convergence of the Newton iteration. !! The closer `rate` is to 0, the faster the Newton iteration is converging. !! the closer `rate` is to 1, the slower the Newton iteration is converging. real(rk), intent(in) :: confac !! Residual scale factor to improve convergence. real(rk), intent(in) :: tolnew !! Tolerance on the norm of the Newton correction in the alternative Newton !! convergence test. integer, intent(in) :: muldel ! @todo: convert to logical !! Flag indicating whether or not to multiply `delta` by `confac`. integer, intent(in) :: maxit !! Maximum number of Newton iterations allowed. integer, intent(out) :: ires !! Error flag from `res`. If `ires < -1`, then `iernew = -1`. integer, intent(out) :: iersl !! Error flag from the linear system solver. See [[dslvk]] for details. !! If `iersl < 0`, then `iernew = -1`. !! If `iersl == 1`, then `iernew = 1`. integer, intent(out) :: iernew !! Error flag for the Newton iteration. !! `0`: the Newton iteration converged. !! `1`: recoverable error inside Newton iteration. !! `-1`: unrecoverable error inside Newton iteration. integer :: m real(rk) :: deltanorm, oldnorm, rate, rhok real(rk), external :: ddwnrm ! @todo: remove this once inside module logical :: converged ! Initialize Newton counter M and accumulation vector E. m = 0 e = zero ! Corrector loop. converged = .false. do iwm(iwloc_nni) = iwm(iwloc_nni) + 1 ! If necessary, multiply residual by convergence factor. if (muldel == 1) then delta = delta*confac end if ! Save residual in SAVR. savr = delta ! Compute a new iterate. Store the correction in DELTA. call dslvk(neq, y, t, ydot, savr, delta, wt, rwm, iwm, & res, ires, psol, iersl, cj, epslin, sqrtn, rsqrtn, rhok, & rpar, ipar) if ((ires /= 0) .or. (iersl /= 0)) exit ! Update Y, E, and YDOT. y = y - delta e = e - delta ydot = ydot - cj*delta ! Test for convergence of the iteration. deltanorm = ddwnrm(neq, delta, wt, rpar, ipar) if (m == 0) then oldnorm = deltanorm if (deltanorm <= tolnew) then converged = .true. exit end if else rate = (deltanorm/oldnorm)**(one/m) if (rate > 0.9_rk) exit s = rate/(one - rate) end if if (s*deltanorm <= epscon) then converged = .true. exit end if ! The corrector has not yet converged. Update M and test whether ! the maximum number of iterations have been tried. m = m + 1 if (m >= maxit) exit ! Evaluate the residual, and go back to do another iteration. ires = 0 call res(t, y, ydot, cj, delta, ires, rpar, ipar) iwm(iwloc_nre) = iwm(iwloc_nre) + 1 if (ires < 0) exit end do if (converged) then iernew = 0 else if ((ires <= -2) .or. (iersl < 0)) then iernew = -1 else iernew = 1 end if end if end subroutine dnsk