Default to pre_apply_bcs=False

firedrake/mg/ufl_utils.py

@@ -199,10 +199,6 @@ def inject_on_restrict(fine, restriction, rscale, injection, coarse):
                     manager.restrict(c, cmapping[c])
                 else:
                     manager.inject(c, cmapping[c])
-            # Apply bcs
-            if cctx.pre_apply_bcs:
-                for bc in cctx._problem.dirichlet_bcs():
-                    bc.apply(cctx._x)
             # When the solution is in the real space
             # PETSc fails to call this hook on coarse levels.
             # As a workaround, we inject into all levels.

firedrake/solving.py

@@ -270,7 +270,7 @@ def _extract_linear_solver_args(*args, **kwargs):
     nullspace_T = kwargs.get("transpose_nullspace", None)
     near_nullspace = kwargs.get("near_nullspace", None)
     options_prefix = kwargs.get("options_prefix", None)
-    pre_apply_bcs = kwargs.get("pre_apply_bcs", True)
+    pre_apply_bcs = kwargs.get("pre_apply_bcs", False)
 
     return P, bcs, solver_parameters, nullspace, nullspace_T, near_nullspace, options_prefix, pre_apply_bcs
 
@@ -326,7 +326,7 @@ def _extract_args(*args, **kwargs):
     solver_parameters = kwargs.get("solver_parameters", {})
     options_prefix = kwargs.get("options_prefix", None)
     restrict = kwargs.get("restrict", False)
-    pre_apply_bcs = kwargs.get("pre_apply_bcs", True)
+    pre_apply_bcs = kwargs.get("pre_apply_bcs", False)
 
     return eq, u, bcs, J, Jp, M, form_compiler_parameters, \
         solver_parameters, nullspace, nullspace_T, near_nullspace, \

firedrake/solving_utils.py

@@ -168,9 +168,6 @@ class _SNESContext(object):
     transfer_manager
         Object that can transfer functions between levels,
         typically a :class:`~.TransferManager`.
-    pre_apply_bcs
-        If `False`, the problem is linearised around the initial guess before
-        imposing the boundary conditions.
 
     The idea here is that the SNES holds a shell DM which contains
     this object as "user context".  When the SNES calls back to the
@@ -188,8 +185,7 @@ def __init__(self, problem,
                  pre_jacobian_callback=None, pre_function_callback=None,
                  post_jacobian_callback=None, post_function_callback=None,
                  options_prefix: str | None = None,
-                 transfer_manager=None,
-                 pre_apply_bcs: bool = True):
+                 transfer_manager=None):
         from firedrake.assemble import get_assembler
 
         if pmat_type is None:
@@ -201,7 +197,6 @@ def __init__(self, problem,
         self.sub_mat_type = sub_mat_type
         self.sub_pmat_type = sub_pmat_type
         self.options_prefix = options_prefix
-        self.pre_apply_bcs = pre_apply_bcs
 
         matfree = mat_type == 'matfree'
         pmatfree = pmat_type == 'matfree'
@@ -251,19 +246,26 @@ def __init__(self, problem,
         self.bcs_Jp = tuple(bc.extract_form('Jp') for bc in problem.bcs)
 
         self._bc_residual = None
-        if not pre_apply_bcs and next(problem.dirichlet_bcs(), None) is not None:
+        if not problem.restrict and next(problem.dirichlet_bcs(), None) is not None:
             # Delayed lifting of DirichletBCs
             self._bc_residual = Function(self._x.function_space())
-            if problem.is_linear:
-                # Drop existing lifting term from the residual
-                assert isinstance(self.F, ufl.BaseForm)
-                self.F = ufl.replace(self.F, {self._x: ufl.zero(self._x.ufl_shape)})
-
-            self.F -= problem.compute_bc_lifting(self.J, self._bc_residual)
+            if problem.is_linear and hasattr(problem, "L"):
+                # So far, F = L - action(J, u), but we need F = L - action(J, u - g).
+                # We could simply subtract action(J, g) as we do for the nonlinear case,
+                # but this doubles the number of integrals to be compiled.
+                # We cannot replace u -> u - g at this point because J or L might depend on u.
+                # The most efficient solution is to reconstruct the lifted residual from scratch.
+                # However, compute_bc_lifting(J, u - g, L) will complain about action not taking
+                # a pure Coefficient/Argument, so we supply a TrialFunction and replace it with u-g afterwards.
+                test, trial = self.J.arguments()
+                Ftrial = problem.compute_bc_lifting(self.J, trial, L=problem.L)
+                self.F = ufl.replace(Ftrial, {trial: self._x - self._bc_residual})
+            else:
+                self.F -= problem.compute_bc_lifting(self.J, self._bc_residual)
 
         self._assemble_residual = get_assembler(self.F, bcs=self.bcs_F,
                                                 form_compiler_parameters=self.fcp,
-                                                zero_bc_nodes=pre_apply_bcs,
+                                                zero_bc_nodes=self._problem.restrict,
                                                 ).assemble
 
         self._jacobian_assembled = False
@@ -289,7 +291,6 @@ def reconstruct(self, problem=None, mat_type=None, pmat_type=None, **kwargs):
             "appctx": self.appctx,
             "options_prefix": self.options_prefix,
             "transfer_manager": self.transfer_manager,
-            "pre_apply_bcs": self.pre_apply_bcs,
         }
         for k, v in default_options.items():
             if kwargs.get(k) is None:
@@ -467,7 +468,7 @@ def form_function(snes, X, F):
         if ctx._pre_function_callback is not None:
             ctx._pre_function_callback(X)
 
-        if not ctx.pre_apply_bcs:
+        if not ctx._problem.restrict:
             # Compute DirichletBC residual
             for bc in ctx._problem.dirichlet_bcs():
                 bc.apply(ctx._bc_residual, u=ctx._x)

firedrake/variational_solver.py

@@ -5,6 +5,7 @@
 from contextlib import ExitStack
 from types import MappingProxyType
 from petsctools import OptionsManager, flatten_parameters
+from warnings import warn
 
 from firedrake import dmhooks, slate, solving, solving_utils, ufl_expr, utils
 from firedrake.petsc import PETSc, DEFAULT_KSP_PARAMETERS, DEFAULT_SNES_PARAMETERS
@@ -90,9 +91,9 @@ def __init__(self, F, u, bcs=None, J=None,
                 bcs = J.bcs
         if bcs and any(isinstance(bc, EquationBC) for bc in bcs):
             restrict = False
-        self.restrict = restrict and bcs
+        self.restrict = restrict and bool(bcs)
 
-        if restrict and bcs:
+        if self.restrict:
             V_res = restricted_function_space(V, extract_subdomain_ids(bcs))
             bcs = [bc.reconstruct(V=V_res, indices=bc._indices) for bc in bcs]
             self.u_restrict = Function(V_res)
@@ -194,7 +195,7 @@ def __init__(self, problem, *, solver_parameters=None,
                  post_jacobian_callback=None,
                  pre_function_callback=None,
                  post_function_callback=None,
-                 pre_apply_bcs=True):
+                 pre_apply_bcs=False):
         r"""
         :arg problem: A :class:`NonlinearVariationalProblem` to solve.
         :kwarg nullspace: an optional :class:`.VectorSpaceBasis` (or
@@ -293,8 +294,10 @@ def update_diffusivity(current_solution):
                                          pre_function_callback=pre_function_callback,
                                          post_jacobian_callback=post_jacobian_callback,
                                          post_function_callback=post_function_callback,
-                                         options_prefix=self.options_prefix,
-                                         pre_apply_bcs=pre_apply_bcs)
+                                         options_prefix=self.options_prefix)
+        if pre_apply_bcs:
+            warn("Setting pre_apply_bcs=True is deprecated.", DeprecationWarning, stacklevel=2)
+        self.pre_apply_bcs = pre_apply_bcs
 
         self.snes = PETSc.SNES().create(comm=problem.dm.comm)
 
@@ -376,7 +379,7 @@ def solve(self, bounds=None):
             # Transfer the initial guess into the RestrictedFunctionSpace
             problem.u_restrict.assign(problem.u)
 
-        if self._ctx.pre_apply_bcs:
+        if self.pre_apply_bcs or problem.restrict:
             for bc in problem.dirichlet_bcs():
                 bc.apply(problem.u_restrict)
 
@@ -440,6 +443,7 @@ def __init__(self, a, L, u, bcs=None, aP=None,
         elif L != 0:
             raise TypeError(f"Provided RHS is a '{type(L).__name__}', not a Form or Slate Tensor")
         F = self.compute_bc_lifting(a, u, L=L)
+        self.L = L
 
         super(LinearVariationalProblem, self).__init__(F, u, bcs=bcs, J=a, Jp=aP,
                                                        form_compiler_parameters=form_compiler_parameters,

tests/firedrake/regression/test_solving_interface.py

@@ -230,7 +230,13 @@ def mesh(request):
     return UnitIntervalMesh(10)
 
 
-def test_solve_cofunction_rhs(mesh):
+@pytest.mark.parametrize("pre_apply_bcs", (False, True, "restrict"))
+def test_solve_cofunction_rhs(mesh, pre_apply_bcs):
+    restrict = False
+    if pre_apply_bcs == "restrict":
+        pre_apply_bcs = True
+        restrict = True
+
     V = FunctionSpace(mesh, "CG", 1)
     x, = SpatialCoordinate(mesh)
 
@@ -245,7 +251,7 @@ def test_solve_cofunction_rhs(mesh):
     Lold = L.copy()
 
     w = Function(V)
-    solve(a == L, w, bcs=bcs)
+    solve(a == L, w, bcs=bcs, pre_apply_bcs=pre_apply_bcs, restrict=restrict)
     assert errornorm(x, w) < 1E-10
     assert np.allclose(L.dat.data, Lold.dat.data)
 
@@ -359,12 +365,12 @@ def test_solve_pre_apply_bcs(mesh, mixed):
     # Raises NaNs if pre_apply_bcs=True
     F = derivative(W, z)
     z.zero()
-    solve(F == 0, z, bcs, pre_apply_bcs=False)
+    solve(F == 0, z, bcs)
     assert errornorm(g, uh) < 1E-10
 
     # Test that pre_apply_bcs=False works with a linear problem
     a = derivative(F, z)
     L = Form([])
     z.zero()
-    solve(a == L, z, bcs, pre_apply_bcs=False)
+    solve(a == L, z, bcs)
     assert errornorm(g, uh) < 1E-10