diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 598e331d9dcaf32efcbb6c91ebf11a1f08717744..cbd881a71b04638cf21014e4a28b9c48801f029c 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -656,6 +656,16 @@ class NavierStokes(_fluid_particle_base):
                     self.statistics[k] = pp_file[k].value
             self.compute_time_averages()
         return None
+    def compute_Reynolds_stress_invariants(
+            self):
+        Rij = self.statistics['R_ij(t)']
+        Rij /= (2*self.statistics['energy(t)'][:, None, None])
+        Rij[:, 0, 0] -= 1./3
+        Rij[:, 1, 1] -= 1./3
+        Rij[:, 2, 2] -= 1./3
+        self.statistics['I2(t)'] = np.sqrt(np.einsum('...ij,...ij', Rij, Rij, optimize = True) / 6)
+        self.statistics['I3(t)'] = np.cbrt(np.einsum('...ij,...jk,...ki', Rij, Rij, Rij, optimize = True) / 6)
+        return None
     def compute_time_averages(self):
         """Compute easy stats.