Allow for specifying perturbation lengthscales in cluster init

docs/cluster.md

@@ -183,14 +183,16 @@ Pertubations are controlled by the following parameters:
 <problem/cluster/init_perturb>
 # for the velocity field
 sigma_v = 0.0      # volume weighted RMS |v| in code velocity; default: 0.0 (disabled)
-k_peak_v = ???     # wavenumber in normalized units where the velocity spectrum peaks. No default value.
+l_peak_v = ???     # lengthscale (in code length units) where the velocity spectrum peaks. No default value.
+k_peak_v = ???     # (alternative to l_peak_v): wavenumber in normalized units where the velocity spectrum peaks. No default value.
 num_modes_v = 40   # (optional) number of wavemodes in spectral space; default: 40
 sol_weight_v = 1.0 # (optional) power in solenoidal (rotational) modes of the perturbation. Range between 0 (fully compressive) and 1.0 (default, fully solenoidal).
 rseed_v = 1        # (optional) integer seed for RNG for wavenumbers and amplitudes
 
 # for the magnetic field
 sigma_b = 0.0      # volume weighted RMS |B| in code magnetic; default: 0.0 (disabled)
-k_peak_b = ???     # wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
+l_peak_b = ???     # lengthscale (in code length units) where the magnetic field spectrum peaks. No default value.
+k_peak_b = ???     # (alternative to l_peak_b): wavenumber in normalized units where the magnetic field spectrum peaks. No default value.
 num_modes_b = 40   # (optional) number of wavemodes in spectral space; default: 40
 rseed_b = 2        # (optional) integer seed for RNG for wavenumbers and amplitudes
 ```

src/pgen/cluster.cpp

@@ -421,7 +421,20 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
   const auto sigma_v = pin->GetOrAddReal("problem/cluster/init_perturb", "sigma_v", 0.0);
   if (sigma_v != 0.0) {
     // peak of init vel perturb
-    auto k_peak_v = pin->GetReal("problem/cluster/init_perturb", "k_peak_v");
+    auto l_peak_v = pin->GetOrAddReal("problem/cluster/init_perturb", "l_peak_v", -1.0);
+    auto k_peak_v = pin->GetOrAddReal("problem/cluster/init_perturb", "k_peak_v", -1.0);
+
+    PARTHENON_REQUIRE_THROWS((l_peak_v > 0.0 && k_peak_v <= 0.0) ||
+                                 (k_peak_v > 0.0 && l_peak_v <= 0.0),
+                             "Setting initial velocity perturbation requires a single "
+                             "length scale by either setting l_peak_v or k_peak_v.");
+    // Set peak wavemode as required by few_modes_fft when not directly given
+    if (l_peak_v > 0) {
+      const auto Lx = pin->GetReal("parthenon/mesh", "x1max") -
+                      pin->GetReal("parthenon/mesh", "x1min");
+      // Note that this assumes a cubic box
+      k_peak_v = Lx / l_peak_v;
+    }
     auto num_modes_v =
         pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_v", 40);
     auto sol_weight_v =
@@ -448,8 +461,20 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *hyd
     PARTHENON_REQUIRE_THROWS(hydro_pkg->Param<Fluid>("fluid") == Fluid::glmmhd,
                              "Requested initial magnetic field perturbations but not "
                              "solving the MHD equations.")
-    // peak of init vel perturb
-    auto k_peak_b = pin->GetReal("problem/cluster/init_perturb", "k_peak_b");
+    // peak of init magnetic field perturb
+    auto l_peak_b = pin->GetOrAddReal("problem/cluster/init_perturb", "l_peak_b", -1.0);
+    auto k_peak_b = pin->GetOrAddReal("problem/cluster/init_perturb", "k_peak_b", -1.0);
+    PARTHENON_REQUIRE_THROWS((l_peak_b > 0.0 && k_peak_b <= 0.0) ||
+                                 (k_peak_b > 0.0 && l_peak_b <= 0.0),
+                             "Setting initial B perturbation requires a single "
+                             "length scale by either setting l_peak_b or k_peak_b.");
+    // Set peak wavemode as required by few_modes_fft when not directly given
+    if (l_peak_b > 0) {
+      const auto Lx = pin->GetReal("parthenon/mesh", "x1max") -
+                      pin->GetReal("parthenon/mesh", "x1min");
+      // Note that this assumes a cubic box
+      k_peak_b = Lx / l_peak_b;
+    }
     auto num_modes_b =
         pin->GetOrAddInteger("problem/cluster/init_perturb", "num_modes_b", 40);
     uint32_t rseed_b = pin->GetOrAddInteger("problem/cluster/init_perturb", "rseed_b", 2);