coarsen EquationBC for LVP

Author: KarsKnook

firedrake/bcs.py

@@ -512,6 +512,8 @@ def __init__(self, *args, bcs=None, J=None, Jp=None, V=None, is_linear=False, Jp
                         raise ValueError("Provided BC RHS is not a linear form")
                     F = ufl_expr.action(J, u) - L
                 self.is_linear = True
+                self.lhs = J
+                self.rhs = L
             # nonlinear
             else:
                 if eq.rhs != 0:

firedrake/mg/ufl_utils.py

@@ -129,15 +129,20 @@ def coarsen_bc(bc, self, coefficient_mapping=None):
 
 @coarsen.register(firedrake.EquationBC)
 def coarsen_equation_bc(ebc, self, coefficient_mapping=None):
-    F = self(ebc._F.f, self, coefficient_mapping=coefficient_mapping)
     J = self(ebc._J.f, self, coefficient_mapping=coefficient_mapping)
     Jp = self(ebc._Jp.f, self, coefficient_mapping=coefficient_mapping)
     u = self(ebc._F.u, self, coefficient_mapping=coefficient_mapping)
     sub_domain = ebc._F.sub_domain
     bcs = [self(bc, self, coefficient_mapping=coefficient_mapping) for bc in ebc.bcs]
     V = self(ebc._F.function_space(), self, coefficient_mapping=coefficient_mapping)
 
-    return type(ebc)(F == 0, u, sub_domain, V=V, bcs=bcs, J=J, Jp=Jp)
+    if ebc.is_linear:
+        lhs = self(ebc.lhs, self, coefficient_mapping=coefficient_mapping)
+        rhs = self(ebc.rhs, self, coefficient_mapping=coefficient_mapping)
+        return type(ebc)(lhs == rhs, u, sub_domain, V=V, bcs=bcs, J=J, Jp=Jp)
+    else:
+        F = self(ebc._F.f, self, coefficient_mapping=coefficient_mapping)
+        return type(ebc)(F == 0, u, sub_domain, V=V, bcs=bcs, J=J, Jp=Jp)
 
 
 @coarsen.register(firedrake.functionspaceimpl.WithGeometryBase)

tests/firedrake/multigrid/test_equation_bc.py

@@ -1,21 +1,21 @@
 from firedrake import *
 
 
-def test_poisson():
+def test_poisson_NLVP():
     mesh = UnitIntervalMesh(10)
     mesh_hierarchy = MeshHierarchy(mesh, 1)
     mesh = mesh_hierarchy[-1]
 
     x = SpatialCoordinate(mesh)
-    u_exact = x[0]*(x[0]-1)
+    u_exact = x[0]**2
     f = -div(grad(u_exact))
 
     V = FunctionSpace(mesh, "CG", 1)
     u = Function(V)
     v = TestFunction(V)
 
     F = inner(grad(u), grad(v))*dx(degree=0) - inner(f, v)*dx
-    bcs = [EquationBC(inner(u, v) * ds == 0, u, "on_boundary", V=V)]
+    bcs = [EquationBC(inner(u - u_exact, v) * ds == 0, u, "on_boundary", V=V)]
     NLVP = NonlinearVariationalProblem(F, u, bcs=bcs)
 
     sp = {"pc_type": "mg"}
@@ -26,6 +26,32 @@ def test_poisson():
     assert errornorm(u_exact, u) < 5e-4
 
 
+def test_poisson_LVP():
+    mesh = UnitIntervalMesh(10)
+    mesh_hierarchy = MeshHierarchy(mesh, 1)
+    mesh = mesh_hierarchy[-1]
+
+    x = SpatialCoordinate(mesh)
+    u_exact = x[0]**2
+    f = -div(grad(u_exact))
+
+    V = FunctionSpace(mesh, "CG", 1)
+    u = TrialFunction(V)
+    u_ = Function(V)
+    v = TestFunction(V)
+
+    a = inner(grad(u), grad(v))*dx(degree=0)
+    L = inner(f, v)*dx
+    bcs = [EquationBC(inner(u, v) * ds == inner(u_exact, v) * ds, u_, "on_boundary", V=V)]
+    LVP = LinearVariationalProblem(a, L, u_, bcs=bcs)
+
+    sp = {"pc_type": "mg"}
+    LVS = LinearVariationalSolver(LVP, solver_parameters=sp)
+    LVS.solve()
+
+    assert errornorm(u_exact, u_) < 5e-4
+
+
 def test_nested_equation_bc():
     mesh = UnitCubeMesh(2, 2, 2)
     mesh_hierarchy = MeshHierarchy(mesh, 2)