Two-body Interactions (QED)

dev/scripts/format.sh

@@ -3,43 +3,43 @@
 verify=false
 
 for arg in "$@"; do
-	case $arg in
-	--verify) verify=true ;;
-	esac
+  case $arg in
+  --verify) verify=true ;;
+  esac
 done
 
 if $verify; then
-	diff_output=""
+  diff_output=""
 
-	if command -v cmake-format &>/dev/null; then
-		while IFS= read -r -d '' f; do
-			if ! diff -q <(cmake-format "$f") "$f" &>/dev/null; then
-				diff_output+="  $f\n"
-			fi
-		done < <(find cmake/ src/ minimal/ tests/ -type f \( -name "*.cmake" -o -name "*.txt" \) -print0)
-	fi
+  if command -v cmake-format &>/dev/null; then
+    while IFS= read -r -d '' f; do
+      if ! diff -q <(cmake-format "$f") "$f" &>/dev/null; then
+        diff_output+="  $f\n"
+      fi
+    done < <(find cmake/ src/ minimal/ tests/ -type f \( -name "*.cmake" -o -name "*.txt" \) -print0)
+  fi
 
-	if command -v clang-format &>/dev/null; then
-		while IFS= read -r -d '' f; do
-			if ! clang-format --style=file --dry-run --Werror "$f" &>/dev/null; then
-				diff_output+="  $f\n"
-			fi
-		done < <(find pgens/ src/ minimal/ tests/ -type f \( -name "*.cpp" -o -name "*.hpp" -o -name "*.h" \) -print0)
-	fi
+  if command -v clang-format &>/dev/null; then
+    while IFS= read -r -d '' f; do
+      if ! clang-format --style=file --dry-run --Werror "$f" &>/dev/null; then
+        diff_output+="  $f\n"
+      fi
+    done < <(find pgens/ examples/ tutorials/ src/ minimal/ tests/ -type f \( -name "*.cpp" -o -name "*.hpp" -o -name "*.h" \) -print0)
+  fi
 
-	if [ -n "$diff_output" ]; then
-		echo "Formatting check failed. The following files need formatting:"
-		printf "$diff_output"
-		exit 1
-	else
-		echo "All files are properly formatted."
-	fi
+  if [ -n "$diff_output" ]; then
+    echo "Formatting check failed. The following files need formatting:"
+    printf '%s' "$diff_output"
+    exit 1
+  else
+    echo "All files are properly formatted."
+  fi
 else
-	if command -v cmake-format &>/dev/null; then
-		find cmake/ src/ minimal/ tests/ -type f -name "*.cmake" -o -name "*.txt" | xargs cmake-format -i
-	fi
+  if command -v cmake-format &>/dev/null; then
+    find cmake/ src/ minimal/ tests/ \( -type f -name "*.cmake" -o -name "*.txt" \) -exec cmake-format -i {} \;
+  fi
 
-	if command -v clang-format &>/dev/null; then
-		find pgens/ src/ minimal/ tests/ -type f -name "*.cpp" -o -name "*.hpp" -o -name "*.h" | xargs clang-format --style=file -i
-	fi
+  if command -v clang-format &>/dev/null; then
+    find pgens/ src/ minimal/ tests/ examples/ pgens/ tutorials/ \( -type f -name "*.cpp" -o -name "*.hpp" -o -name "*.h" \) -exec clang-format --style=file -i {} \;
+  fi
 fi

examples/compton_jones/compton_jones.py

@@ -0,0 +1,63 @@
+import nt2
+import matplotlib.pyplot as plt
+import numpy as np
+
+data = nt2.Data("compton_jones")
+
+photons = data.particles.sel(sp=3).isel(t=-1).load()
+photons = photons[np.sqrt(photons.ux**2 + photons.uy**2 + photons.uz**2) > 0.01]
+
+plt.rcParams["figure.dpi"] = 300
+plt.rcParams["font.family"] = "serif"
+plt.rcParams["mathtext.fontset"] = "stix"
+
+fig = plt.figure(figsize=(9, 4))
+gs = fig.add_gridspec(1, 2, wspace=0.35)
+ax1 = fig.add_subplot(gs[0, 0])
+ax2 = fig.add_subplot(gs[0, 1])
+
+gamma = np.sqrt(1 + data.attrs["setup.electron_4vel"] ** 2)
+e0 = data.attrs["setup.photon_energy"]
+Gamma = 4 * e0 * gamma
+emax = gamma * Gamma / (1 + Gamma)
+
+es = data.spectra.E[1:-1] / emax
+
+dnde = data.spectra.N_3.isel(t=-1)[1:-1]
+dnde /= np.trapezoid(dnde, es)
+ax1.plot(es, dnde)
+
+es = np.linspace(es.values.min(), es.values.max(), 250)
+qs = es / (1 + Gamma * (1 - es))
+dnde_th = (
+    2 * qs * np.log(qs)
+    + (1 + 2 * qs) * (1 - qs)
+    + 0.5 * Gamma**2 * qs**2 / (1 + Gamma * qs) * (1 - qs)
+)
+
+dnde_th /= np.trapezoid(dnde_th, es)
+
+ax1.plot(es, dnde_th, c="k", ls=":")
+ax1.set(
+    xlim=(0, 1),
+    ylim=(0, 4),
+    xlabel=r"$\varepsilon_{\rm ph} / \varepsilon_{\rm max}$",
+    ylabel=r"$dn_{\rm ph}/d\varepsilon_{\rm ph}$",
+)
+
+plt.scatter(
+    photons.ux / emax,
+    photons.uy / emax,
+    s=1,
+    linewidth=0,
+)
+xs = np.linspace(0, 1, 100)
+ys = 2 / gamma * xs
+ax2.plot(xs, ys, c="k", ls="--", lw=0.5)
+ax2.plot(xs, -ys, c="k", ls="--", lw=0.5)
+ax2.set(
+    xlabel=r"$p_{\rm ph}^x / \varepsilon_{\rm max}$",
+    ylabel=r"$p_{\rm ph}^y / \varepsilon_{\rm max}$",
+)
+
+plt.savefig("compton_jones.png", bbox_inches="tight")

examples/compton_jones/compton_jones.toml

@@ -0,0 +1,89 @@
+[simulation]
+  name    = "compton_jones"
+  engine  = "srpic"
+  runtime = 10.0
+
+[grid]
+  resolution = [32, 32]
+  extent     = [[0.0, 1.0], [0.0, 1.0]]
+
+  [grid.metric]
+    metric = "minkowski"
+
+  [grid.boundaries]
+    fields    = [["PERIODIC"], ["PERIODIC"]]
+    particles = [["PERIODIC"], ["PERIODIC"]]
+
+[scales]
+  larmor0    = 1.0
+  skindepth0 = 1.0
+
+[two_body]
+  thomson_optical_depth = 100.0
+
+  [[two_body.interaction]]
+    type      = "compton"
+    group1    = [1]
+    group2    = [2]
+    interval  = 1
+    tile_size = 5
+    recoil1   = false
+    recoil2   = true
+
+[algorithms]
+  current_filters = 0
+
+  [algorithms.deposit]
+    enable = false
+
+  [algorithms.fieldsolver]
+    enable = false
+
+[particles]
+  ppc0           = 5.0
+  clear_interval = 1
+
+  [[particles.species]]
+    label    = "e-"
+    mass     = 1.0
+    charge   = -1.0
+    maxnpart = 1e6
+
+  [[particles.species]]
+    label    = "ph"
+    mass     = 0.0
+    charge   = 0.0
+    maxnpart = 1e6
+
+  [[particles.species]]
+    label    = "ph_out"
+    mass     = 0.0
+    charge   = 0.0
+    maxnpart = 1e7
+    pusher   = "none"
+
+[setup]
+  electron_4vel = 999.9995
+  photon_energy = 1e-2
+
+[output]
+  interval_time = 9.0
+
+  [output.fields]
+    quantities = ["N_1", "N_2", "N_3"]
+
+  [output.particles]
+    species = [1, 2, 3]
+    stride  = 1
+
+  [output.spectra]
+    log_bins = false
+    e_min    = 0
+    e_max    = 1100
+    n_bins   = 100
+
+  [output.stats]
+    quantities = ["T00_1", "T00_2", "T00_3"]
+
+[checkpoint]
+  keep = 0

examples/compton_jones/pgen.hpp

@@ -0,0 +1,148 @@
+#ifndef PROBLEM_GENERATOR_H
+#define PROBLEM_GENERATOR_H
+
+#include "enums.h"
+#include "global.h"
+
+#include "arch/kokkos_aliases.h"
+#include "traits/pgen.h"
+
+#include "archetypes/particle_injector.h"
+#include "framework/domain/metadomain.h"
+
+namespace user {
+  using namespace ntt;
+
+  template <Dimension D>
+  struct DeltaDistribution {
+    const real_t         energy0;
+    bool                 monodirectional;
+    random_number_pool_t random_pool;
+
+    DeltaDistribution(real_t                energy0,
+                      bool                  monodirectional,
+                      random_number_pool_t& random_pool)
+      : energy0 { energy0 }
+      , monodirectional { monodirectional }
+      , random_pool { random_pool } {}
+
+    Inline void operator()(const coord_t<D>&, vec_t<Dim::_3D>& v) const {
+      if (not monodirectional) {
+        auto gen  = random_pool.get_state();
+        auto rnd1 = Random<real_t>(gen);
+        auto rnd2 = Random<real_t>(gen);
+        random_pool.free_state(gen);
+        // random direction
+        const auto phi = static_cast<real_t>(constant::TWO_PI) * rnd1;
+        const auto ct  = 2.0 * rnd2 - 1.0;
+        const auto st  = math::sqrt(1.0 - ct * ct);
+        v[0]           = energy0 * st * math::cos(phi);
+        v[1]           = energy0 * st * math::sin(phi);
+        v[2]           = energy0 * ct;
+      } else {
+        v[0] = energy0;
+        v[1] = 0.0;
+        v[2] = 0.0;
+      }
+    }
+  };
+
+  template <SimEngine::type S, class M>
+  struct PGen {
+
+    static constexpr auto engines {
+      ::traits::pgen::compatible_with<SimEngine::SRPIC> {}
+    };
+    static constexpr auto metrics {
+      ::traits::pgen::compatible_with<Metric::Minkowski> {}
+    };
+    static constexpr auto dimensions {
+      ::traits::pgen::compatible_with<Dim::_1D, Dim::_2D, Dim::_3D> {}
+    };
+
+    const SimulationParams& params;
+    const Metadomain<S, M>& metadomain;
+
+    PGen(const SimulationParams& p, const Metadomain<S, M>& m)
+      : params { p }
+      , metadomain { m } {}
+
+    void InitPrtls(Domain<S, M>& domain) {
+      auto delta_electrons = DeltaDistribution<M::Dim> {
+        params.template get<real_t>("setup.electron_4vel"),
+        true,
+        domain.random_pool()
+      };
+      arch::InjectUniform<S, M, decltype(delta_electrons)>(params,
+                                                           domain,
+                                                           1u,
+                                                           delta_electrons,
+                                                           ONE);
+    }
+
+    void CustomPostStep(timestep_t /*step*/, simtime_t /*time*/, Domain<S, M>& domain) {
+      // copy all photons from species #2 (idx 1) to #3 (idx 2) with an offset
+      const auto offset     = domain.species[2].npart();
+      const auto new_copies = domain.species[1].npart();
+      const auto new_size   = offset + new_copies;
+      const auto from_slice = prtl_slice_t { 0, new_copies };
+      const auto to_slice   = prtl_slice_t { offset, new_size };
+
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].i1, to_slice),
+                        Kokkos::subview(domain.species[1].i1, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].i1_prev, to_slice),
+                        Kokkos::subview(domain.species[1].i1_prev, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].dx1, to_slice),
+                        Kokkos::subview(domain.species[1].dx1, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].dx1_prev, to_slice),
+                        Kokkos::subview(domain.species[1].dx1_prev, from_slice));
+      if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].i2, to_slice),
+                          Kokkos::subview(domain.species[1].i2, from_slice));
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].i2_prev, to_slice),
+                          Kokkos::subview(domain.species[1].i2_prev, from_slice));
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].dx2, to_slice),
+                          Kokkos::subview(domain.species[1].dx2, from_slice));
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].dx2_prev, to_slice),
+                          Kokkos::subview(domain.species[1].dx2_prev, from_slice));
+      }
+      if constexpr (M::Dim == Dim::_3D) {
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].i3, to_slice),
+                          Kokkos::subview(domain.species[1].i3, from_slice));
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].i3_prev, to_slice),
+                          Kokkos::subview(domain.species[1].i3_prev, from_slice));
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].dx3, to_slice),
+                          Kokkos::subview(domain.species[1].dx3, from_slice));
+        Kokkos::deep_copy(Kokkos::subview(domain.species[2].dx3_prev, to_slice),
+                          Kokkos::subview(domain.species[1].dx3_prev, from_slice));
+      }
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].ux1, to_slice),
+                        Kokkos::subview(domain.species[1].ux1, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].ux2, to_slice),
+                        Kokkos::subview(domain.species[1].ux2, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].ux3, to_slice),
+                        Kokkos::subview(domain.species[1].ux3, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].weight, to_slice),
+                        Kokkos::subview(domain.species[1].weight, from_slice));
+      Kokkos::deep_copy(Kokkos::subview(domain.species[2].tag, to_slice),
+                        Kokkos::subview(domain.species[1].tag, from_slice));
+
+      domain.species[1].set_npart(0);
+      domain.species[2].set_npart(new_size);
+
+      auto delta_photons = DeltaDistribution<M::Dim> {
+        params.template get<real_t>("setup.photon_energy"),
+        false,
+        domain.random_pool()
+      };
+      arch::InjectUniform<S, M, decltype(delta_photons)>(params,
+                                                         domain,
+                                                         2u,
+                                                         delta_photons,
+                                                         ONE);
+    }
+  };
+
+} // namespace user
+
+#endif