Merge pull request #1257 from firedrakeproject/pefarrell/la-solve-mg

* pefarrell/la-solve-mg:
  solving: Simplify assembled solve LVP setup
  solving: Allow literal zero for RHS of LVP
  Enable the use of MG etc for the algebraic solve interface. Fixes #1256.

firedrake/solving.py

@@ -23,6 +23,8 @@
 
 import firedrake.linear_solver as ls
 import firedrake.variational_solver as vs
+from firedrake import solving_utils
+from firedrake import dmhooks
 
 
 def solve(*args, **kwargs):
@@ -221,7 +223,19 @@ def _la_solve(A, x, b, **kwargs):
                              near_nullspace=near_nullspace,
                              options_prefix=options_prefix)
 
-    solver.solve(x, b)
+    # linear MG doesn't need RHS, supply zero.
+    lvp = vs.LinearVariationalProblem(a=A.a, L=0, u=x, bcs=bcs)
+    mat_type = solver_parameters.get("mat_type")
+    appctx = solver_parameters.get("appctx", {})
+    ctx = solving_utils._SNESContext(lvp,
+                                     mat_type=mat_type,
+                                     pmat_type=mat_type,
+                                     appctx=appctx,
+                                     options_prefix=options_prefix)
+    dm = solver.ksp.dm
+
+    with dmhooks.appctx(dm, ctx):
+        solver.solve(x, b)
 
 
 def _extract_linear_solver_args(*args, **kwargs):

firedrake/variational_solver.py

@@ -266,12 +266,14 @@ def __init__(self, a, L, u, bcs=None, aP=None,
         # In the linear case, the Jacobian is the equation LHS.
         J = a
         # Jacobian is checked in superclass, but let's check L here.
-        if not isinstance(L, (ufl.Form, slate.slate.TensorBase)):
-            raise TypeError("Provided RHS is a '%s', not a Form or Slate Tensor" % type(L).__name__)
-        if len(L.arguments()) != 1:
-            raise ValueError("Provided RHS is not a linear form")
-
-        F = ufl_expr.action(J, u) - L
+        if L is 0:
+            F = ufl_expr.action(J, u)
+        else:
+            if not isinstance(L, (ufl.Form, slate.slate.TensorBase)):
+                raise TypeError("Provided RHS is a '%s', not a Form or Slate Tensor" % type(L).__name__)
+            if len(L.arguments()) != 1:
+                raise ValueError("Provided RHS is not a linear form")
+            F = ufl_expr.action(J, u) - L
 
         super(LinearVariationalProblem, self).__init__(F, u, bcs, J, aP,
                                                        form_compiler_parameters=form_compiler_parameters)