add vertical_driving_only param

inputs/precipitator_lowres.in

@@ -97,7 +97,7 @@ thermostat_Ki = 10.0      # time constant (inverse cooling times)
 numHist = 64
 h_smooth_heatcool_kpc = 5.0
 
-perturb_sin_drho_over_rho = 0.01
+perturb_sin_drho_over_rho = 0.0
 perturb_kx = 16
 perturb_ky = 16
 perturb_kz = 32
@@ -107,8 +107,8 @@ perturb_exponent = 1    # P(k) ~ k^{-2}
 #sigma_v = 0
 sigma_v = 0.0001   # kpc Myr^{-1}
 k_peak = 2
-num_modes = 40
+num_modes = 12
+xy_modes_only = true
 sol_weight = 1.0
-vertical_weight = 1.0
 rseed = 42
 t_corr = 500.0   # Myr

src/pgen/precipitator.cpp

@@ -287,6 +287,7 @@ void TurbSrcTerm(MeshData<Real> *md, const parthenon::SimTime /*time*/, const Re
   IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
 
   const auto sigma_v = hydro_pkg->Param<Real>("sigma_v");
+  const auto vertical_driving_only = hydro_pkg->Param<bool>("xy_modes_only");
 
   const Real h_smooth = hydro_pkg->Param<Real>("h_smooth_heatcool");
 
@@ -303,9 +304,14 @@ void TurbSrcTerm(MeshData<Real> *md, const parthenon::SimTime /*time*/, const Re
         "normalize_perturb_v", 0, md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
         KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lsum) {
           const auto &coords = cons.GetCoords(b);
-          const Real dv_x = perturb_pack(b, 0, k, j, i);
-          const Real dv_y = perturb_pack(b, 1, k, j, i);
+          Real dv_x = 0;
+          Real dv_y = 0;
+          if (!vertical_driving_only) {
+            dv_x = perturb_pack(b, 0, k, j, i);
+            dv_y = perturb_pack(b, 1, k, j, i);
+          }
           const Real dv_z = perturb_pack(b, 2, k, j, i);
+
           lsum += (SQR(dv_x) + SQR(dv_y) + SQR(dv_z)) * coords.CellVolume(k, j, i);
         },
         v2_sum);
@@ -322,8 +328,12 @@ void TurbSrcTerm(MeshData<Real> *md, const parthenon::SimTime /*time*/, const Re
         "apply_perturb_v", 0, md->NumBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
         KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
           // compute delta_v
-          const Real dv_x = perturb_pack(b, 0, k, j, i) / v_norm;
-          const Real dv_y = perturb_pack(b, 1, k, j, i) / v_norm;
+          Real dv_x = 0;
+          Real dv_y = 0;
+          if (!vertical_driving_only) {
+            dv_x = perturb_pack(b, 0, k, j, i);
+            dv_y = perturb_pack(b, 1, k, j, i);
+          }
           const Real dv_z = perturb_pack(b, 2, k, j, i) / v_norm;
 
           // compute old kinetic energy
@@ -764,6 +774,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
 
   if (sigma_v > 0) {
     auto k_peak_v = pin->GetReal("precipitator/driving", "k_peak");
+    // NOTE: in 2D, there are only 12 modes.
     auto num_modes_v = pin->GetOrAddInteger("precipitator/driving", "num_modes", 40);
     auto sol_weight_v = pin->GetOrAddReal("precipitator/driving", "sol_weight", 1.0);
     uint32_t rseed_v = pin->GetOrAddInteger("precipitator/driving", "rseed", 1);
@@ -774,8 +785,13 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
                        std::vector<int>({3}));
     hydro_pkg->AddField("tmp_perturb", m_perturb);
 
-    auto k_vec_v = utils::few_modes_ft::MakeRandomModes(num_modes_v, k_peak_v, rseed_v);
-
+    auto xy_modes_only =
+        pin->GetOrAddBoolean("precipitator/driving", "xy_modes_only", false);
+    hydro_pkg->AddParam("xy_modes_only", xy_modes_only);
+
+    auto k_vec_v = utils::few_modes_ft::MakeRandomModes(num_modes_v, k_peak_v,
+                                                        xy_modes_only, rseed_v);
+
     auto few_modes_ft = FewModesFT(pin, hydro_pkg, "precipitator_perturb_v", num_modes_v,
                                    k_vec_v, k_peak_v, sol_weight_v, t_corr, rseed_v);
     hydro_pkg->AddParam<>("precipitator/few_modes_ft_v", few_modes_ft);

src/utils/few_modes_ft.cpp

@@ -9,7 +9,6 @@
 
 // C++ headers
 #include <random>
-#include <utility>
 
 // Parthenon headers
 #include "basic_types.hpp"
@@ -21,7 +20,6 @@
 #include "mesh/meshblock_pack.hpp"
 
 // AthenaPK headers
-#include "../main.hpp"
 #include "few_modes_ft.hpp"
 #include "utils/error_checking.hpp"
 
@@ -51,9 +49,9 @@ 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(std::abs(k_vec_host(0, i)) <= static_cast<Real>(gnx1) / 2, "k_vec x1 mode too large");
+    PARTHENON_REQUIRE(std::abs(k_vec_host(1, i)) <= static_cast<Real>(gnx2) / 2, "k_vec x2 mode too large");
+    PARTHENON_REQUIRE(std::abs(k_vec_host(2, i)) <= static_cast<Real>(gnx3) / 2, "k_vec x3 mode too large");
   }
 
   const auto nx1 = pin->GetInteger("parthenon/meshblock", "nx1");
@@ -482,7 +480,7 @@ void FewModesFT::Generate(MeshData<Real> *md, const Real dt,
 
 // Creates a random set of wave vectors with k_mag within k_peak/2 and 2*k_peak
 ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak,
-                                 uint32_t rseed = 31224) {
+                                 uint32_t rseed = 31224, const bool xy_modes_only) {
   auto k_vec = parthenon::ParArray2D<Real>("k_vec", 3, num_modes);
   auto k_vec_h = Kokkos::create_mirror_view_and_copy(parthenon::HostMemSpace(), k_vec);
 
@@ -495,15 +493,22 @@ ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak,
 
   int n_mode = 0;
   int n_attempt = 0;
-  constexpr int max_attempts = 1000000;
-  Real kx1, kx2, kx3, k_mag, ampl;
+  constexpr int max_attempts = 1e6;
+  Real kx1 = NAN;
+  Real kx2 = NAN;
+  Real kx3 = 0.;
+  Real k_mag = NAN;
+  Real ampl = NAN;
+
   bool mode_exists = false;
   while (n_mode < num_modes && n_attempt < max_attempts) {
     n_attempt += 1;
 
     kx1 = dist(rng);
     kx2 = dist(rng);
-    kx3 = dist(rng);
+    if (!xy_modes_only) {
+      kx3 = dist(rng);
+    }
     k_mag = std::sqrt(SQR(kx1) + SQR(kx2) + SQR(kx3));
 
     // Expected amplitude of the spectral function. If this is changed, it also needs to
@@ -528,6 +533,17 @@ ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak,
     k_vec_h(2, n_mode) = kx3;
     n_mode++;
   }
+
+#if 0
+  // print modes
+  for (int n = 0; n < n_mode; ++n) {
+    std::cout << "k_vec_h(0, " << n << ") = " << k_vec_h(0, n) << "\n";
+    std::cout << "k_vec_h(1, " << n << ") = " << k_vec_h(1, n) << "\n";
+    std::cout << "k_vec_h(2, " << n << ") = " << k_vec_h(2, n) << "\n";
+    std::cout << "\n";
+  }
+#endif
+
   PARTHENON_REQUIRE_THROWS(
       n_attempt < max_attempts,
       "Cluster init did not succeed in calculating perturbation modes.")

src/utils/few_modes_ft.hpp

@@ -67,6 +67,6 @@ class FewModesFT {
 };
 
 // Creates a random set of wave vectors with k_mag within k_peak/2 and 2*k_peak
-ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak, uint32_t rseed);
+ParArray2D<Real> MakeRandomModes(const int num_modes, const Real k_peak, uint32_t rseed, const bool xy_modes_only = false);
 
 } // namespace utils::few_modes_ft