tests: check also grad in gauge-invariance test

tests/test_basic.py

@@ -201,3 +201,12 @@ def test_gauge_invariance(omega, use_numdifftools):
         D = numdifftools.Derivative(S_s)
         dS_ds = D(0.0) / 16
         assert_allclose(div_J, dS_ds, rtol=1e-12)
+
+    # The derivative must also match the gradient
+    dS_dPhi = s._core.grad(Phi).reshape(Phi.shape)
+    dPhi = 0 * Phi
+    dPhi[i, j, 0] = 2 * Phi[i, j, 0]
+    dPhi[i, j, 1] = -2 * Phi[i, j, 1]
+    dS_ds2 = (dS_dPhi * dPhi).sum() / 16
+
+    assert_allclose(div_J, dS_ds2, rtol=1e-12)