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(psol_t) | :: | dum1 |
Dummy argument. |
|||
real(kind=rk), | intent(inout) | :: | wt(neq) |
Weights for error criterion. |
||
real(kind=rk), | intent(inout) | :: | rpar(*) |
User real workspace. |
||
integer, | intent(inout) | :: | ipar(*) |
User integer workspace. |
||
real(kind=rk), | intent(in) | :: | dum2(neq) |
Dummy argument. |
||
real(kind=rk), | intent(inout) | :: | delta(neq) |
Real work array. |
||
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) | :: | 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. |
||
real(kind=rk), | intent(inout) | :: | 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(in) | :: | dum6 |
Dummy argument. |
||
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 | ||||
logical, | public | :: | converged |
subroutine dnsd( & t, y, ydot, neq, res, dum1, wt, rpar, ipar, & dum2, delta, e, rwm, iwm, cj, dum3, dum4, dum5, epscon, & s, confac, tolnew, muldel, maxit, ires, dum6, iernew) !! 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(psol_t) :: dum1 !! Dummy argument. real(rk), intent(inout) :: wt(neq) !! Weights for error criterion. real(rk), intent(inout) :: rpar(*) !! User real workspace. integer, intent(inout) :: ipar(*) !! User integer workspace. real(rk), intent(in) :: dum2(neq) !! Dummy argument. real(rk), intent(inout) :: delta(neq) !! Real work array. 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) :: 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. real(rk), intent(inout) :: s !! Convergence factor for the Newton iteration. `s=rate/(1 - rate)`, where !! `rate` is the estimated rate of convergence of the Newton iteration. 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`. !! If `ires < -1`, then `iernew = -1`. integer, intent(in) :: dum6 !! Dummy argument. 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 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 .eq. 1) then delta = delta*confac end if ! Compute a new iterate (back-substitution). ! Store the correction in DELTA. call dslvd(neq, delta, rwm, iwm) ! 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 .eq. 0) then oldnorm = deltanorm if (deltanorm .le. tolnew) then converged = .true. exit end if else rate = (deltanorm/oldnorm)**(one/m) if (rate .gt. 0.9_rk) exit s = rate/(1 - rate) end if if (s*deltanorm .le. 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 .ge. 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 .lt. 0) exit end do if (converged) then iernew = 0 else if (ires .le. -2) then iernew = -1 else iernew = 1 end if end if end subroutine dnsd