Preconditioner Tools for Reaction-Transport Problems. Part II: Block-Grouping in Block-Diagonal Reaction-Based Factor

This module is provided to assist in the generation and solution of preconditioner matrices for problems arising from reaction-transport systems. More specifically, the routines jac_rbgpre, and psol_rbgpre are intended as auxiliary routines for the user-supplied routines jac and psol called by daskr when the Krylov method is selected.

These routines are meant for a DAE system obtained from a system of reaction-transport PDEs, in which some of the PDE variables obey differential equations, and the rest obey algebraic equations [2, section 4]. It is assumed that the right-hand sides of all the equations have the form of a sum of a reaction term $R$ and a transport term $S$, that the transport term is discretized by finite differences, and that in the spatial discretization the PDE variables at each spatial point are kept together. Thus, the DAE system function, in terms of a dependent variable vector $u$, has the form:

$$G(t,u,\dot{u}) = I_d \dot{u} - R(t,u) - S(t,u)$$

where $I_d$ is the identity matrix with zeros in the positions corresponding to the algebraic components and ones in those for the differential components, $R(t,u)$ denotes the reaction terms (spatial coupling absent), and $S(t,u)$ denotes the spatial transport terms.

As shown in Brown et al. [2], two possible preconditioners for such a system are:

• $P_R = c_J I_d - \partial R/ \partial u$, based on the reaction term $R$ alone.
• $P_{SR} = (I - c_J^{-1} \partial S/\partial u) (c_J I_d - \partial R / \partial u)$, the product of two factors (in either order), one being $P_R$ and the other being based on the spatial term $S$ alone.

The routines given here can be used to compute the reaction-based factor $P_R$. More precisely, they provide an approximation $A_R$ to $P_R$. The matrix $P_R$ is block-diagonal, with each block corresponding to one spatial point. In $A_R$, we compute each block by difference quotient approximations, by way of calls to a user-supplied subroutine rblock, that evaluates the reaction terms at a single spatial point. In addition, rather than evaluating a block at every spatial point in the mesh, we use a block-grouping scheme, described in Brown et al. [1]. In this scheme, the mesh points are grouped, as in domain decomposition, and only one block of $\partial R / \partial u)$ is computed for each group; then in solving $A_R x = b$, the inverse of the representative block is applied to all the blocks of unknowns in the group. Block-grouping greatly reduces the storage required for the preconditioner.

The routines given here are specialized to the case of a 2D problem on a rectangular mesh in the x-y plane, and for a block-grouping arrangement that is rectangular (i.e. the Cartesian product of two 1D groupings). However, they can be easily modified for a different problem geometry or a different grouping arrangement. It is also assumed that the PDE variables are ordered so that the differential variables appear first, followed by the algebraic variables.

Usage

To use these routines in conjunction with daskr, the user's calling program should include the following, in addition to setting the other daskr input parameters:

• Set info(12) = 1 to select the Krylov iterative method and info(15) = 1 to indicate that a jac routine exists.

• Call setup_rbgpre to set mesh-related and block-grouping inputs.

• A jac routine, as prescribed by the daskr instructions, which calls jac_rbgpre, and does any other Jacobian-related preprocessing needed for preconditioning. The latter routine takes as argument r0 the current value of the reaction vector $R$. This can be done, for example, by taking r0 to be rpar, and loading rpar with the vector $R$ in the last call to the res routine; in that case, the calling program must declare rpar to have length at least neq. Alternatively, insert a call to rblock within the loop over mesh points inside jac_rbgpre.

• A psol routine, as prescribed by the daskr instructions, which calls psol_rbgpre for the solution of systems $A_R x = b$, and does any other linear system solving required by the preconditioner.

Example

The program example_web demonstrates the use of this preconditioner.

References

1. Peter N. Brown and Alan C. Hindmarsh, Reduced Storage Matrix Methods in Stiff ODE Systems, J. Appl. Math. & Comp., 31 (1989), pp. 40-91.
2. Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold, Using Krylov Methods in the Solution of Large-Scale Differential-Algebraic Systems, SIAM J. Sci. Comput., 15 (1994), pp. 1467-1488.
3. Peter N. Brown, Alan C. Hindmarsh, and Linda R. Petzold, Consistent Initial Condition Calculation for Differential-Algebraic Systems, SIAM Journal on Scientific Computing 19.5 (1998): 1495-1512.