Set the function pointers. This must be called before dvode is called.
The two routines (dewset and dvnorm) relate to the measurement of errors.
either routine can be replaced by a user-supplied version, if desired.
however, since such a replacement may have a major impact on performance,
it should be done only when absolutely necessary, and only with great caution.
(note: the means by which the package version of a routine is
superseded by the user's version may be system-dependent.)
If not changed by the users, the defaults are used:
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(dvode_t), | intent(inout) | :: | me | |||
| procedure(f_func) | :: | f |
the name of the user-supplied subroutine defining the ode system. the system must be put in the first-order form dy/dt = f(t,y), where f is a vector-valued function of the scalar t and the vector y. subroutine f is to compute the function f. it is to have the form subroutine f (neq, t, y, ydot) real(wp) t, y(neq), ydot(neq) where neq, t, and y are input, and the array ydot = f(t,y) is output. y and ydot are arrays of length neq. subroutine f should not alter y(1),...,y(neq). f must be declared external in the calling program. if quantities computed in the f routine are needed externally to dvode, an extra call to f should be made for this purpose, for consistent and accurate results. if only the derivative dy/dt is needed, use dvindy instead. |
|||
| procedure(f_jac), | optional | :: | jac |
the name of the user-supplied routine (miter = 1 or 4) to compute the jacobian matrix, df/dy, as a function of the scalar t and the vector y. it is to have the form f_jac, where neq, t, y, ml, mu, and nrowpd are input and the array pd is to be loaded with partial derivatives (elements of the jacobian matrix) in the output. pd must be given a first dimension of nrowpd. t and y have the same meaning as in subroutine f.
jac need not provide df/dy exactly. a crude approximation (possibly with a smaller bandwidth) will do. in either case, pd is preset to zero by the solver, so that only the nonzero elements need be loaded by jac. each call to jac is preceded by a call to f with the same arguments neq, t, and y. thus to gain some efficiency, intermediate quantities shared by both calculations may be saved in a user common block by f and not recomputed by jac, if desired. also, jac may alter the y array, if desired. jac must be declared external in the calling program. |
||
| procedure(f_dewset), | optional | :: | dewset |
the following subroutine is called just before each internal
integration step, and sets the array of error weights, ewt, as
described under itol/rtol/atol above:
if the user supplies this subroutine, it must return in ewt(i) (i = 1,...,neq) a positive quantity suitable for comparison with errors in y(i). the ewt array returned by dewset is passed to the dvnorm routine (see below.), and also used by dvode in the computation of the optional output imxer, the diagonal jacobian approximation, and the increments for difference quotient jacobians. in the user-supplied version of dewset, it may be desirable to use the current values of derivatives of y. derivatives up to order nq are available from the history array yh, described above under optional output. in dewset, yh is identical to the ycur array, extended to nq + 1 columns with a column length of nyh and scale factors of h**j/factorial(j). on the first call for the problem, given by nst = 0, nq is 1 and h is temporarily set to 1.0. nyh is the initial value of neq. the quantities nq, h, and nst can be obtained by including in dewset the statements: thus, for example, the current value of dy/dt can be obtained as
|
||
| procedure(f_dvnorm), | optional | :: | dvnorm |
the following is a real function routine which computes the weighted
root-mean-square norm of a vector v:
dvnorm is called with n = neq and with w(i) = 1.0/ewt(i), where ewt is as set by subroutine dewset. if the user supplies this function, it should return a non-negative value of dvnorm suitable for use in the error control in dvode. none of the arguments should be altered by dvnorm. for example, a user-supplied dvnorm routine might:
|
subroutine initialize(me,f,jac,dewset,dvnorm) implicit none class(dvode_t),intent(inout) :: me procedure(f_func) :: f !! the name of the user-supplied subroutine defining the !! ode system. the system must be put in the first-order !! form dy/dt = f(t,y), where f is a vector-valued function !! of the scalar t and the vector y. subroutine f is to !! compute the function f. it is to have the form !! subroutine f (neq, t, y, ydot) !! real(wp) t, y(neq), ydot(neq) !! where neq, t, and y are input, and the array ydot = f(t,y) !! is output. y and ydot are arrays of length neq. !! subroutine f should not alter y(1),...,y(neq). !! f must be declared external in the calling program. !! !! if quantities computed in the f routine are needed !! externally to dvode, an extra call to f should be made !! for this purpose, for consistent and accurate results. !! if only the derivative dy/dt is needed, use [[dvindy]] instead. procedure(f_jac),optional :: jac !! the name of the user-supplied routine (miter = 1 or 4) to !! compute the jacobian matrix, df/dy, as a function of !! the scalar t and the vector y. it is to have the form [[f_jac]], !! where neq, t, y, ml, mu, and nrowpd are input and the array !! pd is to be loaded with partial derivatives (elements of the !! jacobian matrix) in the output. pd must be given a first !! dimension of nrowpd. t and y have the same meaning as in !! subroutine f. !! !! * in the full matrix case (miter = 1), ml and mu are !! ignored, and the jacobian is to be loaded into pd in !! columnwise manner, with df(i)/dy(j) loaded into pd(i,j). !! !! * in the band matrix case (miter = 4), the elements !! within the band are to be loaded into pd in columnwise !! manner, with diagonal lines of df/dy loaded into the rows !! of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j). !! ml and mu are the half-bandwidth parameters. (see iwork). !! the locations in pd in the two triangular areas which !! correspond to nonexistent matrix elements can be ignored !! or loaded arbitrarily, as they are overwritten by dvode. !! !! jac need not provide df/dy exactly. a crude !! approximation (possibly with a smaller bandwidth) will do. !! !! in either case, pd is preset to zero by the solver, !! so that only the nonzero elements need be loaded by jac. !! each call to jac is preceded by a call to f with the same !! arguments neq, t, and y. thus to gain some efficiency, !! intermediate quantities shared by both calculations may be !! saved in a user common block by f and not recomputed by jac, !! if desired. also, jac may alter the y array, if desired. !! jac must be declared external in the calling program. procedure(f_dewset),optional :: dewset !! the following subroutine is called just before each internal !! integration step, and sets the array of error weights, ewt, as !! described under itol/rtol/atol above: !! `subroutine dewset (neq, itol, rtol, atol, ycur, ewt)` !! where neq, itol, rtol, and atol are as in the [[dvode]] call sequence, !! ycur contains the current dependent variable vector, and !! ewt is the array of weights set by dewset. !! !! if the user supplies this subroutine, it must return in ewt(i) !! (i = 1,...,neq) a positive quantity suitable for comparison with !! errors in y(i). the ewt array returned by dewset is passed to the !! dvnorm routine (see below.), and also used by dvode in the computation !! of the optional output imxer, the diagonal jacobian approximation, !! and the increments for difference quotient jacobians. !! !! in the user-supplied version of dewset, it may be desirable to use !! the current values of derivatives of y. derivatives up to order nq !! are available from the history array yh, described above under !! optional output. in dewset, yh is identical to the ycur array, !! extended to nq + 1 columns with a column length of nyh and scale !! factors of h**j/factorial(j). on the first call for the problem, !! given by nst = 0, nq is 1 and h is temporarily set to 1.0. !! nyh is the initial value of neq. the quantities nq, h, and nst !! can be obtained by including in dewset the statements: !!```fortran !! type(dvode_data_t) :: sav !! call me%dvsrco(sav,job=1) !! nq = sav%nq !! h = sav%h !! nst = sav%nst !!``` !! thus, for example, the current value of dy/dt can be obtained as !! `ycur(nyh+i)/h (i=1,...,neq)` (and the division by h is !! unnecessary when nst = 0). procedure(f_dvnorm),optional :: dvnorm !! the following is a real function routine which computes the weighted !! root-mean-square norm of a vector v: !! `d = dvnorm (n, v, w)` !! where: !! !! * n = the length of the vector, !! * v = real array of length n containing the vector, !! * w = real array of length n containing weights, !! * d = `sqrt( (1/n) * sum(v(i)*w(i))**2 )`. !! !! dvnorm is called with n = neq and with w(i) = 1.0/ewt(i), where !! ewt is as set by subroutine dewset. !! !! if the user supplies this function, it should return a non-negative !! value of dvnorm suitable for use in the error control in dvode. !! none of the arguments should be altered by dvnorm. !! for example, a user-supplied dvnorm routine might: !! !! * substitute a max-norm of (v(i)*w(i)) for the rms-norm, or !! * ignore some components of v in the norm, with the effect of !! suppressing the error control on those components of y. me%f => f if (present(jac)) then me%jac => jac else me%jac => null() end if if (present(dewset)) then me%dewset => dewset else me%dewset => dewset_default end if if (present(dvnorm)) then me%dvnorm => dvnorm else me%dvnorm => dvnorm_default end if end subroutine initialize