tests: test for gauge invariance

tests/test_basic.py

@@ -3,6 +3,7 @@ import logging
 import pytest
 import numpy as np
 from numpy.testing import assert_allclose
+from scipy.linalg import expm
 
 import usadelndsoc
 import usadelndsoc.solver
@@ -77,3 +78,67 @@ def test_result_S_J0(omega, h, terminals):
     # There are no gradients, hence no currents either
     J = res.J
     assert_allclose(J, 0, atol=1e-5)
+
+
+def test_gauge_invariance():
+    nx = 10
+    ny = 7
+    h = 0.4
+
+    np.random.seed(1)
+
+    s = usadelndsoc.solver.Solver(nx=nx, ny=ny)
+    s.Omega[...] = 0
+    s.U[:, :, 0, :2, :2] = np.random.randn(nx, ny, 2, 2) + 1j * np.random.randn(
+        nx, ny, 2, 2
+    )
+    s.U[:, :, 1, :2, :2] = np.random.randn(nx, ny, 2, 2) + 1j * np.random.randn(
+        nx, ny, 2, 2
+    )
+    s.U[:, :, 0, 2:, 2:] = np.random.randn(nx, ny, 2, 2) + 1j * np.random.randn(
+        nx, ny, 2, 2
+    )
+    s.U[:, :, 1, 2:, 2:] = np.random.randn(nx, ny, 2, 2) + 1j * np.random.randn(
+        nx, ny, 2, 2
+    )
+    s.omega = 0.987
+    s.eta = 0.789
+    s.D = 1.0
+    s.Lx = 10
+    s.Ly = 10
+
+    s._core.init_U()
+
+    Phi = np.random.randn(nx, ny, 2, 2, 2) + 1j * np.random.randn(nx, ny, 2, 2, 2)
+    S0 = s._core.eval(Phi)
+
+    s.reset()
+
+    # Gauge transform
+    i = 3
+    j = 4
+
+    W = np.zeros((4, 4), dtype=complex)
+    W[:2, :2] = np.random.randn(2, 2)
+    W[2:, 2:] = np.random.randn(2, 2)
+    iW = np.linalg.inv(W)
+
+    s.U[i, j, 0] = s.U[i, j, 0] @ iW
+    s.U[i, j, 1] = s.U[i, j, 1] @ iW
+    s.U[i, j, 2] = s.U[i, j, 2] @ iW
+    s.U[i, j, 3] = s.U[i, j, 3] @ iW
+
+    s.U[i + 1, j, 2] = W @ s.U[i + 1, j, 2]
+    s.U[i - 1, j, 0] = W @ s.U[i - 1, j, 0]
+
+    s.U[i, j + 1, 3] = W @ s.U[i, j + 1, 3]
+    s.U[i, j - 1, 1] = W @ s.U[i, j - 1, 1]
+
+    Phi2 = Phi.copy()
+    Phi2[i, j, 0] = W[:2, :2] @ Phi[i, j, 0] @ iW[2:, 2:]
+    Phi2[i, j, 1] = W[2:, 2:] @ Phi[i, j, 1] @ iW[:2, :2]
+
+    S1 = s._core.eval(Phi2)
+
+    # Should be invariant to numerical accuracy
+    assert_allclose(S1, S0, atol=1e-9)