Hypergeometric2F1

.github/workflows/core.yml

@@ -214,6 +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 
           firedrake-clean
           pip list
 

pyop2/compilation.py

@@ -344,7 +344,7 @@ class LinuxGnuCompiler(Compiler):
 
     _cflags = ("-fPIC", "-Wall", "-std=gnu11")
     _cxxflags = ("-fPIC", "-Wall")
-    _ldflags = ("-shared",)
+    _ldflags = ("-shared", "-lgsl", "-lgslcblas")
 
     _optflags = ("-march=native", "-O3", "-ffast-math")
     _debugflags = ("-O0", "-g")

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():
+    sp_hyp2f1 = scipy.special.hyp2f1
+    a = 1
+    b = 1
+    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,
+                       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))
+
+    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))

tsfc/loopy.py

@@ -480,7 +480,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

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)