[spec] Mixed Function Spaces

NOTE: This is, at present, design documentation, not documentation of implemented functionality.

Declaring a Sparsity of a Mixed Space

Sparsities for mixed spaces are declared as a block sparsity, using the following syntax:

# Weird spacing due to square brackets getting interpreted as Markdown
sparsity = Sparsity([ [block_11, block_12],
                      [block_21, block_22] ],
                      [ [dim_11, dim_12],
                      [dim_21, dim_22] ],
                      "sparsity")

Since we special-cased mixed spaces where all the bases are the same by providing a dimension argument to the sparsity, (e.g. for vector fields), we need to now pass a matrix of dimensions.

Declaration of individual blocks

For a single block, a tuple of pairs of maps represents the sparsity:

maps = ((rmap[0], cmap[0]), (rmap[1], cmap[1]), ...)

For a valid sparsity, the following must hold:

r = len(maps)
for i in xrange(r):
    assert rmap[i].iterset == cmap[i].iterset

for i in xrange(r):
    for j in xrange(r):
        assert rmap[i].dataset == rmap[j].dataset

Declaration of dimensions

The dimensions of an individual block is a pair of integers specifying the number of elements in the row- and column-directions. For example, the dimensions of a vector-vector block in 2D would be:

(2,2)

and a vector-pressure block would be:

(2,1)

Matrix instantiation

A matrix for a block sparsity is instantiated similarly to a matrix for a normal sparsity:

m = Mat(sparsity, valuetype, name)

Assembly into a matrix block

Matrix assembly into a blocked matrix is performed by assembling into an individual block that is selected by subscripting the matrix, e.g.:

par_loop(kernel, it_space,
         m[i,j](map_pair, op2.INC))

Fluidity interface

Fluidity does not natively support mixed fields. In order to overcome this, we will define a MixedField class in flop.py. This class will hold a list of its constituent fields, and the MixedElement that corresponds to the product of the elements of its constituents. The use of it to define a mixed matrix is given in the following example:

u = state.vector_fields['Velocity']
h = state.scalar_fields['LayerThickness']

# Dolfin syntax is W=U*H for mixing elements, but this operator
# is too overloaded in flop.py to do this.
w = MixedField(u,h)

v, q = TestFunctions(w)
u, p = TrialFunctions(w)

a_u = inner(v,u)*dx
a_h = q*p*dx
a = a_u + a_h
L_u = inner(v,w)*dx
L_h = inner(q,w)*dx
L = L_u + L_h

solve(a == L, w)

Implementation Plan

For PETSc: each block of the sparsity will be represented by a Mat object. The compound matrix of all the blocks will be represented by a MatNest that contains each of the individual Mats. This requires a small extension to petsc4py to allow the creation of a MatNest object.

As an example, the block matrix with the following sparsity:

 X   X |   X 
   X   | X
 X   X |   X
 ------+----
   X   | X X
 X   X | X X

can be declared as follows:

m11 = PETSc.Mat()
m12 = PETSc.Mat()
m21 = PETSc.Mat()
m22 = PETSc.Mat()

rowptr11 = [ 0, 2, 3, 5 ]
colidx11 = [ 0, 1, 2, 0, 2 ]
val11    = [0]*5
rowptr12 = [ 0, 1, 2, 3 ]
colidx12 = [ 1, 0, 1 ]
val12    = [0]*3
rowptr21 = [ 0, 1, 3 ]
colidx21 = [ 1, 0, 2 ]
val21    = [0]*3
rowptr22 = [ 0, 2, 4 ]
colidx22 = [ 0, 1, 0, 1 ]
val22    = [0]*4

m11.createAIJWithArrays((3, 3), (rowptr11, colidx11, np.asarray(val11,dtype=PETSc.RealType)))
m12.createAIJWithArrays((3, 2), (rowptr12, colidx12, np.asarray(val12,dtype=PETSc.RealType)))
m21.createAIJWithArrays((2, 3), (rowptr21, colidx21, np.asarray(val21,dtype=PETSc.RealType)))
m22.createAIJWithArrays((2, 2), (rowptr22, colidx22, np.asarray(val22,dtype=PETSc.RealType)))

i = PETSc.IS()
i.createGeneral([0,1])

m = PETSc.Mat()
m.createNest(i,i,[m11, m12, m21, m22])

In order to do this in PyOP2, build_sparsity_pattern needs to be called for each individual sub-block with the correct row and column maps for the block, and these can then be passed to createAIJWithArrays.

It is possible to assemble into the individual blocks m11, m12, m21, m22 in the usual way.

API Implementation

In base.py and runtime_base.py, the Sparsity classes will need to be modified to hold row and column maps for each block, as well as the dimensions of entries for each block.

In order to implement the subscripting syntax:

mat[i,j](map, access)

the __getitem__ method of Mat will be overridden, to return a MatBlock object that provides the same methods as a Mat object, but applied to the individual block. The MatBlock object will hold a reference to its parent Mat object and the index of its block, in order that the data in the parent block can be re-used in the execution of the MatBlock's ,methods.

It should be an error to call a Mat object that has a block structure, since it isn't clear what this would mean at present.