Fix FacetSplit.update

firedrake/preconditioners/facet_split.py

@@ -75,27 +75,30 @@ def initialize(self, pc):
             self.needs_zeroing = len(indices) < V.dof_count
             self.subset = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF)
 
+        self.reset_operators = False
         if mat_type != "submatrix":
             from firedrake.assemble import get_assembler
             form_assembler = get_assembler(mixed_operator, bcs=mixed_bcs,
                                            form_compiler_parameters=fcp,
                                            mat_type=mat_type,
                                            options_prefix=options_prefix)
             self.P = form_assembler.allocate()
-            self._assemble_mixed_op = form_assembler.assemble
-            self._assemble_mixed_op(tensor=self.P)
             self.mixed_opmat = self.P.petscmat
-            self.set_nullspaces(pc)
-            self.work_vecs = self.mixed_opmat.createVecs()
+            self._permute_op = partial(form_assembler.assemble, tensor=self.P)
         elif self.subset:
+            self.mixed_opmat = PETSc.Mat()
             global_indices = V.dof_dset.lgmap.apply(self.subset.indices)
             self._global_iperm = PETSc.IS().createGeneral(global_indices, comm=pc.comm)
-            self._permute_op = partial(PETSc.Mat().createSubMatrixVirtual, P, self._global_iperm, self._global_iperm)
-            self.mixed_opmat = self._permute_op()
-            self.set_nullspaces(pc)
-            self.work_vecs = self.mixed_opmat.createVecs()
+            self._permute_op = partial(self.mixed_opmat.createSubMatrixVirtual, P, self._global_iperm)
+            self.reset_operators = True
         else:
             self.mixed_opmat = P
+            self._permute_op = lambda *args: None
+
+        aux = self._permute_op()
+        if aux is not None:
+            self.set_nullspaces(pc)
+            self.work_vecs = self.mixed_opmat.createVecs()
 
         # Internally, we just set up a PC object that the user can configure
         # however from the PETSc command line.  Since PC allows the user to specify
@@ -119,7 +122,7 @@ def initialize(self, pc):
         scpc.setOptionsPrefix(options_prefix)
         scpc.setOperators(A=self.mixed_opmat, P=self.mixed_opmat)
         self.pc = scpc
-        with dmhooks.add_hooks(mixed_dm, self, appctx=self._ctx_ref, save=False):
+        with dmhooks.add_hooks(mixed_dm, self, appctx=self._ctx_ref, save=True):
             scpc.setFromOptions()
 
     def set_nullspaces(self, pc):
@@ -142,14 +145,10 @@ def _restrict_nullspace(nsp):
         Pmat.setTransposeNullSpace(_restrict_nullspace(P.getTransposeNullSpace()))
 
     def update(self, pc):
-        if hasattr(self, "_permute_op"):
-            for mat in self.pc.getOperators():
-                mat.destroy()
-            P = self._permute_op()
-            self.pc.setOperators(A=P, P=P)
-            self.mixed_opmat = P
-        elif hasattr(self, "P"):
-            self._assemble_mixed_op(tensor=self.P)
+        self._permute_op()
+        if self.reset_operators:
+            self.pc.setOperators(A=self.mixed_opmat, P=self.mixed_opmat)
+            self.pc.setFromOptions()
 
     def prolong(self, x, y):
         if x is not y:

firedrake/preconditioners/pmg.py

@@ -156,11 +156,10 @@ def initialize(self, obj):
 
         ppc.incrementTabLevel(1, parent=obj)
         ppc.setFromOptions()
-        ppc.setUp()
         self.ppc = ppc
 
     def update(self, obj):
-        self.ppc.setUp()
+        pass
 
     def view(self, obj, viewer=None):
         if viewer is None:

tests/firedrake/regression/test_facet_split.py

@@ -2,7 +2,9 @@
 from firedrake import *
 
 
-def run_facet_split(quadrilateral, pc_type, refine=2):
+def run_facet_split(quadrilateral, pc_type, refine=2, betas=None):
+    if betas is None:
+        betas = [0]
     if pc_type == "lu":
         parameters = {
             "mat_type": "matfree",
@@ -17,6 +19,7 @@ def run_facet_split(quadrilateral, pc_type, refine=2):
                 "pc_fieldsplit_schur_precondition": "selfp",
                 "fieldsplit_0_pc_type": "jacobi",
                 "fieldsplit_1_pc_type": "lu",
+                "fieldsplit_ksp_type": "preonly",
             },
         }
     elif pc_type == "jacobi":
@@ -37,6 +40,26 @@ def run_facet_split(quadrilateral, pc_type, refine=2):
                 "fieldsplit_1_ksp_rtol": 1E-12,
             },
         }
+    elif pc_type == "fdm":
+        parameters = {
+            "mat_type": "matfree",
+            "ksp_type": "preonly",
+            "pc_type": "python",
+            "pc_python_type": "firedrake.FacetSplitPC",
+            "facet": {
+                "mat_type": "submatrix",
+                "pc_type": "fieldsplit",
+                "pc_fieldsplit_type": "schur",
+                "pc_fieldsplit_schur_fact_type": "full",
+                "pc_fieldsplit_schur_precondition": "a11",
+                "fieldsplit_0_pc_type": "jacobi",
+                "fieldsplit_1_pc_type": "python",
+                "fieldsplit_1_pc_python_type": "firedrake.FDMPC",
+                "fieldsplit_1_fdm_static_condensation": True,
+                "fieldsplit_1_fdm_pc_type": "lu",
+                "fieldsplit_ksp_type": "preonly",
+            },
+        }
 
     r = refine
     variant = "fdm" if quadrilateral else None
@@ -46,15 +69,21 @@ def run_facet_split(quadrilateral, pc_type, refine=2):
     u = TrialFunction(V)
     v = TestFunction(V)
     uh = Function(V)
-
-    a = inner(grad(u), grad(v)) * dx
-    L = inner(Constant(0), v) * dx
     x = SpatialCoordinate(mesh)
     u_exact = 42 * x[1]
     bcs = [DirichletBC(V, Constant(0), 3),
            DirichletBC(V, Constant(42), 4)]
 
-    solve(a == L, uh, bcs=bcs, solver_parameters=parameters)
+    beta = Constant(0)
+    a = inner(grad(u), grad(v)) * dx + inner(u*beta, v) * dx
+    L = inner(u_exact * beta, v) * dx
+    problem = LinearVariationalProblem(a, L, uh, bcs=bcs)
+    solver = LinearVariationalSolver(problem, solver_parameters=parameters)
+    for val in betas:
+        beta.assign(val)
+        uh.assign(0)
+        solver.solve()
+
     return sqrt(assemble(inner(uh - u_exact, uh - u_exact) * dx))
 
 
@@ -68,3 +97,9 @@ def test_facet_split(quadrilateral, pc_type):
 @pytest.mark.parametrize("pc_type", ["lu", "jacobi"])
 def test_facet_split_parallel(pc_type):
     assert run_facet_split(True, pc_type, refine=3) < 1E-10
+
+
+@pytest.mark.parametrize("quadrilateral", [True])
+@pytest.mark.parametrize("pc_type", ["fdm"])
+def test_facet_split_update(quadrilateral, pc_type):
+    assert run_facet_split(quadrilateral, pc_type, betas=[1E4, 0]) < 1E-10