tests: higher accuracy gauge invariance test

pyproject.toml

@@ -29,7 +29,8 @@ dependencies = [
 
 [project.optional-dependencies]
 test = [
-  'pytest >= 6'
+  'pytest >= 6',
+  'numdifftools >= 0.9',
 ]
 
 [tool.pytest.ini_options]

tests/test_basic.py

@@ -5,6 +5,11 @@ import numpy as np
 from numpy.testing import assert_allclose
 from scipy.linalg import expm
 
+try:
+    import numdifftools
+except ImportError:
+    numdifftools = None
+
 import usadelndsoc
 import usadelndsoc.solver
 
@@ -157,18 +162,26 @@ def test_gauge_invariance():
     div_J_y = J[1] - J[3]
     div_J = div_J_x + div_J_y
 
+    print("->", div_J)
+
     # Compute source term (numdiff)
-    Sa = s._core.eval(Phi)
+    def S_s(ds):
+        Phi2 = Phi.copy()
+        dW = expm(ds * t3)
+        idW = np.linalg.inv(dW)
+        Phi2[i, j, 0] = dW[:2, :2] @ Phi[i, j, 0] @ idW[2:, 2:]
+        Phi2[i, j, 1] = dW[2:, 2:] @ Phi[i, j, 1] @ idW[:2, :2]
+        return s._core.eval(Phi2)
 
-    Phi2 = Phi.copy()
     ds = 1e-7
-    dW = expm(ds * t3)
-    idW = np.linalg.inv(dW)
-    Phi2[i, j, 0] = dW[:2, :2] @ Phi[i, j, 0] @ idW[2:, 2:]
-    Phi2[i, j, 1] = dW[2:, 2:] @ Phi[i, j, 1] @ idW[:2, :2]
-
-    Sb = s._core.eval(Phi2)
-    dS_ds = (Sb - Sa) / ds / 16
+    dS_ds = (S_s(ds) - S_s(0)) / ds / 16
 
     # They must match
     assert_allclose(div_J, dS_ds, rtol=1e-6)
+
+    if numdifftools is not None:
+        # More accurate numdiff: the equality should hold to high
+        # numerical accuracy
+        D = numdifftools.Derivative(S_s)
+        dS_ds = D(0.0) / 16
+        assert_allclose(div_J, dS_ds, rtol=1e-12)