diff --git a/bfps/cpp/full_code/NSVE.cpp b/bfps/cpp/full_code/NSVE.cpp
index 9c40968a66f050c318222690a94b273488e46b61..1e24c7af531e7184f75b1f14257d42b822db7a9c 100644
--- a/bfps/cpp/full_code/NSVE.cpp
+++ b/bfps/cpp/full_code/NSVE.cpp
@@ -88,6 +88,16 @@ int NSVE<rnumber>::finalize(void)
     return EXIT_SUCCESS;
 }
 
+/** \brief Compute standard statistics for velocity and vorticity fields.
+ *
+ *  IMPORTANT: at the end of this subroutine, `this->fs->cvelocity` contains
+ *  the Fourier space representation of the velocity field, and
+ *  `this->tmp_vec_field` contains the real space representation of the
+ *  velocity field.
+ *  This behavior is relied upon in the `NSVEparticles` class, so please
+ *  don't break it.
+ */
+
 template <typename rnumber>
 int NSVE<rnumber>::do_stats()
 {
@@ -101,22 +111,23 @@ int NSVE<rnumber>::do_stats()
                 H5P_DEFAULT);
     else
         stat_group = 0;
-    fs->compute_velocity(fs->cvorticity);
-    *tmp_vec_field = fs->cvelocity->get_cdata();
+
+    *tmp_vec_field = fs->cvorticity->get_cdata();
     tmp_vec_field->compute_stats(
             fs->kk,
             stat_group,
-            "velocity",
+            "vorticity",
             fs->iteration / niter_stat,
-            max_velocity_estimate/sqrt(3));
+            max_vorticity_estimate/sqrt(3));
 
-    *tmp_vec_field = fs->cvorticity->get_cdata();
+    fs->compute_velocity(fs->cvorticity);
+    *tmp_vec_field = fs->cvelocity->get_cdata();
     tmp_vec_field->compute_stats(
             fs->kk,
             stat_group,
-            "vorticity",
+            "velocity",
             fs->iteration / niter_stat,
-            max_vorticity_estimate/sqrt(3));
+            max_velocity_estimate/sqrt(3));
 
     if (this->myrank == 0)
         H5Gclose(stat_group);