I’m trying to implement a version of the SIMPLE algorithm in FEniCS (interested in comparing to a Newton solver for steady state). One of the steps involves solving a Poisson-type system of the form C*D^{-1}*B*p = C*u, where C is basically a discretized divergence and B is a gradient. D is the diagonal of a separate A matrix, which is derived from the linearized Navier-Stokes equations.
Obviously the standard approach to something like a Poisson problem would be to construct the bilinear operator like assemble(inner(grad(p), grad(q))*dx), but the problem in this case is the diagonal term. As far as I can tell there’s no way to do that in the usual variational form, since it gets pulled out of the already-assembled matrix A.
My approach was to construct the operators separately and then do a matrix-matrix product:
from fenics import *
mesh = UnitSquareMesh(16, 16)
V = VectorFunctionSpace(mesh, 'P', 2) # Velocity vector space
Q = FunctionSpace(mesh, 'P', 1) # Pressure vector space
bc = DirichletBC(Q, Constant(0), 'near(x[0], 1)')
u, v = TrialFunction(V), TestFunction(V)
p, q = TrialFunction(Q), TestFunction(Q)
B = assemble( -inner(grad(p), v)*dx )
C = assemble( div(u)*q*dx )
L = Matrix( PETScMatrix(
as_backend_type(C).mat().matMult(
as_backend_type(B).mat() ) ) )
bc.apply(L)
From there it’s relatively straightforward to include the D^{-1} term. The problem I’m having is in the final line, which throws an error:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to set given (local) rows to identity matrix.
*** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where: This error was encountered inside PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------
I can’t really make sense of this; assembling B and C with keep_diagonal doesn’t fix it, and the diagonal is dense:
>>> min(abs(np.diag(L.array())))
0.00021701388888888739
I’ve tried using the same order of elements for pressure and velocity rather than Taylor-Hood, but that doesn’t seem to be the issue. Also tried using a mixed FunctionSpace but that doesn’t help either. Does anyone know what’s going wrong here? Is there a problem with composing operators like this?
Also, this is a bit off-topic, but I’m also open to suggestions if SIMPLE is a bad idea. It just seemed like the Newton solver was starting to use a lot of memory in 3D, and a lot of commercial codes use SIMPLE.