DDAPLUS

Another common property of multireaction networks is stiffness, that is, the presence of kinetic steps with widely different rate coefficients. This property was pointed out by Curtiss and Hirschfelder (1952), and has had a major impact on the development of numerical solvers such as BASSE (Petzold 1983) and DDAPLUS of Appendix B. Since stiff equations take added computational effort, there is some incentive to reduce the stiffness of a model at the formulation stage this can be done by substituting Eq. (2.5-2b) or (2.5-3) for some of the reaction or production rate expressions. This strategy replaces some differential equations in the reactor model by algebraic ones to expedite numerical computations. 

The derivatives F r are called the first-order parametric sensitivities of the model. Their direct computation via Newton s method is implemented in Subroutines DDAPLUS (Appendix B) and PDAPLUS. Finite-difference approximations are also provided as options in GREGPLUS to produce the matrix A in either the Gauss-Newton or the full Newton form these approximations are treated in Problems 6.B and 6.C. 

This appendix explains how to use DDAPLUS to solve nonlinear initial-value problems containing ordinary differential equations with or without algebraic equations, or to solve purely algebraic nonlinear equation systems by a damped Newton method. Three detailed examples are given. 

The DDAPLUS algorithm (Caracotsios and Stewart, 1985), updated here, is an extension of the DDASSL (Petzold, 1982) implicit integrator. DDAPLUS solves differential-algebraic equation systems of the form... 

The object code for DDAPLUS and its dependencies is provided in the software package Athena Visual Studio, which is available at www.AthenaVisual.com. 

The code DDAPLUS is normally called via Athena Visual Studio. For very detailed control of the available options, DDAPLUS may be called as follows from a user-provided MAIN program for the given problem ... 

CALL DDAPLUS(t,tout,Nstvar,U,UPRIME,Rtol,Atol,Info,... 

In this section we describe the arguments used in calls to DDAPLUS. 

MAIN For each problem Driver calls DDAPLUS... 

Table B.l Codes the User Provides to Go with DDAPLUS. t, tout...

On the initial entry to DDAPLUS, the initial values Uio just before to should be provided as the first Nstvar elements of U. An initial guess is needed, as in Example B.4, for any element ui whose f-derivative does not appear. If Info(12)> 0, as in Examples B.2 and B.3, the values duio /dOj should be provided as the remaining elements U(i,l+j) of the starting vector U, as in Examples B.2 and B.3. No input value of UPRIME is required, since DDAPLUS computes this vector. On return from DDAPLUS, U and UPRIME will contain the solution and its f-derivative at the current t. 

These are relative and absolute error tolerances set by the user. They can be either scalars or arrays, depending on the user s specification of Info(2) as described below. DDAPLUS uses these inputs to construct a vector of mixed local error tolerances... 

This parameter must be set to zero whenever beginning a problem or when changing any equations. It is reserved for use by DDAPLUS at all other times, except as noted under Idid=—1. 

Set lnfo(l)=0 to indicate that this is an initial call of DDAPLUS, or the start of a new region with changed equations, with to reset to the current t. 

DDAPLUS solves the corrector equations for each t-step by Newton s method, using an iteration matrix G(t 7 to), of order Nstvar, defined by Eq. (B.l-5a). On request, Eq. (B.1-2) is solved similarly with the matrix G(to) defined by Eq. (B.l-5b). If the user does not provide a detailed subroutine Jac (see Section B.6 below), DDAPLUS will approximate G(t) by finite differences. 

Set Info(5)=0 if G(t) is to be computed by DDAPLUS via finite differences. This method is very efficient when G(t) is banded and the option Info(6)=l is used. 

DDAPLUS will perform better if it is told to use band-matrix algorithms w henever the lower and upper bandwddths of G satisfy 2ml - - mu < Nstvar. The matrix G (and 1 if nondiagonal) wall then be stored more... 

The initial stepsize can be chosen by DDAPLUS or by the user. In either case, the sign of the initial step is taken from Rwork(3), which DDAPLUS sets to +0.001 whenever the initial Rwork(3) is zero. 

This input can be used to activate special procedures for maintaining nonnegative solution elements Ui. DDAPLUS will report if it is unable to do this in that event the problem probably needs to be reformulated (for example, by restoring some neglected reverse reaction steps). 

This input value is now ignored, since DDAPLUS provides consistent initial conditions by solving any algebraic rows present in Eq. (B.l-la). The initial derivative vector UPRIME(to) is calculated automatically, so its input value need not be given. 

The calculation of UPRIME will involve the derivative (9/ /5t). 0 if any algebraic equation in the set (B.l-la) has a term depending explicitly on t. In this case, provide a nonzero At in Rwork(44) DDAPLUS will then evaluate dfA/dt)y 0 as a divided difference. Use Rwork(44)=0 otherwise. 

The DDAPLUS sensitivity algorithm uses the matrix function B t) defined in Eqs. (B.l-4a) and (B.l-4b). This function will be approximated by divided differences in DDAPLUS. unless the user sets Info(14) = l and provides the derivatives djijdBj with respect to one or more parameters in Subroutine Bsub. Set Info(14)=0 if you want DDAPLUS to approximate B(t) by differencing. 

This integer, when positive, is called lYstop and indicates that the state or sensitivity variable U(IYstop) must not cross the value Cstop in the current integration. If the index Info(17) is outside the set of variables indicated in Info(15), DDAPLUS will write a diagnostic message and return control to the calling program. 

This integer controls the amount of printout during the integration by DDAPLUS. Set Info(29)=-l to suppress any printout set Info(29) to a zero or positive value to request various levels of solution information from DDAPLUS. 

Rwork is a double precision work array of length Lrw used by DDAPLUS. A value Lrw=5000 suffices for many medium-sized problems. If the input value Lrw is too small, (for example, if Lrw=l), DDAPLUS will promptly terminate and state the additional length required. Some of the contents of Rwork are summarized here the user may access these via the program MAIN. 

Rwork(3) is an input/output register. The t-steps are taken forward if the initial value here is zero or positive, and backward otherwise. This register is initialized by DDAPLUS to 0.001 unless MAIN sets Info(8) = l and jRwork(3) > 0. On return, this register contains the signed step h to be attempted next. 

Rwork(41) contains the relative perturbation step size for finite-difference computation of G(t). The default stepsize is the square root of the machine precision DDAPLUS uses this unless MAIN provides a larger value. Rwork(42) must be set to the desired stopping value Cstop of the state or sensitivity variable U(IYstop). whenever the user inserts an index lYstop into Info (17). 

Iwork(4), when nonzero, contains the index i of the dominant variable in the test failure that stopped the DDAPLUS computation. Info(18) can be used to get related information. 

Rpar is a double precision work array, and Ipar is an integer work array. These are provided so that DDAPLUS can pass special information to the user-provided subroutines. The parameters whose sensitivity coefficients are desired must be stored in the Rpar array in the locations specified by the elements l,.Nspar of the array Ipar. For example, the input Ipar(2)=3 would indicate that the second sensitivity parameter was stored in Rpar(3). The location Ipar(O) is reserved for a user-defined constant, such as the address (later in Ipar) of a list of pointers to other information. The location Rpar(O) may be used to store values of a user-defined continuation parameter. 

The main task of DDAPLUS is to compute the state vector u (and the sensitivity matrix W if requested) at the next output point. The status of the solution upon return to the calling program is described by the parameter Idid, whose possible values are explained below. 

The initial UPRIME and h are ready at t = to, along with U updated from to- whenever algebraic equations are present. This is the normal return when DDAPLUS is called with tout=t. To commence integration, MAIN must call DDAPLUS again with tout moved from t in the direction indicated by the sign of h in Rwork(3). 

---