From de74328dcd05c56e4ed34e9b5d6fb8c3ca15e809 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 16 Apr 2025 14:55:38 +0100 Subject: [PATCH 1/4] Fix FacetSplit.update --- firedrake/preconditioners/facet_split.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/firedrake/preconditioners/facet_split.py b/firedrake/preconditioners/facet_split.py index 2d49f9e908..4a96733c1d 100644 --- a/firedrake/preconditioners/facet_split.py +++ b/firedrake/preconditioners/facet_split.py @@ -88,9 +88,10 @@ def initialize(self, pc): self.set_nullspaces(pc) self.work_vecs = self.mixed_opmat.createVecs() elif self.subset: + self.P = 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._permute_op = partial(self.P.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() @@ -143,11 +144,10 @@ def _restrict_nullspace(nsp): 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 + with dmhooks.add_hooks(self.pc.getDM(), self, appctx=self._ctx_ref, save=True): + self.pc.setFromOptions() elif hasattr(self, "P"): self._assemble_mixed_op(tensor=self.P) From b0fb4dfe396bbaa5b226a23c0c1696f1fd96b284 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 16 Apr 2025 15:55:19 +0100 Subject: [PATCH 2/4] add a test --- firedrake/preconditioners/facet_split.py | 6 +-- .../firedrake/regression/test_facet_split.py | 45 ++++++++++++++++--- 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/firedrake/preconditioners/facet_split.py b/firedrake/preconditioners/facet_split.py index 4a96733c1d..2da82c20d2 100644 --- a/firedrake/preconditioners/facet_split.py +++ b/firedrake/preconditioners/facet_split.py @@ -88,11 +88,11 @@ def initialize(self, pc): self.set_nullspaces(pc) self.work_vecs = self.mixed_opmat.createVecs() elif self.subset: - self.P = PETSc.Mat() + 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(self.P.createSubMatrixVirtual, P, self._global_iperm, self._global_iperm) - self.mixed_opmat = self._permute_op() + self._permute_op = partial(self.mixed_opmat.createSubMatrixVirtual, P, self._global_iperm, self._global_iperm) + self._permute_op() self.set_nullspaces(pc) self.work_vecs = self.mixed_opmat.createVecs() else: diff --git a/tests/firedrake/regression/test_facet_split.py b/tests/firedrake/regression/test_facet_split.py index 3f0748c86a..d66fcea0fc 100644 --- a/tests/firedrake/regression/test_facet_split.py +++ b/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 From 05c3523389c482af6b9034b64644b5aeab1a7e18 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 23 Apr 2025 15:57:20 +0100 Subject: [PATCH 3/4] address review comments --- firedrake/preconditioners/facet_split.py | 28 +++++++++++------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/firedrake/preconditioners/facet_split.py b/firedrake/preconditioners/facet_split.py index 2da82c20d2..09c65a6ac9 100644 --- a/firedrake/preconditioners/facet_split.py +++ b/firedrake/preconditioners/facet_split.py @@ -75,6 +75,7 @@ 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, @@ -82,21 +83,22 @@ def initialize(self, pc): 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(self.mixed_opmat.createSubMatrixVirtual, P, self._global_iperm, self._global_iperm) - self._permute_op() - self.set_nullspaces(pc) - self.work_vecs = self.mixed_opmat.createVecs() + 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 @@ -120,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): @@ -143,13 +145,9 @@ def _restrict_nullspace(nsp): Pmat.setTransposeNullSpace(_restrict_nullspace(P.getTransposeNullSpace())) def update(self, pc): - if hasattr(self, "_permute_op"): - P = self._permute_op() - self.pc.setOperators(A=P, P=P) - with dmhooks.add_hooks(self.pc.getDM(), self, appctx=self._ctx_ref, save=True): - self.pc.setFromOptions() - 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) def prolong(self, x, y): if x is not y: From f6f122c0de1de6cbd699bf4be7365d193214814c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 1 May 2025 23:13:49 +0100 Subject: [PATCH 4/4] WIP call setFromOptions as a workaround --- firedrake/preconditioners/facet_split.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/firedrake/preconditioners/facet_split.py b/firedrake/preconditioners/facet_split.py index 09c65a6ac9..bbeefb7cef 100644 --- a/firedrake/preconditioners/facet_split.py +++ b/firedrake/preconditioners/facet_split.py @@ -122,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=True): + with dmhooks.add_hooks(mixed_dm, self, appctx=self._ctx_ref, save=False): scpc.setFromOptions() def set_nullspaces(self, pc): @@ -149,6 +149,10 @@ def update(self, pc): if self.reset_operators: self.pc.setOperators(A=self.mixed_opmat, P=self.mixed_opmat) + dm = self.pc.getDM() + with dmhooks.add_hooks(dm, self, appctx=self._ctx_ref): + self.pc.setFromOptions() + def prolong(self, x, y): if x is not y: if self.needs_zeroing: