Test MTW tetrahedron (#4884)

Author: pbrubeck

tests/firedrake/regression/test_mtw.py

@@ -1,20 +1,28 @@
 from firedrake import *
+import pytest
 import numpy as np
 
 
-convergence_orders = lambda x: np.log2(np.array(x)[:-1] / np.array(x)[1:])
-
-
-def test_mtw():
+@pytest.fixture(params=(2, 3))
+def mh(request):
+    dim = request.param
     N_base = 2
-    msh = UnitSquareMesh(N_base, N_base)
-    mh = MeshHierarchy(msh, 5)
+    if dim == 2:
+        refine = 3
+        msh = UnitSquareMesh(N_base, N_base)
+    elif dim == 3:
+        refine = 2
+        msh = UnitCubeMesh(N_base, N_base, N_base)
+    else:
+        raise ValueError("Unexpected dimension")
+    mh = MeshHierarchy(msh, refine)
 
     V = FunctionSpace(msh, msh.coordinates.ufl_element())
     eps = Constant(1 / 2**(N_base-1))
-    x, y = SpatialCoordinate(msh)
+    x, y, *z = SpatialCoordinate(msh)
+
     new = Function(V).interpolate(as_vector([x + eps*sin(2*pi*x)*sin(2*pi*y),
-                                             y - eps*sin(2*pi*x)*sin(2*pi*y)]))
+                                             y - eps*sin(2*pi*x)*sin(2*pi*y), *z]))
 
     # And propagate to refined meshes
     coords = [new]
@@ -26,63 +34,86 @@ def test_mtw():
 
     for msh, coord in zip(mh, coords):
         msh.coordinates.assign(coord)
-
-    params = {"snes_type": "newtonls",
-              "snes_linesearch_type": "basic",
-              "snes_monitor": None,
-              "mat_type": "aij",
-              "snes_max_it": 10,
-              "snes_lag_jacobian": -2,
-              "snes_lag_preconditioner": -2,
-              "ksp_type": "preonly",
-              "pc_type": "lu",
-              "pc_factor_shift_type": "inblocks",
-              "snes_rtol": 1e-16,
-              "snes_atol": 1e-25}
-
+    return mh
+
+
+def mesh_sizes(mh):
+    mesh_size = []
+    for msh in mh:
+        DG0 = FunctionSpace(msh, "DG", 0)
+        h = Function(DG0).interpolate(CellDiameter(msh))
+        with h.dat.vec as hvec:
+            _, maxh = hvec.max()
+        mesh_size.append(maxh)
+    return mesh_size
+
+
+def convergence_orders(error, h):
+    return np.diff(np.log2(error)) / np.diff(np.log2(h))
+
+
+def test_mtw_darcy_convergence(mh):
+    sp = {
+        "ksp_monitor": None,
+        "mat_type": "matfree",
+        "pmat_type": "nest",
+        "ksp_type": "minres",
+        "ksp_norm_type": "preconditioned",
+        "pc_type": "fieldsplit",
+        "pc_fieldsplit_type": "additive",
+        "fieldsplit_ksp_type": "preonly",
+        "fieldsplit_0_pc_type": "lu",
+        "fieldsplit_0_pc_factor_mat_solver_type": "mumps",
+        "fieldsplit_1_pc_type": "jacobi",
+    }
+    gamma = Constant(1E4)
     l2_u = []
     l2_p = []
     for msh in mh[1:]:
-        x, y = SpatialCoordinate(msh)
+        x, y, *z = SpatialCoordinate(msh)
         pex = sin(pi * x) * sin(2 * pi * y)
+        if z:
+            pex *= sin(pi*z[0])
+
         uex = -grad(pex)
         f = div(uex)
 
-        V = FunctionSpace(msh, "MTW", 3)
-        W = FunctionSpace(msh, "DG", 0)
-        Z = V * W
+        V = FunctionSpace(msh, "MTW", 1)
+        Q = FunctionSpace(msh, "DG", 0)
+        Z = V * Q
 
-        up = Function(Z)
-        u, p = split(up)
-        v, w = TestFunctions(Z)
+        u, p = TrialFunctions(Z)
+        v, q = TestFunctions(Z)
+
+        a = -inner(u, v) * dx + inner(p, div(v)) * dx + inner(div(u), q) * dx
+        L = inner(f, q) * dx
+
+        Jp = inner(u, v)*dx + inner(div(u), gamma*div(v))*dx + inner(p/gamma, q)*dx
 
-        F = (inner(u, v) * dx - inner(p, div(v)) * dx
-             + inner(div(u), w) * dx - inner(f, w) * dx)
+        up = Function(Z)
 
-        solve(F == 0, up, solver_parameters=params)
+        solve(a == L, up, Jp=Jp, solver_parameters=sp)
 
         u, p = up.subfunctions
         l2_u.append(errornorm(uex, u))
         l2_p.append(errornorm(pex, p))
 
-    assert min(convergence_orders(l2_u)) > 1.8
-    assert min(convergence_orders(l2_p)) > 0.8
+    h = mesh_sizes(mh[1:])
+    assert min(convergence_orders(l2_u, h)) > 1.75
+    assert min(convergence_orders(l2_p, h)) > 0.8
 
 
 def test_mtw_interior_facet():
     mesh = UnitSquareMesh(4, 4)
-    V = FunctionSpace(mesh, mesh.coordinates.ufl_element())
     eps = Constant(0.5 / 2**3)
 
     x, y = SpatialCoordinate(mesh)
-    new = Function(V).interpolate(as_vector([x + eps*sin(2*pi*x)*sin(2*pi*y),
-                                             y - eps*sin(2*pi*x)*sin(2*pi*y)]))
-    mesh = Mesh(new)
+    mesh.coordinates.interpolate(as_vector([x + eps*sin(2*pi*x)*sin(2*pi*y),
+                                            y - eps*sin(2*pi*x)*sin(2*pi*y)]))
 
-    V = FunctionSpace(mesh, 'Mardal-Tai-Winther', 3)
+    V = FunctionSpace(mesh, 'Mardal-Tai-Winther', 1)
 
-    x, y = SpatialCoordinate(mesh)
-    uh = project(as_vector((x+y, 2*x-y)), V)
+    uh = Function(V).interpolate(as_vector((x+y, 2*x-y)))
 
     volume = assemble(div(uh)*dx)