From ad9fe2c762253b2317f5d24d68ae56e19b944b4e Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Wed, 15 Jan 2025 16:26:02 +0000
Subject: [PATCH] Test high order CR and variants (#3969)

---
 tests/firedrake/regression/test_helmholtz.py  |  2 +-
 .../test_helmholtz_crouzeix_raviart.py        | 59 +++++++++++++++++++
 .../regression/test_helmholtz_serendipity.py  |  7 +--
 3 files changed, 62 insertions(+), 6 deletions(-)
 create mode 100644 tests/firedrake/regression/test_helmholtz_crouzeix_raviart.py

diff --git a/tests/firedrake/regression/test_helmholtz.py b/tests/firedrake/regression/test_helmholtz.py
index 9097840275..6eeba6e71a 100644
--- a/tests/firedrake/regression/test_helmholtz.py
+++ b/tests/firedrake/regression/test_helmholtz.py
@@ -49,7 +49,7 @@ def helmholtz(r, quadrilateral=False, degree=2, mesh=None):
     assemble(L)
     sol = Function(V)
     solve(a == L, sol, solver_parameters={'ksp_type': 'cg'})
-    # Analytical solution
+    # Error norm
     return sqrt(assemble(inner(sol - expect, sol - expect) * dx)), sol, expect
 
 
diff --git a/tests/firedrake/regression/test_helmholtz_crouzeix_raviart.py b/tests/firedrake/regression/test_helmholtz_crouzeix_raviart.py
new file mode 100644
index 0000000000..9c8faef8a3
--- /dev/null
+++ b/tests/firedrake/regression/test_helmholtz_crouzeix_raviart.py
@@ -0,0 +1,59 @@
+"""This demo program solves Helmholtz's equation
+
+  - div grad u(x, y) + u(x,y) = f(x, y)
+
+on the unit square with source f given by
+
+  f(x, y) = (1.0 + 8.0*pi**2)*cos(x[0]*2*pi)*cos(x[1]*2*pi)
+
+and the analytical solution
+
+  u(x, y) = cos(x[0]*2*pi)*cos(x[1]*2*pi)
+"""
+
+import numpy as np
+import pytest
+
+from firedrake import *
+
+
+def helmholtz(r, quadrilateral=False, degree=1, variant=None, mesh=None):
+    # Create mesh and define function space
+    if mesh is None:
+        mesh = UnitSquareMesh(2 ** r, 2 ** r, quadrilateral=quadrilateral)
+    V = FunctionSpace(mesh, "CR", degree, variant=variant)
+
+    x, y = SpatialCoordinate(mesh)
+
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    uex = cos(x*pi*2)*cos(y*pi*2)
+    f = -div(grad(uex)) + uex
+
+    a = (inner(grad(u), grad(v)) + inner(u, v))*dx
+    L = inner(f, v)*dx(degree=12)
+
+    params = {"snes_type": "ksponly",
+              "ksp_type": "preonly",
+              "pc_type": "lu"}
+
+    # Compute solution
+    sol = Function(V)
+    solve(a == L, sol, solver_parameters=params)
+    # Error norm
+    return sqrt(assemble(dot(sol - uex, sol - uex) * dx)), sol, uex
+
+
+@pytest.mark.parametrize(('testcase', 'convrate'),
+                         [((1, (4, 6)), 1.9),
+                          ((3, (2, 4)), 3.9),
+                          ((5, (2, 4)), 5.7)])
+@pytest.mark.parametrize("variant", ("point", "integral"))
+def test_firedrake_helmholtz_scalar_convergence(variant, testcase, convrate):
+    degree, (start, end) = testcase
+    l2err = np.zeros(end - start)
+    for ii in [i + start for i in range(len(l2err))]:
+        l2err[ii - start] = helmholtz(ii, degree=degree, variant=variant)[0]
+    assert (np.array([np.log2(l2err[i]/l2err[i+1]) for i in range(len(l2err)-1)]) > convrate).all()
diff --git a/tests/firedrake/regression/test_helmholtz_serendipity.py b/tests/firedrake/regression/test_helmholtz_serendipity.py
index 77578d3c86..937b02899b 100644
--- a/tests/firedrake/regression/test_helmholtz_serendipity.py
+++ b/tests/firedrake/regression/test_helmholtz_serendipity.py
@@ -35,7 +35,7 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
     uex = cos(x*pi*2)*cos(y*pi*2)
     f = -div(grad(uex)) + uex
 
-    a = (inner(grad(u), grad(v)) + inner(u, v))*dx(degree=12)
+    a = (inner(grad(u), grad(v)) + inner(u, v))*dx
     L = inner(f, v)*dx(degree=12)
 
     params = {"snes_type": "ksponly",
@@ -45,10 +45,7 @@ def helmholtz(r, quadrilateral=True, degree=2, mesh=None):
     # Compute solution
     sol = Function(V)
     solve(a == L, sol, solver_parameters=params)
-
-    # Analytical solution
-    f = Function(V)
-    f.project(cos(x*pi*2)*cos(y*pi*2))
+    # Error norm
     return sqrt(assemble(dot(sol - uex, sol - uex) * dx)), sol, uex