diff --git a/lib/gpt/algorithms/base.py b/lib/gpt/algorithms/base.py
index 023ab5ac..2602525c 100644
--- a/lib/gpt/algorithms/base.py
+++ b/lib/gpt/algorithms/base.py
@@ -16,7 +16,8 @@
 #    with this program; if not, write to the Free Software Foundation, Inc.,
 #    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #
-import gpt, sys
+import gpt, sys, os
+from datetime import datetime
 
 
 class base:
@@ -69,15 +70,34 @@ class base_iterative(base):
     def __init__(self, name=None):
         super().__init__(name)
         self.verbose_convergence = gpt.default.is_verbose(self.name + "_convergence")
+        self.verbose_log_convergence = gpt.default.is_verbose(self.name + "_log_convergence")
         self.converged = None
 
+    def get_log_file(self):
+        if not self.verbose_log_convergence or gpt.rank() != 0:
+            return None
+
+        self.log_directory = "log/" + datetime.today().strftime("%Y-%m-%d")
+        if not os.path.exists(self.log_directory):
+            os.makedirs(self.log_directory, exist_ok=True)
+
+        log_filename = f"{self.log_directory}/{self.name}." + datetime.now().strftime("%H-%M-%f")
+        gpt.message(f"Convergence of {self.name} saved in {log_filename}")
+        return open(log_filename, "wt")
+
     def log_convergence(self, iteration, value, target=None):
         if (isinstance(iteration, int) and iteration == 0) or (
             isinstance(iteration, tuple) and all([x == 0 for x in iteration])
         ):
             self.history = []
+            self.log_file = self.get_log_file()
 
         self.history.append(value)
+        if self.log_file is not None:
+            if target is None:
+                self.log_file.write(f"{iteration} {value}\n")
+            else:
+                self.log_file.write(f"{iteration} {value:e} {target:e}\n")
 
         if target is not None:
             self.converged = value <= target
diff --git a/lib/gpt/core/group/__init__.py b/lib/gpt/core/group/__init__.py
index 730e84aa..c6d561dc 100644
--- a/lib/gpt/core/group/__init__.py
+++ b/lib/gpt/core/group/__init__.py
@@ -23,6 +23,7 @@
     compose,
     projected_convert,
     inverse,
+    invariant_distance,
 )
 from gpt.core.group.differentiable_functional import differentiable_functional
 from gpt.core.group.diffeomorphism import diffeomorphism, local_diffeomorphism
diff --git a/lib/gpt/core/group/operation.py b/lib/gpt/core/group/operation.py
index f53add57..660da66a 100644
--- a/lib/gpt/core/group/operation.py
+++ b/lib/gpt/core/group/operation.py
@@ -23,6 +23,12 @@ def defect(field):
     return field.otype.defect(field)
 
 
+def invariant_distance(field, field_prime):
+    field = g(field)
+    field_prime = g(field_prime)
+    return field.otype.invariant_distance(field, field_prime)
+
+
 def cartesian(field):
     if isinstance(field, list):
         return [cartesian(f) for f in field]
diff --git a/lib/gpt/core/object_type/su_n.py b/lib/gpt/core/object_type/su_n.py
index 0005c513..34e47f02 100644
--- a/lib/gpt/core/object_type/su_n.py
+++ b/lib/gpt/core/object_type/su_n.py
@@ -88,6 +88,11 @@ def __init__(self, Nc, Ndim, name):
     def cartesian(self):
         return self
 
+    def invariant_distance(self, a, b):
+        c = gpt(a - b)
+        ds2 = 2.0 * gpt.sum(gpt.trace(c * c)).real
+        return ds2**0.5
+
     def defect(self, A):
         err2 = gpt.norm2(gpt.adj(A) - A)
         err2 += gpt.norm2(gpt.trace(A))
@@ -129,6 +134,13 @@ def __init__(self, Nc, Ndim, name):
     def compose(self, a, b):
         return a * b
 
+    def invariant_distance(self, a, b):
+        d = gpt(a - b)
+        v = gpt(a + b)
+        c = self.infinitesimal_to_cartesian(v, d)
+        ds2 = 2.0 * gpt.sum(gpt.trace(c * c)).real
+        return ds2**0.5
+
     def defect(self, U):
         I = gpt.identity(U)
         I_s = gpt.identity(gpt.complex(U.grid))
diff --git a/lib/gpt/qcd/gauge/smear/__init__.py b/lib/gpt/qcd/gauge/smear/__init__.py
index 50fad848..acf0af6c 100644
--- a/lib/gpt/qcd/gauge/smear/__init__.py
+++ b/lib/gpt/qcd/gauge/smear/__init__.py
@@ -18,5 +18,5 @@
 #
 from gpt.qcd.gauge.smear.stout import stout, differentiable_stout
 from gpt.qcd.gauge.smear.local_stout import local_stout
-from gpt.qcd.gauge.smear.wilson_flow import wilson_flow
+from gpt.qcd.gauge.smear.wilson_flow import wilson_flow, gradient_flow
 from gpt.qcd.gauge.smear.differentiable import differentiable_field_transformation
diff --git a/lib/gpt/qcd/gauge/smear/differentiable.py b/lib/gpt/qcd/gauge/smear/differentiable.py
index b236a7fa..c6de2ef5 100644
--- a/lib/gpt/qcd/gauge/smear/differentiable.py
+++ b/lib/gpt/qcd/gauge/smear/differentiable.py
@@ -18,16 +18,24 @@
 #
 import gpt as g
 from gpt.core.group import diffeomorphism, differentiable_functional
+from gpt.ad import reverse as rad
+
+
+def assert_compatible(a, b, tag=""):
+    if type(a) is not type(b):
+        raise Exception(f"Incompatible types: {type(a)} and {type(b)}{tag}")
+    if isinstance(a, rad.node_base):
+        assert_compatible(a.value, b.value, " root " + str(type(a)))
 
 
 class dft_diffeomorphism(diffeomorphism):
     def __init__(self, U, ft):
-        rad = g.ad.reverse
         self.ft = ft
-        self.aU = [rad.node(u.new()) for u in U]
+        self.aU = [rad.node(u) for u in U]
         self.aUft = ft(self.aU)
 
     def __call__(self, fields):
+        # ft needs to be callable with a node or a lattice
         res = self.ft(fields)
         return [g(x) for x in res]
 
@@ -37,6 +45,7 @@ def jacobian(self, fields, fields_prime, dfields):
         assert len(dfields) == N
         aU_prime = [g(2j * dfields[mu] * fields_prime[mu]) for mu in range(N)]
         for mu in range(N):
+            assert_compatible(self.aU[mu].value, fields[mu])
             self.aU[mu].value = fields[mu]
         gradient = [None] * N
         for mu in range(N):
@@ -69,19 +78,19 @@ def adjoint_jacobian(self, fields, dfields):
 
 
 class dft_action_log_det_jacobian(differentiable_functional):
-    def __init__(self, U, ft, dfm, inverter_force, inverter_action):
+    def __init__(self, U, ft, dfm, dfm_node, inverter_force, inverter_action):
         self.dfm = dfm
+        self.dfm_node = dfm_node
         self.inverter_force = inverter_force
         self.inverter_action = inverter_action
         self.N = len(U)
         mom = [g.group.cartesian(u) for u in U]
-        rad = g.ad.reverse
 
         _U = [rad.node(g.copy(u)) for u in U]
         _left = [rad.node(g.copy(u), with_gradient=False) for u in mom]
         _right = [rad.node(g.copy(u), with_gradient=False) for u in mom]
-        _Up = dfm(_U)
-        J_right = dfm.jacobian(_U, _Up, _right)
+        _Up = dfm_node(_U)
+        J_right = dfm_node.jacobian(_U, _Up, _right)
 
         act = None
         for mu in range(self.N):
@@ -167,6 +176,7 @@ def __init__(self, U, ft, inverter_force, inverter_action, optimizer):
         self.ft = ft
         self.U = U
         self.dfm = dft_diffeomorphism(self.U, self.ft)
+        self.dfm_node = dft_diffeomorphism([rad.node(u) for u in self.U], self.ft)
         self.inverter_force = inverter_force
         self.inverter_action = inverter_action
         self.optimizer = optimizer
@@ -175,7 +185,6 @@ def diffeomorphism(self):
         return self.dfm
 
     def inverse(self, Uft):
-        rad = g.ad.reverse
         aU = [rad.node(g.copy(u)) for u in Uft]
         aUft_target = [rad.node(u, with_gradient=False) for u in Uft]
         aUft = self.ft(aU)
@@ -186,5 +195,5 @@ def inverse(self, Uft):
 
     def action_log_det_jacobian(self):
         return dft_action_log_det_jacobian(
-            self.U, self.ft, self.dfm, self.inverter_force, self.inverter_action
+            self.U, self.ft, self.dfm, self.dfm_node, self.inverter_force, self.inverter_action
         )
diff --git a/lib/gpt/qcd/gauge/smear/stout.py b/lib/gpt/qcd/gauge/smear/stout.py
index 840be117..d9cf9e82 100644
--- a/lib/gpt/qcd/gauge/smear/stout.py
+++ b/lib/gpt/qcd/gauge/smear/stout.py
@@ -223,7 +223,7 @@ def __call__(self, aU):
 
         U_prime = []
         for mu in range(nd):
-            U_mu_prime = (
+            U_mu_prime = g(
                 g.matrix.exp(g.qcd.gauge.project.traceless_anti_hermitian(C[mu] * g.adj(aU[mu])))
                 * aU[mu]
             )
diff --git a/lib/gpt/qcd/gauge/smear/wilson_flow.py b/lib/gpt/qcd/gauge/smear/wilson_flow.py
index 47c5c231..10257715 100644
--- a/lib/gpt/qcd/gauge/smear/wilson_flow.py
+++ b/lib/gpt/qcd/gauge/smear/wilson_flow.py
@@ -19,8 +19,12 @@
 import gpt as g
 
 
-def wilson_flow(U, epsilon):
-    action = g.qcd.gauge.action.wilson(2.0 * U[0].otype.shape[0])
+def gradient_flow(U, epsilon, action):
     return g.algorithms.integrator.runge_kutta_4(
         U, lambda Up: [g(-u) for u in action.gradient(Up, Up)], epsilon
     )
+
+
+def wilson_flow(U, epsilon):
+    action = g.qcd.gauge.action.wilson(2.0 * U[0].otype.shape[0])
+    return gradient_flow(U, epsilon, action)
diff --git a/tests/qcd/gauge.py b/tests/qcd/gauge.py
index 17609dc4..9fa0dfca 100755
--- a/tests/qcd/gauge.py
+++ b/tests/qcd/gauge.py
@@ -20,6 +20,28 @@
 # reference plaquette
 P = g.qcd.gauge.plaquette(U)
 
+# invariant distance (test axioms)
+ds = g.group.invariant_distance(U[0], U[1])
+dsp = g.group.invariant_distance(g.group.compose(V, U[0]), g.group.compose(V, U[1]))
+assert abs(ds / dsp - 1) < 1e-12
+
+ds = g.group.invariant_distance(U[0], U[1])
+dsp = g.group.invariant_distance(g.group.compose(U[0], V), g.group.compose(U[1], V))
+assert abs(ds / dsp - 1) < 1e-12
+
+ds = g.group.invariant_distance(U[0], U[0])
+assert ds < 1e-12
+
+ds = g.group.invariant_distance(U[0], U[1])
+dsp = g.group.invariant_distance(U[1], U[0])
+assert abs(ds / dsp - 1) < 1e-12
+
+mom = [g.group.cartesian(u) for u in U]
+rng.element(mom)
+ds = g.group.invariant_distance(mom[0], mom[1])
+dsp = g.group.invariant_distance(g.group.compose(mom[2], mom[0]), g.group.compose(mom[2], mom[1]))
+assert abs(ds / dsp - 1) < 1e-12
+
 # import cgpt, sys
 
 # smr = []