From c2c4541f8d9747eba316072697e4fce8e0987827 Mon Sep 17 00:00:00 2001 From: Florian Wechsung Date: Mon, 26 Aug 2019 17:52:36 +0100 Subject: [PATCH 1/5] hdiv-l2 is working for bfs2d but not for ldc2d... strange --- alfi/driver.py | 9 +++-- alfi/solver.py | 97 ++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 100 insertions(+), 6 deletions(-) diff --git a/alfi/driver.py b/alfi/driver.py index bca6ac1..b4cd2d6 100644 --- a/alfi/driver.py +++ b/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, diff --git a/alfi/solver.py b/alfi/solver.py index 245abfd..c92010f 100644 --- a/alfi/solver.py +++ b/alfi/solver.py @@ -200,6 +200,7 @@ def visprolong(u): self.z_last = z.copy(deepcopy=True) + self.bcs = bcs F = self.residual() """ Stabilisation """ @@ -239,7 +240,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 @@ -247,7 +247,8 @@ def visprolong(u): nullspace=nsp, options_prefix="ns_", appctx=appctx) self.transfers = self.get_transfers() - self.solver.set_transfer_operators(*self.transfers) + if self.transfers is not None and len(self.transfers) > 0: + self.solver.set_transfer_operators(*self.transfers) self.check_nograddiv_residual = True if self.check_nograddiv_residual: self.message(GREEN % "Checking residual without grad-div term") @@ -512,7 +513,8 @@ def setup_adjoint(self, J): options_prefix="ns_adj", appctx=self.appctx) - solver.set_transfer_operators(*self.transfers) + if self.transfers is not None and len(self.transfers) > 0: + solver.set_transfer_operators(*self.transfers) self.solver_adjoint = solver def load_balance(self, mesh): @@ -645,3 +647,92 @@ 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 = FacetArea(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) + prolongation = EmbeddedDGTransfer(V.ufl_element()) + return [dmhooks.transfer_operators(V, prolong=prolongation.prolong, inject=prolongation.inject, restrict=prolongation.restrict)] + + 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) + elep = FiniteElement("Discontinuous Lagrange L2", 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) + elep = FiniteElement("Discontinuous Lagrange L2", mesh.ufl_cell(), k-1) + V = FunctionSpace(mesh, eleu) + Q = FunctionSpace(mesh, elep) + return MixedFunctionSpace([V, Q]) From 970fd3d4734894952b2eb40e62f69a6db9309972 Mon Sep 17 00:00:00 2001 From: Florian Wechsung Date: Mon, 7 Oct 2019 18:11:30 +0100 Subject: [PATCH 2/5] wip --- alfi/solver.py | 5 +++-- examples/Makefile | 15 +++++++++------ examples/bfs2d/bfs2d.py | 10 ++++++---- examples/ldc2d/ldc2d.py | 14 ++++++++++---- examples/mms.py | 21 ++++++++++++--------- examples/mmsldc2d/mmsldc2d.py | 16 ++++++++++------ 6 files changed, 50 insertions(+), 31 deletions(-) diff --git a/alfi/solver.py b/alfi/solver.py index c92010f..8ac41f6 100644 --- a/alfi/solver.py +++ b/alfi/solver.py @@ -319,6 +319,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", @@ -722,7 +723,7 @@ class RTSolver(HdivSolver): def function_space(self, mesh, k): eleu = FiniteElement("RT", mesh.ufl_cell(), k) - elep = FiniteElement("Discontinuous Lagrange L2", mesh.ufl_cell(), k-1) + elep = FiniteElement("Discontinuous Lagrange", mesh.ufl_cell(), k-1) V = FunctionSpace(mesh, eleu) Q = FunctionSpace(mesh, elep) return MixedFunctionSpace([V, Q]) @@ -732,7 +733,7 @@ class BDMSolver(HdivSolver): def function_space(self, mesh, k): eleu = FiniteElement("BDM", mesh.ufl_cell(), k) - elep = FiniteElement("Discontinuous Lagrange L2", mesh.ufl_cell(), k-1) + elep = FiniteElement("Discontinuous Lagrange", mesh.ufl_cell(), k-1) V = FunctionSpace(mesh, eleu) Q = FunctionSpace(mesh, elep) return MixedFunctionSpace([V, Q]) diff --git a/examples/Makefile b/examples/Makefile index 9f35932..180c03a 100644 --- a/examples/Makefile +++ b/examples/Makefile @@ -1,4 +1,4 @@ -exec=mpirun -n 12 python3 +exec=mpiexec -n 3 python3 pre: mkdir -p logs/ @@ -9,11 +9,11 @@ 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 @@ -31,3 +31,6 @@ 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 bdm --solver-type allu --mh uniform --patch star --gamma 1e3 2>&1 | tee logs/mms2dhdiv.log diff --git a/examples/bfs2d/bfs2d.py b/examples/bfs2d/bfs2d.py index 26d779a..4abf566 100644 --- a/examples/bfs2d/bfs2d.py +++ b/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))) diff --git a/examples/ldc2d/ldc2d.py b/examples/ldc2d/ldc2d.py index 83a4f75..31ea29c 100644 --- a/examples/ldc2d/ldc2d.py +++ b/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) diff --git a/examples/mms.py b/examples/mms.py index 939a8bb..20089c2 100644 --- a/examples/mms.py +++ b/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 = [0, 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) diff --git a/examples/mmsldc2d/mmsldc2d.py b/examples/mmsldc2d/mmsldc2d.py index 278428d..ffc157c 100644 --- a/examples/mmsldc2d/mmsldc2d.py +++ b/examples/mmsldc2d/mmsldc2d.py @@ -26,15 +26,17 @@ def mesh(self, distribution_parameters): return base def bcs(self, Z): - bcs = [DirichletBC(Z.sub(0), self.actual_solution(Z)[0], [4]), - DirichletBC(Z.sub(0), Constant((0., 0.)), [1, 2, 3])] + bcs = [DirichletBC(Z.sub(0), self.actual_solution(Z)[0], 4), + 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 def interpolate_initial_guess(self, z): - w_expr = self.actual_solution(z)[0] - z.sub(0).interpolate(w_expr) + w_expr = self.actual_solution(z.function_space().sub(0))[0] + z.sub(0).project(w_expr) def char_length(self): return 2.0 @@ -59,7 +61,8 @@ def actual_solution(self, Z): G1 = -24 * y**5+8*y**3-4*y u = 8 * f * dg v = -8 * df * g - p = (8./self.Re) * (F * dddg + df*dg) + 64 * F2 * (g*ddg-dg**2) + # p = (8./self.Re) * (F * dddg + df*dg) + 64 * F2 * (g*ddg-dg**2) + p = (8.) * (F * dddg + df*dg) + 64 * F2 * (g*ddg-dg**2) u = replace(u, {X: 0.5 * X}) v = replace(v, {X: 0.5 * X}) p = replace(p, {X: 0.5 * X}) @@ -73,7 +76,8 @@ def actual_solution(self, Z): def rhs(self, Z): (u, p) = self.actual_solution(Z) - nu = self.char_length() * self.char_velocity() / self.Re + nu = self.char_length() * self.char_velocity() / conditional(self.Re>0, self.Re, 1) f1 = -nu * div(2*sym(grad(u))) + dot(grad(u), u) + grad(p) + # f1 = -nu * div(2*sym(grad(u))) + grad(p) f2 = -div(u) return f1, f2 From d25feebbccc57db0ab55051e6f16883e1d23b95b Mon Sep 17 00:00:00 2001 From: Fabian Laakmann Date: Mon, 13 Jul 2020 16:13:08 +0100 Subject: [PATCH 3/5] replaced EmbeddedDGTransfer with new transfer syntax --- alfi/solver.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/alfi/solver.py b/alfi/solver.py index dbd4682..1d5b6e4 100644 --- a/alfi/solver.py +++ b/alfi/solver.py @@ -720,8 +720,9 @@ def b(u, q): def get_transfers(self): V = self.Z.sub(0) - prolongation = EmbeddedDGTransfer(V.ufl_element()) - return [dmhooks.transfer_operators(V, prolong=prolongation.prolong, inject=prolongation.inject, restrict=prolongation.restrict)] + 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" From 90eb6b1d66409b4b527d835743971487f62fc431 Mon Sep 17 00:00:00 2001 From: Fabian Laakmann Date: Tue, 20 Oct 2020 21:27:25 +0100 Subject: [PATCH 4/5] update branch + make mms in 3d run through. I see reduced convergence rates for BDM and no convergence for RT --- alfi/bary.py | 34 ++++++++++++++++++++-------------- alfi/solver.py | 4 ++-- examples/Makefile | 9 +++++++-- examples/mms.py | 2 +- examples/mmsldc3d/mmsldc3d.py | 8 ++++++-- 5 files changed, 36 insertions(+), 21 deletions(-) diff --git a/alfi/bary.py b/alfi/bary.py index e9a6195..301d332 100644 --- a/alfi/bary.py +++ b/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") diff --git a/alfi/solver.py b/alfi/solver.py index 1d5b6e4..4abbd93 100644 --- a/alfi/solver.py +++ b/alfi/solver.py @@ -736,7 +736,7 @@ def distribution_parameters(self): class RTSolver(HdivSolver): def function_space(self, mesh, k): - eleu = FiniteElement("RT", mesh.ufl_cell(), 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) @@ -746,7 +746,7 @@ def function_space(self, mesh, k): class BDMSolver(HdivSolver): def function_space(self, mesh, k): - eleu = FiniteElement("BDM", mesh.ufl_cell(), 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) diff --git a/examples/Makefile b/examples/Makefile index 180c03a..706b15a 100644 --- a/examples/Makefile +++ b/examples/Makefile @@ -14,7 +14,7 @@ iters2dsv: pre $(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 @@ -33,4 +33,9 @@ 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 1e4 2>&1 | tee logs/idealalgamma1e4.log || 1 mms2dhdiv: pre - $(exec) mms.py --dim 2 --nref 2 --k 2 --discretisation bdm --solver-type allu --mh uniform --patch star --gamma 1e3 2>&1 | tee logs/mms2dhdiv.log + $(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 diff --git a/examples/mms.py b/examples/mms.py index 20089c2..f9499f4 100644 --- a/examples/mms.py +++ b/examples/mms.py @@ -22,7 +22,7 @@ else: raise NotImplementedError -res = [0, 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] = {} diff --git a/examples/mmsldc3d/mmsldc3d.py b/examples/mmsldc3d/mmsldc3d.py index 2319bd6..70021cf 100644 --- a/examples/mmsldc3d/mmsldc3d.py +++ b/examples/mmsldc3d/mmsldc3d.py @@ -22,8 +22,12 @@ def mesh(self, distribution_parameters): return base def bcs(self, Z): - bcs = [DirichletBC(Z.sub(0), self.actual_solution(Z)[0], [4, 5, 6]), - DirichletBC(Z.sub(0), Constant((0., 0., 0.)), [1, 2, 3])] + bcs = [DirichletBC(Z.sub(0), self.actual_solution(Z)[0], 4), + DirichletBC(Z.sub(0), self.actual_solution(Z)[0], 5), + DirichletBC(Z.sub(0), self.actual_solution(Z)[0], 6), + DirichletBC(Z.sub(0), Constant((0., 0., 0.)), 1), + DirichletBC(Z.sub(0), Constant((0., 0., 0.)), 2), + DirichletBC(Z.sub(0), Constant((0., 0., 0.)), 3)] return bcs def has_nullspace(self): return True From c322a5b3f9ca6dafd491d7b6fc5179f2fc3a658d Mon Sep 17 00:00:00 2001 From: Fabian Laakmann Date: Wed, 4 Nov 2020 15:17:53 +0000 Subject: [PATCH 5/5] we need CellDiameter instead of FacetArea in the penalty term to make 3d working --- alfi/solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/alfi/solver.py b/alfi/solver.py index 4abbd93..0a0e77b 100644 --- a/alfi/solver.py +++ b/alfi/solver.py @@ -672,7 +672,7 @@ 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 = FacetArea(self.z.ufl_domain()) + h = CellDiameter(self.z.ufl_domain()) n = FacetNormal(self.z.ufl_domain()) ulast = split(self.z_last)[0]