"TwoPhaseFlow app Parameters.py, we need to force the AddedMass PETSc configuration to always use the iterative solver configuration. The superlu/superlu_dist configuration seems to work for me reliably, but since the system is singular (constant null space), I think we generally need to be using an iterative solver configured to be robust and scalable for pure Neumann Poisson problems. " - @cekees