Fixed data race condition in few modes IFT (no practical implication) (#172)

* Remove incorrect enforced symmetry with race condition and add explicit checks

* Remove complex conjugate modes from sample input

* Make REQUIRE a WARN for backwards compatibility

* Remove Hermitian modes from gen fmturb modes

* Add Changelog

* Update reference data given new sample vectors

* use ref data created on CI machine

* Remove phase factor from hermitian symmetric modes and exclude nyquist ones

* Python formatting

* Adjust tolerance

Author: pgrete

CHANGELOG.md

@@ -4,13 +4,12 @@
 
 ### General notes
 
-### Fixed (not changing behavior/API/variables/...)
-
 ### Added (new features/APIs/variables/...)
 
 ### Changed (changing behavior/API/variables/...)
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 172]](https://github.com/parthenon-hpc-lab/athenapk/pull/172) Fixed data race condition in few modes IFT (no practical implication)
 
 ### Infrastructure
 - [[PR 168]](https://github.com/parthenon-hpc-lab/athenapk/pull/168) Document agentic coding guidelines (and add PR template)
@@ -44,14 +43,14 @@ To enable, set `do_coalesced_comms=true` in the `<parthenon/mesh>` block of the
 - Input parameters can now be [automatically documented](https://github.com/parthenon-hpc-lab/parthenon/pull/1283)
 by adding an optional string as last argument to any `ParameterInput` `Get` or `GetOrAdd` call.
 
-### Fixed (not changing behavior/API/variables/...)
-- [[PR 163]](https://github.com/parthenon-hpc-lab/athenapk/pull/163) Add normalization by volume for relative B field divergence in history file
-
 ### Added (new features/APIs/variables/...)
 - [[PR 162]](https://github.com/parthenon-hpc-lab/athenapk/pull/162) Add pgen for cloud shattering setup
 - [[PR 158]](https://github.com/parthenon-hpc-lab/athenapk/pull/158) Update particle id handling (now automated `uint64`). Extend particle history lookback in turbulence pgen and include in turbulence test
 - [[PR 157]](https://github.com/parthenon-hpc-lab/athenapk/pull/157) Support injection of blobs with density/temp contrast in turbulence simulations
 
+### Changed (changing behavior/API/variables/...)
+- [[PR 163]](https://github.com/parthenon-hpc-lab/athenapk/pull/163) Add normalization by volume for relative B field divergence in history file
+
 ### Fixed (not changing behavior/API/variables/...)
 - [[PR 160]](https://github.com/parthenon-hpc-lab/athenapk/pull/160) Backport HLLD degeneracy check from Athena++
 

inputs/generate_fmturb_modes.py

@@ -26,16 +26,26 @@
 
 all_vec = []
 
-for i in range(k_high):
-    for j in range(-k_high, k_high):
-        for k in range(-k_high, k_high):
+for i in range(k_high + 1):
+    for j in range(-k_high, k_high + 1):
+        for k in range(-k_high, k_high + 1):
 
             k_mag = np.sqrt(i**2 + j**2 + k**2)
             if k_mag > k_high or k_mag < k_low:
                 continue
             # this is the spectral shape of the implemented forcing function
             if (k_mag / k_peak) ** 2.0 * (2.0 - (k_mag / k_peak) ** 2.0) < 0:
                 continue
+
+            # skip Hermitian partner
+            if i == 0:
+                skip = False
+                for vec in all_vec:
+                    if j == -vec[1] and k == -vec[2]:
+                        skip = True
+                if skip:
+                    continue
+
             all_vec.append((i, j, k))
 
 # select all wavevectors

inputs/turbulence.in

@@ -72,7 +72,7 @@ corr_time    = 1.0      # autocorrelation time of the OU forcing process
 rseed        = 20190729 # random seed of the OU forcing process
 sol_weight   = 1.0      # solenoidal weight of the acceleration field
 accel_rms    = 0.5      # root mean square value of the acceleration field
-num_modes    = 30       # number of wavemodes
+num_modes    = 20       # number of wavemodes
 
 <modes>
 k_1_0	= +2
@@ -93,9 +93,9 @@ k_5_2	= -1
 k_6_0	= +1
 k_6_1	= -1
 k_6_2	= -2
-k_7_0	= +0
-k_7_1	= +0
-k_7_2	= -2
+k_7_0	= +1
+k_7_1	= +2
+k_7_2	= -1
 k_8_0	= +1
 k_8_1	= +0
 k_8_2	= -1
@@ -117,9 +117,9 @@ k_13_2	= +1
 k_14_0	= +1
 k_14_1	= -1
 k_14_2	= +0
-k_15_0	= +0
-k_15_1	= -1
-k_15_2	= -1
+k_15_0	= +1
+k_15_1	= +1
+k_15_2	= +2
 k_16_0	= +1
 k_16_1	= +0
 k_16_2	= +1
@@ -135,34 +135,3 @@ k_19_2	= +1
 k_20_0	= +2
 k_20_1	= +1
 k_20_2	= -1
-k_21_0	= +0
-k_21_1	= -1
-k_21_2	= -2
-k_22_0	= +2
-k_22_1	= -1
-k_22_2	= +1
-k_23_0	= +0
-k_23_1	= +1
-k_23_2	= +1
-k_24_0	= +1
-k_24_1	= -2
-k_24_2	= +1
-k_25_0	= +1
-k_25_1	= -2
-k_25_2	= +0
-k_26_0	= +1
-k_26_1	= +2
-k_26_2	= +0
-k_27_0	= +1
-k_27_1	= -2
-k_27_2	= -1
-k_28_0	= +2
-k_28_1	= -1
-k_28_2	= -1
-k_29_0	= +1
-k_29_1	= +2
-k_29_2	= -1
-k_30_0	= +1
-k_30_1	= +1
-k_30_2	= +2
-

src/pgen/turbulence.cpp

@@ -19,6 +19,7 @@
 #include "globals.hpp"
 #include "interface/metadata.hpp"
 #include "kokkos_abstraction.hpp"
+#include "kokkos_types.hpp"
 #include "mesh/mesh.hpp"
 #include <iomanip>
 #include <ios>
@@ -359,7 +360,7 @@ void ProblemInitTracerData(ParameterInput *pin, parthenon::StateDescriptor *trac
     PARTHENON_THROW(
         "This line should not be reached as invalid n_lookback are caught above.");
   }
-  auto dncycles_d =
+  parthenon::ParArray1D<int> dncycles_d =
       Kokkos::create_mirror_view_and_copy(parthenon::DevMemSpace(), dncycles_h);
   tracer_pkg->AddParam("turbulence/n_lookback", n_lookback);
   tracer_pkg->AddParam("turbulence/dncycles_h", dncycles_h);

src/utils/few_modes_ft.cpp

@@ -9,6 +9,7 @@
 
 // C++ headers
 #include <random>
+#include <string>
 
 // Parthenon headers
 #include "basic_types.hpp"
@@ -47,9 +48,50 @@ FewModesFT::FewModesFT(parthenon::ParameterInput *pin, parthenon::StateDescripto
   // lambda cannot live in the constructor of an object.
   auto k_vec_host = k_vec.GetHostMirrorAndCopy();
   for (int i = 0; i < num_modes; i++) {
-    PARTHENON_REQUIRE(std::abs(k_vec_host(0, i)) <= gnx1 / 2, "k_vec x1 mode too large");
-    PARTHENON_REQUIRE(std::abs(k_vec_host(1, i)) <= gnx2 / 2, "k_vec x2 mode too large");
-    PARTHENON_REQUIRE(std::abs(k_vec_host(2, i)) <= gnx3 / 2, "k_vec x3 mode too large");
+    PARTHENON_REQUIRE(
+        k_vec_host(0, i) >= 0,
+        "k_x mode must be positive as explicit Hermitian symmetry is assumed.");
+    PARTHENON_REQUIRE(
+        !(k_vec_host(0, i) == 0 && k_vec_host(1, i) == 0 && k_vec_host(2, i) == 0),
+        "Forcing normalization is handled separately so do not include the 0 mode.");
+    PARTHENON_REQUIRE(std::abs(k_vec_host(0, i)) < gnx1 / 2, "k_vec x1 mode too large");
+    PARTHENON_REQUIRE(std::abs(k_vec_host(1, i)) < gnx2 / 2, "k_vec x2 mode too large");
+    PARTHENON_REQUIRE(std::abs(k_vec_host(2, i)) < gnx3 / 2, "k_vec x3 mode too large");
+  }
+
+  // Ensure that the provided k_vec do not include their conjugate and/or itself again as
+  // otherwise that mode would be counted twice
+  for (int i = 0; i < num_modes; i++) {
+    for (int j = 0; j < num_modes; j++) {
+      if (i == j) {
+        continue;
+      }
+
+      // The sample input file shipped for a couple of years
+      // unforunately came with complex conjugate modes...
+      // So in order to not crash old/existing sims upon restart we allow this and just
+      // issue a warning
+      if (k_vec_host(0, i) == 0 && k_vec_host(0, i) == k_vec_host(0, j) &&
+          k_vec_host(1, i) == -k_vec_host(1, j) &&
+          k_vec_host(2, i) == -k_vec_host(2, j)) {
+        PARTHENON_WARN(
+            "The given set of k_vec include complex conjugate partners for mode " +
+            std::to_string(i) + " with components " + std::to_string(k_vec_host(0, i)) +
+            " " + std::to_string(k_vec_host(1, i)) + " " +
+            std::to_string(k_vec_host(2, i)) +
+            "."
+            "In theory, this results in these modes being counted double. In practice, "
+            "this will have little to no effect as the normalization is done separate "
+            "for the real space field rather than the spectral field. However, if this "
+            "is a fresh simulation (and not a restarted one with given/fixed k_vec) it "
+            "is recommended to update the set of k_vec in the input file.");
+      }
+      PARTHENON_REQUIRE_THROWS(!(k_vec_host(0, i) == k_vec_host(0, j) &&
+                                 k_vec_host(1, i) == k_vec_host(1, j) &&
+                                 k_vec_host(2, i) == k_vec_host(2, j)),
+                               "The given set of k_vec include identical modes. "
+                               "Please provide a set without dublicates.");
+    }
   }
 
   const auto nx1 = pin->GetInteger("parthenon/meshblock", "nx1");
@@ -161,12 +203,7 @@ void FewModesFT::SetPhases(MeshBlock *pmb, ParameterInput *pin) {
 
         for (int m = 0; m < num_modes; m++) {
           w_kx = k_vec(0, m) * 2. * M_PI / static_cast<Real>(gnx1);
-          // adjust phase factor to Complex->Real IFT: u_hat*(k) = u_hat(-k)
-          if (k_vec(0, m) == 0.0) {
-            phase = 0.5 * Kokkos::exp(I * w_kx * gi);
-          } else {
-            phase = Kokkos::exp(I * w_kx * gi);
-          }
+          phase = Kokkos::exp(I * w_kx * gi);
           phases_i(i, m, 0) = phase.real();
           phases_i(i, m, 1) = phase.imag();
         }
@@ -255,19 +292,6 @@ void FewModesFT::Generate(MeshData<Real> *md, const Real dt,
             Complex(tmp * norm * random_num(n, m, 0), tmp * norm * random_num(n, m, 1));
       });
 
-  // enforce symmetry of complex to real transform
-  pmb->par_for(
-      "forcing: enforce symmetry", 0, 2, 0, num_modes - 1,
-      KOKKOS_LAMBDA(const int n, const int m) {
-        if (k_vec(0, m) == 0.) {
-          for (int m2 = 0; m2 < m; m2++) {
-            if (k_vec(1, m) == -k_vec(1, m2) && k_vec(2, m) == -k_vec(2, m2))
-              var_hat_new(n, m) =
-                  Complex(var_hat_new(n, m2).real(), -var_hat_new(n, m2).imag());
-          }
-        }
-      });
-
   const auto sol_weight = sol_weight_;
   if (sol_weight_ >= 0.0) {
     // project
@@ -372,6 +396,7 @@ ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak,
   constexpr int max_attempts = 1000000;
   Real kx1, kx2, kx3, k_mag, ampl;
   bool mode_exists = false;
+  bool hermitian_exists = false;
   while (n_mode < num_modes && n_attempt < max_attempts) {
     n_attempt += 1;
 
@@ -393,8 +418,20 @@ ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak,
       }
     }
 
+    // Check is Hermitian symmetric partner already exist
+    hermitian_exists = false;
+    if (kx1 == 0) {
+      for (int n_mode_exsist = 0; n_mode_exsist < n_mode; n_mode_exsist++) {
+        if (k_vec_h(0, n_mode_exsist) == kx1 && -k_vec_h(1, n_mode_exsist) == kx2 &&
+            -k_vec_h(2, n_mode_exsist) == kx3) {
+          hermitian_exists = true;
+        }
+      }
+    }
+
     // kx1 < 0.0 because we use a explicit symmetric Complex to Real transform
-    if (ampl < 0 || k_mag < k_low || k_mag > k_high || mode_exists || kx1 < 0.0) {
+    if (ampl < 0 || k_mag < k_low || k_mag > k_high || mode_exists || hermitian_exists ||
+        kx1 < 0.0) {
       continue;
     }
     k_vec_h(0, n_mode) = kx1;
@@ -409,4 +446,4 @@ ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak,
 
   return k_vec;
 }
-} // namespace utils::few_modes_ft
+} // namespace utils::few_modes_ft

tst/regression/test_suites/turbulence/ref_data.pkl

@@ -45,17 +45,37 @@ def Analyse(self, parameters):
         data = np.genfromtxt(data_filename)
 
         # Check Ms
-        if not (data[-1, -3] > 0.45 and data[-1, -3] < 0.50):
+        if data[-1, -3] != 4.70153e-01:
             print(f"ERROR: Mismatch in Ms={data[-1, -3]}")
             success = False
 
         # Check Ma
-        if not (data[-1, -2] > 12.8 and data[-1, -2] < 13.6):
+        if data[-1, -2] != 1.37253e01:
             print(f"ERROR: Mismatch in Ma={data[-1, -2]}")
             success = False
 
         # Loading the data
         data = phdf.phdf(f"{parameters.output_path}/parthenon.restart.final.rhdf")
+
+        components = data.GetComponents(data.Info["ComponentNames"], flatten=False)
+        density_sum = components["cons_density"].sum()
+        try:
+            np.testing.assert_array_max_ulp(density_sum, 64**3, maxulp=2)
+        except AssertionError as ar:
+            print(
+                f"TEST FAIL: incorrect density sum\n"
+                f"Got {density_sum} and expected {64**3}."
+            )
+            print(ar)
+            success = False
+        energy_sum = components["cons_total_energy_density"].sum()
+        if energy_sum != 2621560462.960244:
+            print(
+                f"TEST FAIL: incorrect energy sum\n"
+                f"Got {energy_sum} and expected {2621560462.960244}."
+            )
+            success = False
+
         tracers = data.GetSwarm("tracers")
         ids = tracers.id
 
@@ -76,13 +96,11 @@ def Analyse(self, parameters):
         # all_var_data = {}
         # for var in tracers.variables:
         #    var_data = tracers.Get(var)
-        #    all_var_data[var] = var_data[order]
         #    if len(var_data.shape) > 1:
         #        out_data = var_data[np.arange(var_data.shape[0])[:, None], order]
         #    else:
         #        out_data = var_data[order]
         #    all_var_data[var] = out_data
-
         # with open("ref_data.pkl", "wb") as outfile:
         #    pickle.dump(all_var_data, outfile)
 
@@ -119,7 +137,7 @@ def Analyse(self, parameters):
                         )
                     else:
                         np.testing.assert_allclose(
-                            var_data_sorted, ref_data[var], rtol=4e-8, strict=True
+                            var_data_sorted, ref_data[var], rtol=8.1e-7, strict=True
                         )
 
                 except AssertionError as ar: