daskr_rbdpre Module

Preconditioner Tools for Reaction-Transport Problems. Part I: Block-Diagonal Reaction-Based Factor without Grouping

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_rbdpre, and psol_rbdpre 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 and a transport term , 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 , has the form:

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

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

  • , based on the reaction term alone.
  • , the product of two factors (in either order), one being and the other being based on the spatial term alone.

The routines given here can be used to compute the reaction-based factor . More precisely, they provide an approximation to . The matrix is block-diagonal, with each block corresponding to one spatial point. In , 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. has one such block for each spatial point in the mesh. For a more economical approximation, see daskr_rbgpre which uses block-grouping in .

The routines given here are specialized to the case of a 2D problem on a rectangular mesh in the x-y plane. However, they can be easily modified for a different problem geometry. It is also assumed that the PDE variables are ordered so that the differential variables appear first, followed by the algebraic variables.

Usage

Warning

This module uses module variables to share the mesh and block-group parameters among routines and is, thus, thread unsafe.

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_rbdpre to set mesh-related inputs.

  • A jac routine, as prescribed by the daskr instructions, which calls jac_rbdpre, and does any other Jacobian-related preprocessing needed for preconditioning. The latter routine takes as argument r0 the current value of the reaction vector . This can be done, for example, by taking r0 to be rpar, and loading rpar with the vector 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_rbdpre.

  • A psol routine, as prescribed by the daskr instructions, which calls psol_rbdpre for the solution of systems , 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.

Uses

  • module~~daskr_rbdpre~~UsesGraph module~daskr_rbdpre daskr_rbdpre module~daskr_kinds daskr_kinds module~daskr_rbdpre->module~daskr_kinds iso_fortran_env iso_fortran_env module~daskr_kinds->iso_fortran_env

Used by

  • module~~daskr_rbdpre~~UsedByGraph module~daskr_rbdpre daskr_rbdpre proc~jac~2 jac proc~jac~2->module~daskr_rbdpre proc~psol psol proc~psol->module~daskr_rbdpre program~example_web example_web program~example_web->module~daskr_rbdpre

Subroutines

public subroutine setup_rbdpre(mx, my, ns, nsd, lid, iwork)

This routine sets the mesh parameters needed to use the routines jac_rbdpre and psol_rbdpre, assuming a 2D rectangular problem. Additionally, it loads the lengths lrwp and liwp in array iwork.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: mx

Number of mesh points in x-direction.

integer, intent(in) :: my

Number of mesh points in y-direction.

integer, intent(in) :: ns

Number of PDE variables, the size of each block in the block-diagonal preconditioner matrix .

integer, intent(in) :: nsd

Number of differential PDE variables. In the DAE system, the first nsd variables at each spatial point have time derivatives, and the remaining (ns - nsd) do not.

integer, intent(in) :: lid

Flag indicating whether to load the ID array in iwork. If lid > 0, set the ID array in iwork, indicating which components are differential and which are algebraic. This value is required if either info(11) = 1 or info(16) = 1, in which case set lid = 40 or lid = 40 + neq, depending on the value of the constraint option info(10). Otherwise, set lid = 0.

integer, intent(inout) :: iwork(*)

Integer work array.

public subroutine jac_rbdpre(t, u, r0, rblock, r1, rewt, cj, bd, ipbd, ierr)

This routine generates and preprocesses a block-diagonal preconditioner matrix, based on the part of the Jacobian corresponding to the reaction terms of the problem. It generates a matrix of the form . It calls dgefa from LINPACK to do the LU decomposition of each diagonal block. The computation of the diagonal blocks uses the mesh information in the module variables. One block per spatial point is computed. The Jacobian elements are generated by difference quotients. This routine calls a user routine rblock to compute the block (jx,jy) of .

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: t

Current independent variable.

real(kind=rk), intent(inout) :: u(*)

Current dependent variables.

real(kind=rk), intent(in) :: r0(*)

Current value of the vector .

private subroutine rblock(t, jx, jy, uxy, rxy)

User routine that computes a single block of .

Arguments
Type IntentOptional Attributes Name
real(kind=rk), intent(in) :: t

Current independent variable.

integer, intent(in) :: jx

Spatial index in x-direction.

integer, intent(in) :: jy

Spatial index in y-direction.

real(kind=rk), intent(in) :: uxy(*)

Block of ns dependent variables at spatial point (jx,jy).

real(kind=rk), intent(out) :: rxy(*)

Block (jx,jy) of .

real(kind=rk), intent(inout) :: r1(*)

Real work space available to this subroutine.

real(kind=rk), intent(in) :: rewt(*)

Reciprocal error weights.

real(kind=rk), intent(in) :: cj

Scalar used in forming the system Jacobian.

real(kind=rk), intent(out) :: bd(*)

LU factors of the diagonal blocks.

integer, intent(out) :: ipbd(*)

Pivots for the LU factorizations.

integer, intent(inout) :: ierr

Error flag. 0: no error occurred; k > 0: a zero pivot was found at the k-th stage in one of the LU factorizations.

public pure subroutine psol_rbdpre(b, bd, ipbd)

This routine solves a linear system , using the LU factors of the diagonal blocks computed in jac_rbdpre and mesh parameters passed as module variables. The solution is carried out by dgesl from LINPACK.

Arguments

Type IntentOptional Attributes Name
real(kind=rk), intent(inout) :: b(*)

Right-hand side vector on entry and solution vector on return.

real(kind=rk), intent(in) :: bd(*)

LU factors of the diagonal blocks. bd corresponds to the segment rwp of rwork.

integer, intent(in) :: ipbd(*)

Pivots for the LU factorizations. ipbd corresponds to the segment iwp of iwork.