BUG: Fix EquationBC + matfree + fieldsplit (#4734)

* BUG: Fix EquationBC + matfree + fieldsplit

Author: pbrubeck

firedrake/assemble.py

@@ -1248,6 +1248,9 @@ def _apply_bc(self, tensor, bc, u=None):
                 bc.apply(r, u=u)
         elif isinstance(bc, EquationBCSplit):
             bc.zero(tensor)
+            if isinstance(bc.f, ufl.ZeroBaseForm) or bc.f.empty():
+                # form is empty, do nothing
+                return
             OneFormAssembler(bc.f, bcs=bc.bcs,
                              form_compiler_parameters=self._form_compiler_params,
                              needs_zeroing=False,

firedrake/bcs.py

@@ -641,9 +641,6 @@ def reconstruct(self, field=None, V=None, subu=None, u=None, row_field=None, col
             rank = len(self.f.arguments())
             splitter = ExtractSubBlock()
             form = splitter.split(self.f, argument_indices=(row_field, col_field)[:rank])
-            if isinstance(form, ufl.ZeroBaseForm) or form.empty():
-                # form is empty, do nothing
-                return
             if u is not None:
                 form = firedrake.replace(form, {self.u: u})
         if action_x is not None:

tests/firedrake/equation_bcs/test_equation_bcs.py

@@ -174,8 +174,8 @@ def linear_poisson_mixed(solver_parameters, mesh_num, porder):
     u1 = cos(2 * pi * y) / 2
     n = FacetNormal(mesh)
 
-    a = (inner(sigma, tau) + inner(u, div(tau)) + inner(div(sigma), v)) * dx
-    L = inner(u1, dot(tau, n)) * ds(1) + inner(f, v) * dx
+    a = inner(sigma, tau)*dx + inner(u, div(tau))*dx + inner(div(sigma), v + div(tau))*dx
+    L = inner(u1, dot(tau, n)) * ds(1) + inner(f, v + div(tau)) * dx
 
     g = as_vector([-2 * pi * sin(2 * pi * x + pi / 3) * cos(2 * pi * y), -2 * pi * cos(2 * pi * x + pi / 3) * sin(2 * pi * y)])
 
@@ -331,24 +331,24 @@ def test_EquationBC_mixedpoisson_matfree_fieldsplit():
     eq_type = "linear"
     porder = 0
     # Mixed poisson with EquationBCs
-    # matfree with fieldsplit pc
+    # matfree with fieldsplit schur direct solver
 
     solver_parameters = {'mat_type': mat_type,
-                         'ksp_type': 'fgmres',
-                         'ksp_atol': 1e-11,
-                         'ksp_max_it': 200,
+                         'ksp_type': 'preonly',
                          'pc_type': 'fieldsplit',
                          'pc_fieldsplit_type': 'schur',
                          'pc_fieldsplit_schur_fact_type': 'full',
-                         'fieldsplit_0_ksp_type': 'cg',
+                         'fieldsplit_0_ksp_type': 'preonly',
                          'fieldsplit_0_pc_type': 'python',
                          'fieldsplit_0_pc_python_type': 'firedrake.AssembledPC',
                          'fieldsplit_0_assembled_pc_type': 'lu',
-                         'fieldsplit_1_ksp_type': 'cg',
-                         'fieldsplit_1_pc_use_amat': True,
+                         'fieldsplit_0_assembled_pc_factor_mat_solver_type': 'mumps',
+                         'fieldsplit_1_ksp_type': 'gmres',
                          'fieldsplit_1_pc_type': 'python',
                          'fieldsplit_1_pc_python_type': 'firedrake.MassInvPC',
-                         'fieldsplit_1_Mp_pc_type': 'icc'}
+                         'fieldsplit_1_Mp_pc_type': 'lu',
+                         'fieldsplit_1_Mp_pc_factor_mat_solver_type': 'mumps',
+                         }
     err = []
     mesh_sizes = [16, 32]
     if eq_type == "linear":