Add initial pertubation for cluster pgen

README.md

@@ -2,9 +2,9 @@
 
 AthenaPK: a performance portable version based on [Athena++](https://github.com/PrincetonUniversity/athena),  [Parthenon](https://github.com/parthenon-hpc-lab/parthenon) and [Kokkos](https://github.com/kokkos/kokkos).
 
-## Current state of the code
+## Overview
 
-For this reason, it is highly recommended to only use AthenaPK with the Kokkos and Parthenon versions that are provided by the submodules (see [building](#building)) and to build everything (AthenaPK, Parthenon, and Kokkos) together from source.
+It is highly recommended to only use AthenaPK with the Kokkos and Parthenon versions that are provided by the submodules (see [building](#building)) and to build everything (AthenaPK, Parthenon, and Kokkos) together from source.
 Neither other versions or nor using preinstalled Parthenon/Kokkos libraries have been tested.
 
 Current features include
@@ -76,19 +76,22 @@ Obtain all (AthenaPK, Parthenon, and Kokkos) sources
 Most of the general build instructions and options for Parthenon (see [here](https://parthenon-hpc-lab.github.io/parthenon/develop/src/building.html)) also apply to AthenaPK.
 The following examples are a few standard cases.
 
-Most simple configuration (only CPU, no MPI, no HDF5).
+Most simple configuration (only CPU, no MPI).
 The `Kokkos_ARCH_...` parameter should be adjusted to match the target machine where AthenaPK will be executed.
 A full list of architecture keywords is available on the [Kokkos wiki](https://kokkos.github.io/kokkos-core-wiki/keywords.html#architecture-keywords).
 
-
-    # configure with enabling Broadwell architecture (AVX2) instructions
-    cmake -S. -Bbuild-host -DKokkos_ARCH_BDW=ON -DPARTHENON_DISABLE_MPI=ON -DPARTHENON_DISABLE_HDF5=ON
+    # configure with enabling Intel Broadwell or similar architecture (AVX2) instructions
+    cmake -S. -Bbuild-host -DKokkos_ARCH_BDW=ON -DPARTHENON_DISABLE_MPI=ON
     # now build with
     cd build-host && make
     # or alternatively
     cmake --build build-host
 
-An Intel Skylake system (AVX512 instructions) with NVidia Volta V100 GPUs and with MPI and HDF5 enabled (the latter is the default option, so they don't need to be specified)
+If `cmake` has troubling finding the HDF5 library (which is required for writing analysis outputs or
+restartings simulation) an additional hint to the location of the library can be provided via
+`-DHDF5_ROOT=/path/to/local/hdf5` on the first `cmake` command for configuration.
+
+An Intel Skylake system (AVX512 instructions) with NVidia Volta V100 GPUs and with MPI enabled (the latter is the default option, so they don't need to be specified)
 
     cmake -S. -Bbuild-gpu -DKokkos_ARCH_SKX=ON -DKokkos_ENABLE_CUDA=ON -DKokkos_ARCH_VOLTA70=ON
     # now build with

docs/cluster.md

@@ -162,6 +162,39 @@ r_sampling = 4.0
 Specifically, the resolution of the 1D profile for each meshblock is either
 `min(dx,dy,dz)/r_sampling` or `r_k/r_sampling`, whichever is smaller.
 
+## Initial perturbations
+
+Initial perturbations for both the velocity field and the magnetic field are
+supported.
+
+*Note* that these intial perturbation are currently incompatible with
+other initial conditions that modify the velocity or magnetic field, e.g.,
+an initial magnetic dipole or a uniform field.
+This restriction could be lifted if required/desired but the normalization
+of the fields would need to be adjusted.
+
+In general, the perturbations will be seeded by an inverse parabolic
+spectral profile centered on a peak wavenumber (in normalized units, i.e.,
+`k_peak = 2` mean half the box size) with a finite number of modes (default 40)
+randomly chosen between `k_peak/2` and `2*k_peak`.
+
+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.
+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.
+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
+```
+
 ## AGN Triggering
 
 If AGN triggering is enabled, at the end of each time step, a mass accretion

src/CMakeLists.txt

@@ -17,6 +17,7 @@ add_executable(
         hydro/srcterms/tabular_cooling.cpp
         refinement/gradient.cpp
         refinement/other.cpp
+        utils/few_modes_ft.cpp
 )
 
 add_subdirectory(pgen)

src/eos/adiabatic_glmmhd.hpp

@@ -24,8 +24,10 @@ using parthenon::Real;
 class AdiabaticGLMMHDEOS : public EquationOfState {
  public:
   AdiabaticGLMMHDEOS(Real pressure_floor, Real density_floor, Real internal_e_floor,
-                     Real gamma)
-      : EquationOfState(pressure_floor, density_floor, internal_e_floor), gamma_{gamma} {}
+                     Real velocity_ceiling, Real internal_e_ceiling, Real gamma)
+      : EquationOfState(pressure_floor, density_floor, internal_e_floor, velocity_ceiling,
+                        internal_e_ceiling),
+        gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
 
@@ -66,6 +68,9 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
     auto pressure_floor_ = GetPressureFloor();
     auto e_floor_ = GetInternalEFloor();
 
+    auto velocity_ceiling_ = GetVelocityCeiling();
+    auto e_ceiling_ = GetInternalECeiling();
+
     Real &u_d = cons(IDN, k, j, i);
     Real &u_m1 = cons(IM1, k, j, i);
     Real &u_m2 = cons(IM2, k, j, i);
@@ -109,6 +114,23 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
     Real e_B = 0.5 * (SQR(u_b1) + SQR(u_b2) + SQR(u_b3));
     w_p = gm1 * (u_e - e_k - e_B);
 
+    // apply velocity ceiling. By default ceiling is std::numeric_limits<Real>::infinity()
+    const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz);
+    if (w_v2 > SQR(velocity_ceiling_)) {
+      const Real w_v = sqrt(w_v2);
+      w_vx *= velocity_ceiling_ / w_v;
+      w_vy *= velocity_ceiling_ / w_v;
+      w_vz *= velocity_ceiling_ / w_v;
+
+      u_m1 *= velocity_ceiling_ / w_v;
+      u_m2 *= velocity_ceiling_ / w_v;
+      u_m3 *= velocity_ceiling_ / w_v;
+
+      Real e_k_new = 0.5 * u_d * SQR(velocity_ceiling_);
+      u_e -= e_k - e_k_new;
+      e_k = e_k_new;
+    }
+
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative pressure is encountered.
     PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
@@ -130,6 +152,14 @@ class AdiabaticGLMMHDEOS : public EquationOfState {
       w_p = eff_pressure_floor;
     }
 
+    // temperature (internal energy) based pressure ceiling
+    const Real eff_pressure_ceiling = gm1 * u_d * e_ceiling_;
+    if (w_p > eff_pressure_ceiling) {
+      // apply temperature ceiling, correct total energy
+      u_e = (u_d * e_ceiling_) + e_k + e_B;
+      w_p = eff_pressure_ceiling;
+    }
+
     // Convert passive scalars
     for (auto n = nhydro; n < nhydro + nscalars; ++n) {
       prim(n, k, j, i) = cons(n, k, j, i) * di;

src/eos/adiabatic_hydro.hpp

@@ -25,8 +25,10 @@ using parthenon::Real;
 class AdiabaticHydroEOS : public EquationOfState {
  public:
   AdiabaticHydroEOS(Real pressure_floor, Real density_floor, Real internal_e_floor,
-                    Real gamma)
-      : EquationOfState(pressure_floor, density_floor, internal_e_floor), gamma_{gamma} {}
+                    Real velocity_ceiling, Real internal_e_ceiling, Real gamma)
+      : EquationOfState(pressure_floor, density_floor, internal_e_floor, velocity_ceiling,
+                        internal_e_ceiling),
+        gamma_{gamma} {}
 
   void ConservedToPrimitive(MeshData<Real> *md) const override;
 
@@ -55,6 +57,9 @@ class AdiabaticHydroEOS : public EquationOfState {
     auto pressure_floor_ = GetPressureFloor();
     auto e_floor_ = GetInternalEFloor();
 
+    auto velocity_ceiling_ = GetVelocityCeiling();
+    auto e_ceiling_ = GetInternalECeiling();
+
     Real &u_d = cons(IDN, k, j, i);
     Real &u_m1 = cons(IM1, k, j, i);
     Real &u_m2 = cons(IM2, k, j, i);
@@ -84,6 +89,23 @@ class AdiabaticHydroEOS : public EquationOfState {
     Real e_k = 0.5 * di * (SQR(u_m1) + SQR(u_m2) + SQR(u_m3));
     w_p = gm1 * (u_e - e_k);
 
+    // apply velocity ceiling. By default ceiling is std::numeric_limits<Real>::infinity()
+    const Real w_v2 = SQR(w_vx) + SQR(w_vy) + SQR(w_vz);
+    if (w_v2 > SQR(velocity_ceiling_)) {
+      const Real w_v = sqrt(w_v2);
+      w_vx *= velocity_ceiling_ / w_v;
+      w_vy *= velocity_ceiling_ / w_v;
+      w_vz *= velocity_ceiling_ / w_v;
+
+      u_m1 *= velocity_ceiling_ / w_v;
+      u_m2 *= velocity_ceiling_ / w_v;
+      u_m3 *= velocity_ceiling_ / w_v;
+
+      Real e_k_new = 0.5 * u_d * SQR(velocity_ceiling_);
+      u_e -= e_k - e_k_new;
+      e_k = e_k_new;
+    }
+
     // Let's apply floors explicitly, i.e., by default floor will be disabled (<=0)
     // and the code will fail if a negative pressure is encountered.
     PARTHENON_REQUIRE(w_p > 0.0 || pressure_floor_ > 0.0 || e_floor_ > 0.0,
@@ -105,6 +127,14 @@ class AdiabaticHydroEOS : public EquationOfState {
       w_p = eff_pressure_floor;
     }
 
+    // temperature (internal energy) based pressure ceiling
+    const Real eff_pressure_ceiling = gm1 * u_d * e_ceiling_;
+    if (w_p > eff_pressure_ceiling) {
+      // apply temperature ceiling, correct total energy
+      u_e = (u_d * e_ceiling_) + e_k;
+      w_p = eff_pressure_ceiling;
+    }
+
     // Convert passive scalars
     for (auto n = nhydro; n < nhydro + nscalars; ++n) {
       prim(n, k, j, i) = cons(n, k, j, i) * di;

src/eos/eos.hpp

@@ -32,9 +32,11 @@ using parthenon::Real;
 
 class EquationOfState {
  public:
-  EquationOfState(Real pressure_floor, Real density_floor, Real internal_e_floor)
+  EquationOfState(Real pressure_floor, Real density_floor, Real internal_e_floor,
+                  Real velocity_ceiling, Real internal_e_ceiling)
       : pressure_floor_(pressure_floor), density_floor_(density_floor),
-        internal_e_floor_(internal_e_floor) {}
+        internal_e_floor_(internal_e_floor), velocity_ceiling_(velocity_ceiling),
+        internal_e_ceiling_(internal_e_ceiling) {}
   virtual void ConservedToPrimitive(MeshData<Real> *md) const = 0;
 
   KOKKOS_INLINE_FUNCTION
@@ -47,8 +49,15 @@ class EquationOfState {
   KOKKOS_INLINE_FUNCTION
   Real GetInternalEFloor() const { return internal_e_floor_; }
 
+  KOKKOS_INLINE_FUNCTION
+  Real GetVelocityCeiling() const { return velocity_ceiling_; }
+
+  KOKKOS_INLINE_FUNCTION
+  Real GetInternalECeiling() const { return internal_e_ceiling_; }
+
  private:
   Real pressure_floor_, density_floor_, internal_e_floor_;
+  Real velocity_ceiling_, internal_e_ceiling_;
 };
 
 #endif // EOS_EOS_HPP_

src/hydro/hydro.cpp

@@ -385,10 +385,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   auto first_order_flux_correct =
       pin->GetOrAddBoolean("hydro", "first_order_flux_correct", false);
-  if (first_order_flux_correct && integrator != Integrator::vl2) {
-    PARTHENON_FAIL("Please use 'vl2' integrator with first order flux correction. Other "
-                   "integrators have not been tested.")
-  }
   pkg->AddParam<>("first_order_flux_correct", first_order_flux_correct);
   if (first_order_flux_correct) {
     if (fluid == Fluid::euler) {
@@ -439,6 +435,23 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
       efloor = Tfloor / mbar_over_kb / (gamma - 1.0);
     }
 
+    // By default disable ceilings by setting to infinity
+    Real vceil =
+        pin->GetOrAddReal("hydro", "vceil", std::numeric_limits<Real>::infinity());
+    Real Tceil =
+        pin->GetOrAddReal("hydro", "Tceil", std::numeric_limits<Real>::infinity());
+    Real eceil = Tceil;
+    if (eceil < std::numeric_limits<Real>::infinity()) {
+      if (!pkg->AllParams().hasKey("mbar_over_kb")) {
+        PARTHENON_FAIL("Temperature ceiling requires units and gas composition. "
+                       "Either set a 'units' block and the 'hydro/He_mass_fraction' in "
+                       "input file or use a pressure floor "
+                       "(defined code units) instead.");
+      }
+      auto mbar_over_kb = pkg->Param<Real>("mbar_over_kb");
+      eceil = Tceil / mbar_over_kb / (gamma - 1.0);
+    }
+
     auto conduction = Conduction::none;
     auto conduction_str = pin->GetOrAddString("diffusion", "conduction", "none");
     if (conduction_str == "spitzer") {
@@ -472,12 +485,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     pkg->AddParam<>("conduction", conduction);
 
     if (fluid == Fluid::euler) {
-      AdiabaticHydroEOS eos(pfloor, dfloor, efloor, gamma);
+      AdiabaticHydroEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma);
       pkg->AddParam<>("eos", eos);
       pkg->FillDerivedMesh = ConsToPrim<AdiabaticHydroEOS>;
       pkg->EstimateTimestepMesh = EstimateTimestep<Fluid::euler>;
     } else if (fluid == Fluid::glmmhd) {
-      AdiabaticGLMMHDEOS eos(pfloor, dfloor, efloor, gamma);
+      AdiabaticGLMMHDEOS eos(pfloor, dfloor, efloor, vceil, eceil, gamma);
       pkg->AddParam<>("eos", eos);
       pkg->FillDerivedMesh = ConsToPrim<AdiabaticGLMMHDEOS>;
       pkg->EstimateTimestepMesh = EstimateTimestep<Fluid::glmmhd>;
@@ -975,6 +988,7 @@ TaskStatus FirstOrderFluxCorrect(MeshData<Real> *u0_data, MeshData<Real> *u1_dat
 
   std::vector<parthenon::MetadataFlag> flags_ind({Metadata::Independent});
   auto u0_cons_pack = u0_data->PackVariablesAndFluxes(flags_ind);
+  auto const &u0_prim_pack = u0_data->PackVariables(std::vector<std::string>{"prim"});
   auto u1_cons_pack = u1_data->PackVariablesAndFluxes(flags_ind);
   auto pkg = pmb->packages.Get("Hydro");
 
@@ -987,14 +1001,6 @@ TaskStatus FirstOrderFluxCorrect(MeshData<Real> *u0_data, MeshData<Real> *u1_dat
   if (fluid == Fluid::glmmhd) {
     c_h = pkg->Param<Real>("c_h");
   }
-  // Using "u1_prim" as "u0_prim" here because all current integrators start with copying
-  // the initial state to the "u0" register, see conditional for `stage == 1` in the
-  // hydro_driver where normally only "cons" is copied but in case for flux correction
-  // "prim", too. This means both during stage 1 and during stage 2 `u1` holds the
-  // original data at the beginning of the timestep. For flux correction we want to make a
-  // full (dt) low order update using the original data and thus use the "prim" data from
-  // u1 here.
-  auto const &u0_prim_pack = u1_data->PackVariables(std::vector<std::string>{"prim"});
 
   const int ndim = pmb->pmy_mesh->ndim;
 

src/main.cpp

@@ -85,7 +85,7 @@ int main(int argc, char *argv[]) {
     Hydro::ProblemInitPackageData = rand_blast::ProblemInitPackageData;
     Hydro::ProblemSourceFirstOrder = rand_blast::RandomBlasts;
   } else if (problem == "cluster") {
-    pman.app_input->ProblemGenerator = cluster::ProblemGenerator;
+    pman.app_input->MeshProblemGenerator = cluster::ProblemGenerator;
     pman.app_input->MeshBlockUserWorkBeforeOutput = cluster::UserWorkBeforeOutput;
     Hydro::ProblemInitPackageData = cluster::ProblemInitPackageData;
     Hydro::ProblemSourceUnsplit = cluster::ClusterSrcTerm;
@@ -119,10 +119,6 @@ int main(int argc, char *argv[]) {
     // This line actually runs the simulation
     driver.Execute();
   }
-  // very ugly cleanup...
-  if (problem == "turbulence") {
-    turbulence::Cleanup();
-  }
 
   // call MPI_Finalize and Kokkos::finalize if necessary
   pman.ParthenonFinalize();