tests

Author: ksagiyam

tests/firedrake/submesh/test_submesh_solve.py

@@ -463,12 +463,12 @@ def test_submesh_solve_cell_cell_equation_bc(nref, degree, simplex):
     assert sqrt(assemble(inner(sol[1] - x * y, sol[1] - x * y) * dx_inner)) < 1.e-12
 
 
-def _test_submesh_solve_quad_triangle_convergence(nref, degree):
+def _test_submesh_solve_quad_triangle_poisson(nref, degree):
     dim = 2
     label_ext = 1
     label_interf = 2
     distribution_parameters_noop = {
-        "partition": False,
+        "partition": True,
         "overlap_type": (DistributedMeshOverlapType.NONE, 0),
     }
     mesh = Mesh(os.path.join(cwd, "..", "meshes", "mixed_cell_unit_square.msh"), distribution_parameters=distribution_parameters_noop)
@@ -523,20 +523,20 @@ def _test_submesh_solve_quad_triangle_convergence(nref, degree):
     bc_q = DirichletBC(V.sub(1), g_q, label_ext)
     solve(a == L, sol, bcs=[bc_q])
     sol_t, sol_q = split(sol)
-    L2error_t = assemble(inner(sol_t - g_t, sol_t - g_t) * dx_t)
-    L2error_q = assemble(inner(sol_q - g_q, sol_q - g_q) * dx_q)
-    H1error_t = assemble((inner(sol_t - g_t, sol_t - g_t) + inner(grad(sol_t - g_t), grad(sol_t - g_t)))* dx_t)
-    H1error_q = assemble((inner(sol_q - g_q, sol_q - g_q) + inner(grad(sol_q - g_q), grad(sol_q - g_q)))* dx_q)
-    return sqrt(L2error_t + L2error_q), sqrt(H1error_t + H1error_q)
+    L2Error_t = assemble(inner(sol_t - g_t, sol_t - g_t) * dx_t)
+    L2Error_q = assemble(inner(sol_q - g_q, sol_q - g_q) * dx_q)
+    H1Error_t = L2Error_t + assemble(inner(grad(sol_t - g_t), grad(sol_t - g_t))* dx_t)
+    H1Error_q = L2Error_q + assemble(inner(grad(sol_q - g_q), grad(sol_q - g_q))* dx_q)
+    return sqrt(L2Error_t + L2Error_q), sqrt(H1Error_t + H1Error_q)
 
 
 @pytest.mark.parallel(nprocs=8)
-def test_submesh_solve_quad_triangle_convergence():
+def test_submesh_solve_quad_triangle_poisson_convergence():
     for degree in range(1, 5):
         L2Errors = []
         H1Errors = []
         for nref in range(4):
-            L2Error, H1Error = _test_submesh_solve_quad_triangle_convergence(nref, degree)
+            L2Error, H1Error = _test_submesh_solve_quad_triangle_poisson(nref, degree)
             L2Errors.append(L2Error)
             H1Errors.append(H1Error)
         L2Errors = [np.log2(c) - np.log2(f) for c, f in zip(L2Errors[:-1], L2Errors[1:])]
@@ -545,25 +545,26 @@ def test_submesh_solve_quad_triangle_convergence():
         assert (np.array(H1Errors) > (degree) * 0.995).all()
 
 
-@pytest.mark.parallel(nprocs=8)
-@pytest.mark.parametrize('simplex', [True, False])
-@pytest.mark.parametrize('direction', [0, 1, 2])
-def test_submesh_solve_hex_quad_convergence(simplex, direction):
+def _test_submesh_solve_3d_2d_poisson(simplex, direction, nref, degree):
     dim = 3
+    distribution_parameters_noop = {
+        "partition": True,
+        "overlap_type": (DistributedMeshOverlapType.NONE, 0),
+    }
     distribution_parameters={
         "overlap_type": (DistributedMeshOverlapType.RIDGE, 1),
     }
     if simplex:
-        nref = 3
-        mesh = BoxMesh(2 ** nref, 2 ** nref, 2 ** nref, 1., 1., 1., hexahedral=False, distribution_parameters=distribution_parameters)
+        nref_simplex = 3
+        mesh = BoxMesh(2 ** nref_simplex, 2 ** nref_simplex, 2 ** nref_simplex, 1., 1., 1., hexahedral=False, distribution_parameters=distribution_parameters_noop)
         xyz = SpatialCoordinate(mesh)
         DG0 = FunctionSpace(mesh, "DG", 0)
         c1 = Function(DG0).interpolate(conditional(xyz[direction] < 0.318310, 1, 0))
         c2 = Function(DG0).interpolate(conditional(xyz[direction] > 0.318310, 1, 0))
         mesh = RelabeledMesh(mesh, [c1, c2], [1, 2])
         family = "P"
     else:
-        mesh = Mesh(join(cwd, "..", "meshes", "cube_hex.msh"), distribution_parameters=distribution_parameters)
+        mesh = Mesh(join(cwd, "..", "meshes", "cube_hex.msh"), distribution_parameters=distribution_parameters_noop)
         xyz = SpatialCoordinate(mesh)
         DG0 = FunctionSpace(mesh, "DQ", 0)
         c1 = Function(DG0).interpolate(conditional(xyz[direction] < 0.318310, 1, 0))
@@ -577,6 +578,13 @@ def test_submesh_solve_hex_quad_convergence(simplex, direction):
         f6 = Function(HDivTrace0).interpolate(conditional(xyz[2] < .999, 1, 0))
         mesh = RelabeledMesh(mesh, [c1, c2, f1, f2, f3, f4, f5, f6], [1, 2, 1, 2, 3, 4, 5, 6])
         family = "Q"
+    plex = mesh.topology_dm
+    for _ in range(nref):
+        plex = plex.refine()
+    plex.removeLabel("pyop2_core")
+    plex.removeLabel("pyop2_owned")
+    plex.removeLabel("pyop2_ghost")
+    mesh = Mesh(plex, distribution_parameters=distribution_parameters)
     mesh1 = Submesh(mesh, dim, 1)
     x1, y1, z1 = SpatialCoordinate(mesh1)
     mesh2 = Submesh(mesh, dim, 2)
@@ -586,17 +594,23 @@ def test_submesh_solve_hex_quad_convergence(simplex, direction):
     dx1 = Measure("dx", mesh1)
     dx2 = Measure("dx", mesh2)
     dx12 = Measure("dx", mesh12)
-    dx12_ds1 = Measure("dx", mesh12, intersect_measures=(Measure("ds", mesh1),))
-    ds1_dx12 = Measure("ds", mesh1, intersect_measures=(Measure("dx", mesh12),))
+    ds1_ds2 = Measure("ds", mesh1, intersect_measures=(Measure("ds", mesh2),))
     ds2_dx12 = Measure("ds", mesh2, intersect_measures=(Measure("dx", mesh12),))
+    dx12_ds1_ds2 = Measure(
+        "dx", mesh12,
+        intersect_measures=(
+            Measure("ds", mesh1),
+            Measure("ds", mesh2),
+        )
+    )
     # Check sanity.
     vol1 = assemble(Constant(1) * dx1)
     vol2 = assemble(Constant(1) * dx2)
-    assert abs(vol1 + vol2 - 1.) < 1.e-14
+    assert abs(vol1 + vol2 - 1.) < 1.e-13
     # Solve Poisson problem.
-    V1 = FunctionSpace(mesh1, family, 3)
-    V12 = FunctionSpace(mesh12, family, 5)
-    V2 = FunctionSpace(mesh2, family, 4)
+    V1 = FunctionSpace(mesh1, family, degree)
+    V12 = FunctionSpace(mesh12, family, degree)
+    V2 = FunctionSpace(mesh2, family, degree)
     V = V1 * V12 * V2
     u = TrialFunction(V)
     v = TestFunction(V)
@@ -607,22 +621,70 @@ def test_submesh_solve_hex_quad_convergence(simplex, direction):
     f1 = 12 * pi**2 * g1
     f2 = 12 * pi**2 * g2
     n1 = FacetNormal(mesh1)
+    n2 = FacetNormal(mesh2)
+    h = 0.1 / 2**nref  # roughly
     a = (
-        inner(grad(u1), grad(v1)) * dx1 - inner(u12, v1) * ds1_dx12(label_interf) +
-        inner(u12, v12) * dx12 +
-        inner(grad(u2), grad(v2)) * dx2 + inner(u12, v2) * ds2_dx12(label_interf)
+        inner(grad(u1), grad(v1)) * dx1 +
+        inner(grad(u2), grad(v2)) * dx2
+        - inner(
+            u12,
+            (v1 - v2)
+        ) * dx12_ds1_ds2
+        - inner(
+            (u1 * n1 + u2 * n2),
+            (grad(v1) + grad(v2)) / 2
+        ) * dx12_ds1_ds2
+        + 100 / h * inner(u1 - u2, v1 - v2) * ds1_ds2(label_interf)  # Can also use dx12_ds1_ds2.
+        + inner(
+            (dot(grad(u1), n1) - dot(grad(u2), n1)) / 2 - u12,
+            v12
+        ) * dx12_ds1_ds2
     )
     L = (
         inner(f1, v1) * dx1 +
-        inner(dot(grad(g1), n1), v12) * dx12_ds1 +
         inner(f2, v2) * dx2
     )
     sol = Function(V)
-    bc1 = DirichletBC(V.sub(0), g1, (2 * direction + 1,))
-    bc2 = DirichletBC(V.sub(2), g2, (2 * direction + 2,))
-    solve(a == L, sol, bcs=[bc1, bc2])
+    bc1 = DirichletBC(V.sub(0), g1, [i for i in range(1, 7) if i != 2 * direction + 2])
+    bc2 = DirichletBC(V.sub(2), g2, [i for i in range(1, 7) if i != 2 * direction + 1])
+    solver_parameters = {
+        "mat_type": "matfree",
+        "ksp_type": "gmres",
+        "ksp_rtol": 1.e-14,
+        "pc_type": "jacobi",
+    }
+    solve(a == L, sol, bcs=[bc1, bc2], solver_parameters=solver_parameters)
     sol1, sol12, sol2 = split(sol)
-    error1 = assemble(inner(sol1 - g1, sol1 - g1) * dx1)
-    error2 = assemble(inner(sol2 - g2, sol2 - g2) * dx2)
-    assert sqrt(error1) < 4.e-4
-    assert sqrt(error2) < 4.e-5
+    L2Error1 = assemble(inner(sol1 - g1, sol1 - g1) * dx1)
+    L2Error2 = assemble(inner(sol2 - g2, sol2 - g2) * dx2)
+    H1Error1 = L2Error1 + assemble(inner(grad(sol1 - g1), grad(sol1 - g1))* dx1)
+    H1Error2 = L2Error2 + assemble(inner(grad(sol2 - g2), grad(sol2 - g2))* dx2)
+    return sqrt(L2Error1 + L2Error2), sqrt(H1Error2 + H1Error2)
+
+
+@pytest.mark.parallel(nprocs=6)
+@pytest.mark.parametrize('simplex', [True, False])
+@pytest.mark.parametrize('direction', [0, 1, 2])
+def test_submesh_solve_3d_2d_poisson_sanity(simplex, direction):
+    nref = 0
+    degree = 4
+    L2Error, H1Error = _test_submesh_solve_3d_2d_poisson(simplex, direction, nref, degree)
+    assert L2Error < 2.e-4
+    assert H1Error < 3.e-2
+
+
+@pytest.mark.parallel(nprocs=8)
+@pytest.mark.parametrize('simplex', [False])
+@pytest.mark.parametrize('direction', [0])
+@pytest.mark.parametrize('degree', [3])
+def test_submesh_solve_3d_2d_poisson_convergence(simplex, direction, degree):
+    L2Errors = []
+    H1Errors = []
+    for nref in range(2):
+        L2Error, H1Error = _test_submesh_solve_3d_2d_poisson(simplex, direction, nref, degree)
+        L2Errors.append(L2Error)
+        H1Errors.append(H1Error)
+    L2Errors = [np.log2(c) - np.log2(f) for c, f in zip(L2Errors[:-1], L2Errors[1:])]
+    H1Errors = [np.log2(c) - np.log2(f) for c, f in zip(H1Errors[:-1], H1Errors[1:])]
+    assert (np.array(L2Errors) > (degree + 1) * 0.98).all()
+    assert (np.array(H1Errors) > (degree) * 0.98).all()