Merge branch 'forrestglines/parthenon-pr930' into pgrete/turb-next

docs/cluster.md

@@ -449,7 +449,7 @@ The parameters $\alpha$ and $\ell$ can be changed with
 ```
 <problem/cluster/magnetic_tower>
 potential_type = li
-alpha = 20
+li_alpha = 20
 l_scale = 0.001
 ```
 When injected as a fraction of 
@@ -498,9 +498,9 @@ with all other components being zero.
 ```
 <problem/cluster/magnetic_tower>
 potential_type = donut
-l_scale = 0.0005   # in code length
-offset = 0.001     # in code length
-thickness = 0.0001 # in code length
+l_scale = 0.0005         # in code length
+donut_offset = 0.001     # in code length
+donut_thickness = 0.0001 # in code length
 ```
 
 It is recommended to match the donut thickness to the thickness of the kinetic jet launching region

docs/pgen.md

@@ -1,11 +1,73 @@
-# Problem specific callbacks
+# Problem generators
 
-## Source functions
+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
+generators by any user via pull requests.
+
+## Addding a new problem generator
+
+In general, four small steps are required:
+
+### 1. Add a new source file to the `src/pgen` folder
+
+The file shoud includes at least the `void ProblemGenerator(Mesh *pmesh, ParameterInput *pin, MeshData<Real> *md)`
+function, which is used to initialize the data at the beginning of the simulation.
+Alternatively, the `MeshBlock` version (`void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin)`) can be used that
+operates on a single block at a time rather than on a collection of blocks -- though this is not recommended from a performance point of view.
+
+The function (and all other functions in that file) should be encapsulated in their own namespace (that, ideally, is named
+identical to the file itself) so that the functions name are unique across different problem generators.
+
+Tip: Do not write the problem generator file from scratch but mix and match from existing ones,
+e.g., start with a simple one like [orszag_tang.cpp](../src/pgen/orszag_tang.cpp).
+
+### 2. Add new function(s) to `pgen.hpp`
+
+All callback functions, i.e., at least the `ProblemGenerator` plus additional optional ones (see step 5 below),
+need to be added to [pgen.hpp](../src/pgen/pgen.hpp).
+Again, just follow the existing pattern in the file and add the new function declarations with the appropriate namespace.
+
+### 3. Add callbacks to `main.cpp`
+
+All problem specific callback functions need to be enrolled in the [`main.cpp`](../src/main.cpp) file.
+The selection (via the input file) is controlled by the `problem_id` string in the `<job>` block.
+Again, for consistency it is recommended to pick a string that matches the namespace and problem generator file.
+
+### 4. Ensure new problem generator is compiled
+
+Add the new source file to [src/pgen/CMakeLists.txt](../src/pgen/CMakeLists.txt) so that it will be compiled
+along all other problem generators.
+
+### 5. (Optional) Add more additional callback
+
+In addition to the `ProblemGenerator` that initializes the data, other callback functions exists
+that allow to modify the behavior of AthenaPK on a problem specific basis.
+See [Callback functions](#Callback-functions)] below for available options.
+
+### 6. (Optional but most likely required) Write an input file
+
+In theory, one can hardcode all paramters in the source file (like in the
+[orszag_tang.cpp](../src/pgen/orszag_tang.cpp) problem generator) but it
+prevents modification of the problem setup once the binary is compiled.
+
+The more common usecase is to create an input file that contains a problem specific
+input block.
+The convention here is to have a block named `<problem/NAME>` where `NAME` is the name
+of the problem generator (or namespace used).
+For example, the Sod shock tube problem generator processes the input file with lines like
+```
+  Real rho_l = pin->GetOrAddReal("problem/sod", "rho_l", 1.0);
+```
+to set the density on the left hand side of the domain.
+
+## Callback functions
+
+### Source functions
 
 Additional (physical) source terms (e.g., the ones typically on the right hand side of
 system of equations) can be added in various ways from an algorithmic point of view.
 
-### Unsplit
+#### Unsplit
 
 Unsplit sources are added at each stage of the (multi-stage) integration after the
 (conserved) fluxes are calculated.
@@ -23,7 +85,7 @@ by implementing a function with that signature and assigning it to the
 
 Note, there is no requirement to call the function `MyUnsplitSource`.
 
-### Split first order (generally not recommended)
+#### Split first order (generally not recommended)
 
 If for special circumstances sources need to be added in a fully operator split way,
 i.e., the source function is only called *after* the full hydro (or MHD) integration,
@@ -38,12 +100,27 @@ Note, as the name suggests, this method is only first order accurate (while ther
 is no requirement to call the function `MyFirstOrderSource`).
 
 
-### Split second order (Strang splitting)
+#### Split second order (Strang splitting)
 
-Not implemented yet.
+Strang splitting achieves second order by interleaving the main hydro/MHD update
+with a source update.
+In practice, AthenaPK calls Strang split source with `0.5*dt` before the first stage of
+each integrator and with `0.5*dt` after the last stage of the integrator.
+Note that for consistency, a Strang split source terms should update both conserved
+and primitive variables in all zones (i.e., also in the ghost zones as those
+are currently not updated after calling a source function).
+
+Users can enroll a custom source function
+```c++
+void MyStrangSplitSource(MeshData<Real> *md, const Real beta_dt);
+```
+by implementing a function with that signature and assigning it to the
+`Hydro::ProblemSourceStrangSplit` callback (currently in `main.cpp`).
+
+Note, there is no requirement to call the function `MyStrangSplitSource`.
 
 
-## Timestep restrictions
+### Timestep restrictions
 
 If additional problem specific physics are implemented (e.g., through the source
 functions above) that require a custom timestep, it can be added via the
@@ -55,4 +132,60 @@ Real ProblemEstimateTimestep(MeshData<Real> *md);
 The return value is expected to be the minimum value over all blocks in the
 contained in the `MeshData` container, cf., the hydro/mhd `EstimateTimestep` function.
 Note, that the hyperbolic CFL value is currently applied after the function call, i.e.,
-it is also applied to the problem specific function.
+it is also applied to the problem specific function.
+
+### Additional initialization on startup (adding variables/fields, custom history output, ...)
+
+Sometimes a problem generator requires a more general processing of the input file
+and/or needs to make certain global adjustments (e.g., adding a custom callback function
+for the history output or adding a new field).
+This can be achieved via modifying the AthenaPK "package" (currently called
+`parthenon::StateDescriptor`) through the following function.
+```c++
+void ProblemInitPackageData(ParameterInput *pin, parthenon::StateDescriptor *pkg)
+```
+
+For example, the [src/pgen/turbulence.cpp](../[src/pgen/turbulence.cpp]) problem generator
+add an additional field (here for an acceleration field) by calling
+```c++
+  Metadata m({Metadata::Cell, Metadata::Derived, Metadata::OneCopy},
+             std::vector<int>({3}));
+  pkg->AddField("acc", m);
+```
+in the `ProblemInitPackageData`.
+
+### Additional initialization on mesh creation/remeshing/load balancing
+
+For some problem generators it is required to initialize data on "new" blocks.
+These new blocks can, for example, be created during mesh refinement
+(or derefinement) and this data is typically not active data (like
+conserved or primitive variables as those are handled automatically)
+but more general data.
+Once example is the phase information in the turbulence driver that
+does not vary over time but spatially (and therefore at a per block level).
+
+The appropriate callback to enroll is
+```c++
+void InitMeshBlockUserData(MeshBlock *pmb, ParameterInput *pin)
+```
+
+### UserWorkAfterLoop
+
+If additional work is required once the main loop finishes (i.e., once the
+simulation reaches its final time) computations can be done inside the
+```c++
+void UserWorkAfterLoop(Mesh *mesh, ParameterInput *pin, parthenon::SimTime &tm)
+```
+callback function.
+This is, for example, done in the linear wave problem generator to calculate the
+error norms for convergence tests.
+
+### MeshBlockUserWorkBeforeOutput
+
+Sometimes it is desirable to further process data before an output file is written.
+For this the
+```c++
+void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin)
+```
+callback is available, that is called once every time right before a data output
+(hdf5 or restart) is being written.

inputs/cluster/cluster.in

@@ -182,7 +182,7 @@ kinetic_jet_offset = 0.05
 
 <problem/cluster/magnetic_tower>
 potential_type = li
-alpha = 20
+li_alpha = 20
 l_scale = 0.001
 initial_field = 0.12431560000204142
 #NOTE: Change to this line to disable initial magnetic tower

inputs/cluster/magnetic_tower.in

@@ -110,7 +110,7 @@ thermal_fraction = 0
 
 <problem/cluster/magnetic_tower>
 potential_type = li
-alpha = 20
+li_alpha = 20
 l_scale = 0.01
 initial_field = 0.12431560000204142
 fixed_field_rate = 12.431560000204144

src/hydro/hydro.cpp

@@ -143,6 +143,14 @@ TaskStatus AddUnsplitSources(MeshData<Real> *md, const SimTime &tm, const Real b
   if (hydro_pkg->Param<Fluid>("fluid") == Fluid::glmmhd) {
     hydro_pkg->Param<GLMMHD::SourceFun_t>("glmmhd_source")(md, beta_dt);
   }
+  const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
+
+  if (enable_cooling == Cooling::tabular) {
+    const TabularCooling &tabular_cooling =
+        hydro_pkg->Param<TabularCooling>("tabular_cooling");
+
+    tabular_cooling.SrcTerm(md, beta_dt);
+  }
   if (ProblemSourceUnsplit != nullptr) {
     ProblemSourceUnsplit(md, tm, beta_dt);
   }
@@ -153,32 +161,13 @@ TaskStatus AddUnsplitSources(MeshData<Real> *md, const SimTime &tm, const Real b
 TaskStatus AddSplitSourcesFirstOrder(MeshData<Real> *md, const SimTime &tm) {
   auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
 
-  const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
-
-  if (enable_cooling == Cooling::tabular) {
-    const TabularCooling &tabular_cooling =
-        hydro_pkg->Param<TabularCooling>("tabular_cooling");
-
-    tabular_cooling.SrcTerm(md, tm.dt);
-  }
   if (ProblemSourceFirstOrder != nullptr) {
     ProblemSourceFirstOrder(md, tm, tm.dt);
   }
   return TaskStatus::complete;
 }
 
 TaskStatus AddSplitSourcesStrang(MeshData<Real> *md, const SimTime &tm) {
-  // auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
-
-  // const auto &enable_cooling = hydro_pkg->Param<Cooling>("enable_cooling");
-
-  // if (enable_cooling == Cooling::tabular) {
-  //   const TabularCooling &tabular_cooling =
-  //       hydro_pkg->Param<TabularCooling>("tabular_cooling");
-
-  //   tabular_cooling.SrcTerm(md, 0.5 * tm.dt);
-  // }
-
   if (ProblemSourceStrangSplit != nullptr) {
     ProblemSourceStrangSplit(md, tm, tm.dt);
   }
@@ -788,8 +777,8 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
   int il, iu, jl, ju, kl, ku;
   jl = jb.s, ju = jb.e, kl = kb.s, ku = kb.e;
   // TODO(pgrete): are these looop limits are likely too large for 2nd order
-  if (pmb->block_size.nx2 > 1) {
-    if (pmb->block_size.nx3 == 1) // 2D
+  if (pmb->block_size.nx(X2DIR) > 1) {
+    if (pmb->block_size.nx(X3DIR) == 1) // 2D
       jl = jb.s - 1, ju = jb.e + 1, kl = kb.s, ku = kb.e;
     else // 3D
       jl = jb.s - 1, ju = jb.e + 1, kl = kb.s - 1, ku = kb.e + 1;
@@ -861,7 +850,7 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshData<Real>> &md) {
         parthenon::ScratchPad2D<Real>::shmem_size(num_scratch_vars, nx1) * 3;
     // set the loop limits
     il = ib.s - 1, iu = ib.e + 1, kl = kb.s, ku = kb.e;
-    if (pmb->block_size.nx3 == 1) // 2D
+    if (pmb->block_size.nx(X3DIR) == 1) // 2D
       kl = kb.s, ku = kb.e;
     else // 3D
       kl = kb.s - 1, ku = kb.e + 1;

src/main.cpp

@@ -88,7 +88,8 @@ int main(int argc, char *argv[]) {
     pman.app_input->MeshProblemGenerator = cluster::ProblemGenerator;
     pman.app_input->MeshBlockUserWorkBeforeOutput = cluster::UserWorkBeforeOutput;
     Hydro::ProblemInitPackageData = cluster::ProblemInitPackageData;
-    Hydro::ProblemSourceUnsplit = cluster::ClusterSrcTerm;
+    Hydro::ProblemSourceUnsplit = cluster::ClusterUnsplitSrcTerm;
+    Hydro::ProblemSourceFirstOrder = cluster::ClusterSplitSrcTerm;
     Hydro::ProblemEstimateTimestep = cluster::ClusterEstimateTimestep;
   } else if (problem == "sod") {
     pman.app_input->ProblemGenerator = sod::ProblemGenerator;

src/pgen/CMakeLists.txt

@@ -8,6 +8,8 @@ target_sources(athenaPK PRIVATE
     cluster.cpp
     cluster/agn_feedback.cpp
     cluster/agn_triggering.cpp
+    cluster/cluster_clips.cpp
+    cluster/cluster_reductions.cpp
     cluster/hydrostatic_equilibrium_sphere.cpp
     cluster/magnetic_tower.cpp
     cluster/snia_feedback.cpp