Add support for tracer particles (#102)

* Add Phoebus tracers

* Fix/adapt tracers to AthenaPK

* Add tracer support to turbulence pgen

* Fix swarm interface

* Fix particle init

* Call FillTracer every cycle

* Added brief instructions for new particle field

* Update tracers.cpp

filled in values for s_i

* Simplify/dedup code

* Add lookback time

* Update turbulence.cpp

* Update turbulence.cpp

* Update turbulence.cpp

small edits in the par_reduce in the computation of the correlations

* Update tracers.cpp

* Update turbulence.cpp

* Update turbulence.cpp

* Use MeshData in FillTracers to allow reduction

* Update FillTracers internals to loop over blocks logic

* Calc correlations

* Dump corr to file

* Import interpolation from Phoebus

* Adjust interface

* User interpolation for FillTracers

* Use 2nd order time integrator for particles

* Update tracers.cpp

fixed a bug in computing the correlations

* Fix csv header

* Add s and sdot to csv output

* Fix interface

* Use upstream interfaces

* Fix swarm interface

* Add simple particle advection test

* update GH artifact

* Prevent artifact upload conflict

* Move particle seeding to tracer pkg

* Update particle seeding logic

* Add tracer seed and init callback function, and doc

* Separate default and pget FillTracer

* add plotting script

* add example tracer input

* add tracer example input

* FillTracers after seeding

* Formatting, copyright, minor review comments

* Add xdfm swarm comment to sample input file

* Add changelog and README

---------

Co-authored-by: mbruggen <mbrueggen@hs.uni-hamburg.de>
Co-authored-by: evan1022 <evan.scannapieco@asu.edu>
Co-authored-by: Ben Wibking <ben@wibking.com>

Author: 4 people

CHANGELOG.md

@@ -10,6 +10,7 @@ In between a subtle bug was introduced that resulted in inconsistent divergence
 refinement.
 
 ### Added (new features/APIs/variables/...)
+- [[PR102]](https://github.com/parthenon-hpc-lab/athenapk/pull/102) Add support for tracer particles
 - [[PR 89]](https://github.com/parthenon-hpc-lab/athenapk/pull/89) Add viscosity and resistivity
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
 

README.md

@@ -22,6 +22,7 @@ Current features include
     - unsplit
     - operator-split, second-order RKL2 supertimestepping
   - optically thin cooling based on tabulated cooling tables with either Townsend 2009 exact integration or operator-split subcycling
+  - tracer particles
 - static and adaptive mesh refinement
 - problem generators for
   - linear waves
@@ -84,7 +85,7 @@ The following examples are a few standard cases.
 
 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).
+A full list of architecture keywords is available on the [Kokkos wiki](https://kokkos.org/kokkos-core-wiki/get-started/configuration-guide.html#keywords-arch).
 
     # configure with enabling Intel Broadwell or similar architecture (AVX2) instructions
     cmake -S. -Bbuild-host -DKokkos_ARCH_BDW=ON -DPARTHENON_DISABLE_MPI=ON

docs/img/tracer_example.png

@@ -323,3 +323,88 @@ d_log_temp_tol = 1e-8              # Tolerance in cooling table between subseque
 Finally, a more comprehensive descriptions of various conventions used for cooling
 functions (and the chosen implementation in AthenaPK) can be found in [Cooling Notes](cooling_notes.md)
 and a notebook comparing various cooling tables (and their conversion) in [cooling/cooling.ipynb](cooling/cooling.ipynb).
+
+### Particles
+
+#### Tracers
+
+Tracer particles can be enabled by setting `tracers/enabled=true` in the input file:
+
+```
+<tracers>
+enabled = true
+
+initial_seed_method  = random_per_block    # alternative: user
+initial_num_tracers_per_cell = 0.125
+### Optional arguments
+#initial_rng_seed = INTEGER
+```
+
+Two seeding methods are currently supported:
+
+- `initial_seed_method=random_per_block`
+  - seeds particles at random positions in each block
+  - NOTE: the random number generator seed uses the unique block id. Therefore, simulations with the mesh decomposition (mesh and meshblock sizes) are identical independent of the number of MPI ranks used, but if the meshblock size is changed for given mesh (and thus the total number of blocks) the initial state will be different.
+  - `initial_num_tracers_per_cell` determines the number of seeded particles per cell
+  - `initial_rng_seed` (optional) is used as seed in addition to the block id.
+- `initial_seed_method=user`
+  - Calls a problem specific callback function (`ProblemSeedInitialTracers`), see [tracer callback documenation](https://github.com/parthenon-hpc-lab/athenapk/blob/main/docs/pgen.md#tracers).
+
+By default, swarm fields are written only to restart files.
+If they are required for "standard" output files (like single precision `hdf5`),
+they need to be added manually to the output block, e.g., (bottom two lines)
+```
+<parthenon/output2>
+file_type  = hdf5       # Binary data dump
+variables   = prim   # variables to be output
+dt         = 0.1        # time increment between outputs
+id         = prim
+single_precision_output = true
+
+swarms = tracers
+tracers_variables = id, x, y, z, rho
+#write_swarm_xdmf=true  # uncomment to create an xdmf output (e.g., for Paraview or Visit)
+```
+
+Tracers can be read/processed by Paraview or Visit (via the xdmf file)
+or by the `phdf` package shipped with the Parthenon submodule.
+A sample plotting script for the latter might look like
+```python
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+sys.path.insert(
+    1,
+    "/path/to/athenapk/external/parthenon"
+    + "/scripts/python/packages/parthenon_tools/parthenon_tools",
+)
+
+try:
+    import phdf
+except ModuleNotFoundError:
+    print("Couldn't find module to read Parthenon hdf5 files.")
+
+data = phdf.phdf("../run-turbulence/parthenon.prim.00010.phdf")
+tracers = data.GetSwarm("tracers")
+xs = tracers.x
+ys = tracers.y
+zs = tracers.z
+ids = tracers.Get("id")
+
+fig = plt.figure()
+ax = fig.add_subplot(projection='3d')
+
+ax.scatter(xs, ys, zs,s=1,c=ids)
+```
+resulting in the following image:
+
+![image](img/tracer_example.png)
+
+Following "restrictions" apply to the current tracer implementation:
+- Only 3D simulations.
+- Only one advection method (RK2/Heun's method).
+- All primitive fields (`rho`, `pressure`, `vel_x`, `vel_y`, `vel_z`, `B_x`, `B_y`, and `B_x`) are traced by default (independent of whether they're needed or not) in addition to the position (`x`, `y`, `z`) and id (`id`) fields.
+- Default tracer values (such as the primitive fields) are only updated right before writing an output file.
+- Ids are only unique if tracers are seeded at the beginning at the simulations and no new tracer particles are added dynamically while the simulation is running.
+
+Please get in touch, if you interested in running simulations that require lifting one (or more) of those restrictions.

docs/input.md

@@ -189,3 +189,29 @@ 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.
+
+### Particles
+
+#### Tracers
+
+In order to allow custom seeding of tracer particles (per problem generator), the
+`ProblemSeedInitialTracers` function can be defined.
+```c++
+void ProblemSeedInitialTracers(Mesh *pmesh, ParameterInput *pin, parthenon::SimTime &tm);
+```
+It is executed once when a simulation is started the very first time (i.e., it
+is not called on subsequent restarts).
+See the standard tracer seeding implementation[`SeedInitialTracers`](https://github.com/parthenon-hpc-lab/athenapk/blob/main/src/tracers/tracers.cpp) for reference on how to deposit particles.
+
+
+To allow adding additional fields, the
+```c++
+void ProblemInitTracerData(ParameterInput * pin, parthenon::StateDescriptor *tracer_pkg);```
+callback is available.
+It is called at the end of the tracer package initialization and can be defined at the per-problem-generator level, see, e.g., the turbulence driver as an example.
+
+Similarly, a callback is available to fill those new fields (or more) via
+```c++
+TaskStatus ProblemFillTracers(MeshData<Real> *md, parthenon::SimTime &tm, const Real dt);
+```
+It is called right after the default `FillTracers` task in the driver.

docs/pgen.md

@@ -0,0 +1,178 @@
+<comment>
+problem   = Turbulence with parabolic forcing profile (using few driving modes)
+reference = https://github.com/parthenon-hpc-lab/athenapk/blob/main/docs/turbulence.md
+
+<job>
+problem_id = turbulence
+
+<parthenon/output1>
+file_type  = hst        # History data dump
+dt         = 0.05       # time increment between outputs
+
+<parthenon/output2>
+file_type  = hdf5       # Binary data dump
+variables  = prim,acc   # variables to be output
+dt         = 0.1        # time increment between outputs
+id         = prim
+single_precision_output = true
+swarms = tracers
+tracers_variables = id, x, y, z, rho, s, sdot
+#write_swarm_xdmf=true  # uncomment to create an xdmf output (e.g., for Paraview or Visit)
+
+<parthenon/output3>
+file_type  = rst        # Binary data dump
+dt         = 0.1        # time increment between outputs
+id         = restart
+
+<parthenon/time>
+cfl = 0.3               # The Courant, Friedrichs, & Lewy (CFL) Number
+nlim       = 100000     # cycle limit
+tlim       = 15.0        # time limit
+integrator  = rk2       # time integration algorithm
+ncycle_out_mesh = -500  # print mesh structure every 500 cyles and on refinement
+
+<parthenon/mesh>
+nghost     = 3
+refinement = none
+nx1        = 64       # Number of zones in X1-direction
+x1min      = 0.0        # minimum value of X1
+x1max      = 1.0        # maximum value of X1
+ix1_bc     = periodic   # inner-X1 boundary flag
+ox1_bc     = periodic   # outer-X1 boundary flag
+
+nx2        = 64       # Number of zones in X2-direction
+x2min      = 0.0        # minimum value of X2
+x2max      = 1.0        # maximum value of X2
+ix2_bc     = periodic   # inner-X2 boundary flag
+ox2_bc     = periodic   # outer-X2 boundary flag
+
+nx3        = 64       # Number of zones in X3-direction
+x3min      = 0.0        # minimum value of X3
+x3max      = 1.0        # maximum value of X3
+ix3_bc     = periodic   # inner-X3 boundary flag
+ox3_bc     = periodic   # outer-X3 boundary flag
+
+pack_size = -1          # pack all blocks in a single pack (currently required)
+
+<parthenon/meshblock>
+nx1        = 32
+nx2        = 32
+nx3        = 32
+
+<hydro>
+fluid = euler
+eos = adiabatic
+riemann = hllc
+reconstruction = ppm
+gamma           = 1.66666         # gamma = C_p/C_v
+first_order_flux_correct = true
+
+<problem/turbulence>
+rho0         = 1.0      #4.0 # initial mean density
+p0           = 0.002101 #0.015625 # initial mean pressure
+b0           = 0.01     # initial magnetic field (x-direction)
+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   = 0.3      # solenoidal weight of the acceleration field
+accel_rms    = 0.45      # root mean square value of the acceleration field
+num_modes    = 30       # number of wavemodes
+
+<tracers>
+enabled = true
+initial_seed_method  = random_per_block    # alternative: user
+initial_num_tracers_per_cell = 0.125
+initial_rng_seed = 42   # optional
+
+<modes>
+k_1_0    = +2
+k_1_1    = -1
+k_1_2    = +0
+k_2_0    = +1
+k_2_1    = +0
+k_2_2    = +2
+k_3_0    = +1
+k_3_1    = +1
+k_3_2    = -1
+k_4_0    = +2
+k_4_1    = +0
+k_4_2    = +1
+k_5_0    = +0
+k_5_1    = +0
+k_5_2    = -1
+k_6_0    = +1
+k_6_1    = -1
+k_6_2    = -2
+k_7_0    = +0
+k_7_1    = +0
+k_7_2    = -2
+k_8_0    = +1
+k_8_1    = +0
+k_8_2    = -1
+k_9_0    = +0
+k_9_1    = +2
+k_9_2    = +1
+k_10_0    = +0
+k_10_1    = -1
+k_10_2    = +2
+k_11_0    = +0
+k_11_1    = +0
+k_11_2    = +2
+k_12_0    = +0
+k_12_1    = +2
+k_12_2    = -1
+k_13_0    = +2
+k_13_1    = +1
+k_13_2    = +1
+k_14_0    = +1
+k_14_1    = -1
+k_14_2    = +0
+k_15_0    = +0
+k_15_1    = -1
+k_15_2    = -1
+k_16_0    = +1
+k_16_1    = +0
+k_16_2    = +1
+k_17_0    = +0
+k_17_1    = -1
+k_17_2    = +1
+k_18_0    = +0
+k_18_1    = +1
+k_18_2    = +0
+k_19_0    = +1
+k_19_1    = -1
+k_19_2    = +1
+k_20_0    = +2
+k_20_1    = +1
+k_20_2    = -1
+k_21_0    = +0
+k_21_1    = -1
+k_21_2    = -2
+k_22_0    = +2
+k_22_1    = -1
+k_22_2    = +1
+k_23_0    = +0
+k_23_1    = +1
+k_23_2    = +1
+k_24_0    = +1
+k_24_1    = -2
+k_24_2    = +1
+k_25_0    = +1
+k_25_1    = -2
+k_25_2    = +0
+k_26_0    = +1
+k_26_1    = +2
+k_26_2    = +0
+k_27_0    = +1
+k_27_1    = -2
+k_27_2    = -1
+k_28_0    = +2
+k_28_1    = -1
+k_28_2    = -1
+k_29_0    = +1
+k_29_1    = +2
+k_29_2    = -1
+k_30_0    = +1
+k_30_1    = +1
+k_30_2    = +2

inputs/turb_with_tracers.in

@@ -1,6 +1,6 @@
 <comment>
 problem   = Turbulence with parabolic forcing profile (using few driving modes)
-reference = to be written, see https://gitlab.com/pgrete/kathena/wikis/turbulence for now
+reference = https://github.com/parthenon-hpc-lab/athenapk/blob/main/docs/turbulence.md
 
 <job>
 problem_id = turbulence

inputs/turbulence.in

@@ -0,0 +1,40 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from pathlib import Path
+import sys
+
+sys.path.insert(
+    1,
+    "./external/parthenon" + "/scripts/python/packages/parthenon_tools/parthenon_tools",
+)
+
+try:
+    import phdf
+except ModuleNotFoundError:
+    print("Couldn't find module to read Parthenon hdf5 files.")
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("filenames", nargs="*")
+    args = parser.parse_args()
+
+    for filename in args.filenames:
+        data = phdf.phdf(filename)
+        tracers = data.GetSwarm("tracers")
+        xs = tracers.x
+        ys = tracers.y
+        zs = tracers.z
+        ids = tracers.Get("id")
+
+        fig = plt.figure()
+        ax = fig.add_subplot(projection="3d")
+
+        ax.scatter(xs, ys, zs, s=1, c=ids)
+
+        filepath = Path(filename)
+        figfile = str(filepath.stem) + ".png"
+        print(f"saving to {figfile}...")
+        plt.savefig(figfile)
+        plt.close()

scripts/plot_particles.py

@@ -21,7 +21,10 @@ add_executable(
         hydro/srcterms/tabular_cooling.cpp
         refinement/gradient.cpp
         refinement/other.cpp
+        tracers/tracers.cpp
+        tracers/tracers.hpp
         utils/few_modes_ft.cpp
+        utils/few_modes_ft.hpp
 )
 
 add_subdirectory(pgen)

src/CMakeLists.txt

@@ -25,6 +25,7 @@
 #include "../recon/weno3_simple.hpp"
 #include "../recon/wenoz_simple.hpp"
 #include "../refinement/refinement.hpp"
+#include "../tracers/tracers.hpp"
 #include "../units.hpp"
 #include "defs.hpp"
 #include "diffusion/diffusion.hpp"
@@ -54,6 +55,7 @@ using parthenon::HistoryOutputVar;
 parthenon::Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   parthenon::Packages_t packages;
   packages.Add(Hydro::Initialize(pin.get()));
+  packages.Add(Tracers::Initialize(pin.get()));
   return packages;
 }