Navier-Stokes: efficient assembly of the convection term exploiting block structure

Question

Dear all,

I'd like to speed up the assembly of the "fully coupled" linearized
Navier-Stokes equations discretized with inf-sup stable mixed elements (e.g., P2P1).
For a linearization of type $\vec u^n\cdot\nabla\vec u^{n+1}$ the unknown velocity
components $u^{n+1}_i$ are decoupled. It is sufficient to assemble the
convection term for just one velocity component and then build the complete
matrix by copying the matrix of the scalar into the corresponding blocks.

This is explained in the FEniCS book on p. 428.

I spent some time trying to do this, an "easy way" did not occurr to me.
My naive approach was:

• build the scalar convection matrix Ac
• extract the values and indices with (indptr, col, val) = Ac.getValuesCSR()
• shift the indptr and the column index col in order to match the DOFs of velocity components in the
  mixed function space, i.e., W.sub(0).sub(i).dofmap().dofs() $\forall i=1,...,d$
• write the extracted values into the complete matrix (suppose the full matrix
  was assembled once in the beginning) with
  A.setValuesCSR(indptr_shift, col_shift, val)

However, there are several problems.

1. DOF ordering: not trivial without parameters["reorder_dofs_serial"] = False since the ordering is not the same for all velocity components
2. In parallel this approach won't work, at least not without fiddling with
   MPI, because the matrices are distributed across the processes...

Does anyone know of a way to implement this? I'd be very glad about any hints.
I wonder if there are other DOLFIN or PETSc functions that would make this
easier..

Thanks!
David

edit: an example to make it clearer:

from dolfin import *

N = 16
mesh = UnitCubeMesh(N, N, N)

P1 = FiniteElement('P', mesh.ufl_cell(), 1)
P2 = VectorElement('P', mesh.ufl_cell(), 2)
W = FunctionSpace(mesh, P2*P1)
bc = DirichletBC(W.sub(0), Constant([0]*3), "on_boundary")

u0 = project(Constant([1]*3), W.sub(0).collapse(), solver_type='cg')

# define trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# Build complete convection matrix for u = (u1, u2, u3)
a = dot(grad(u)*u0, v)*dx
K = PETScMatrix()
assemble(a, tensor=K)

# Build convection matrix for just one component "uc"
Vc = FunctionSpace(mesh, FiniteElement('P', mesh.ufl_cell(), 2))
uc = TrialFunction(Vc)
vc = TestFunction(Vc)
ac = dot(u0, nabla_grad(uc))*vc*dx

C = PETScMatrix()
assemble(ac, tensor=C)

# In theory, K = diag((C, C, C))
(ai, aj, av) = C.mat().getValuesCSR()

for i in range(3):
    iset = W.sub(0).sub(i).dofmap().dofs()
    # some computations to determine what ai, aj mean with respect to the index
    # set of the i'th velocity sub space
    ...
    K.mat().setValuesCSR(ai_new, ai_new, av)

K.mat().assemble()

Answer 1

Hello,

Regarding the 1) dof ordering, you could do something like this to work around it....

dofs_is_list = []
    for i in range(numberOfFunctionSpaces):
        dofs_is_list.append([self._ME.sub(k).dofmap().dofs()])

After that you get a list of each subspace's DOF. From that,

# for all dofs
for j in range(len(dofs_is_list)):
    index = 0
    for dofs_i in dofs_is_list[j]:
         # Extract subfunction value and do whatever you want
         U_i = u_termsVectorList[j][dofs_i]
         # do things
         ... 
         ...
         # Store it back
         u_termsVectorList[j][dofs_i] = U_i

where, u_termsVectorList is the MIXEDFUNCTIONSPACE[INDEX].vector() function...

Hope this helps...

Comments

Just added a small example that I'll refer to. The DOF ordering problem is the following: I want to store 3 $n \times n$ matrices $C$ at the right locations inside the $(3n + m) \times (3n + m)$ matrix $K$. $C$ has a defined order, but the orderings of "K.sub(0)", "K.sub(1)" differ (shorthand for sub matrix). If the ordering of C corresponds to "K.sub(0)", the procedure above will work. But "K.sub(1)" has a different node numbering than C... (and K.sub(0)) so I would have to rearrange C first or assemble it using the same ordering.
You see?