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 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:
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. 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 is computed for each group; then in solving
, 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.
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_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 . 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_rbgpre.
A psol
routine, as prescribed by the daskr instructions, which calls psol_rbgpre
for the solution of systems , and does any other linear system solving required
by the preconditioner.
The program example_web demonstrates the use of this preconditioner.
This routine sets the mesh and block-grouping parameters needed to use the routines
jac_rbgpre and psol_rbgpre, assuming a 2D rectangular problem with uniform
rectangular block-grouping. Additionally, it loads the lengths lrwp
and liwp
in array
iwork
.
Type | Intent | Optional | 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 |
||
integer, | intent(in) | :: | nxg |
Number of groups in x-direction. |
||
integer, | intent(in) | :: | nyg |
Number of groups in y-direction. |
||
integer, | intent(in) | :: | lid |
Flag indicating whether to load the ID array in |
||
integer, | intent(inout) | :: | iwork(*) |
Integer work array. |
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, using
block-grouping. 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 group is computed. The Jacobian elements are generated by difference quotients.
This routine calls a user routine rblock
to compute the block (jx,jy)
of .
Type | Intent | Optional | 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
|
||||||||||||||||||||||||||||||||||||||||||||||||
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.
|
This routine solves a linear system , using the LU factors of the diagonal blocks computed in jac_rbgpre and mesh parameters passed as module variables. The solution is carried out by dgesl from LINPACK.
Type | Intent | Optional | 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. |
||
integer, | intent(in) | :: | ipbd(*) |
Pivots for the LU factorizations. |