Degraded SNES performance with static condensation and mixed Dirichlet/Neumann BCs

[lindsayad]

In the attached code we (together with @rhebergens and @bensworth) are solving a nonlinear vector advection-diffusion equation with HDG. If we use a direct solver on the full linearized system (no static condensation), run command python3 minimalHDG.py --sol lu, then SNES converges in 10 iterations. If we use firedrake.SCPC and a direct solver on the statically condensed linear problem, run command python3 minimalHDG.py --sol scpclu, then SNES converges in 27 iterations. Any ideas why this could be? We only noticed this for problems that have both Dirichlet and Neumann boundary conditions. Also, if we impose the Dirichlet boundary conditions weakly in the formulation instead of using DirichletBC, we noticed that the SNES count is the same whether we apply a direct solver on the full system or on the statically condensed system.

from firedrake import *
from firedrake.petsc import PETSc
import argparse

parser = argparse.ArgumentParser(add_help=False)
parser.add_argument("--sol", choices=["lu", "scpclu"], required=True)
args, _ = parser.parse_known_args()
solver = args.sol

nu = Constant(0.01)

mesh = RectangleMesh(20, 20, 10, 2)
V = VectorFunctionSpace(mesh, "DG", 2)
VT = VectorFunctionSpace(mesh, "DGT", 2)
Z = V*VT

def elBoundary(g):
    return g('-')*dS + g('+')*dS + g*ds

def ad_hdg(nu, mesh):
    
    n, h = FacetNormal(mesh), CellSize(mesh)

    solution = Function(Z)
    u, uT = split(solution)
    v, vT = TestFunctions(Z)

    x, y = SpatialCoordinate(mesh)
    uin = as_vector([4*(2-y)*(y-1)*(y>1), 0])
    bcs = [DirichletBC(Z.sub(1), uin, 1),
           DirichletBC(Z.sub(1), Constant((0.0, 0.0)), (3,4))]

    F = (2*nu*inner(sym(grad(u)), sym(grad(v)))-inner(outer(u, u), grad(v)))*dx
    F += elBoundary((80*nu/h)*inner(vT-v, uT-u)
                    - 2*nu*inner(sym(grad(u))*n, v-vT)
                    - 2*nu*inner(sym(grad(v))*n, u-uT)
                    + inner(u-uT,n)*inner(v-vT,n)
                    + 0.5*dot(u,n)*dot(u+uT, v-vT)
                    + 0.5*abs(dot(u,n))*dot(u-uT, v-vT))
    F += dot(uT,n)*dot(uT, vT)*ds(2)

    return (solution, F, bcs)

if solver == "lu":
    params = {"snes_type": "newtonls",
              "snes_monitor": None,
              "ksp_type": "preonly",
              "pc_type": "lu",
              "pc_factor_mat_solver_type": "mumps"}
else:
    params = {"snes_type": "newtonls",
              "snes_monitor": None,
              "ksp_type": "preonly",
              "mat_type": "matfree",
              "pmat_type": "matfree",
              "pc_type": "python",
              "pc_python_type": "firedrake.SCPC",
              "pc_sc_eliminate_fields": "0",
              "condensed_field": {"ksp_type": "preonly",
                                  "pc_type": "lu",
                                  "pc_factor_mat_solver_type": "mumps"}}

solution, F, bcs = ad_hdg(nu, mesh)
solution.assign(0)
problem = NonlinearVariationalProblem(F, solution, bcs=bcs)
solver = NonlinearVariationalSolver(problem, solver_parameters=params)
solver.solve()