From cb3ce18f90bcb174c471342be0bf26648a639086 Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Thu, 20 Jun 2019 15:46:20 +0100 Subject: [PATCH 01/11] testing for l2pullbacks --- tests/regression/test_assemble.py | 25 +++++++++++- tests/regression/test_assemble_inverse.py | 14 +++++++ tests/regression/test_function_spaces.py | 13 ++++-- tests/regression/test_helmholtz_mixed.py | 6 ++- tests/regression/test_interpolate.py | 40 +++++++++++++++++++ tests/regression/test_stokes_hdiv_parallel.py | 6 ++- 6 files changed, 96 insertions(+), 8 deletions(-) diff --git a/tests/regression/test_assemble.py b/tests/regression/test_assemble.py index 61f1dab79e..ed72352839 100644 --- a/tests/regression/test_assemble.py +++ b/tests/regression/test_assemble.py @@ -12,7 +12,9 @@ def mesh(): 'cg1cg1', 'cg1cg1[0]', 'cg1cg1[1]', 'cg1vcg1[0]', 'cg1vcg1[1]', 'cg1dg0', 'cg1dg0[0]', 'cg1dg0[1]', - 'cg2dg1', 'cg2dg1[0]', 'cg2dg1[1]']) + 'cg2dg1', 'cg2dg1[0]', 'cg2dg1[1]', + 'cg1dg0l2', 'cg1dg0l2[0]', 'cg1dg0l2[1]', + 'cg2dg1l2', 'cg2dg1l2[0]', 'cg2dg1l2[1]']) def fs(request, mesh): cg1 = FunctionSpace(mesh, "CG", 1) cg2 = FunctionSpace(mesh, "CG", 2) @@ -20,6 +22,8 @@ def fs(request, mesh): tcg1 = TensorFunctionSpace(mesh, "CG", 1) dg0 = FunctionSpace(mesh, "DG", 0) dg1 = FunctionSpace(mesh, "DG", 1) + dg0l2 = FunctionSpace(mesh, "DG L2", 0) + dg1l2 = FunctionSpace(mesh, "DG L2", 1) return {'cg1': cg1, 'vcg1': vcg1, 'tcg1': tcg1, @@ -34,7 +38,13 @@ def fs(request, mesh): 'cg1dg0[1]': (cg1*dg0)[1], 'cg2dg1': cg2*dg1, 'cg2dg1[0]': (cg2*dg1)[0], - 'cg2dg1[1]': (cg2*dg1)[1]}[request.param] + 'cg2dg1[1]': (cg2*dg1)[1], + 'cg1dg0l2': cg1*dg0l2, + 'cg1dg0l2[0]': (cg1*dg0l2)[0], + 'cg1dg0l2[1]': (cg1*dg0l2)[1], + 'cg2dg1l2': cg2*dg1l2, + 'cg2dg1l2[0]': (cg2*dg1l2)[0], + 'cg2dg1l2[1]': (cg2*dg1l2)[1]}[request.param] @pytest.fixture @@ -122,3 +132,14 @@ def test_assemble_mat_with_tensor(mesh): M = assemble(Constant(2)*a, M) # Make sure we get the result of the last assembly assert np.allclose(M.M.values, 2*assemble(a).M.values, rtol=1e-14) + +def test_assemble_mat_with_tensor_l2(mesh): + V = FunctionSpace(mesh, "DG L2", 0) + u = TestFunction(V) + v = TrialFunction(V) + a = u*v*dx + M = assemble(a) + # Assemble a different form into M + M = assemble(Constant(2)*a, M) + # Make sure we get the result of the last assembly + assert np.allclose(M.M.values, 2*assemble(a).M.values, rtol=1e-14) diff --git a/tests/regression/test_assemble_inverse.py b/tests/regression/test_assemble_inverse.py index 66fa8664d6..d13cd7e5b0 100644 --- a/tests/regression/test_assemble_inverse.py +++ b/tests/regression/test_assemble_inverse.py @@ -16,3 +16,17 @@ def test_assemble_inverse(degree): eye = np.dot(m_forward.M.values, m_inverse.M.values) assert ((eye - np.eye(fs.node_count)).round(12) == 0).all() + +@pytest.mark.parametrize("degree", range(4)) +def test_assemble_inverse_l2(degree): + m = UnitSquareMesh(2, 1) + fs = FunctionSpace(m, "DG L2", degree) + u = TrialFunction(fs) + v = TestFunction(fs) + + m_forward = assemble(u*v*dx) + m_inverse = assemble(u*v*dx, inverse=True) + + eye = np.dot(m_forward.M.values, m_inverse.M.values) + + assert ((eye - np.eye(fs.node_count)).round(12) == 0).all() diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 3d90a42306..705c991102 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -26,13 +26,18 @@ def cg2(mesh): def dg0(mesh): return VectorFunctionSpace(mesh, "DG", 0) +@pytest.fixture(scope="module") +def dg0l2(mesh): + return VectorFunctionSpace(mesh, "DG L2", 0) -@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0']) +@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0', 'cg1dg0l2', 'cg2dg0l2']) def fs(request, cg1, cg2, dg0): return {'cg1cg1': cg1*cg1, 'cg1cg2': cg1*cg2, 'cg1dg0': cg1*dg0, - 'cg2dg0': cg2*dg0}[request.param] + 'cg2dg0': cg2*dg0, + 'cg1dg0l2': cg1*dg0l2, + 'cg2dg0l2': cg2*dg0l2}[request.param] def test_function_space_cached(mesh): @@ -105,7 +110,9 @@ def test_function_space_collapse(cg1): TensorProductElement(FiniteElement("CG", triangle, 1), FiniteElement("CG", interval, 3)), MixedElement(FiniteElement("CG", triangle, 1), - FiniteElement("DG", triangle, 2))], + FiniteElement("DG", triangle, 2)), + MixedElement(FiniteElement("CG", triangle, 1), + FiniteElement("DG L2", triangle, 2))], ids=["FE", "Enriched", "TPE", "Mixed"]) def test_validation(modifier, element): with pytest.raises(ValueError): diff --git a/tests/regression/test_helmholtz_mixed.py b/tests/regression/test_helmholtz_mixed.py index d9472bfe7b..e145ef32f6 100644 --- a/tests/regression/test_helmholtz_mixed.py +++ b/tests/regression/test_helmholtz_mixed.py @@ -48,7 +48,11 @@ def helmholtz_mixed(r, V1, V2, action=False): [(('RT', 1), ('DG', 0), 1.9, False), (('BDM', 1), ('DG', 0), 1.89, False), (('BDM', 1), ('DG', 0), 1.89, True), - (('BDFM', 2), ('DG', 1), 1.9, False)]) + (('BDFM', 2), ('DG', 1), 1.9, False), + (('RT', 1), ('DG L2', 0), 1.9, False), + (('BDM', 1), ('DG L2', 0), 1.89, False), + (('BDM', 1), ('DG L2', 0), 1.89, True), + (('BDFM', 2), ('DG L2', 1), 1.9, False)]) def test_firedrake_helmholtz(V1, V2, threshold, action): import numpy as np diff = np.array([helmholtz_mixed(i, V1, V2) for i in range(3, 6)]) diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index 90571aed82..13242e7877 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -122,6 +122,21 @@ def test_cell_orientation(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 +def test_cell_orientation_l2(): + m = UnitCubedSphereMesh(2) + x = SpatialCoordinate(m) + m.init_cell_orientations(x) + x = m.coordinates + U = FunctionSpace(m, 'RTCF', 1) + V = VectorFunctionSpace(m, 'DQ L2', 1) + + f = project(as_tensor([x[1], -x[0], 0.0]), U) + g = interpolate(f, V) + + # g shall be close to: + h = project(f, V) + + assert abs(g.dat.data - h.dat.data).max() < 1e-2 def test_cellvolume(): m = UnitSquareMesh(2, 2) @@ -131,6 +146,13 @@ def test_cellvolume(): assert np.allclose(f.dat.data_ro, 0.125) +def test_cellvolume_l2(): + m = UnitSquareMesh(2, 2) + V = FunctionSpace(m, 'DG L2', 0) + + f = interpolate(CellVolume(m), V) + + assert np.allclose(f.dat.data_ro, 0.125) def test_cellvolume_higher_order_coords(): m = UnitTriangleMesh() @@ -151,6 +173,24 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) +def test_cellvolume_higher_order_coords_l2(): + m = UnitTriangleMesh() + V = VectorFunctionSpace(m, 'P', 3) + f = Function(V) + f.interpolate(m.coordinates) + + # Warp mesh so that the bottom triangle line is: + # x(x - 1)(x + a) with a = 19/12.0 + def warp(x): + return x * (x - 1)*(x + 19/12.0) + + f.dat.data[1, 1] = warp(1.0/3.0) + f.dat.data[2, 1] = warp(2.0/3.0) + + mesh = Mesh(f) + g = interpolate(CellVolume(mesh), FunctionSpace(mesh, 'DG L2', 0)) + + assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) def test_mixed(): m = UnitTriangleMesh() diff --git a/tests/regression/test_stokes_hdiv_parallel.py b/tests/regression/test_stokes_hdiv_parallel.py index 223f260118..9a2f0ff6bd 100644 --- a/tests/regression/test_stokes_hdiv_parallel.py +++ b/tests/regression/test_stokes_hdiv_parallel.py @@ -9,8 +9,10 @@ def mat_type(request): @pytest.fixture(params=[(("RT", 3), ("DG", 2)), - (("BDM", 2), ("DG", 1))], - ids=["RT3-DG2", "BDM2-DG1"]) + (("BDM", 2), ("DG", 1)), + (("RT", 3), ("DG L2", 2)), + (("BDM", 2), ("DG L2", 1))], + ids=["RT3-DG2", "BDM2-DG1", "RT3-DG2L2", "BDM2-DG1L2"]) def element_pair(request): return request.param From 781e230ef5bf2fd4ec5635cc175a2a7b9238b2ab Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Thu, 20 Jun 2019 15:49:41 +0100 Subject: [PATCH 02/11] DROP BEFORE MERGE --- Jenkinsfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Jenkinsfile b/Jenkinsfile index 7684a36f23..9f75ae1a09 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -26,7 +26,7 @@ pipeline { sh 'mkdir tmp' dir('tmp') { timestamps { - sh '../scripts/firedrake-install --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)' + sh '../scripts/firedrake-install --package-branch ufl l2pullbacks --package-branch tsfc l2pullbacks --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)' } } } From dfdc1010cdf28593b1ffe16897ad8b05f6491ffb Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Thu, 20 Jun 2019 16:39:39 +0100 Subject: [PATCH 03/11] flake8 --- tests/regression/test_assemble.py | 1 + tests/regression/test_assemble_inverse.py | 1 + tests/regression/test_function_spaces.py | 4 +++- tests/regression/test_interpolate.py | 6 ++++++ 4 files changed, 11 insertions(+), 1 deletion(-) diff --git a/tests/regression/test_assemble.py b/tests/regression/test_assemble.py index ed72352839..728e3d82c0 100644 --- a/tests/regression/test_assemble.py +++ b/tests/regression/test_assemble.py @@ -133,6 +133,7 @@ def test_assemble_mat_with_tensor(mesh): # Make sure we get the result of the last assembly assert np.allclose(M.M.values, 2*assemble(a).M.values, rtol=1e-14) + def test_assemble_mat_with_tensor_l2(mesh): V = FunctionSpace(mesh, "DG L2", 0) u = TestFunction(V) diff --git a/tests/regression/test_assemble_inverse.py b/tests/regression/test_assemble_inverse.py index d13cd7e5b0..9f410217cc 100644 --- a/tests/regression/test_assemble_inverse.py +++ b/tests/regression/test_assemble_inverse.py @@ -17,6 +17,7 @@ def test_assemble_inverse(degree): assert ((eye - np.eye(fs.node_count)).round(12) == 0).all() + @pytest.mark.parametrize("degree", range(4)) def test_assemble_inverse_l2(degree): m = UnitSquareMesh(2, 1) diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 705c991102..433ca2071c 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -26,11 +26,13 @@ def cg2(mesh): def dg0(mesh): return VectorFunctionSpace(mesh, "DG", 0) + @pytest.fixture(scope="module") def dg0l2(mesh): return VectorFunctionSpace(mesh, "DG L2", 0) -@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0', 'cg1dg0l2', 'cg2dg0l2']) + +@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0', 'cg1dg0l2', 'cg2dg0l2']) def fs(request, cg1, cg2, dg0): return {'cg1cg1': cg1*cg1, 'cg1cg2': cg1*cg2, diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index 13242e7877..ef95f7159b 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -122,6 +122,7 @@ def test_cell_orientation(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 + def test_cell_orientation_l2(): m = UnitCubedSphereMesh(2) x = SpatialCoordinate(m) @@ -138,6 +139,7 @@ def test_cell_orientation_l2(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 + def test_cellvolume(): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'DG', 0) @@ -146,6 +148,7 @@ def test_cellvolume(): assert np.allclose(f.dat.data_ro, 0.125) + def test_cellvolume_l2(): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'DG L2', 0) @@ -154,6 +157,7 @@ def test_cellvolume_l2(): assert np.allclose(f.dat.data_ro, 0.125) + def test_cellvolume_higher_order_coords(): m = UnitTriangleMesh() V = VectorFunctionSpace(m, 'P', 3) @@ -173,6 +177,7 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) + def test_cellvolume_higher_order_coords_l2(): m = UnitTriangleMesh() V = VectorFunctionSpace(m, 'P', 3) @@ -192,6 +197,7 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) + def test_mixed(): m = UnitTriangleMesh() x = m.coordinates From 2cd121776c1cc2ca4712cbb8d221575edba2e963 Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Fri, 21 Jun 2019 13:39:04 +0100 Subject: [PATCH 04/11] More testing changes --- tests/regression/test_interpolate.py | 46 ---------------------------- 1 file changed, 46 deletions(-) diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index ef95f7159b..90571aed82 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -123,23 +123,6 @@ def test_cell_orientation(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 -def test_cell_orientation_l2(): - m = UnitCubedSphereMesh(2) - x = SpatialCoordinate(m) - m.init_cell_orientations(x) - x = m.coordinates - U = FunctionSpace(m, 'RTCF', 1) - V = VectorFunctionSpace(m, 'DQ L2', 1) - - f = project(as_tensor([x[1], -x[0], 0.0]), U) - g = interpolate(f, V) - - # g shall be close to: - h = project(f, V) - - assert abs(g.dat.data - h.dat.data).max() < 1e-2 - - def test_cellvolume(): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'DG', 0) @@ -149,15 +132,6 @@ def test_cellvolume(): assert np.allclose(f.dat.data_ro, 0.125) -def test_cellvolume_l2(): - m = UnitSquareMesh(2, 2) - V = FunctionSpace(m, 'DG L2', 0) - - f = interpolate(CellVolume(m), V) - - assert np.allclose(f.dat.data_ro, 0.125) - - def test_cellvolume_higher_order_coords(): m = UnitTriangleMesh() V = VectorFunctionSpace(m, 'P', 3) @@ -178,26 +152,6 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) -def test_cellvolume_higher_order_coords_l2(): - m = UnitTriangleMesh() - V = VectorFunctionSpace(m, 'P', 3) - f = Function(V) - f.interpolate(m.coordinates) - - # Warp mesh so that the bottom triangle line is: - # x(x - 1)(x + a) with a = 19/12.0 - def warp(x): - return x * (x - 1)*(x + 19/12.0) - - f.dat.data[1, 1] = warp(1.0/3.0) - f.dat.data[2, 1] = warp(2.0/3.0) - - mesh = Mesh(f) - g = interpolate(CellVolume(mesh), FunctionSpace(mesh, 'DG L2', 0)) - - assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) - - def test_mixed(): m = UnitTriangleMesh() x = m.coordinates From 4948a6a3d6bdd3a4763791c836ecb722d83c2c1f Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Mon, 24 Jun 2019 15:33:15 +0100 Subject: [PATCH 05/11] blank l2 pullback testing file --- tests/regression/test_l2pullbacks.py | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 tests/regression/test_l2pullbacks.py diff --git a/tests/regression/test_l2pullbacks.py b/tests/regression/test_l2pullbacks.py new file mode 100644 index 0000000000..6e8a6e7cdf --- /dev/null +++ b/tests/regression/test_l2pullbacks.py @@ -0,0 +1,7 @@ + + +# assembly/projections/function spaces/etc. + +# laplace/helmholtz mixed on plane and sphere + +# facet integrals- interior and exterior From 7df156490e1a61e9691fb2654a6e99e6dde4f94b Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Mon, 24 Jun 2019 15:33:17 +0100 Subject: [PATCH 06/11] Revert "More testing changes" This reverts commit 2cd121776c1cc2ca4712cbb8d221575edba2e963. --- tests/regression/test_interpolate.py | 46 ++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index 90571aed82..ef95f7159b 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -123,6 +123,23 @@ def test_cell_orientation(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 +def test_cell_orientation_l2(): + m = UnitCubedSphereMesh(2) + x = SpatialCoordinate(m) + m.init_cell_orientations(x) + x = m.coordinates + U = FunctionSpace(m, 'RTCF', 1) + V = VectorFunctionSpace(m, 'DQ L2', 1) + + f = project(as_tensor([x[1], -x[0], 0.0]), U) + g = interpolate(f, V) + + # g shall be close to: + h = project(f, V) + + assert abs(g.dat.data - h.dat.data).max() < 1e-2 + + def test_cellvolume(): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'DG', 0) @@ -132,6 +149,15 @@ def test_cellvolume(): assert np.allclose(f.dat.data_ro, 0.125) +def test_cellvolume_l2(): + m = UnitSquareMesh(2, 2) + V = FunctionSpace(m, 'DG L2', 0) + + f = interpolate(CellVolume(m), V) + + assert np.allclose(f.dat.data_ro, 0.125) + + def test_cellvolume_higher_order_coords(): m = UnitTriangleMesh() V = VectorFunctionSpace(m, 'P', 3) @@ -152,6 +178,26 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) +def test_cellvolume_higher_order_coords_l2(): + m = UnitTriangleMesh() + V = VectorFunctionSpace(m, 'P', 3) + f = Function(V) + f.interpolate(m.coordinates) + + # Warp mesh so that the bottom triangle line is: + # x(x - 1)(x + a) with a = 19/12.0 + def warp(x): + return x * (x - 1)*(x + 19/12.0) + + f.dat.data[1, 1] = warp(1.0/3.0) + f.dat.data[2, 1] = warp(2.0/3.0) + + mesh = Mesh(f) + g = interpolate(CellVolume(mesh), FunctionSpace(mesh, 'DG L2', 0)) + + assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) + + def test_mixed(): m = UnitTriangleMesh() x = m.coordinates From 50dc415176e079e7d3ef07751acabcb85115f46f Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Mon, 24 Jun 2019 15:34:23 +0100 Subject: [PATCH 07/11] revert changed testing files --- tests/regression/test_assemble.py | 1 - tests/regression/test_assemble_inverse.py | 1 - tests/regression/test_function_spaces.py | 4 +--- tests/regression/test_interpolate.py | 6 ------ 4 files changed, 1 insertion(+), 11 deletions(-) diff --git a/tests/regression/test_assemble.py b/tests/regression/test_assemble.py index 728e3d82c0..ed72352839 100644 --- a/tests/regression/test_assemble.py +++ b/tests/regression/test_assemble.py @@ -133,7 +133,6 @@ def test_assemble_mat_with_tensor(mesh): # Make sure we get the result of the last assembly assert np.allclose(M.M.values, 2*assemble(a).M.values, rtol=1e-14) - def test_assemble_mat_with_tensor_l2(mesh): V = FunctionSpace(mesh, "DG L2", 0) u = TestFunction(V) diff --git a/tests/regression/test_assemble_inverse.py b/tests/regression/test_assemble_inverse.py index 9f410217cc..d13cd7e5b0 100644 --- a/tests/regression/test_assemble_inverse.py +++ b/tests/regression/test_assemble_inverse.py @@ -17,7 +17,6 @@ def test_assemble_inverse(degree): assert ((eye - np.eye(fs.node_count)).round(12) == 0).all() - @pytest.mark.parametrize("degree", range(4)) def test_assemble_inverse_l2(degree): m = UnitSquareMesh(2, 1) diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 433ca2071c..705c991102 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -26,13 +26,11 @@ def cg2(mesh): def dg0(mesh): return VectorFunctionSpace(mesh, "DG", 0) - @pytest.fixture(scope="module") def dg0l2(mesh): return VectorFunctionSpace(mesh, "DG L2", 0) - -@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0', 'cg1dg0l2', 'cg2dg0l2']) +@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0', 'cg1dg0l2', 'cg2dg0l2']) def fs(request, cg1, cg2, dg0): return {'cg1cg1': cg1*cg1, 'cg1cg2': cg1*cg2, diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index ef95f7159b..13242e7877 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -122,7 +122,6 @@ def test_cell_orientation(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 - def test_cell_orientation_l2(): m = UnitCubedSphereMesh(2) x = SpatialCoordinate(m) @@ -139,7 +138,6 @@ def test_cell_orientation_l2(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 - def test_cellvolume(): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'DG', 0) @@ -148,7 +146,6 @@ def test_cellvolume(): assert np.allclose(f.dat.data_ro, 0.125) - def test_cellvolume_l2(): m = UnitSquareMesh(2, 2) V = FunctionSpace(m, 'DG L2', 0) @@ -157,7 +154,6 @@ def test_cellvolume_l2(): assert np.allclose(f.dat.data_ro, 0.125) - def test_cellvolume_higher_order_coords(): m = UnitTriangleMesh() V = VectorFunctionSpace(m, 'P', 3) @@ -177,7 +173,6 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) - def test_cellvolume_higher_order_coords_l2(): m = UnitTriangleMesh() V = VectorFunctionSpace(m, 'P', 3) @@ -197,7 +192,6 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) - def test_mixed(): m = UnitTriangleMesh() x = m.coordinates From 8dbf47748c8c377ff88ff1c86a654cd1bbecb30d Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Mon, 24 Jun 2019 15:34:27 +0100 Subject: [PATCH 08/11] Revert "testing for l2pullbacks" This reverts commit cb3ce18f90bcb174c471342be0bf26648a639086. --- tests/regression/test_assemble.py | 25 +----------- tests/regression/test_assemble_inverse.py | 14 ------- tests/regression/test_function_spaces.py | 13 ++---- tests/regression/test_helmholtz_mixed.py | 6 +-- tests/regression/test_interpolate.py | 40 ------------------- tests/regression/test_stokes_hdiv_parallel.py | 6 +-- 6 files changed, 8 insertions(+), 96 deletions(-) diff --git a/tests/regression/test_assemble.py b/tests/regression/test_assemble.py index ed72352839..61f1dab79e 100644 --- a/tests/regression/test_assemble.py +++ b/tests/regression/test_assemble.py @@ -12,9 +12,7 @@ def mesh(): 'cg1cg1', 'cg1cg1[0]', 'cg1cg1[1]', 'cg1vcg1[0]', 'cg1vcg1[1]', 'cg1dg0', 'cg1dg0[0]', 'cg1dg0[1]', - 'cg2dg1', 'cg2dg1[0]', 'cg2dg1[1]', - 'cg1dg0l2', 'cg1dg0l2[0]', 'cg1dg0l2[1]', - 'cg2dg1l2', 'cg2dg1l2[0]', 'cg2dg1l2[1]']) + 'cg2dg1', 'cg2dg1[0]', 'cg2dg1[1]']) def fs(request, mesh): cg1 = FunctionSpace(mesh, "CG", 1) cg2 = FunctionSpace(mesh, "CG", 2) @@ -22,8 +20,6 @@ def fs(request, mesh): tcg1 = TensorFunctionSpace(mesh, "CG", 1) dg0 = FunctionSpace(mesh, "DG", 0) dg1 = FunctionSpace(mesh, "DG", 1) - dg0l2 = FunctionSpace(mesh, "DG L2", 0) - dg1l2 = FunctionSpace(mesh, "DG L2", 1) return {'cg1': cg1, 'vcg1': vcg1, 'tcg1': tcg1, @@ -38,13 +34,7 @@ def fs(request, mesh): 'cg1dg0[1]': (cg1*dg0)[1], 'cg2dg1': cg2*dg1, 'cg2dg1[0]': (cg2*dg1)[0], - 'cg2dg1[1]': (cg2*dg1)[1], - 'cg1dg0l2': cg1*dg0l2, - 'cg1dg0l2[0]': (cg1*dg0l2)[0], - 'cg1dg0l2[1]': (cg1*dg0l2)[1], - 'cg2dg1l2': cg2*dg1l2, - 'cg2dg1l2[0]': (cg2*dg1l2)[0], - 'cg2dg1l2[1]': (cg2*dg1l2)[1]}[request.param] + 'cg2dg1[1]': (cg2*dg1)[1]}[request.param] @pytest.fixture @@ -132,14 +122,3 @@ def test_assemble_mat_with_tensor(mesh): M = assemble(Constant(2)*a, M) # Make sure we get the result of the last assembly assert np.allclose(M.M.values, 2*assemble(a).M.values, rtol=1e-14) - -def test_assemble_mat_with_tensor_l2(mesh): - V = FunctionSpace(mesh, "DG L2", 0) - u = TestFunction(V) - v = TrialFunction(V) - a = u*v*dx - M = assemble(a) - # Assemble a different form into M - M = assemble(Constant(2)*a, M) - # Make sure we get the result of the last assembly - assert np.allclose(M.M.values, 2*assemble(a).M.values, rtol=1e-14) diff --git a/tests/regression/test_assemble_inverse.py b/tests/regression/test_assemble_inverse.py index d13cd7e5b0..66fa8664d6 100644 --- a/tests/regression/test_assemble_inverse.py +++ b/tests/regression/test_assemble_inverse.py @@ -16,17 +16,3 @@ def test_assemble_inverse(degree): eye = np.dot(m_forward.M.values, m_inverse.M.values) assert ((eye - np.eye(fs.node_count)).round(12) == 0).all() - -@pytest.mark.parametrize("degree", range(4)) -def test_assemble_inverse_l2(degree): - m = UnitSquareMesh(2, 1) - fs = FunctionSpace(m, "DG L2", degree) - u = TrialFunction(fs) - v = TestFunction(fs) - - m_forward = assemble(u*v*dx) - m_inverse = assemble(u*v*dx, inverse=True) - - eye = np.dot(m_forward.M.values, m_inverse.M.values) - - assert ((eye - np.eye(fs.node_count)).round(12) == 0).all() diff --git a/tests/regression/test_function_spaces.py b/tests/regression/test_function_spaces.py index 705c991102..3d90a42306 100644 --- a/tests/regression/test_function_spaces.py +++ b/tests/regression/test_function_spaces.py @@ -26,18 +26,13 @@ def cg2(mesh): def dg0(mesh): return VectorFunctionSpace(mesh, "DG", 0) -@pytest.fixture(scope="module") -def dg0l2(mesh): - return VectorFunctionSpace(mesh, "DG L2", 0) -@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0', 'cg1dg0l2', 'cg2dg0l2']) +@pytest.fixture(scope='module', params=['cg1cg1', 'cg1cg2', 'cg1dg0', 'cg2dg0']) def fs(request, cg1, cg2, dg0): return {'cg1cg1': cg1*cg1, 'cg1cg2': cg1*cg2, 'cg1dg0': cg1*dg0, - 'cg2dg0': cg2*dg0, - 'cg1dg0l2': cg1*dg0l2, - 'cg2dg0l2': cg2*dg0l2}[request.param] + 'cg2dg0': cg2*dg0}[request.param] def test_function_space_cached(mesh): @@ -110,9 +105,7 @@ def test_function_space_collapse(cg1): TensorProductElement(FiniteElement("CG", triangle, 1), FiniteElement("CG", interval, 3)), MixedElement(FiniteElement("CG", triangle, 1), - FiniteElement("DG", triangle, 2)), - MixedElement(FiniteElement("CG", triangle, 1), - FiniteElement("DG L2", triangle, 2))], + FiniteElement("DG", triangle, 2))], ids=["FE", "Enriched", "TPE", "Mixed"]) def test_validation(modifier, element): with pytest.raises(ValueError): diff --git a/tests/regression/test_helmholtz_mixed.py b/tests/regression/test_helmholtz_mixed.py index e145ef32f6..d9472bfe7b 100644 --- a/tests/regression/test_helmholtz_mixed.py +++ b/tests/regression/test_helmholtz_mixed.py @@ -48,11 +48,7 @@ def helmholtz_mixed(r, V1, V2, action=False): [(('RT', 1), ('DG', 0), 1.9, False), (('BDM', 1), ('DG', 0), 1.89, False), (('BDM', 1), ('DG', 0), 1.89, True), - (('BDFM', 2), ('DG', 1), 1.9, False), - (('RT', 1), ('DG L2', 0), 1.9, False), - (('BDM', 1), ('DG L2', 0), 1.89, False), - (('BDM', 1), ('DG L2', 0), 1.89, True), - (('BDFM', 2), ('DG L2', 1), 1.9, False)]) + (('BDFM', 2), ('DG', 1), 1.9, False)]) def test_firedrake_helmholtz(V1, V2, threshold, action): import numpy as np diff = np.array([helmholtz_mixed(i, V1, V2) for i in range(3, 6)]) diff --git a/tests/regression/test_interpolate.py b/tests/regression/test_interpolate.py index 13242e7877..90571aed82 100644 --- a/tests/regression/test_interpolate.py +++ b/tests/regression/test_interpolate.py @@ -122,21 +122,6 @@ def test_cell_orientation(): assert abs(g.dat.data - h.dat.data).max() < 1e-2 -def test_cell_orientation_l2(): - m = UnitCubedSphereMesh(2) - x = SpatialCoordinate(m) - m.init_cell_orientations(x) - x = m.coordinates - U = FunctionSpace(m, 'RTCF', 1) - V = VectorFunctionSpace(m, 'DQ L2', 1) - - f = project(as_tensor([x[1], -x[0], 0.0]), U) - g = interpolate(f, V) - - # g shall be close to: - h = project(f, V) - - assert abs(g.dat.data - h.dat.data).max() < 1e-2 def test_cellvolume(): m = UnitSquareMesh(2, 2) @@ -146,13 +131,6 @@ def test_cellvolume(): assert np.allclose(f.dat.data_ro, 0.125) -def test_cellvolume_l2(): - m = UnitSquareMesh(2, 2) - V = FunctionSpace(m, 'DG L2', 0) - - f = interpolate(CellVolume(m), V) - - assert np.allclose(f.dat.data_ro, 0.125) def test_cellvolume_higher_order_coords(): m = UnitTriangleMesh() @@ -173,24 +151,6 @@ def warp(x): assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) -def test_cellvolume_higher_order_coords_l2(): - m = UnitTriangleMesh() - V = VectorFunctionSpace(m, 'P', 3) - f = Function(V) - f.interpolate(m.coordinates) - - # Warp mesh so that the bottom triangle line is: - # x(x - 1)(x + a) with a = 19/12.0 - def warp(x): - return x * (x - 1)*(x + 19/12.0) - - f.dat.data[1, 1] = warp(1.0/3.0) - f.dat.data[2, 1] = warp(2.0/3.0) - - mesh = Mesh(f) - g = interpolate(CellVolume(mesh), FunctionSpace(mesh, 'DG L2', 0)) - - assert np.allclose(g.dat.data_ro, 0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0)) def test_mixed(): m = UnitTriangleMesh() diff --git a/tests/regression/test_stokes_hdiv_parallel.py b/tests/regression/test_stokes_hdiv_parallel.py index 9a2f0ff6bd..223f260118 100644 --- a/tests/regression/test_stokes_hdiv_parallel.py +++ b/tests/regression/test_stokes_hdiv_parallel.py @@ -9,10 +9,8 @@ def mat_type(request): @pytest.fixture(params=[(("RT", 3), ("DG", 2)), - (("BDM", 2), ("DG", 1)), - (("RT", 3), ("DG L2", 2)), - (("BDM", 2), ("DG L2", 1))], - ids=["RT3-DG2", "BDM2-DG1", "RT3-DG2L2", "BDM2-DG1L2"]) + (("BDM", 2), ("DG", 1))], + ids=["RT3-DG2", "BDM2-DG1"]) def element_pair(request): return request.param From 74ec42fb7184e28ef24347f03ab316c0fa1615b6 Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Wed, 26 Jun 2019 17:10:21 +0100 Subject: [PATCH 09/11] proper testing --- tests/regression/test_l2pullbacks.py | 239 ++++++++++++++++++++++++++- 1 file changed, 237 insertions(+), 2 deletions(-) diff --git a/tests/regression/test_l2pullbacks.py b/tests/regression/test_l2pullbacks.py index 6e8a6e7cdf..3c92d45bc5 100644 --- a/tests/regression/test_l2pullbacks.py +++ b/tests/regression/test_l2pullbacks.py @@ -1,7 +1,242 @@ +import pytest +import numpy as np +from firedrake import * -# assembly/projections/function spaces/etc. +def get_complex(family, hdegree, vdegree=None): -# laplace/helmholtz mixed on plane and sphere + if family == 'RT': + l2elem = FiniteElement('DG L2', triangle, hdegree-1) + hdivelem = FiniteElement('RT', triangle, hdegree) + elif family == 'BDM': + l2elem = FiniteElement('DG L2', triangle, hdegree-1) + hdivelem = FiniteElement('BDM', triangle, hdegree) + elif family == 'BDFM': + l2elem = FiniteElement('DG L2', triangle, hdegree-1) + hdivelem = FiniteElement('BDFM', triangle, hdegree) + elif family == 'RTCF': + l2elem = FiniteElement('DQ L2', quadrilateral, hdegree-1) + hdivelem = FiniteElement('RTCF', quadrilateral, hdegree) + elif family == 'ext': + vdegree = vdegree or hdegree + h1_horiz = FiniteElement("CG", "interval", hdegree) + l2_horiz = FiniteElement("DG L2", "interval", hdegree-1) + h1_vert = FiniteElement("CG", "interval", vdegree) + l2_vert = FiniteElement("DG L2", "interval", vdegree-1) + hdivelem = HDiv(TensorProductElement(h1_horiz, l2_vert)) + HDiv(TensorProductElement(l2_horiz, h1_vert)) + l2elem = TensorProductElement(l2_horiz, l2_vert) + return hdivelem, l2elem + + +def get_mesh(r, meshtype, meshd=None): + + meshd = meshd or 1 + if meshtype == 'spherical-quad': + mesh = UnitCubedSphereMesh(refinement_level=r, degree=meshd) + x = SpatialCoordinate(mesh) + mesh.init_cell_orientations(x) + elif meshtype == 'spherical-tri': + mesh = UnitIcosahedralSphereMesh(refinement_level=r, degree=meshd) + x = SpatialCoordinate(mesh) + mesh.init_cell_orientations(x) + elif meshtype == 'planar-quad': + mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=True) + elif meshtype == 'planar-tri': + mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=False) + elif meshtype == 'extruded': + basemesh = UnitIntervalMesh(2**r) + mesh = ExtrudedMesh(basemesh, 2**r) + return mesh + + +def helmholtz_mixed(r, meshtype, family, hdegree, vdegree=None, meshd=None, useaction=False): + + mesh = get_mesh(r, meshtype, meshd=meshd) + hdivelem, l2elem = get_complex(family, hdegree, vdegree=vdegree) + + V1 = FunctionSpace(mesh, hdivelem, name="V") + V2 = FunctionSpace(mesh, l2elem, name="P") + W = V1 * V2 + + # Define variational problem + lmbda = 1 + u, p = TrialFunctions(W) + v, q = TestFunctions(W) + f = Function(V2) + + x = SpatialCoordinate(mesh) + if meshtype in ['planar-quad', 'planar-tri', 'extruded']: + f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype in ['spherical-quad', 'spherical-tri']: + f.project(x[0]*x[1]*x[2]) + a = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx + L = f*q*dx + + # Compute solution + sol = Function(W) + + if useaction: + system = action(a, sol) - L == 0 + else: + system = a == L + + # Block system is: + # V Ct + # Ch P + # Eliminate V by forming a schur complement + solve(system, sol, solver_parameters={'pc_type': 'fieldsplit', + 'pc_fieldsplit_type': 'schur', + 'ksp_type': 'cg', + 'pc_fieldsplit_schur_fact_type': 'FULL', + 'fieldsplit_V_ksp_type': 'cg', + 'fieldsplit_P_ksp_type': 'cg'}) + + # Analytical solution + + if meshtype in ['planar-quad', 'planar-tri', 'extruded']: + f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) + return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx)) + elif meshtype in ['spherical-quad', 'spherical-tri']: + _, u = sol.split() + f.project(x[0]*x[1]*x[2]/13.0) + return errornorm(f, u, degree_rise=0) + + +# helmholtz mixed on plane + +@pytest.mark.parametrize(('family', 'degree', 'celltype', 'action', 'threshold'), + [('RT', 1, 'tri', False, 1.9), + ('RT', 2, 'tri', False, 2.9), + ('BDM', 1, 'tri', False, 1.87), + ('BDM', 2, 'tri', False, 2.9), + ('BDFM', 2, 'tri', False, 2.9), + ('RTCF', 1, 'quad', False, 1.9), + ('RTCF', 2, 'quad', False, 2.9), + ('RT', 2, 'tri', True, 2.9), + ('BDM', 2, 'tri', True, 2.9)]) +def test_firedrake_helmholtz_mixed_l2pullbacks(family, degree, celltype, action, threshold): + diff = np.array([helmholtz_mixed(r, 'planar-' + celltype, family, degree, useaction=action) for r in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz mixed on sphere + +@pytest.mark.parametrize(('family', 'degree', 'celltype', 'md', 'action', 'threshold'), + [('RT', 1, 'tri', 1, False, 1.9), + ('RT', 2, 'tri', 1, False, 1.9), + ('RT', 2, 'tri', 2, False, 2.9), + ('BDM', 1, 'tri', 1, False, 1.88), + ('BDM', 2, 'tri', 1, False, 1.9), + ('BDM', 2, 'tri', 2, False, 2.9), + ('BDFM', 2, 'tri', 1, False, 1.9), + ('BDFM', 2, 'tri', 2, False, 2.9), + ('RTCF', 1, 'quad', 1, False, 1.67), + ('RTCF', 2, 'quad', 1, False, 1.9), + ('RTCF', 2, 'quad', 2, False, 2.9), + ('RT', 2, 'tri', 1, True, 1.9), + ('BDM', 2, 'tri', 1, True, 1.9)]) +def test_firedrake_helmholtz_mixed_l2pullbacks_sphere(family, degree, celltype, md, action, threshold): + diff = np.array([helmholtz_mixed(r, 'spherical-' + celltype, family, degree, meshd=md, useaction=action) for r in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz mixed on extruded + +@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold', 'action'), + [(1, 1, 1.9, False), + (2, 1, 1.9, False), + (1, 2, 1.9, False), + (2, 2, 2.9, False), + (1, 1, 1.9, True)]) +def test_firedrake_helmholtz_mixed_l2pullbacks_extruded(hdegree, vdegree, threshold, action): + diff = np.array([helmholtz_mixed(i, 'extruded', 'ext', hdegree, vdegree=vdegree, useaction=action) for i in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() # facet integrals- interior and exterior +# tested using an interior penalty formulation for the Helmholtz using DG/DQ L2 elements + + +def laplacian_IP(r, degree, meshd, meshtype): + + dIF = dS + dEF = ds + if meshtype == 'extruded': + dIF = dS_v + dS_h + dEF = ds_t + ds_b + ds_v + hdivelem, l2elem = get_complex('ext', degree+1, vdegree=degree+1) + elif meshtype in ['spherical-quad', 'planar-quad']: + hdivelem, l2elem = get_complex('RTCF', degree+1) + elif meshtype in ['spherical-tri', 'planar-tri']: + hdivelem, l2elem = get_complex('RT', degree+1) + + mesh = get_mesh(r, meshtype, meshd=meshd) + + V = FunctionSpace(mesh, l2elem) + + # Define variational problem + u = TrialFunction(V) + v = TestFunction(V) + f = Function(V) + + x = SpatialCoordinate(mesh) + if meshtype in ['planar-quad', 'planar-tri', 'extruded']: + f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype in ['spherical-quad', 'spherical-tri']: + f.project(x[0]*x[1]*x[2]) + + FA = FacetArea(mesh) + CV = CellVolume(mesh) + ddx = CV/FA + ddx_avg = (ddx('+') + ddx('-'))/2. + alpha = Constant(4. * degree * (degree + 1.)) + gamma = Constant(8. * degree * (degree + 1.)) + penalty_int = alpha / ddx_avg + penalty_ext = gamma / ddx + + n = FacetNormal(mesh) + aV = (inner(grad(u), grad(v)) + inner(u, v)) * dx # volume term + aIF = (inner(jump(u, n), jump(v, n)) * penalty_int - inner(avg(grad(u)), jump(v, n)) - inner(avg(grad(v)), jump(u, n))) * dIF # interior facet term + aEF = (inner(u, v) * penalty_ext - inner(grad(u), v*n) - inner(grad(v), u*n)) * dEF # exterior facet term + a = aV + aEF + aIF + L = f*v*dx + + # Compute solution + sol = Function(V) + + solve(a == L, sol, solver_parameters={'pc_type': 'ilu', + 'ksp_type': 'lgmres'}) + + # Analytical solution + if meshtype in ['planar-quad', 'planar-tri', 'extruded']: + f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype in ['spherical-quad', 'spherical-tri']: + f.project(x[0]*x[1]*x[2]/13.0) + return sqrt(assemble(dot(sol - f, sol - f) * dx)) + + +@pytest.mark.parametrize(('degree', 'meshd', 'meshtype', 'threshold'), + [(1, 1, 'planar-quad', 1.8), + (2, 1, 'planar-quad', 2.9), + (1, 1, 'planar-tri', 1.5), + (2, 1, 'planar-tri', 2.9), + (1, 1, 'extruded', 1.8), + (2, 1, 'extruded', 2.9), + (2, 2, 'spherical-quad', 2.9), + (3, 3, 'spherical-quad', 3.9), + (2, 2, 'spherical-tri', 2.9), + (3, 3, 'spherical-tri', 3.9)]) +def test_firedrake_laplacian_IP_l2pullbacks(degree, meshd, meshtype, threshold): + diff = np.array([laplacian_IP(i, degree, meshd, meshtype) for i in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() From b52bafe00831e39cc9258367730badf42ef407ad Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Wed, 26 Jun 2019 17:46:00 +0100 Subject: [PATCH 10/11] DROP BEFORE MERGE --- Jenkinsfile | 2 +- requirements-git.txt | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Jenkinsfile b/Jenkinsfile index 9f75ae1a09..fc8446c4e4 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -26,7 +26,7 @@ pipeline { sh 'mkdir tmp' dir('tmp') { timestamps { - sh '../scripts/firedrake-install --package-branch ufl l2pullbacks --package-branch tsfc l2pullbacks --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)' + sh '../scripts/firedrake-install --package-branch ufl mse --package-branch tsfc mse --package-branch fiat mse --package-branch FInAT mse --disable-ssh --minimal-petsc --slepc --documentation-dependencies --install thetis --install gusto --install icepack --install pyadjoint --no-package-manager || (cat firedrake-install.log && /bin/false)' } } } diff --git a/requirements-git.txt b/requirements-git.txt index d935e8dcd8..96651ca2f9 100644 --- a/requirements-git.txt +++ b/requirements-git.txt @@ -1,5 +1,5 @@ -git+https://github.com/firedrakeproject/ufl.git#egg=ufl -git+https://github.com/firedrakeproject/fiat.git#egg=fiat -git+https://github.com/FInAT/FInAT.git#egg=finat -git+https://github.com/firedrakeproject/tsfc.git#egg=tsfc +git+https://github.com/firedrakeproject/ufl.git@mse#egg=ufl +git+https://github.com/firedrakeproject/fiat.git@mse#egg=fiat +git+https://github.com/FInAT/FInAT.git@mse#egg=finat +git+https://github.com/firedrakeproject/tsfc.git@mse#egg=tsfc git+https://github.com/OP2/PyOP2.git#egg=pyop2 From eb8e3d286c0f1f0e216712467fd6a0da5f7f622b Mon Sep 17 00:00:00 2001 From: Chris Eldred Date: Thu, 27 Jun 2019 18:55:46 +0100 Subject: [PATCH 11/11] testing for mse and dualmse spaces --- tests/regression/test_dualmse.py | 365 +++++++++++++++++++++++++++++++ tests/regression/test_mse.py | 291 ++++++++++++++++++++++++ 2 files changed, 656 insertions(+) create mode 100644 tests/regression/test_dualmse.py create mode 100644 tests/regression/test_mse.py diff --git a/tests/regression/test_dualmse.py b/tests/regression/test_dualmse.py new file mode 100644 index 0000000000..333757c6bb --- /dev/null +++ b/tests/regression/test_dualmse.py @@ -0,0 +1,365 @@ +import pytest +import numpy as np +from firedrake import * + + +def get_complex(meshtype, hdegree, vdegree=None): + + if meshtype == 'extruded': + vdegree = vdegree or hdegree + h1_horiz = FiniteElement("EGL", "interval", hdegree) + l2_horiz = FiniteElement("EGL-Edge L2", "interval", hdegree-1) + h1_vert = FiniteElement("EGL", "interval", vdegree) + l2_vert = FiniteElement("EGL-Edge L2", "interval", vdegree-1) + hcurlelem = HCurl(TensorProductElement(h1_horiz, l2_vert)) + HCurl(TensorProductElement(l2_horiz, h1_vert)) + hdivelem = HDiv(TensorProductElement(h1_horiz, l2_vert)) + HDiv(TensorProductElement(l2_horiz, h1_vert)) + l2elem = TensorProductElement(l2_horiz, l2_vert) + h1elem = TensorProductElement(h1_horiz, h1_vert) + else: + l2elem = FiniteElement('DQ L2', quadrilateral, hdegree-1, variant='dualmse') + hdivelem = FiniteElement('RTCF', quadrilateral, hdegree, variant='dualmse') + hcurlelem = FiniteElement('RTCE', quadrilateral, hdegree, variant='dualmse') + h1elem = FiniteElement('Q', quadrilateral, hdegree, variant='dualmse') + + return h1elem, hcurlelem, hdivelem, l2elem + + +def get_mesh(r, meshtype, meshd=None): + + meshd = meshd or 1 + if meshtype == 'spherical': + mesh = UnitCubedSphereMesh(refinement_level=r, degree=meshd) + x = SpatialCoordinate(mesh) + mesh.init_cell_orientations(x) + elif meshtype == 'planar': + mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=True) + elif meshtype == 'extruded': + basemesh = UnitIntervalMesh(2**r) + mesh = ExtrudedMesh(basemesh, 2**r) + return mesh + + +def helmholtz(r, meshtype, hdegree, vdegree=None, meshd=None): + + mesh = get_mesh(r, meshtype, meshd=meshd) + h1elem, _, _, _ = get_complex(meshtype, hdegree, vdegree=vdegree) + + V = FunctionSpace(mesh, h1elem) + + # Define variational problem + lmbda = 1 + u = TrialFunction(V) + v = TestFunction(V) + f = Function(V) + x = SpatialCoordinate(mesh) + if meshtype in ['planar', 'extruded']: + f.project((1+8*pi*pi)*cos(x[0]*pi*2)*cos(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]) + a = (dot(grad(v), grad(u)) + lmbda * v * u) * dx + L = f * v * dx + + # Compute solution + assemble(a) + assemble(L) + sol = Function(V) + solve(a == L, sol, solver_parameters={'ksp_type': 'cg'}) + + # Analytical solution + if meshtype in ['planar', 'extruded']: + f.project(cos(x[0]*pi*2)*cos(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]/13.0) + return sqrt(assemble(dot(sol - f, sol - f) * dx)) + +# helmholtz on plane + + +@pytest.mark.parametrize(('degree', 'threshold'), + [(2, 2.9), + (3, 3.9)]) +def test_firedrake_helmholtz_dualmse(degree, threshold): + diff = np.array([helmholtz(r, 'planar', degree) for r in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz on sphere + +@pytest.mark.parametrize(('degree', 'md', 'threshold'), + [(2, 1, 1.9), + (2, 2, 2.9), + (3, 1, 1.87), + (3, 2, 2.9), + (3, 3, 3.9)]) +def test_firedrake_helmholtz_dualmse_sphere(degree, md, threshold): + diff = np.array([helmholtz(r, 'spherical', degree, meshd=md) for r in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz on extruded + +@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold'), + [(2, 2, 2.9), + (2, 3, 2.9), + (3, 2, 2.9), + (3, 3, 3.9)]) +def test_firedrake_helmholtz_dualmse_extruded(hdegree, vdegree, threshold): + diff = np.array([helmholtz(i, 'extruded', hdegree, vdegree=vdegree) for i in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +def helmholtz_mixed(r, meshtype, hdegree, vdegree=None, meshd=None, useaction=False): + + mesh = get_mesh(r, meshtype, meshd=meshd) + _, _, hdivelem, l2elem = get_complex(meshtype, hdegree, vdegree=vdegree) + + V1 = FunctionSpace(mesh, hdivelem, name="V") + V2 = FunctionSpace(mesh, l2elem, name="P") + W = V1 * V2 + + # Define variational problem + lmbda = 1 + u, p = TrialFunctions(W) + v, q = TestFunctions(W) + f = Function(V2) + + x = SpatialCoordinate(mesh) + if meshtype in ['planar', 'extruded']: + f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]) + a = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx + L = f*q*dx + + # Compute solution + sol = Function(W) + + if useaction: + system = action(a, sol) - L == 0 + else: + system = a == L + + # Block system is: + # V Ct + # Ch P + # Eliminate V by forming a schur complement + solve(system, sol, solver_parameters={'pc_type': 'fieldsplit', + 'pc_fieldsplit_type': 'schur', + 'ksp_type': 'cg', + 'pc_fieldsplit_schur_fact_type': 'FULL', + 'fieldsplit_V_ksp_type': 'cg', + 'fieldsplit_P_ksp_type': 'cg'}) + + if meshtype in ['planar', 'extruded']: + f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) + return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx)) + elif meshtype == 'spherical': + _, u = sol.split() + f.project(x[0]*x[1]*x[2]/13.0) + return errornorm(f, u, degree_rise=0) + + +# helmholtz mixed on plane + +@pytest.mark.parametrize(('degree', 'action', 'threshold'), + [(2, False, 2.9), + (3, False, 3.9), + (2, True, 2.9)]) +def test_firedrake_helmholtz_mixed_dualmse(degree, action, threshold): + diff = np.array([helmholtz_mixed(r, 'planar', degree, useaction=action) for r in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz mixed on sphere + +@pytest.mark.parametrize(('degree', 'md', 'action', 'threshold'), + [(2, 1, False, 1.9), + (2, 2, False, 2.9), + (3, 1, False, 1.87), + (3, 2, False, 2.9), + (3, 3, False, 3.9), + (2, 2, True, 2.9)]) +def test_firedrake_helmholtz_mixed_dualmse_sphere(degree, md, action, threshold): + diff = np.array([helmholtz_mixed(r, 'spherical', degree, meshd=md, useaction=action) for r in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz mixed on extruded + +@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold', 'action'), + [(2, 2, 2.9, False), + (2, 3, 2.9, False), + (3, 2, 2.9, False), + (3, 3, 3.9, False), + (2, 2, 2.9, True)]) +def test_firedrake_helmholtz_mixed_dualmse_extruded(hdegree, vdegree, threshold, action): + diff = np.array([helmholtz_mixed(i, 'extruded', hdegree, vdegree=vdegree, useaction=action) for i in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# facet integrals- interior and exterior +# tested using an interior penalty formulation for the Helmholtz using DG/DQ L2 elements + + +def laplacian_IP(r, degree, meshd, meshtype): + + dIF = dS + dEF = ds + if meshtype == 'extruded': + dIF = dS_v + dS_h + dEF = ds_t + ds_b + ds_v + + mesh = get_mesh(r, meshtype, meshd=meshd) + _, _, _, l2elem = get_complex(meshtype, degree+1, vdegree=degree+1) + + V = FunctionSpace(mesh, l2elem) + + # Define variational problem + u = TrialFunction(V) + v = TestFunction(V) + f = Function(V) + + x = SpatialCoordinate(mesh) + if meshtype in ['planar', 'extruded']: + f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]) + + FA = FacetArea(mesh) + CV = CellVolume(mesh) + ddx = CV/FA + ddx_avg = (ddx('+') + ddx('-'))/2. + alpha = Constant(4. * degree * (degree + 1.)) + gamma = Constant(8. * degree * (degree + 1.)) + penalty_int = alpha / ddx_avg + penalty_ext = gamma / ddx + + n = FacetNormal(mesh) + aV = (inner(grad(u), grad(v)) + inner(u, v)) * dx # volume term + aIF = (inner(jump(u, n), jump(v, n)) * penalty_int - inner(avg(grad(u)), jump(v, n)) - inner(avg(grad(v)), jump(u, n))) * dIF # interior facet term + aEF = (inner(u, v) * penalty_ext - inner(grad(u), v*n) - inner(grad(v), u*n)) * dEF # exterior facet term + a = aV + aEF + aIF + L = f*v*dx + + # Compute solution + sol = Function(V) + + solve(a == L, sol, solver_parameters={'pc_type': 'ilu', + 'ksp_type': 'lgmres'}) + + # Analytical solution + if meshtype in ['planar', 'extruded']: + f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]/13.0) + return sqrt(assemble(dot(sol - f, sol - f) * dx)) + + +@pytest.mark.parametrize(('degree', 'meshd', 'meshtype', 'threshold'), + [(2, 1, 'planar', 2.9), + (3, 1, 'planar', 3.8), + (2, 1, 'extruded', 2.9), + (3, 1, 'extruded', 3.8), + (2, 2, 'spherical', 2.9), + (3, 3, 'spherical', 3.9)]) +def test_firedrake_laplacian_IP_dualmse(degree, meshd, meshtype, threshold): + diff = np.array([laplacian_IP(i, degree, meshd, meshtype) for i in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +def vector_laplace(r, meshtype, hdegree, vdegree=None): + vdegree = vdegree or hdegree + + mesh = get_mesh(r, meshtype) + h1elem, hcurlelem, _, _ = get_complex(meshtype, hdegree, vdegree=vdegree) + h1elem_high, _, _, _ = get_complex(meshtype, hdegree+1, vdegree=vdegree+1) + + # spaces for calculation + V0 = FunctionSpace(mesh, h1elem) + V1 = FunctionSpace(mesh, hcurlelem) + V = V0*V1 + + # spaces to store 'analytic' functions + W0 = FunctionSpace(mesh, h1elem_high) + W1 = VectorFunctionSpace(mesh, h1elem_high) + + # constants + k = 1.0 + l = 2.0 + + xs = SpatialCoordinate(mesh) + f_expr = as_vector([pi*pi*(k*k + l*l)*sin(k*pi*xs[0])*cos(l*pi*xs[1]), pi*pi*(k*k + l*l)*cos(k*pi*xs[0])*sin(l*pi*xs[1])]) + exact_s_expr = -(k+l)*pi*cos(k*pi*xs[0])*cos(l*pi*xs[1]) + exact_u_expr = as_vector([sin(k*pi*xs[0])*cos(l*pi*xs[1]), cos(k*pi*xs[0])*sin(l*pi*xs[1])]) + + f = Function(W1).project(f_expr) + exact_s = Function(W0).project(exact_s_expr) + exact_u = Function(W1).project(exact_u_expr) + + sigma, u = TrialFunctions(V) + tau, v = TestFunctions(V) + a = (sigma*tau - dot(u, grad(tau)) + dot(grad(sigma), v) + dot(curl(u), curl(v)))*dx + L = dot(f, v)*dx + + out = Function(V) + + # preconditioner for H1 x H(curl) + aP = (dot(grad(sigma), grad(tau)) + sigma*tau + dot(curl(u), curl(v)) + dot(u, v))*dx + + solve(a == L, out, Jp=aP, + solver_parameters={'pc_type': 'fieldsplit', + 'pc_fieldsplit_type': 'additive', + 'fieldsplit_0_pc_type': 'lu', + 'fieldsplit_1_pc_type': 'lu', + 'ksp_monitor': None}) + + out_s, out_u = out.split() + + return (sqrt(assemble(dot(out_u - exact_u, out_u - exact_u)*dx)), + sqrt(assemble((out_s - exact_s)*(out_s - exact_s)*dx))) + + +@pytest.mark.parametrize(('degree', 'threshold'), + [(2, 1.9), + (3, 2.9), + (2, 1.9)]) +def test_firedrake_helmholtz_vector_dualmse(degree, threshold): + diff = np.array([vector_laplace(r, 'planar', degree) for r in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold'), + [(2, 2, 1.9), + (2, 3, 1.9), + (3, 2, 1.9), + (3, 3, 2.9)]) +def test_firedrake_helmholtz_vector_dualmse_extruded(hdegree, vdegree, threshold): + diff = np.array([vector_laplace(i, 'extruded', hdegree, vdegree=vdegree) for i in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1, :] / diff[1:, :]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() diff --git a/tests/regression/test_mse.py b/tests/regression/test_mse.py new file mode 100644 index 0000000000..2776686094 --- /dev/null +++ b/tests/regression/test_mse.py @@ -0,0 +1,291 @@ +import pytest +import numpy as np +from firedrake import * + + +def get_complex(meshtype, hdegree, vdegree=None): + + if meshtype == 'extruded': + vdegree = vdegree or hdegree + h1_horiz = FiniteElement("GLL", "interval", hdegree) + l2_horiz = FiniteElement("GLL-Edge L2", "interval", hdegree-1) + h1_vert = FiniteElement("GLL", "interval", vdegree) + l2_vert = FiniteElement("GLL-Edge L2", "interval", vdegree-1) + hcurlelem = HCurl(TensorProductElement(h1_horiz, l2_vert)) + HCurl(TensorProductElement(l2_horiz, h1_vert)) + hdivelem = HDiv(TensorProductElement(h1_horiz, l2_vert)) + HDiv(TensorProductElement(l2_horiz, h1_vert)) + l2elem = TensorProductElement(l2_horiz, l2_vert) + h1elem = TensorProductElement(h1_horiz, h1_vert) + else: + l2elem = FiniteElement('DQ L2', quadrilateral, hdegree-1, variant='mse') + hdivelem = FiniteElement('RTCF', quadrilateral, hdegree, variant='mse') + hcurlelem = FiniteElement('RTCE', quadrilateral, hdegree, variant='mse') + h1elem = FiniteElement('Q', quadrilateral, hdegree, variant='mse') + + return h1elem, hcurlelem, hdivelem, l2elem + + +def get_mesh(r, meshtype, meshd=None): + + meshd = meshd or 1 + if meshtype == 'spherical': + mesh = UnitCubedSphereMesh(refinement_level=r, degree=meshd) + x = SpatialCoordinate(mesh) + mesh.init_cell_orientations(x) + elif meshtype == 'planar': + mesh = UnitSquareMesh(2**r, 2**r, quadrilateral=True) + elif meshtype == 'extruded': + basemesh = UnitIntervalMesh(2**r) + mesh = ExtrudedMesh(basemesh, 2**r) + return mesh + + +def helmholtz_mixed(r, meshtype, hdegree, vdegree=None, meshd=None, useaction=False): + + mesh = get_mesh(r, meshtype, meshd=meshd) + _, _, hdivelem, l2elem = get_complex(meshtype, hdegree, vdegree=vdegree) + + V1 = FunctionSpace(mesh, hdivelem, name="V") + V2 = FunctionSpace(mesh, l2elem, name="P") + W = V1 * V2 + + # Define variational problem + lmbda = 1 + u, p = TrialFunctions(W) + v, q = TestFunctions(W) + f = Function(V2) + + x = SpatialCoordinate(mesh) + if meshtype in ['planar', 'extruded']: + f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]) + a = (p*q - q*div(u) + lmbda*inner(v, u) + div(v)*p) * dx + L = f*q*dx + + # Compute solution + sol = Function(W) + + if useaction: + system = action(a, sol) - L == 0 + else: + system = a == L + + # Block system is: + # V Ct + # Ch P + # Eliminate V by forming a schur complement + solve(system, sol, solver_parameters={'pc_type': 'fieldsplit', + 'pc_fieldsplit_type': 'schur', + 'ksp_type': 'cg', + 'pc_fieldsplit_schur_fact_type': 'FULL', + 'fieldsplit_V_ksp_type': 'cg', + 'fieldsplit_P_ksp_type': 'cg'}) + + if meshtype in ['planar', 'extruded']: + f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) + return sqrt(assemble(dot(sol[2] - f, sol[2] - f) * dx)) + elif meshtype == 'spherical': + _, u = sol.split() + f.project(x[0]*x[1]*x[2]/13.0) + return errornorm(f, u, degree_rise=0) + + +# helmholtz mixed on plane + +@pytest.mark.parametrize(('degree', 'action', 'threshold'), + [(1, False, 1.9), + (2, False, 2.9), + (3, False, 3.9), + (2, True, 2.9)]) +def test_firedrake_helmholtz_mixed_mse(degree, action, threshold): + diff = np.array([helmholtz_mixed(r, 'planar', degree, useaction=action) for r in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz mixed on sphere + +@pytest.mark.parametrize(('degree', 'md', 'action', 'threshold'), + [(1, 1, False, 1.67), + (2, 1, False, 1.9), + (2, 2, False, 2.9), + (3, 1, False, 1.87), + (3, 2, False, 2.9), + (3, 3, False, 3.9), + (2, 2, True, 2.9)]) +def test_firedrake_helmholtz_mixed_mse_sphere(degree, md, action, threshold): + diff = np.array([helmholtz_mixed(r, 'spherical', degree, meshd=md, useaction=action) for r in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# helmholtz mixed on extruded + +@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold', 'action'), + [(1, 1, 1.9, False), + (2, 1, 1.9, False), + (1, 2, 1.9, False), + (2, 2, 2.9, False), + (1, 1, 1.9, True)]) +def test_firedrake_helmholtz_mixed_mse_extruded(hdegree, vdegree, threshold, action): + diff = np.array([helmholtz_mixed(i, 'extruded', hdegree, vdegree=vdegree, useaction=action) for i in range(3, 6)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +# facet integrals- interior and exterior +# tested using an interior penalty formulation for the Helmholtz using DG/DQ L2 elements + + +def laplacian_IP(r, degree, meshd, meshtype): + + dIF = dS + dEF = ds + if meshtype == 'extruded': + dIF = dS_v + dS_h + dEF = ds_t + ds_b + ds_v + + mesh = get_mesh(r, meshtype, meshd=meshd) + _, _, _, l2elem = get_complex(meshtype, degree+1, vdegree=degree+1) + + V = FunctionSpace(mesh, l2elem) + + # Define variational problem + u = TrialFunction(V) + v = TestFunction(V) + f = Function(V) + + x = SpatialCoordinate(mesh) + if meshtype in ['planar', 'extruded']: + f.project((1+8*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]) + + FA = FacetArea(mesh) + CV = CellVolume(mesh) + ddx = CV/FA + ddx_avg = (ddx('+') + ddx('-'))/2. + alpha = Constant(4. * degree * (degree + 1.)) + gamma = Constant(8. * degree * (degree + 1.)) + penalty_int = alpha / ddx_avg + penalty_ext = gamma / ddx + + n = FacetNormal(mesh) + aV = (inner(grad(u), grad(v)) + inner(u, v)) * dx # volume term + aIF = (inner(jump(u, n), jump(v, n)) * penalty_int - inner(avg(grad(u)), jump(v, n)) - inner(avg(grad(v)), jump(u, n))) * dIF # interior facet term + aEF = (inner(u, v) * penalty_ext - inner(grad(u), v*n) - inner(grad(v), u*n)) * dEF # exterior facet term + a = aV + aEF + aIF + L = f*v*dx + + # Compute solution + sol = Function(V) + + solve(a == L, sol, solver_parameters={'pc_type': 'ilu', + 'ksp_type': 'lgmres'}) + + # Analytical solution + if meshtype in ['planar', 'extruded']: + f.project(sin(x[0]*pi*2)*sin(x[1]*pi*2)) + elif meshtype == 'spherical': + f.project(x[0]*x[1]*x[2]/13.0) + return sqrt(assemble(dot(sol - f, sol - f) * dx)) + + +@pytest.mark.parametrize(('degree', 'meshd', 'meshtype', 'threshold'), + [(1, 1, 'planar', 1.8), + (2, 1, 'planar', 2.9), + (1, 1, 'extruded', 1.8), + (2, 1, 'extruded', 2.9), + (2, 2, 'spherical', 2.9), + (3, 3, 'spherical', 3.9)]) +def test_firedrake_laplacian_IP_mse(degree, meshd, meshtype, threshold): + diff = np.array([laplacian_IP(i, degree, meshd, meshtype) for i in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +def vector_laplace(r, meshtype, hdegree, vdegree=None): + vdegree = vdegree or hdegree + + mesh = get_mesh(r, meshtype) + h1elem, hcurlelem, _, _ = get_complex(meshtype, hdegree, vdegree=vdegree) + h1elem_high, _, _, _ = get_complex(meshtype, hdegree+1, vdegree=vdegree+1) + + # spaces for calculation + V0 = FunctionSpace(mesh, h1elem) + V1 = FunctionSpace(mesh, hcurlelem) + V = V0*V1 + + # spaces to store 'analytic' functions + W0 = FunctionSpace(mesh, h1elem_high) + W1 = VectorFunctionSpace(mesh, h1elem_high) + + # constants + k = 1.0 + l = 2.0 + + xs = SpatialCoordinate(mesh) + f_expr = as_vector([pi*pi*(k*k + l*l)*sin(k*pi*xs[0])*cos(l*pi*xs[1]), pi*pi*(k*k + l*l)*cos(k*pi*xs[0])*sin(l*pi*xs[1])]) + exact_s_expr = -(k+l)*pi*cos(k*pi*xs[0])*cos(l*pi*xs[1]) + exact_u_expr = as_vector([sin(k*pi*xs[0])*cos(l*pi*xs[1]), cos(k*pi*xs[0])*sin(l*pi*xs[1])]) + + f = Function(W1).project(f_expr) + exact_s = Function(W0).project(exact_s_expr) + exact_u = Function(W1).project(exact_u_expr) + + sigma, u = TrialFunctions(V) + tau, v = TestFunctions(V) + a = (sigma*tau - dot(u, grad(tau)) + dot(grad(sigma), v) + dot(curl(u), curl(v)))*dx + L = dot(f, v)*dx + + out = Function(V) + + # preconditioner for H1 x H(curl) + aP = (dot(grad(sigma), grad(tau)) + sigma*tau + dot(curl(u), curl(v)) + dot(u, v))*dx + + solve(a == L, out, Jp=aP, + solver_parameters={'pc_type': 'fieldsplit', + 'pc_fieldsplit_type': 'additive', + 'fieldsplit_0_pc_type': 'lu', + 'fieldsplit_1_pc_type': 'lu', + 'ksp_monitor': None}) + + out_s, out_u = out.split() + + return (sqrt(assemble(dot(out_u - exact_u, out_u - exact_u)*dx)), + sqrt(assemble((out_s - exact_s)*(out_s - exact_s)*dx))) + + +@pytest.mark.parametrize(('degree', 'threshold'), + [(1, 0.9), + (2, 1.9), + (3, 2.9), + (2, 0.9)]) +def test_firedrake_helmholtz_vector_mse(degree, threshold): + diff = np.array([vector_laplace(r, 'planar', degree) for r in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1] / diff[1:]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all() + + +@pytest.mark.parametrize(('hdegree', 'vdegree', 'threshold'), + [(1, 1, 0.9), + (2, 1, 0.9), + (1, 2, 0.9), + (2, 2, 1.9), + (1, 1, 0.9)]) +def test_firedrake_helmholtz_vector_mse_extruded(hdegree, vdegree, threshold): + diff = np.array([vector_laplace(i, 'extruded', hdegree, vdegree=vdegree) for i in range(2, 5)]) + print("l2 error norms:", diff) + conv = np.log2(diff[:-1, :] / diff[1:, :]) + print("convergence order:", conv) + assert (np.array(conv) > threshold).all()