Fw/hdiv

alfi/bary.py

@@ -4,6 +4,7 @@
 from mpi4py import MPI
 from fractions import Fraction
 from firedrake.cython import mgimpl as impl
+from firedrake.cython.dmcommon import FACE_SETS_LABEL
 from firedrake.petsc import *
 import numpy
 
@@ -474,7 +475,7 @@ def cdm_coords(p):
 
 def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, callbacks=None, reorder=None,
                       refinements_per_level=1):
-    cdm = mesh._plex
+    cdm = mesh._topology_dm
     cdm.setRefinementUniform(True)
     dms = []
     if mesh.comm.size > 1 and mesh._grown_halos:
@@ -509,7 +510,10 @@ def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, cal
         # Remove vertex (and edge) points from labels on exterior
         # facets.  Interior facets will be relabeled in Mesh
         # construction below.
-        impl.filter_exterior_facet_labels(rdm)
+        impl.filter_labels(rdm, rdm.getHeightStratum(1),
+                           "exterior_facets", "boundary_faces",
+                           FACE_SETS_LABEL)
+
         rdm.removeLabel("pyop2_core")
         rdm.removeLabel("pyop2_owned")
         rdm.removeLabel("pyop2_ghost")
@@ -526,10 +530,12 @@ def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, cal
             scale = mesh._radius / np.linalg.norm(coords, axis=1).reshape(-1, 1)
             coords *= scale
 
-    barydms = (bary(mesh._plex), ) + tuple(bary(dm) for dm in dms)
+    barydms = (bary(mesh._topology_dm), ) + tuple(bary(dm) for dm in dms)
 
     for bdm in barydms:
-        impl.filter_exterior_facet_labels(bdm)
+        impl.filter_labels(bdm, bdm.getHeightStratum(1),
+                           "exterior_facets", "boundary_faces",
+                           FACE_SETS_LABEL)
 
     barymeshes = [firedrake.Mesh(dm, dim=mesh.ufl_cell().geometric_dimension(),
                                  distribution_parameters=distribution_parameters,
@@ -543,10 +549,10 @@ def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, cal
 
     lgmaps = []
     for i, m in enumerate(meshes):
-        no = impl.create_lgmap(m._plex)
+        no = impl.create_lgmap(m._topology_dm)
         m.init()
-        o = impl.create_lgmap(m._plex)
-        m._plex.setRefineLevel(i)
+        o = impl.create_lgmap(m._topology_dm)
+        m._topology_dm.setRefineLevel(i)
         lgmaps.append((no, o))
 
     coarse_to_fine_cells = []
@@ -559,10 +565,10 @@ def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, cal
 
     lgmaps = []
     for i, m in enumerate(barymeshes):
-        no = impl.create_lgmap(m._plex)
+        no = impl.create_lgmap(m._topology_dm)
         m.init()
-        o = impl.create_lgmap(m._plex)
-        m._plex.setRefineLevel(i)
+        o = impl.create_lgmap(m._topology_dm)
+        m._topology_dm.setRefineLevel(i)
         lgmaps.append((no, o))
 
     d = mesh.topological_dimension()
@@ -574,8 +580,8 @@ def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, cal
                zip(lgmaps[:-1], lgmaps[1:]),
                coarse_to_fine_cells):
 
-        cdm = coarseu._plex
-        fdm = fineu._plex
+        cdm = coarseu._topology_dm
+        fdm = fineu._topology_dm
         _, cn2o = impl.get_entity_renumbering(cdm, coarseu._cell_numbering, "cell")
         _, fn2o = impl.get_entity_renumbering(fdm, fineu._cell_numbering, "cell")
         plex_uniform_coarse_to_fine = numpy.empty_like(uniform_coarse_to_fine)
@@ -595,8 +601,8 @@ def BaryMeshHierarchy(mesh, refinement_levels, distribution_parameters=None, cal
                 bary_cells.extend(list(range(fc*(d+1), (fc+1)*(d+1))))
             plex_coarse_bary_to_fine_bary[c] = bary_cells
 
-        cdm = coarse._plex
-        fdm = fine._plex
+        cdm = coarse._topology_dm
+        fdm = fine._topology_dm
 
         co2n, _ = impl.get_entity_renumbering(cdm, coarse._cell_numbering, "cell")
         fo2n, _ = impl.get_entity_renumbering(fdm, fine._cell_numbering, "cell")

alfi/driver.py

@@ -1,4 +1,4 @@
-from alfi.solver import ConstantPressureSolver, ScottVogeliusSolver
+from alfi.solver import ConstantPressureSolver, ScottVogeliusSolver, RTSolver, BDMSolver
 from mpi4py import MPI
 from firedrake.petsc import PETSc
 from firedrake import *
@@ -26,7 +26,7 @@ def get_default_parser():
     parser.add_argument("--stabilisation-type", type=str, default=None,
                         choices=["none", "burman", "gls", "supg"])
     parser.add_argument("--discretisation", type=str, required=True,
-                        choices=["pkp0", "sv"])
+                        choices=["pkp0", "sv", "rt", "bdm"])
     parser.add_argument("--gamma", type=float, default=1e4)
     parser.add_argument("--clear", dest="clear", default=False,
                         action="store_true")
@@ -50,7 +50,10 @@ def get_default_parser():
 
 def get_solver(args, problem, hierarchy_callback=None):
     solver_t = {"pkp0": ConstantPressureSolver,
-                "sv": ScottVogeliusSolver}[args.discretisation]
+                "sv": ScottVogeliusSolver,
+                "rt": RTSolver,
+                "bdm": BDMSolver}[args.discretisation]
+
     solver = solver_t(
         problem,
         solver_type=args.solver_type,

alfi/solver.py

@@ -199,6 +199,7 @@ def visprolong(u):
         self.z_last = z.copy(deepcopy=True)
 
 
+        self.bcs = bcs
         F = self.residual()
 
         """ Stabilisation """
@@ -238,7 +239,6 @@ def visprolong(u):
 
         appctx = {"nu": self.nu, "gamma": self.gamma}
         problem = NonlinearVariationalProblem(F, z, bcs=bcs)
-        self.bcs = bcs
         self.params = params
         self.nsp = nsp
         self.appctx = appctx
@@ -317,6 +317,7 @@ def get_parameters(self):
             "ksp_convergence_test": "skip",
             "pc_type": "python",
             "pc_python_type": "firedrake.PatchPC",
+            # "pc_python_type": "matpatch.MatPatch",
             "patch_pc_patch_save_operators": True,
             "patch_pc_patch_partition_of_unity": False,
             "patch_pc_patch_local_type": "multiplicative" if multiplicative else "additive",
@@ -660,3 +661,93 @@ def configure_patch_solver(self, opts):
 
     def distribution_parameters(self):
         return {"partition": True, "overlap_type": (DistributedMeshOverlapType.VERTEX, 2)}
+
+
+class HdivSolver(NavierStokesSolver):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def residual(self):
+        u, p = split(self.z)
+        v, q = TestFunctions(self.Z)
+        sigma = Constant(10) * self.Z.sub(0).ufl_element().degree()**2
+        h = CellDiameter(self.z.ufl_domain())
+        n = FacetNormal(self.z.ufl_domain())
+        ulast = split(self.z_last)[0]
+
+        def a(u, v):
+            return self.nu * inner(2*sym(grad(u)), grad(v))*dx \
+                - self.nu * inner(avg(2*sym(grad(u))), 2*avg(outer(v, n))) * dS \
+                - self.nu * inner(avg(2*sym(grad(v))), 2*avg(outer(u, n))) * dS \
+                + self.nu * sigma/avg(h) * inner(2*avg(outer(u,n)),2*avg(outer(v,n))) * dS
+
+        def a_bc(u, v, bid, g): 
+            nu = self.nu
+            return -inner(outer(v,n),2*nu*sym(grad(u)))*ds(bid) \
+                -inner(outer(u-g,n),2*nu*sym(grad(v)))*ds(bid) \
+                +nu*(sigma/h)*inner(v,u-g)*ds(bid)
+
+        def c(u, v):
+            uflux_int = 0.5*(dot(u, n) + abs(dot(u, n)))*u
+            return -self.advect * dot(u ,div(outer(v,u)))*dx \
+                +self.advect * dot(v('+')-v('-'), uflux_int('+')-uflux_int('-'))*dS
+
+        def c_bc(u, v, bid, g):
+            if g is None:
+                uflux_ext = 0.5*(inner(u,n)+abs(inner(u,n)))*u
+            else:
+                uflux_ext = 0.5*(inner(u,n)+abs(inner(u,n)))*u + 0.5*(inner(u,n)-abs(inner(u,n)))*g
+            return self.advect * dot(v,uflux_ext)*ds(bid)
+
+        def b(u, q):
+            return div(u) * q * dx
+
+        F = a(u, v) + c(u, v) - b(u, q) - b(v, p) + self.gamma * inner(div(u), div(v)) * dx
+
+        exterior_markers = list(self.mesh.exterior_facets.unique_markers)
+        for bc in self.bcs:
+            if "DG" in str(bc._function_space):
+                continue
+            g = bc.function_arg
+            bid = bc.sub_domain
+            exterior_markers.remove(bid)
+            F += a_bc(u, v, bid, g)
+            F += c_bc(u, v, bid, g)
+        for bid in exterior_markers:
+            F += c_bc(u, v, bid, None)
+        return F
+
+    def get_transfers(self):
+        V = self.Z.sub(0)
+        dgtransfer = DGInjection()
+        transfers = {VectorElement("DG", V.mesh().ufl_cell(), V.ufl_element().degree()): (dgtransfer.prolong, restrict, dgtransfer.inject)}
+        return transfers
+
+    def configure_patch_solver(self, opts):
+        opts["patch_pc_patch_sub_mat_type"] = "seqdense"
+        opts["patch_sub_pc_factor_mat_solver_type"] = "petsc"
+        opts["pc_python_type"] = "matpatch.MatPatch"
+
+    def distribution_parameters(self):
+        return {"partition": True, "overlap_type": (DistributedMeshOverlapType.VERTEX, 1)}
+
+
+class RTSolver(HdivSolver):
+
+    def function_space(self, mesh, k):
+        eleu = FiniteElement("RT", mesh.ufl_cell(), k, variant="integral")
+        elep = FiniteElement("Discontinuous Lagrange", mesh.ufl_cell(), k-1)
+        V = FunctionSpace(mesh, eleu)
+        Q = FunctionSpace(mesh, elep)
+        return MixedFunctionSpace([V, Q])
+
+
+class BDMSolver(HdivSolver):
+
+    def function_space(self, mesh, k):
+        eleu = FiniteElement("BDM", mesh.ufl_cell(), k, variant="integral")
+        elep = FiniteElement("Discontinuous Lagrange", mesh.ufl_cell(), k-1)
+        V = FunctionSpace(mesh, eleu)
+        Q = FunctionSpace(mesh, elep)
+        return MixedFunctionSpace([V, Q])

examples/Makefile

@@ -1,4 +1,4 @@
-exec=mpirun -n 12 python3
+exec=mpiexec -n 3 python3
 pre:
 	mkdir -p logs/
 
@@ -9,12 +9,12 @@ iters2dpkp0: pre
 	$(exec) iters.py --problem bfs2d --nref-start 1 --nref-end 3 --mesh coarse03.msh --discretisation pkp0 --mh uniform --k 2 --solver-type almg --stabilisation-type supg --patch star --smoothing 6 --restriction 2>&1 | tee logs/iters2dbfs03addrestrict6.log
 
 iters2dsv: pre
-	$(exec) iters.py --nref-start 1 --nref-end 4 --problem ldc2d --k 2 --baseN 10          --solver-type almg --discretisation sv   --mh bary    --stabilisation-type burman --patch macro  --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvldcrestrict5em3k2.log
-	$(exec) iters.py --nref-start 1 --nref-end 4 --problem ldc2d --k 3 --baseN 10          --solver-type almg --discretisation sv   --mh bary    --stabilisation-type burman --patch macro  --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvldcrestrict5em3k3.log
-	$(exec) iters.py --nref-start 1 --nref-end 4 --problem ldc2d --k 2 --baseN 16          --solver-type almg --discretisation pkp0 --mh uniform --stabilisation-type supg   --patch star   --smoothing 6 --restriction 2>&1 | tee logs/iters2dsvldccompp2p0restrict.log
-	$(exec) iters.py --nref-start 1 --nref-end 4 --problem bfs2d --k 2 --mesh coarse09.msh --solver-type almg --discretisation sv --mh bary --stabilisation-type burman --patch macro --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvbfs09restrict5em3k2.log
-	$(exec) iters.py --nref-start 1 --nref-end 3 --problem bfs2d --k 3 --mesh coarse09.msh --solver-type almg --discretisation sv --mh bary --stabilisation-type burman --patch macro --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvbfs09restrict5em3k3.log
-	
+	$(exec) iters.py --nref-start 1 --nref-end 1 --problem ldc2d --k 2 --baseN 10          --solver-type almg --discretisation sv   --mh bary    --stabilisation-type burman --patch macro  --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvldcrestrict5em3k2.log
+	$(exec) iters.py --nref-start 1 --nref-end 1 --problem ldc2d --k 3 --baseN 10          --solver-type almg --discretisation sv   --mh bary    --stabilisation-type burman --patch macro  --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvldcrestrict5em3k3.log
+	$(exec) iters.py --nref-start 1 --nref-end 1 --problem ldc2d --k 2 --baseN 16          --solver-type almg --discretisation pkp0 --mh uniform --stabilisation-type supg   --patch star   --smoothing 6 --restriction 2>&1 | tee logs/iters2dsvldccompp2p0restrict.log
+	$(exec) iters.py --nref-start 1 --nref-end 1 --problem bfs2d --k 2 --mesh coarse09.msh --solver-type almg --discretisation sv --mh bary --stabilisation-type burman --patch macro --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvbfs09restrict5em3k2.log
+	$(exec) iters.py --nref-start 1 --nref-end 1 --problem bfs2d --k 3 --mesh coarse09.msh --solver-type almg --discretisation sv --mh bary --stabilisation-type burman --patch macro --smoothing 6 --restriction --stabilisation-weight 5e-3 2>&1 | tee logs/iters2dsvbfs09restrict5em3k3.log
+
 
 mms3dpkp0: pre
 	$(exec) mms.py --dim 3 --k 1 --nref 5 --discretisation pkp0 --solver-type almg --mh uniform --patch star --gamma 1e4 --stabilisation-type supg --stabilisation--weight 0.05 --baseN 4 2>&1 | tee logs/mms3dpkp0highgamma.log
@@ -31,3 +31,11 @@ idealal:
 	$(exec) iters.py --baseN 16 --nref-start 2 --nref-end 2 --problem ldc2d --k 2 --solver-type allu --discretisation pkp0 --mh uniform --stabilisation-type supg --smoothing 6 --restriction --time --re-max 10000 --gamma 1e2 2>&1 | tee logs/idealalgamma1e2.log || 1
 	$(exec) iters.py --baseN 16 --nref-start 2 --nref-end 2 --problem ldc2d --k 2 --solver-type allu --discretisation pkp0 --mh uniform --stabilisation-type supg --smoothing 6 --restriction --time --re-max 10000 --gamma 1e3 2>&1 | tee logs/idealalgamma1e3.log || 1
 	$(exec) iters.py --baseN 16 --nref-start 2 --nref-end 2 --problem ldc2d --k 2 --solver-type allu --discretisation pkp0 --mh uniform --stabilisation-type supg --smoothing 6 --restriction --time --re-max 10000 --gamma 1e4 2>&1 | tee logs/idealalgamma1e4.log || 1
+
+mms2dhdiv: pre
+	$(exec) mms.py --dim 2 --nref 2 --k 2 --discretisation rt --solver-type allu --mh uniform --patch star --gamma 1e3 2>&1 | tee logs/mms2dhdivrt.log
+	$(exec) mms.py --dim 2 --nref 2 --k 2 --discretisation bdm --solver-type allu --mh uniform --patch star --gamma 1e3 2>&1 | tee logs/mms2dhdivbdm.log
+
+mms3dhdiv: pre
+	$(exec) mms.py --dim 3 --nref 2 --k 2 --discretisation rt --solver-type allu --mh uniform --patch star --gamma 1e3 2>&1 | tee logs/mms2dhdivrt.log
+	$(exec) mms.py --dim 3 --nref 2 --k 2 --discretisation bdm --solver-type allu --mh uniform --patch star --gamma 1e3 2>&1 | tee logs/mms2dhdivbdm.log

examples/bfs2d/bfs2d.py

@@ -18,7 +18,7 @@ def mesh(self, distribution_parameters):
     @staticmethod
     def poiseuille_flow(domain):
         (x, y) = SpatialCoordinate(domain)
-        return as_vector([4 * (2-y)*(y-1)*(y>1), 0])
+        return as_vector([conditional(ge(y, 1), 4 * (2-y)*(y-1), 0), 0])
 
     def bcs(self, Z):
         bcs = [DirichletBC(Z.sub(0), self.poiseuille_flow(Z.mesh()), 1),
@@ -43,9 +43,11 @@ def relaxation_direction(self): return "0+:1-"
     problem = TwoDimBackwardsFacingStepProblem(args.mesh)
     solver = get_solver(args, problem)
 
-    start = 250
-    end = 10000
+    start = 500
+    end = 1000
     step = 250
     res = [0, 1, 10, 100] + list(range(start, end+step, step))
-    res = [1, 10, 100] # + list(range(start, end+step, step))
+    res = [0, 1, 3, 10, 20, 100, 150, 200, 250, 300, 350, 400, 500] + list(range(start, end+step, step))
     results = run_solver(solver, res, args)
+    u = solver.z.split()[0]
+    print("div(u)",norm(div(u)))

examples/ldc2d/ldc2d.py

@@ -21,7 +21,9 @@ def mesh(self, distribution_parameters):
 
     def bcs(self, Z):
         bcs = [DirichletBC(Z.sub(0), self.driver(Z.ufl_domain()), 4),
-               DirichletBC(Z.sub(0), Constant((0., 0.)), [1, 2, 3])]
+               DirichletBC(Z.sub(0), Constant((0, 0)), 1),
+               DirichletBC(Z.sub(0), Constant((0, 0)), 2),
+               DirichletBC(Z.sub(0), Constant((0, 0)), 3)]
         return bcs
 
     def has_nullspace(self): return True
@@ -48,10 +50,14 @@ def relaxation_direction(self): return "0+:1-"
     args, _ = parser.parse_known_args()
     problem = TwoDimLidDrivenCavityProblem(args.baseN, args.diagonal)
     solver = get_solver(args, problem)
-
-    start = 250
+    start = 500
     end = 10000
     step = 250
     res = [0, 1, 10, 100] + list(range(start, end+step, step))
-    res = [1, 10, 50, 100, 150, 200]# + list(range(start, end+step, step))
+    res = [0, 1, 10, 50, 100, 150, 199, 200]# + list(range(start, end+step, step))
+    res = [0, 1, 0]
     results = run_solver(solver, res, args)
+    u = solver.z.split()[0]
+    print("div(u)",norm(div(u)))
+    divu = Function(solver.Z.sub(1)).project(div(u))
+    File("divu.pvd").write(divu)

examples/mms.py

@@ -5,8 +5,8 @@
 import os
 from firedrake import *
 import numpy as np
-import inflect
-eng = inflect.engine()
+# import inflect
+# eng = inflect.engine()
 
 convergence_orders = lambda x: np.log2(np.array(x)[:-1] / np.array(x)[1:])
 
@@ -22,7 +22,7 @@
 else:
     raise NotImplementedError
 
-res = [1, 9, 10, 50, 90, 100, 400, 500, 900, 1000]
+res = [1]#, 9, 10, 50, 90, 100, 400, 500, 900, 1000]
 results = {}
 for re in res:
     results[re] = {}
@@ -49,14 +49,16 @@
         u, p = z.split()
         Z = z.function_space()
 
-        # uviz = solver.visprolong(u)
-        # (u_, p_) = problem.actual_solution(uviz.function_space())
+        uviz = solver.visprolong(u)
+        (u_, p_) = problem.actual_solution(uviz.function_space())
+        warning("velocity error: %f" % norm(u-u_))
         # File("output/u-re-%i-nref-%i.pvd" % (re, nref)).write(uviz.interpolate(uviz))
-        # File("output/uerr-re-%i-nref-%i.pvd" % (re, nref)).write(uviz.interpolate(uviz-u_))
+        # File("output/uerr-re-%i-nref-%i.pvd" % (re, nref)).write(uviz.project(uviz-u_))
         # File("output/uex-re-%i-nref-%i.pvd" % (re, nref)).write(uviz.interpolate(u_))
         (u_, p_) = problem.actual_solution(Z)
         # File("output/perr-re-%i-nref-%i.pvd" % (re, nref)).write(Function(Z.sub(1)).interpolate(p-p_))
         veldiv = norm(div(u))
+        veldivex = norm(div(u_))
         pressureintegral = assemble(p_ * dx)
         uerr = norm(u_-u)
         ugraderr = norm(grad(u_-u))
@@ -74,9 +76,9 @@
         results[re]["divergence"].append(veldiv)
         if comm.rank == 0:
             print("|div(u_h)| = ", veldiv)
-            print("p_exact * dx = ", pressureintegral)
-            print("p_approx * dx = ", pintegral)
-
+            print("|div(u)| = ", veldivex)
+            print("p_h * dx = ", pintegral)
+            print("p * dx = ", pressureintegral)
 if comm.rank == 0:
     for re in res:
         print("Results for Re =", re)
@@ -86,6 +88,7 @@
         print("convergence orders:", convergence_orders(results[re]["pressure"]))
     print("gamma =", args.gamma)
     print("h =", hs)
+    import sys; sys.exit()
 
     for re in [10, 100, 500, 1000]:
         print("%%Re = %i" % re)