Lawrence,

Thanks for your help. Unfortunately we do not have a pure PETSc example, although we possibly could attempt to modify Matt's DMPlex examples since all we need are CG1 finite elements.

Our problems are relatively small so the direct solver is enough for our results, but if we do encounter bigger problems we will give this another go.

Thanks again,
Justin

On Tue, Nov 29, 2016 at 5:28 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,

On 29/11/16 08:38, Lawrence Mitchell wrote:
>
>
> On 29 Nov 2016, at 08:31, Justin Chang <jychang48@gmail.com
> <mailto:jychang48@gmail.com>> wrote:
>
>> Okay, switched over to the firedrake list
>
> I wonder. Normally this appears if you didn't get the right splittings
> for a three or more field problem.
> I'm on my phone right now so can't debug. I'll have a look a bit later.

There are a few issues here.

1. The VIRS active set method creates the Jacobian by calling
MatGetSubMatrix with arbitrary index sets. MatNest matrices don't
support extracting these.  I see you were setting a parameter to turn
these off, but we changed the way you do this recently.

Either set:

parameters["default_matrix_type"] = "aij"

or set it per-matrix with

assemble(..., mat_type="aij")

2. With that out the way, we run into a bug in PETSc (which Barry has
a branch fixing) that PCReset on a fieldsplit PC didn't clear
everything out.

3. Cherry-picking that fix, we then run into a further issue in that
the DM I provide you doesn't know anything about the currently active
set, and therefore spits back ISes that define the unconstrained
fields, so we can no longer split correctly.

In short, it looks like maybe when I said this was a firedrake
problem, it's not?

To me it looks like there's currently no support in PETSc for using
fieldsplit PCs with the VIRS if you're using a DM to define the
splits, but I could be mistaken.  There's special case code to handle
interpolation and restriction, but I think one needs more
special-casing for DMCreateFieldDecomposition and friends.

Do you have a pure PETSc example that breaks?  I guess we could
support this stuff if I knew what to do.

Lawrence