Fix KHI pgen. Add turbulence pgen doc.

CHANGELOG.md

@@ -7,6 +7,7 @@
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 119]](https://github.com/parthenon-hpc-lab/athenapk/pull/119) Fixed Athena++ paper test case for KHI pgen. Added turbulence pgen doc.
 - [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions.
 - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
 

docs/README.md

@@ -0,0 +1,36 @@
+# AthenaPK documentation
+
+Note that we're aware the rendering of equations in markdown on GitHub for the documentation
+is currently not great.
+Eventually, the docs will transition to Read the Docs (or similar).
+
+## Overview
+
+The documentation currently includes
+
+- [Configuring solvers in the input file](input.md)
+  - [An extended overview of various conventions for optically thin cooling implementations](cooling_notes.md)
+  - [Notebooks to calculate cooling tables from literature](cooling)
+- [Brief notes on developing code for AthenaPK](development.md)
+- [How to add a custom/user problem generator](pgen.md)
+- Detailed descriptions of more complex problem generators
+  - [Galaxy Cluster and Cluster-like Problem Setup](cluster.md)
+  - [Driven turbulence](turbulence.md)
+
+## Tutorial
+
+An AthenaPK tutorial was given as part of the **Towards exascale-ready astrophysics**
+workshop https://indico3-jsc.fz-juelich.de/event/169/ taking place 25-27 Sep 2024 online.
+
+The material is currently located at https://github.com/pgrete/athenapk_tutorial
+
+While the instructions for building the code are specific to the workshop environment
+the tutorial itself should translate directly to other environments/systems.
+
+## Parthenon documenation
+
+Many paramters/options are directly controlled through the Parthenon framework
+(both with regard to building and in the input file).
+
+While the [Parthenon documenation](https://parthenon-hpc-lab.github.io/parthenon) is
+more geared towards developers it also contains useful information for users.

docs/pgen.md

@@ -1,4 +1,4 @@
-# Problem generators
+# Custom problem generators
 
 Different simulation setups in AthenaPK are controlled via the so-called problem generators.
 New problem generators can easily be added and we are happy to accept and merge contibuted problem

docs/turbulence.md

@@ -0,0 +1,100 @@
+# Driven turbulence simulations
+
+The turbulence problem generator uses explicit inverse Fourier transformations (iFTs)
+on each meshblock in order to reduce communication during the iFT.
+Thus, it is only efficient if comparatively few modes are used (say < 100).
+
+Quite generally, driven turbulence simulations start from uniform initial conditions
+(uniform density and pressure, some initial magnetic field configuration in case of an
+MHD setup, and the fluid at rest) and reach a state of stationary, isotropic (or anisotropic
+depending on the strength of the background magnetic field) turbulence after one to few
+large eddy turnover times (again depending on the background magnetic field strength).
+The large eddy turnover time is usually defined as `T = L/U` with `L` being the scale
+of the largest eddies and `U` the root mean square Mach number in the stationary regime.
+
+The current implementation uses the following forcing spectrum
+`(k/k_peak)^2 * (2 - (k/k_peak)^2)`.
+Here, `k_peak` is the peak wavenumber of the forcing spectrum. It is related the scales of the largest eddies as
+`L = 1/k_f` given that a box size of 1 is currently assumed/hardcoded.
+
+## Problem setup
+
+An example parameter file can be found in `inputs/turbulence.in`.
+
+A typical setup contains the following blocks in the input file:
+
+```
+<job>
+problem_id = turbulence
+
+<problem/turbulence>
+rho0         = 1.0      # initial mean density
+p0           = 1.0      # initial mean pressure
+b0           = 0.01     # initial magnetic field strength
+b_config     = 0        # 0 - net flux; 1 - no net flux uniform B; 2 - non net flux sin B; 4 - field loop
+kpeak        = 2.0      # characteristic wavenumber
+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
+
+<modes>
+k_1_0	= +2
+k_1_1	= -1
+k_1_2	= +0
+k_2_0	= +1
+...
+```
+
+Following parameters can be changed to control both the initial state 
+
+- `rho0` initial mean density
+- `p0` initial mean thermal pressure
+- `b0` initial mean magnetic field strength
+- `b_config`
+  - `0`: net flux case (uniform B_x)
+  - `1`: no net flux case (uniform B_x with changing sign in each half of the box)
+  - `2`: no net flux with initial sinosoidal B_x field
+  - `3`: deprecated
+  - `4`: closed field loop/cylinder in the box (in x-y plane) located at
+    - `x0=0.5` (default)
+    - `y0=0.5` (default)
+    - `z0=0.5` (default)
+    - and radius `loop_rad=0.25`
+
+as well as the driving field
+
+- `kpeak` peak wavenumber of the forcing spectrum. Make sure to update the wavemodes to match `kpeak`, see below.
+- `corr_time` autocorrelation time of the acceleration field (in code units).
+Using delta-in-time forcing, i.e., a very low value, is discouraged, see [Grete et al. 2018 ApJL](https://iopscience.iop.org/article/10.3847/2041-8213/aac0f5).
+- `rseed` random seed for the OU process. Only change for new simulation, but keep unchanged for restarting simulations.
+- `sol_weight` solenoidal weight of the acceleration field. `1.0` is purely solenoidal/rotational and `0.0` is purely dilatational/compressive. Any value between `0.0` and `1.0` is possible. The parameter is related to the resulting rotational power in the 3D acceleration field as
+`1. - ((1-sol_weight)^2/(1-2*sol_weight+3*sol_weight^2))`, see eq (9) in [Federrath et al. 2010 A&A](
+https://doi.org/10.1051/0004-6361/200912437).
+- `accel_rms` root mean square value of the acceleration (controls the "strength" of the forcing field)
+- `num_modes` number of wavemodes that are specified in the `<modes>` section of the parameter file.
+In order to generate a full set of modes run the `inputs/generate_fmturb_modes.py` script and replace
+the corresponding parts of the parameter file with the output of the script.
+Within the script, the top three variables (`k_peak`, `k_high`, and `k_low`) need to be adjusted in
+order to generate a complete set (i.e., all) of wavemodes.
+Alternatively, wavemodes can be chosen/defined manually, e.g., if not all wavemodes are desired or
+only individual modes should be forced.
+
+## Typical results
+
+The results shown here are obtained from running simulations with the parameters given in the next section.
+
+### High level temporal evolution
+![image](img/turb_evol.png)
+
+### Power spectra
+![image](img/turb_spec.png)
+
+### Consistency of acceleration field
+As each meshblock does a full iFT of all modes the following slices from a run with 8 meshblocks 
+illustrate that there is no discontinuities at the meshblock boundary.
+
+Plot shows x-, y-, and z-acceleration (in rows top to bottom) slices in the x-, y-, and z-direction (in columns from left to right).
+
+![image](img/turb_acc.png)

src/main.cpp

@@ -81,7 +81,7 @@ int main(int argc, char *argv[]) {
     pman.app_input->ProblemGenerator = field_loop::ProblemGenerator;
     Hydro::ProblemInitPackageData = field_loop::ProblemInitPackageData;
   } else if (problem == "kh") {
-    pman.app_input->ProblemGenerator = kh::ProblemGenerator;
+    pman.app_input->MeshProblemGenerator = kh::ProblemGenerator;
   } else if (problem == "rand_blast") {
     pman.app_input->ProblemGenerator = rand_blast::ProblemGenerator;
     Hydro::ProblemInitPackageData = rand_blast::ProblemInitPackageData;

src/pgen/kh.cpp

@@ -32,61 +32,32 @@
 
 // AthenaPK headers
 #include "../main.hpp"
+#include "utils/error_checking.hpp"
 
 namespace kh {
 using namespace parthenon::driver::prelude;
 
 //----------------------------------------------------------------------------------------
-//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin)
+//! \fn void Mesh::ProblemGenerator(Mesh *pm, ParameterInput *pin, MeshData<Real> *md)
 //  \brief Problem Generator for the Kelvin-Helmholtz test
 
-void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
+void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md) {
   auto vflow = pin->GetReal("problem/kh", "vflow");
   auto iprob = pin->GetInteger("problem/kh", "iprob");
+  // Get pointer to first block (always exists) for common data like loop bounds
+  auto pmb = md->GetBlockData(0)->GetBlockPointer();
   auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
   auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
   auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);
   auto gam = pin->GetReal("hydro", "gamma");
   auto gm1 = (gam - 1.0);
 
-  // initialize conserved variables
-  auto &rc = pmb->meshblock_data.Get();
-  auto &u_dev = rc->Get("cons").data;
-  auto &coords = pmb->coords;
-  // initializing on host
-  auto u = u_dev.GetHostMirrorAndCopy();
+  // Initialize conserved variables
+  // Get a MeshBlockPack on device with all conserved variables
+  const auto &cons = md->PackVariables(std::vector<std::string>{"cons"});
+  const auto num_blocks = md->NumBlocks();
 
-  std::mt19937 gen(pmb->gid); // Standard mersenne_twister_engine seeded with gid
-  std::uniform_real_distribution<Real> ran(-0.5, 0.5);
-
-  //--- iprob=1.  Uniform stream with density ratio "drat" located in region -1/4<y<1/4
-  // moving at (-vflow) seperated by two slip-surfaces from background medium with d=1
-  // moving at (+vflow), random perturbations.  This is the classic, unresolved K-H test.
-
-  if (iprob == 1) {
-    // Read problem parameters
-    Real drat = pin->GetReal("problem/kh", "drat");
-    Real amp = pin->GetReal("problem/kh", "amp");
-    for (int k = kb.s; k <= kb.e; k++) {
-      for (int j = jb.s; j <= jb.e; j++) {
-        for (int i = ib.s; i <= ib.e; i++) {
-          u(IDN, k, j, i) = 1.0;
-          u(IM1, k, j, i) = vflow + amp * ran(gen);
-          u(IM2, k, j, i) = amp * ran(gen);
-          u(IM3, k, j, i) = 0.0;
-          if (std::abs(coords.Xc<2>(j)) < 0.25) {
-            u(IDN, k, j, i) = drat;
-            u(IM1, k, j, i) = -drat * (vflow + amp * ran(gen));
-            u(IM2, k, j, i) = drat * amp * ran(gen);
-          }
-          // Pressure scaled to give a sound speed of 1 with gamma=1.4
-          u(IEN, k, j, i) =
-              2.5 / gm1 +
-              0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i);
-        }
-      }
-    }
-  }
+  //--- iprob=1. This was the classic, unresolved K-H test.
 
   //--- iprob=2. Uniform density medium moving at +/-vflow seperated by a single shear
   // layer with tanh() profile at y=0 with a single mode perturbation, reflecting BCs at
@@ -97,9 +68,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
     Real amp = pin->GetReal("problem/kh", "amp");
     Real a = 0.02;
     Real sigma = 0.2;
-    for (int k = kb.s; k <= kb.e; k++) {
-      for (int j = jb.s; j <= jb.e; j++) {
-        for (int i = ib.s; i <= ib.e; i++) {
+    pmb->par_for(
+        "KHI: iprob2", 0, num_blocks - 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) {
+          const auto &u = cons(b);
+          const auto &coords = cons.GetCoords(b);
           u(IDN, k, j, i) = 1.0;
           u(IM1, k, j, i) = vflow * std::tanh((coords.Xc<2>(j)) / a);
           u(IM2, k, j, i) = amp * std::cos(2.0 * M_PI * coords.Xc<1>(i)) *
@@ -108,23 +81,22 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
           u(IEN, k, j, i) =
               1.0 / gm1 +
               0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i);
-        }
-      }
-    }
-  }
+        });
 
-  //--- iprob=3.  Test in SR paper (Beckwith & Stone, ApJS 193, 6, 2011).  Gives two
-  // resolved shear layers with tanh() profiles for velocity and density located at
-  // y = +/- 0.5, density one in middle and 0.01 elsewhere, single mode perturbation.
+  } else if (iprob == 3) {
+    //--- iprob=3.  Test in SR paper (Beckwith & Stone, ApJS 193, 6, 2011).  Gives two
+    // resolved shear layers with tanh() profiles for velocity and density located at
+    // y = +/- 0.5, density one in middle and 0.01 elsewhere, single mode perturbation.
 
-  if (iprob == 3) {
     // Read/set problem parameters
     Real amp = pin->GetReal("problem/kh", "amp");
     Real a = 0.01;
     Real sigma = 0.1;
-    for (int k = kb.s; k <= kb.e; k++) {
-      for (int j = jb.s; j <= jb.e; j++) {
-        for (int i = ib.s; i <= ib.e; i++) {
+    pmb->par_for(
+        "KHI: iprob3", 0, num_blocks - 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) {
+          const auto &u = cons(b);
+          const auto &coords = cons.GetCoords(b);
           u(IDN, k, j, i) =
               0.505 + 0.495 * std::tanh((std::abs(coords.Xc<2>(j)) - 0.5) / a);
           u(IM1, k, j, i) = vflow * std::tanh((std::abs(coords.Xc<2>(j)) - 0.5) / a);
@@ -139,17 +111,15 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
           u(IEN, k, j, i) =
               1.0 / gm1 +
               0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i);
-        }
-      }
-    }
-  }
+        });
 
-  //--- iprob=4.  "Lecoanet" test, resolved shear layers with tanh() profiles for velocity
-  // and density located at z1=0.5, z2=1.5 two-mode perturbation for fully periodic BCs
+  } else if (iprob == 4) {
+    //--- iprob=4.  "Lecoanet" test, resolved shear layers with tanh() profiles for
+    // velocity
+    // and density located at z1=0.5, z2=1.5 two-mode perturbation for fully periodic BCs
 
-  // To promote symmetry of FP errors about midplanes, rescale z' = z - 1. ; x' = x - 0.5
-  // so that domain x1 = [-0.5, 0.5] and x2 = [-1.0, 1.0] is centered about origin
-  if (iprob == 4) {
+    // To promote symmetry of FP errors about midplanes, rescale z' = z - 1. ; x' = x -
+    // 0.5 so that domain x1 = [-0.5, 0.5] and x2 = [-1.0, 1.0] is centered about origin
     // Read/set problem parameters
     Real amp = pin->GetReal("problem/kh", "amp");
     // unstratified problem is the default
@@ -164,9 +134,11 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
     Real z1 = -0.5; // z1' = z1 - 1.0
     Real z2 = 0.5;  // z2' = z2 - 1.0
 
-    for (int k = kb.s; k <= kb.e; k++) {
-      for (int j = jb.s; j <= jb.e; j++) {
-        for (int i = ib.s; i <= ib.e; i++) {
+    pmb->par_for(
+        "KHI: iprob4", 0, num_blocks - 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) {
+          const auto &u = cons(b);
+          const auto &coords = cons.GetCoords(b);
           // Lecoanet (2015) equation 8a)
           Real dens = 1.0 + 0.5 * drho_rho0 *
                                 (std::tanh((coords.Xc<2>(j) - z1) / a) -
@@ -233,40 +205,39 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
               P0 / gm1 +
               0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i)) + SQR(u(IM3, k, j, i))) /
                   u(IDN, k, j, i);
-        }
-      }
-    }
-    // copy initialized vars to device
-    u_dev.DeepCopy(u);
-  }
+        });
 
-  //--- iprob=5. Uniform stream with density ratio "drat" located in region -1/4<y<1/4
-  // moving at (-vflow) seperated by two resolved slip-surfaces from background medium
-  // with d=1 moving at (+vflow), with m=2 perturbation, for the AMR test.
+  } else if (iprob == 5) {
+    //--- iprob=5. Uniform stream with density ratio "drat" located in region -1/4<y<1/4
+    // moving at (-vflow) seperated by two resolved slip-surfaces from background medium
+    // with d=1 moving at (+vflow), with m=2 perturbation, for the AMR test.
+    // Note that the code below should match what's in the Athena++ method paper (but it
+    // does not match anymore what's implemented in Athena++ itself).
 
-  if (iprob == 5) {
     // Read problem parameters
     Real a = pin->GetReal("problem/kh", "a");
     Real sigma = pin->GetReal("problem/kh", "sigma");
     Real drat = pin->GetReal("problem/kh", "drat");
     Real amp = pin->GetReal("problem/kh", "amp");
-    for (int k = kb.s; k <= kb.e; k++) {
-      for (int j = jb.s; j <= jb.e; j++) {
-        for (int i = ib.s; i <= ib.e; i++) {
+    pmb->par_for(
+        "KHI: iprob5", 0, num_blocks - 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) {
+          const auto &u = cons(b);
+          const auto &coords = cons.GetCoords(b);
           Real w = (std::tanh((std::abs(coords.Xc<2>(j)) - 0.25) / a) + 1.0) * 0.5;
           u(IDN, k, j, i) = w + (1.0 - w) * drat;
-          u(IM1, k, j, i) = w * vflow - (1.0 - w) * vflow * drat;
+          u(IM1, k, j, i) = u(IDN, k, j, i) * vflow * (w - 0.5);
           u(IM2, k, j, i) =
-              u(IDN, k, j, i) * amp * std::sin(2.0 * 2.0 * M_PI * coords.Xc<1>(i)) *
+              u(IDN, k, j, i) * amp * std::cos(2.0 * 2.0 * M_PI * coords.Xc<1>(i)) *
               std::exp(-SQR(std::abs(coords.Xc<2>(j)) - 0.25) / (sigma * sigma));
           u(IM3, k, j, i) = 0.0;
           // Pressure scaled to give a sound speed of 1 with gamma=1.4
           u(IEN, k, j, i) =
               2.5 / gm1 +
-              0.25 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i);
-        }
-      }
-    }
+              0.5 * (SQR(u(IM1, k, j, i)) + SQR(u(IM2, k, j, i))) / u(IDN, k, j, i);
+        });
+  } else {
+    PARTHENON_FAIL("Unknow iprob for KHI pgen.")
   }
 }
 

src/pgen/pgen.hpp

@@ -84,7 +84,7 @@ void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg
 namespace kh {
 using namespace parthenon::driver::prelude;
 
-void ProblemGenerator(MeshBlock *pmb, parthenon::ParameterInput *pin);
+void ProblemGenerator(Mesh *pm, parthenon::ParameterInput *pin, MeshData<Real> *md);
 } // namespace kh
 namespace rand_blast {
 using namespace parthenon::driver::prelude;