From bddcba870d19acb0132c4e38d6b29a8f042808bb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 4 Mar 2026 11:18:43 +0000 Subject: [PATCH 1/9] Hypergeometric2F1 --- tsfc/loopy.py | 8 +++++++- tsfc/ufl2gem.py | 3 +++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/tsfc/loopy.py b/tsfc/loopy.py index 6826f0b672..32b34d2c5c 100644 --- a/tsfc/loopy.py +++ b/tsfc/loopy.py @@ -449,7 +449,13 @@ def _expression_power(expr, ctx): @_expression.register(gem.MathFunction) def _expression_mathfunction(expr, ctx): - if expr.name.startswith('cyl_bessel_'): + if expr.name == "hyp2f1": + assert isinstance(ctx.target, lp.target.c.CWithGNULibcTarget) + # Generate right functions calls to gnulibc hypergeometric function + name = "hyperg_2F1" + return p.Variable(name)(*(expression(c, ctx) for c in expr.children)) + + elif expr.name.startswith('cyl_bessel_'): # Bessel functions if is_complex(ctx.scalar_type): raise NotImplementedError("Bessel functions for complex numbers: " diff --git a/tsfc/ufl2gem.py b/tsfc/ufl2gem.py index e283fb1414..19e4967d17 100644 --- a/tsfc/ufl2gem.py +++ b/tsfc/ufl2gem.py @@ -81,6 +81,9 @@ def math_function(self, o, expr): def atan2(self, o, y, x): return MathFunction("atan2", y, x) + def hypergeometric2_f1(self, o, a, b, c, arg): + return MathFunction(o._name, a, b, c, arg) + def bessel_i(self, o, nu, arg): return MathFunction(o._name, nu, arg) From 2864ea3ea56c85b22c5c799a7455af3f46300ef8 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 4 Mar 2026 11:19:13 +0000 Subject: [PATCH 2/9] test --- .../test_hypergeometric_function.py | 31 +++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 tests/firedrake/regression/test_hypergeometric_function.py diff --git a/tests/firedrake/regression/test_hypergeometric_function.py b/tests/firedrake/regression/test_hypergeometric_function.py new file mode 100644 index 0000000000..0ee14b5e0a --- /dev/null +++ b/tests/firedrake/regression/test_hypergeometric_function.py @@ -0,0 +1,31 @@ +from firedrake import * +import numpy as np +import scipy +import pytest + + +@pytest.mark.skipif(utils.complex_mode, reason="Complex bessel functions are not implemented.") +def test_hypergeometric_function(): + + a = 1 + b = 2 + c = 1 + mesh = UnitDiskMesh(3) + V = FunctionSpace(mesh, "CG", 1) + + x, y = SpatialCoordinate(mesh) + + expr = x / 10 + uexact = hyp2f1(a, b, c, expr) + assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, + scipy.hypf1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + + expr = sqrt(x**2+y**2) / 10 + uexact = hyp2f1(a, b, c, expr) + assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, + scipy.hypf1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + + expr = sqrt(x*x + y*y) / 10 + uexact = hyp2f1(a, b, c, expr) + assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, + scipy.hypf1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) From 52678a3eb02656bf7d445b02fd651ac813c46ecb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 4 Mar 2026 11:21:12 +0000 Subject: [PATCH 3/9] DROP BEFORE MERGE --- .github/workflows/core.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index c0682518a4..0a41ae171e 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -214,6 +214,7 @@ jobs: --extra-index-url https://download.pytorch.org/whl/cpu \ "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]" + pip install -I git+https://github.com/firedrakeproject/ufl.git@firedrakeproject:pbrubeck/hypergeometric firedrake-clean pip list From 1b736b7e0f37b57711b1e3b00d662e7467ec769f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 4 Mar 2026 17:19:34 +0000 Subject: [PATCH 4/9] WIP --- pyop2/compilation.py | 2 +- .../regression/test_hypergeometric_function.py | 10 +++++----- tsfc/loopy.py | 13 ++++++------- 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/pyop2/compilation.py b/pyop2/compilation.py index b37f004d30..a271bae635 100644 --- a/pyop2/compilation.py +++ b/pyop2/compilation.py @@ -344,7 +344,7 @@ class LinuxGnuCompiler(Compiler): _cflags = ("-fPIC", "-Wall", "-std=gnu11") _cxxflags = ("-fPIC", "-Wall") - _ldflags = ("-shared",) + _ldflags = ("-shared", "-lgsl") _optflags = ("-march=native", "-O3", "-ffast-math") _debugflags = ("-O0", "-g") diff --git a/tests/firedrake/regression/test_hypergeometric_function.py b/tests/firedrake/regression/test_hypergeometric_function.py index 0ee14b5e0a..f18cf597d5 100644 --- a/tests/firedrake/regression/test_hypergeometric_function.py +++ b/tests/firedrake/regression/test_hypergeometric_function.py @@ -6,9 +6,9 @@ @pytest.mark.skipif(utils.complex_mode, reason="Complex bessel functions are not implemented.") def test_hypergeometric_function(): - + sp_hyp2f1 = scipy.special.hyp2f1 a = 1 - b = 2 + b = 1 c = 1 mesh = UnitDiskMesh(3) V = FunctionSpace(mesh, "CG", 1) @@ -18,14 +18,14 @@ def test_hypergeometric_function(): expr = x / 10 uexact = hyp2f1(a, b, c, expr) assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, - scipy.hypf1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + sp_hyp2f1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) expr = sqrt(x**2+y**2) / 10 uexact = hyp2f1(a, b, c, expr) assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, - scipy.hypf1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + sp_hyp2f1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) expr = sqrt(x*x + y*y) / 10 uexact = hyp2f1(a, b, c, expr) assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, - scipy.hypf1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + sp_hyp2f1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) diff --git a/tsfc/loopy.py b/tsfc/loopy.py index 32b34d2c5c..5bb6e645a0 100644 --- a/tsfc/loopy.py +++ b/tsfc/loopy.py @@ -449,13 +449,8 @@ def _expression_power(expr, ctx): @_expression.register(gem.MathFunction) def _expression_mathfunction(expr, ctx): - if expr.name == "hyp2f1": - assert isinstance(ctx.target, lp.target.c.CWithGNULibcTarget) - # Generate right functions calls to gnulibc hypergeometric function - name = "hyperg_2F1" - return p.Variable(name)(*(expression(c, ctx) for c in expr.children)) - elif expr.name.startswith('cyl_bessel_'): + if expr.name.startswith('cyl_bessel_'): # Bessel functions if is_complex(ctx.scalar_type): raise NotImplementedError("Bessel functions for complex numbers: " @@ -486,7 +481,11 @@ def _expression_mathfunction(expr, ctx): else: return p.Variable(f"{name}n")(nu_, arg_) else: - if expr.name == "ln": + if expr.name == "hyp2f1": + assert isinstance(ctx.target, lp.target.c.CWithGNULibcTarget) + # Generate right functions calls to gsl hypergeometric function + name = "hyperg_2F1" + elif expr.name == "ln": name = "log" else: name = expr.name From c4fb8b59a932051c4b17a7616550c6d0dc3141dc Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 4 Mar 2026 17:34:48 +0000 Subject: [PATCH 5/9] DROP BEFORE MERGE --- .github/workflows/core.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index 0a41ae171e..69783c8d1a 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -215,6 +215,7 @@ jobs: "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]" pip install -I git+https://github.com/firedrakeproject/ufl.git@firedrakeproject:pbrubeck/hypergeometric + pip install -I git@github.com:pbrubeck/loopy.git@pbrubeck:pbrubeck/hypergeometric firedrake-clean pip list From cbd9b4b2b1b2355db7f8f17e74219d16020897a1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 5 Mar 2026 08:56:22 +0000 Subject: [PATCH 6/9] Add another compiler flag --- pyop2/compilation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyop2/compilation.py b/pyop2/compilation.py index a271bae635..a64ee6ce1a 100644 --- a/pyop2/compilation.py +++ b/pyop2/compilation.py @@ -344,7 +344,7 @@ class LinuxGnuCompiler(Compiler): _cflags = ("-fPIC", "-Wall", "-std=gnu11") _cxxflags = ("-fPIC", "-Wall") - _ldflags = ("-shared", "-lgsl") + _ldflags = ("-shared", "-lgsl", "-lgslcblas") _optflags = ("-march=native", "-O3", "-ffast-math") _debugflags = ("-O0", "-g") From 3176605d0bf5040e312ed522066ef5c85b382c42 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 5 Mar 2026 09:00:15 +0000 Subject: [PATCH 7/9] empty line --- tsfc/loopy.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tsfc/loopy.py b/tsfc/loopy.py index 5bb6e645a0..88c7e2403e 100644 --- a/tsfc/loopy.py +++ b/tsfc/loopy.py @@ -449,7 +449,6 @@ def _expression_power(expr, ctx): @_expression.register(gem.MathFunction) def _expression_mathfunction(expr, ctx): - if expr.name.startswith('cyl_bessel_'): # Bessel functions if is_complex(ctx.scalar_type): From 2f2ba9a3c66d459382fb1b7c9142a5ff2f6559b0 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 5 Mar 2026 09:06:51 +0000 Subject: [PATCH 8/9] DROP BEFORE MERGE --- .github/workflows/core.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/core.yml b/.github/workflows/core.yml index 69783c8d1a..b1cf6132b6 100644 --- a/.github/workflows/core.yml +++ b/.github/workflows/core.yml @@ -214,8 +214,8 @@ jobs: --extra-index-url https://download.pytorch.org/whl/cpu \ "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]" - pip install -I git+https://github.com/firedrakeproject/ufl.git@firedrakeproject:pbrubeck/hypergeometric - pip install -I git@github.com:pbrubeck/loopy.git@pbrubeck:pbrubeck/hypergeometric + pip install -I git+https://github.com/firedrakeproject/ufl.git@pbrubeck/hypergeometric + pip install -I git+https://github.com/pbrubeck/loopy.git@pbrubeck/hypergeometric firedrake-clean pip list From d69227a93e95524f42ff1c31229c0bd87370144f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 5 Mar 2026 10:35:35 +0000 Subject: [PATCH 9/9] fix libraries --- pyop2/compilation.py | 2 +- pyop2/global_kernel.py | 2 +- .../test_hypergeometric_function.py | 32 +++++++------------ 3 files changed, 14 insertions(+), 22 deletions(-) diff --git a/pyop2/compilation.py b/pyop2/compilation.py index a64ee6ce1a..b37f004d30 100644 --- a/pyop2/compilation.py +++ b/pyop2/compilation.py @@ -344,7 +344,7 @@ class LinuxGnuCompiler(Compiler): _cflags = ("-fPIC", "-Wall", "-std=gnu11") _cxxflags = ("-fPIC", "-Wall") - _ldflags = ("-shared", "-lgsl", "-lgslcblas") + _ldflags = ("-shared",) _optflags = ("-march=native", "-O3", "-ffast-math") _debugflags = ("-O0", "-g") diff --git a/pyop2/global_kernel.py b/pyop2/global_kernel.py index 4dd404b1f3..c9c294aab3 100644 --- a/pyop2/global_kernel.py +++ b/pyop2/global_kernel.py @@ -396,7 +396,7 @@ def _cppargs(self): def _ldargs(self): ldargs = [f"-L{d}/lib" for d in get_petsc_dir()] ldargs.extend(f"-Wl,-rpath,{d}/lib" for d in get_petsc_dir()) - ldargs.extend(["-lpetsc", "-lm"]) + ldargs.extend(["-lpetsc", "-lm", "-lgsl", "-lgslcblas"]) ldargs.extend(self.local_kernel.ldargs) return tuple(ldargs) diff --git a/tests/firedrake/regression/test_hypergeometric_function.py b/tests/firedrake/regression/test_hypergeometric_function.py index f18cf597d5..20c8a7e205 100644 --- a/tests/firedrake/regression/test_hypergeometric_function.py +++ b/tests/firedrake/regression/test_hypergeometric_function.py @@ -1,31 +1,23 @@ from firedrake import * +from scipy.special import hyp2f1 as sp_hyp2f1 import numpy as np -import scipy import pytest -@pytest.mark.skipif(utils.complex_mode, reason="Complex bessel functions are not implemented.") def test_hypergeometric_function(): - sp_hyp2f1 = scipy.special.hyp2f1 - a = 1 - b = 1 - c = 1 mesh = UnitDiskMesh(3) V = FunctionSpace(mesh, "CG", 1) - + u = Function(V) + w = Function(V) x, y = SpatialCoordinate(mesh) - expr = x / 10 - uexact = hyp2f1(a, b, c, expr) - assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, - sp_hyp2f1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) - - expr = sqrt(x**2+y**2) / 10 - uexact = hyp2f1(a, b, c, expr) - assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, - sp_hyp2f1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + expressions = [((1/3, 1/2, 1), Constant(0.9999) * x), + ((1/2, 1/2, 1), Constant(0.9999) * sqrt(x**2+y**2)), + ((-1/2, 1/2, 1), Constant(0.4999) * sqrt(x*x + y*y))] - expr = sqrt(x*x + y*y) / 10 - uexact = hyp2f1(a, b, c, expr) - assert np.allclose(assemble(Function(V).interpolate(uexact)).dat.data, - sp_hyp2f1(a, b, c, assemble(Function(V).interpolate(expr)).dat.data)) + for (a, b, c), expr in expressions: + u.interpolate(hyp2f1(a, b, c, expr)) + result = u.dat.data + w.interpolate(expr) + expect = sp_hyp2f1(a, b, c, w.dat.data) + assert np.allclose(result, expect)