Fieldsplit: update MatNest Jacobian

Author: pbrubeck

firedrake/dmhooks.py

@@ -363,43 +363,36 @@ def create_subdm(dm, fields, *args, **kwargs):
     :arg fields: The fields in the new sub-DM.
     """
     W = get_function_space(dm)
-    ctx = get_appctx(dm)
-    coarsen = get_ctx_coarsener(dm)
-    parent = get_parent(dm)
     if len(fields) == 1:
         # Subspace is just a single FunctionSpace.
         idx, = fields
-        subdm = W[idx].dm
+        subspace = W[idx]
         iset = W._ises[idx]
-        add_hook(parent, setup=partial(push_parent, subdm, parent), teardown=partial(pop_parent, subdm, parent),
-                 call_setup=True)
-
-        if ctx is not None:
-            ctx, = ctx.split([(idx, )])
-            add_hook(parent, setup=partial(push_appctx, subdm, ctx), teardown=partial(pop_appctx, subdm, ctx),
-                     call_setup=True)
-            add_hook(parent, setup=partial(push_ctx_coarsener, subdm, coarsen), teardown=partial(pop_ctx_coarsener, subdm, coarsen),
-                     call_setup=True)
-        return iset, subdm
     else:
         # Need to build an MFS for the subspace
         subspace = firedrake.MixedFunctionSpace([W[f] for f in fields])
 
-        add_hook(parent, setup=partial(push_parent, subspace.dm, parent), teardown=partial(pop_parent, subspace.dm, parent),
-                 call_setup=True)
         # Index set mapping from W into subspace.
-        iset = PETSc.IS().createGeneral(numpy.concatenate([W._ises[f].indices
-                                                           for f in fields]),
+        iset = PETSc.IS().createGeneral(numpy.concatenate([W.dof_dset.field_ises[f].indices for f in fields]),
                                         comm=W._comm)
-        if ctx is not None:
-            ctx, = ctx.split([fields])
-            add_hook(parent, setup=partial(push_appctx, subspace.dm, ctx),
-                     teardown=partial(pop_appctx, subspace.dm, ctx),
-                     call_setup=True)
-            add_hook(parent, setup=partial(push_ctx_coarsener, subspace.dm, coarsen),
-                     teardown=partial(pop_ctx_coarsener, subspace.dm, coarsen),
-                     call_setup=True)
-        return iset, subspace.dm
+
+    subdm = subspace.dm
+    parent = get_parent(dm)
+    add_hook(parent, setup=partial(push_parent, subdm, parent),
+             teardown=partial(pop_parent, subdm, parent),
+             call_setup=True)
+
+    ctx = get_appctx(dm)
+    coarsen = get_ctx_coarsener(dm)
+    if ctx is not None:
+        ctx, = ctx.split([fields])
+        add_hook(parent, setup=partial(push_appctx, subdm, ctx),
+                 teardown=partial(pop_appctx, subdm, ctx),
+                 call_setup=True)
+        add_hook(parent, setup=partial(push_ctx_coarsener, subdm, coarsen),
+                 teardown=partial(pop_ctx_coarsener, subdm, coarsen),
+                 call_setup=True)
+    return iset, subdm
 
 
 @PETSc.Log.EventDecorator()

firedrake/formmanipulation.py

@@ -1,4 +1,3 @@
-
 import numpy
 import collections
 
@@ -161,6 +160,7 @@ def cofunction(self, o):
             return Cofunction(W, val=MixedDat(o.dat[i] for i in indices))
 
     def matrix(self, o):
+        from firedrake.bcs import DirichletBC, EquationBC
         ises = []
         args = []
         for a in o.arguments():
@@ -180,8 +180,22 @@ def matrix(self, o):
             args.append(asplit)
 
         submat = o.petscmat.createSubMatrix(*ises)
-        bcs = ()
-        return AssembledMatrix(tuple(args), bcs, submat)
+        bcs = []
+        spaces = [a.function_space() for a in o.arguments()]
+        for bc in o.bcs:
+            W = bc.function_space()
+            W = W.parent or W
+
+            number = spaces.index(W)
+            V = args[number].function_space()
+            field = self.blocks[number]
+            if isinstance(bc, DirichletBC):
+                sub_bc = bc.reconstruct(field=field, V=V, g=bc.function_arg)
+            elif isinstance(bc, EquationBC):
+                raise NotImplementedError("Please get in touch if you need this")
+            if sub_bc is not None:
+                bcs.append(sub_bc)
+        return AssembledMatrix(tuple(args), tuple(bcs), submat)
 
     def zero_base_form(self, o):
         return ZeroBaseForm(tuple(map(self, o.arguments())))

firedrake/matrix.py

@@ -170,9 +170,7 @@ def __init__(self, a, bcs, mat_type, *args, **kwargs):
         self.mat_type = mat_type
 
     def assemble(self):
-        raise NotImplementedError("API compatibility to apply bcs after 'assemble(a)'\
-                                  has been removed.  Use 'assemble(a, bcs=bcs)', which\
-                                  now returns an assembled matrix.")
+        self.M.assemble()
 
 
 class ImplicitMatrix(MatrixBase):
@@ -250,3 +248,9 @@ def __init__(self, a, bcs, petscmat, *args, **kwargs):
 
     def mat(self):
         return self.petscmat
+
+    def assemble(self):
+        # Bump petsc matrix state by assembling it.
+        # Ensures that if the matrix changed, the preconditioner is
+        # updated if necessary.
+        self.petscmat.assemble()

firedrake/solving_utils.py

@@ -407,23 +407,33 @@ def split(self, fields):
                 F += problem.compute_bc_lifting(J, subu)
             else:
                 F = replace(F, {problem.u_restrict: u})
+
             if problem.Jp is not None:
                 Jp = splitter.split(problem.Jp, argument_indices=(field, field))
                 Jp = replace(Jp, {problem.u_restrict: u})
             else:
                 Jp = None
-            bcs = []
-            for bc in problem.bcs:
-                if isinstance(bc, DirichletBC):
-                    bc_temp = bc.reconstruct(field=field, V=V, g=bc.function_arg, sub_domain=bc.sub_domain)
-                elif isinstance(bc, EquationBC):
-                    bc_temp = bc.reconstruct(V, subu, u, field, problem.is_linear)
-                if bc_temp is not None:
-                    bcs.append(bc_temp)
+
+            if isinstance(J, MatrixBase) and J.has_bcs:
+                bcs = None
+            else:
+                bcs = []
+                for bc in problem.bcs:
+                    if isinstance(bc, DirichletBC):
+                        bc_temp = bc.reconstruct(field=field, V=V, g=bc.function_arg)
+                    elif isinstance(bc, EquationBC):
+                        bc_temp = bc.reconstruct(V, subu, u, field, problem.is_linear)
+                    if bc_temp is not None:
+                        bcs.append(bc_temp)
+
             new_problem = NLVP(F, subu, bcs=bcs, J=J, Jp=Jp, is_linear=problem.is_linear,
                                form_compiler_parameters=problem.form_compiler_parameters)
             new_problem._constant_jacobian = problem._constant_jacobian
-            splits.append(type(self)(new_problem, mat_type=self.mat_type, pmat_type=self.pmat_type,
+            splits.append(type(self)(new_problem,
+                                     mat_type=self.mat_type,
+                                     pmat_type=self.pmat_type,
+                                     sub_mat_type=self.sub_mat_type,
+                                     sub_pmat_type=self.sub_pmat_type,
                                      appctx=self.appctx,
                                      transfer_manager=self.transfer_manager,
                                      pre_apply_bcs=self.pre_apply_bcs))
@@ -504,6 +514,15 @@ def form_jacobian(snes, X, J, P):
         ctx.set_nullspace(ctx._nullspace_T, ises, transpose=True, near=False)
         ctx.set_nullspace(ctx._near_nullspace, ises, transpose=False, near=True)
 
+        # Bump petsc matrix state of each split by assembling them.
+        # Ensures that if the matrix changed, the preconditioner is
+        # updated if necessary.
+        for field, splits in ctx._splits.items():
+            for subctx in splits:
+                subctx._jac.assemble()
+                if subctx.Jp is not None:
+                    subctx._pjac.assemble()
+
     @staticmethod
     def compute_operators(ksp, J, P):
         r"""Form the Jacobian for this problem

tests/firedrake/regression/test_nested_fieldsplit_solves.py

@@ -152,6 +152,40 @@ def test_nested_fieldsplit_solve_parallel(W, A, b, expect):
     assert norm(f) < 1e-11
 
 
+def test_nonlinear_fielsplit():
+    mesh = UnitIntervalMesh(1)
+    V = FunctionSpace(mesh, "DG", 0)
+    Z = V * V * V
+
+    u = Function(Z)
+    u0, u1, u2 = split(u)
+    v0, v1, v2 = TestFunctions(Z)
+
+    F = inner(u0, v0) * dx
+    F += inner(0.5*u1**2 + u1, v1) * dx
+    F += inner(u2, v2) * dx
+    u.subfunctions[1].assign(Constant(1))
+
+    sp = {
+        "mat_type": "nest",
+        "snes_max_it": 10,
+        "ksp_type": "fgmres",
+        "pc_type": "fieldsplit",
+        "pc_fieldsplit_type": "additive",
+        "pc_fieldsplit_0_fields": "0",
+        "pc_fieldsplit_1_fields": "1,2",
+        "fieldsplit_1_ksp_view_eigenvalues": None,
+        "fieldsplit": {
+            "ksp_type": "gmres",
+            "pc_type": "jacobi",
+        },
+    }
+    J = derivative(F, u)
+    solver = NonlinearVariationalSolver(NonlinearVariationalProblem(F, u), solver_parameters=sp)
+    solver.solve()
+    assert np.allclose(solver.snes.ksp.pc.getFieldSplitSubKSP()[1].computeEigenvalues(), 1)
+
+
 def test_matrix_types(W):
     a = inner(TrialFunction(W), TestFunction(W))*dx