diff --git a/.gitignore b/.gitignore index cd508a596..85331517d 100644 --- a/.gitignore +++ b/.gitignore @@ -65,6 +65,7 @@ ignore-* tombi/ tidy/ .claude +.understand-anything # devenv .devenv* diff --git a/CODEGUIDE.md b/CODEGUIDE.md index d69cf03b2..4e95ca1ee 100644 --- a/CODEGUIDE.md +++ b/CODEGUIDE.md @@ -16,7 +16,7 @@ entity │ ├── styling.cmake # styling functions │ └── tests.cmake # root cmake for tests ├── dev # developer-specific tools -│ ├── nix # nix-shells +│ ├── nix # nix-shell & devenv environments │ ├── runners # dockerfiles for github runners on different architectures │ ├── scripts # developer-specific scripts │ ├── Dockerfile.common # parent docker environment for development diff --git a/dev/nix/adios2.nix b/dev/nix/adios2.nix index f7114a017..d3c0011ac 100644 --- a/dev/nix/adios2.nix +++ b/dev/nix/adios2.nix @@ -39,6 +39,7 @@ stdenv.mkDerivation { ]; propagatedBuildInputs = [ + pkgs.libgcc pkgs.gcc15 ] ++ ( diff --git a/dev/nix/devenv.lock b/dev/nix/devenv.lock index 9e6f7f3b4..9c8f8637c 100644 --- a/dev/nix/devenv.lock +++ b/dev/nix/devenv.lock @@ -42,4 +42,4 @@ }, "root": "root", "version": 7 -} \ No newline at end of file +} diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix index 2af11c164..ed428ab82 100644 --- a/dev/nix/kokkos.nix +++ b/dev/nix/kokkos.nix @@ -27,6 +27,7 @@ let ]; "NONE" = [ pkgs.clang-tools + pkgs.libgcc pkgs.gcc15 ]; }; diff --git a/dev/scripts/format.sh b/dev/scripts/format.sh index e12bef2ec..ca4aab0f4 100755 --- a/dev/scripts/format.sh +++ b/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 diff --git a/entity.schema.json b/entity.schema.json index 6ca79d81d..2900c5f90 100644 --- a/entity.schema.json +++ b/entity.schema.json @@ -3,7 +3,7 @@ "$id": "https://entity-toolkit.github.io/schema/entity.schema.json", "title": "entity", "description": "Entity simulation input file", - "$comment": "GENERATOR CONVENTION -- this schema is the single source of truth for `input.template.toml`. Standard keywords carry the machine-checkable part (type/enum/minimum/maximum/items/prefixItems/minItems/maxItems/pattern/default/required/deprecated). The `x-entity` object carries the parts JSON Schema cannot express, all verbatim from the template comments: `type` = the literal `@type:` annotation (use it for the comment line whenever present, else derive from the standard keywords); `default` = the literal `@default:` annotation (use it whenever present, else format the standard `default`); `notes` = ordered `@note:` lines; `examples` = ordered `@example:` lines; `enum` = an illustrative, NON-exhaustive value list that must not be validated as a real `enum`; `deprecated` = the `@deprecated:` text. `x-entity.inferred` on an object lists quantities the code derives rather than reads; they are NOT valid input keys (so they are absent from `properties` and rejected by `additionalProperties: false`) and should be emitted as an `@inferred:` comment block after that table's own keys and before its sub-tables. Property order in `properties` is the emission order. Objects with `additionalProperties: true` are free-form. `$ref` is used once, for `$defs/colormap`.", + "$comment": "GENERATOR CONVENTION -- this schema is the single source of truth for `input.default.toml`. Standard keywords carry the machine-checkable part (type/enum/minimum/maximum/items/prefixItems/minItems/maxItems/pattern/default/required/deprecated). The `x-entity` object carries the parts JSON Schema cannot express, all verbatim from the template comments: `type` = the literal `@type:` annotation (use it for the comment line whenever present, else derive from the standard keywords); `default` = the literal `@default:` annotation (use it whenever present, else format the standard `default`); `notes` = ordered `@note:` lines; `examples` = ordered `@example:` lines; `enum` = an illustrative, NON-exhaustive value list that must not be validated as a real `enum`; `deprecated` = the `@deprecated:` text. `x-entity.inferred` on an object lists quantities the code derives rather than reads; they are NOT valid input keys (so they are absent from `properties` and rejected by `additionalProperties: false`) and should be emitted as an `@inferred:` comment block after that table's own keys and before its sub-tables. Property order in `properties` is the emission order. Objects with `additionalProperties: true` are free-form. `$ref` is used once, for `$defs/colormap`.", "type": "object", "additionalProperties": false, "required": [ @@ -844,6 +844,100 @@ } } }, + "two_body": { + "description": "Parameters for two-body interactions", + "type": "object", + "additionalProperties": false, + "properties": { + "thomson_optical_depth": { + "description": "Nominal Thomson optical depth: `tau = n0 * sigma_T * 1` over a distance of 1 in physical units (n0 = nominal density)", + "type": "number", + "exclusiveMinimum": 0.0, + "default": 1.0 + }, + "interaction": { + "description": "Specific interaction parameters", + "type": "array", + "minItems": 1, + "x-entity": { + "array_of_tables": true + }, + "items": { + "type": "object", + "additionalProperties": false, + "required": [ + "type", + "group1" + ], + "properties": { + "type": { + "description": "Type of the two-body interaction", + "anyOf": [ + { + "enum": [ + "Compton" + ] + }, + { + "type": "string", + "pattern": "(?i)^(Compton)$" + } + ], + "x-entity": { + "type": "string" + } + }, + "group1": { + "description": "First group of species indices participating in the interaction", + "type": "array", + "minItems": 1, + "items": { + "type": "integer", + "minimum": 1 + } + }, + "group2": { + "description": "Second group of species indices participating in the interaction", + "type": "array", + "minItems": 0, + "default": [], + "x-entity": { + "notes": [ + "For interactions between particles of the same group, leave `group2` empty" + ] + }, + "items": { + "type": "integer", + "minimum": 1 + } + }, + "interval": { + "description": "Interval in timesteps between checking for the interaction", + "type": "integer", + "minimum": 1, + "default": 1 + }, + "tile_size": { + "description": "Size of interaction tile in cells", + "type": "integer", + "minimum": 1, + "default": 4 + }, + "recoil1": { + "description": "Whether to apply recoil to the first group of particles", + "type": "boolean", + "default": true + }, + "recoil2": { + "description": "Whether to apply recoil to the second group of particles", + "type": "boolean", + "default": true + } + } + } + } + } + }, "algorithms": { "description": "Algorithm and solver tuning", "type": "object", @@ -1650,7 +1744,7 @@ "e_min": { "description": "Minimum energy for the spectra output", "type": "number", - "exclusiveMinimum": 0.0, + "minimum": 0.0, "default": 0.001, "x-entity": { "type": "float", diff --git a/examples/compton_jones/compton_jones.py b/examples/compton_jones/compton_jones.py new file mode 100644 index 000000000..505b777bf --- /dev/null +++ b/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") diff --git a/examples/compton_jones/compton_jones.toml b/examples/compton_jones/compton_jones.toml new file mode 100644 index 000000000..8230a04c5 --- /dev/null +++ b/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 diff --git a/examples/compton_jones/pgen.hpp b/examples/compton_jones/pgen.hpp new file mode 100644 index 000000000..f2da19802 --- /dev/null +++ b/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 + 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&, vec_t& v) const { + if (not monodirectional) { + auto gen = random_pool.get_state(); + auto rnd1 = Random(gen); + auto rnd2 = Random(gen); + random_pool.free_state(gen); + // random direction + const auto phi = static_cast(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 + struct PGen { + + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + + const SimulationParams& params; + const Metadomain& metadomain; + + PGen(const SimulationParams& p, const Metadomain& m) + : params { p } + , metadomain { m } {} + + void InitPrtls(Domain& domain) { + auto delta_electrons = DeltaDistribution { + params.template get("setup.electron_4vel"), + true, + domain.random_pool() + }; + arch::InjectUniform(params, + domain, + 1u, + delta_electrons, + ONE); + } + + void CustomPostStep(timestep_t /*step*/, simtime_t /*time*/, Domain& 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 { + params.template get("setup.photon_energy"), + false, + domain.random_pool() + }; + arch::InjectUniform(params, + domain, + 2u, + delta_photons, + ONE); + } + }; + +} // namespace user + +#endif diff --git a/examples/compton_kompaneets/compton_kompaneets.py b/examples/compton_kompaneets/compton_kompaneets.py new file mode 100644 index 000000000..b507964c5 --- /dev/null +++ b/examples/compton_kompaneets/compton_kompaneets.py @@ -0,0 +1,66 @@ +import nt2 +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +data = nt2.Data("compton_kompaneets") +stats = pd.read_csv("compton_kompaneets/compton_kompaneets_stats.csv") +stats.columns = stats.columns.str.strip() + +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.3) +ax1 = fig.add_subplot(gs[0, 0]) + +tC = 1 / ( + data.attrs["setup.temperature"] * data.attrs["two_body.thomson_optical_depth"] +) + +tvals = len(data.spectra.t.values) + +nphot = data.spectra.N_3.isel(t=-1).sum().values[()] +for ti, t in enumerate([0, 0.5 * tC, tC, 3 * tC]): + ax1.plot( + data.spectra.E.coarsen(E=2).mean().values / data.attrs["setup.temperature"], + data.spectra.N_3.sel(t=t, method="nearest").coarsen(E=2).mean().values, + c=plt.get_cmap("viridis")(ti / 3), + label=f"$y={t / tC:.1f}$", + lw=1, + ) + +es = data.spectra.E.values / data.attrs["setup.temperature"] +dndes = es**2 * np.exp(-es) +dndes /= np.sum(dndes) +dndes *= nphot +ax1.plot( + es, + dndes, + c="k", + ls=":", + label=r"$\propto \varepsilon_{\rm ph}^2 e^{-\varepsilon_{\rm ph} / T_\pm}$", +) + +ax1.set( + yscale="log", + ylim=(1, 1e5), + xlim=(0, 10), + xlabel=r"$\varepsilon / T_\pm$", + ylabel=r"$dn_{\rm ph}/d\varepsilon$", +) +ax1.legend() + +ax2 = fig.add_subplot(gs[0, 1]) +ax2.plot(stats["time"] / tC, stats["T00_3"], c="C0") +ax2.set(ylabel=r"total photon energy", xlabel=r"$y\equiv t/t_C$") +ax2.yaxis.label.set_color("C0") +ax2.tick_params(axis="y", labelcolor="C0") +ax2twin = ax2.twinx() +ax2twin.plot(data.spectra.t.values / tC, data.spectra.N_3.sum("E"), c="C2") +ax2twin.set(ylabel=r"photon number") +ax2twin.yaxis.label.set_color("C2") +ax2twin.tick_params(axis="y", labelcolor="C2") + +plt.savefig("compton_kompaneets_plot.png", bbox_inches="tight") diff --git a/examples/compton_kompaneets/compton_kompaneets.toml b/examples/compton_kompaneets/compton_kompaneets.toml new file mode 100644 index 000000000..68631ed0f --- /dev/null +++ b/examples/compton_kompaneets/compton_kompaneets.toml @@ -0,0 +1,87 @@ +[simulation] + name = "compton_kompaneets" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [128, 128] + 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 = 60.0 + + [[two_body.interaction]] + type = "compton" + group1 = [1, 2] + group2 = [3] + interval = 1 + tile_size = 5 + recoil1 = false + recoil2 = true + +[algorithms] + current_filters = 0 + + [algorithms.deposit] + enable = false + + [algorithms.fieldsolver] + enable = false + +[particles] + ppc0 = 2.0 + clear_interval = 1 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e6 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e6 + + [[particles.species]] + label = "ph" + mass = 0.0 + charge = 0.0 + maxnpart = 1e6 + +[setup] + temperature = 0.01 + photon_energy = 1e-3 + +[output] + interval_time = 0.1 + + [output.fields] + enable = false + + [output.particles] + enable = false + + [output.spectra] + log_bins = false + e_min = 0.0 + e_max = 1.0 + num_energy_bins = 500 + + [output.stats] + quantities = ["T00_3"] + +[checkpoint] + keep = 0 diff --git a/examples/compton_kompaneets/pgen.hpp b/examples/compton_kompaneets/pgen.hpp new file mode 100644 index 000000000..0b0158ada --- /dev/null +++ b/examples/compton_kompaneets/pgen.hpp @@ -0,0 +1,74 @@ +#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 "archetypes/utils.h" +#include "framework/domain/metadomain.h" + +namespace user { + using namespace ntt; + + template + struct DeltaDistribution { + const real_t photon_energy0; + random_number_pool_t random_pool; + + DeltaDistribution(real_t photon_energy0, random_number_pool_t& random_pool) + : photon_energy0 { photon_energy0 } + , random_pool { random_pool } {} + + Inline void operator()(const coord_t&, vec_t& v) const { + auto gen = random_pool.get_state(); + auto rnd1 = Random(gen); + auto rnd2 = Random(gen); + random_pool.free_state(gen); + // random direction + const auto phi = static_cast(constant::TWO_PI) * rnd1; + const auto ct = 2.0 * rnd2 - 1.0; + const auto st = math::sqrt(1.0 - ct * ct); + v[0] = photon_energy0 * st * math::cos(phi); + v[1] = photon_energy0 * st * math::sin(phi); + v[2] = photon_energy0 * ct; + } + }; + + template + struct PGen { + + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + + const SimulationParams& params; + const Metadomain& metadomain; + + PGen(const SimulationParams& p, const Metadomain& m) + : params { p } + , metadomain { m } {} + + void InitPrtls(Domain& domain) { + const auto temperature = params.template get("setup.temperature"); + arch::InjectUniformMaxwellian(params, domain, ONE, temperature, { 1u, 2u }); + + auto delta = DeltaDistribution { params.template get( + "setup.photon_energy"), + domain.random_pool() }; + arch::InjectUniform(params, domain, 3u, delta, ONE); + } + }; + +} // namespace user + +#endif diff --git a/input.default.toml b/input.default.toml index 17daf3819..b9e060543 100644 --- a/input.default.toml +++ b/input.default.toml @@ -4,12 +4,12 @@ # @required # @type: string # @note: The name is used for the output files - name = "" + name = "" # Simulation engine to use # @required # @type: string # @enum: "SRPIC", "GRPIC" - engine = "SRPIC" + engine = "SRPIC" # Max runtime in physical (code) units # @required # @type: float [> 0] @@ -21,7 +21,7 @@ # Number of domains # @type: int # @default: 1 [no MPI]; MPI_SIZE [MPI] - number = 1 + number = 1 # Decomposition of the domain (for MPI) in each of the directions # @type: array [size 1 :->: 3] # @default: [-1, -1, -1] @@ -40,11 +40,11 @@ # Enable dynamic load balancing # @type: bool # @default: false - enable = false + enable = false # Run the rebalancer every `interval` timesteps (0 disables) # @type: int # @default: 0 - interval = 0 + interval = 0 # Dimensions along which load is redistributed (1 = x1, 2 = x2, 3 = x3) # @type: array [subset of {1, 2, 3}] # @default: [1] @@ -54,13 +54,13 @@ # particle count is below this fraction # @type: float # @default: 0.1 - tolerance = 0.1 + tolerance = 0.1 # Maximum cell-shift per interior boundary per event; clamped at compile # time to N_GHOSTS so the migrating field strip is already cached in the # rank's ghost zone. # @type: int # @default: N_GHOSTS - max_shift = 0 + max_shift = 0 # Parameters specific to grid geometry [grid] @@ -77,7 +77,7 @@ # are set automatically # @note: For cartesian geometry, cell aspect ratio has to be 1: `dx=dy=dz` # @example: [[0.0, 1.0], [-1.0, 1.0]] - extent = [[0.0, 0.0]] + extent = [[0.0, 0.0]] # @inferred: # - dim @@ -93,7 +93,7 @@ # @type: string # @enum: "Minkowski", "Spherical", "QSpherical", "Kerr_Schild", # "QKerr_Schild", "Kerr_Schild_0" - metric = "Minkowski" + metric = "Minkowski" # `r0` paramter for the QSpherical metric `x1 = log(r-r0)` # @type: float [-inf -> rmin] # @default: 0.0 @@ -103,11 +103,11 @@ # (pi-2*x2)*(pi-x2)/pi^2` # @type: float [-1 :->: 1] # @default: 0.0 - qsph_h = 0.0 + qsph_h = 0.0 # Spin parameter for the Kerr Schild metric # @type: float [0 :-> 1] # @default: 0.0 - ks_a = 0.0 + ks_a = 0.0 # @inferred: # - coord @@ -139,7 +139,7 @@ # @note: In GR, the horizon boundary is set automatically (only specify bc # @ rmax): [["MATCH"]] # @example: [["CUSTOM", "MATCH"]] (for 2D spherical `[[rmin, rmax]]`) - fields = [["PERIODIC"]] + fields = [["PERIODIC"]] # Boundary conditions for particles # @required # @type: array> [size 1 :->: 3] @@ -187,18 +187,18 @@ temperature = 0.0 # Peak number density of the atmosphere at base in units of `n0` # @type: float - density = 0.0 + density = 0.0 # Pressure scale-height in physical units # @type: float - height = 1.0 + height = 1.0 # Species indices of particles that populate the atmosphere # @type: array [size 2] - species = [1, 1] + species = [1, 1] # Distance from the edge to which the gravity is imposed in physical units # @type: float # @default: 0.0 # @note: 0.0 means no limit - ds = 0.0 + ds = 0.0 # @inferred: # - g @@ -213,7 +213,7 @@ # Fiducial larmor radius # @required # @type: float [> 0.0] - larmor0 = 1.0 + larmor0 = 1.0 # Fiducial plasma skin depth # @required # @type: float [> 0.0] @@ -287,7 +287,7 @@ # c^2` in fiducial magnetic field `B0` # @type: float [> 1.0] # @default: 10.0 - gamma_qed = 10.0 + gamma_qed = 10.0 # Minimum photon energy for synchrotron emission (units of `m0 c^2`) # @type: float [> 0.0] # @default: 1e-3 @@ -295,11 +295,11 @@ # Weights for the emitted synchrotron photons # @type: float [> 0.0] # @default: 1.0 - photon_weight = 1.0 + photon_weight = 1.0 # Index of species for the emitted photon # @required # @type: ushort [> 0] - photon_species = 1 + photon_species = 1 # @inferred: # - nominal_probability @@ -324,7 +324,7 @@ # `m0 c^2` in fiducial magnetic field `B0` # @type: float [> 1.0] # @default: 10.0 - gamma_qed = 10.0 + gamma_qed = 10.0 # Minimum photon energy for inverse Compton emission (units of `m0 c^2`) # @type: float [> 0.0] # @default: 1e-3 @@ -332,11 +332,11 @@ # Weights for the emitted inverse Compton photons # @type: float [> 0.0] # @default: 1.0 - photon_weight = 1.0 + photon_weight = 1.0 # Index of species for the emitted photon # @required # @type: ushort [> 0] - photon_species = 1 + photon_species = 1 # @inferred: # - nominal_probability @@ -354,6 +354,48 @@ # @from: `.gamma_qed` # @value: `(1 / gamma_qed)^2` +# Parameters for two-body interactions +[two_body] + # Nominal Thomson optical depth: `tau = n0 * sigma_T * 1` over a distance of 1 + # in physical units (n0 = nominal density) + # @type: number + # @default: 1.0 + thomson_optical_depth = 1.0 + + # Specific interaction parameters + [[two_body.interaction]] + # Type of the two-body interaction + # @required + # @type: string + # @enum: "Compton" + type = "Compton" + # First group of species indices participating in the interaction + # @required + # @type: array + group1 = [1] + # Second group of species indices participating in the interaction + # @type: array + # @default: [] + # @note: For interactions between particles of the same group, leave + # `group2` empty + group2 = [] + # Interval in timesteps between checking for the interaction + # @type: integer + # @default: 1 + interval = 1 + # Size of interaction tile in cells + # @type: integer + # @default: 4 + tile_size = 4 + # Whether to apply recoil to the first group of particles + # @type: boolean + # @default: true + recoil1 = true + # Whether to apply recoil to the second group of particles + # @type: boolean + # @default: true + recoil2 = true + # Algorithm and solver tuning [algorithms] # Number of current smoothing passes @@ -367,7 +409,7 @@ # @type: float [0.0 -> 1.0] # @default: 0.95 # @note: CFL number determines the timestep duration - CFL = 0.95 + CFL = 0.95 # Correction factor for the speed of light used in field solver # @type: float # @default: 1.0 @@ -385,12 +427,12 @@ # Enable the current deposition # @type: bool # @default: true - enable = true + enable = true # Tiled-deposit work-group (team) size # @type: uint [>= 0] # @default: 0 # @deprecated: removed in 1.6+, use `tiled_deposit_team_size` instead - team_policy_team_size = 0 + team_policy_team_size = 0 # Tiled-deposit work-group (team) size # @type: uint [>= 0] # @default: 0 @@ -412,7 +454,7 @@ # Stepsize for numerical differentiation in GR pusher # @type: float [> 0.0] # @default: 1e-6 - pusher_eps = 1e-6 + pusher_eps = 1e-6 # Number of iterations for the Newton-Raphson method in GR pusher # @type: ushort [> 0] # @default: 10 @@ -428,7 +470,7 @@ # @type: float # @default: 0.0 # @note: When `larmor_max` == 0, the limit is disabled - larmor_max = 0.0 + larmor_max = 0.0 # Stencil coefficients for the field solver [notation as in Blinne+ (2018)] # @note: Standard Yee solver: `delta_i = beta_ij = 0.0` @@ -436,7 +478,7 @@ # Enable the fieldsolver # @type: bool # @default: true - enable = true + enable = true # delta_x coefficient (for `F_{i +/- 3/2, j, k}`) # @type: float # @default: 0.0 @@ -487,16 +529,16 @@ # Fiducial number of particles per cell # @required # @type: float [> 0.0] - ppc0 = 1.0 + ppc0 = 1.0 # Toggle for using particle weights # @type: bool # @default: false - use_weights = false + use_weights = false # Timesteps between particle re-sorting by tags (removing dead particles) # @type: uint # @default: 100 # @note: Set to 0 to disable re-sorting - clear_interval = 100 + clear_interval = 100 # Timesteps between spatial sorting of particles (for better cache # performance) # @type: uint @@ -517,41 +559,41 @@ # @default: "s" # @note: `` is the index of the species in the list starting from 1 # @example: "e-" - label = "s" + label = "s" # Mass of the species (in units of fiducial mass) # @required # @type: float [>= 0.0] - mass = 0.0 + mass = 0.0 # Charge of the species (in units of fiducial charge) # @required # @type: float - charge = 0.0 + charge = 0.0 # Maximum number of particles per task # @required # @type: uint [> 0] # @note: Read as a float, so exponential notation is fine (e.g. `1e8`) - maxnpart = 1.0 + maxnpart = 1.0 # Pusher algorithm for the species # @type: string # @default: "Boris" [massive]; "Photon" [massless] # @enum: "Boris", "Vay", "Boris,GCA", "Vay,GCA", "Photon", "None" - pusher = "Boris" + pusher = "Boris" # Number of additional real-valued variables (payloads) for each particle of # the given species # @type: ushort # @default: 0 - n_payloads_real = 0 + n_payloads_real = 0 # Number of additional integer-valued variables (payloads) for each particle # of the given species # @type: ushort # @default: 0 # @note: If tracking is enabled, one or two extra integer payloads are # reserved (depending on whether MPI is enabled) - n_payloads_int = 0 + n_payloads_int = 0 # Enable tracking of particles using indices for the given species # @type: bool # @default: false - tracking = false + tracking = false # Radiation reaction to use for the species # @type: string # @default: "None" @@ -559,7 +601,7 @@ # @note: Can also be coma-separated combination, e.g., # "Synchrotron,Compton" # @note: Relevant radiation.drag parameters should also be provided - radiative_drag = "None" + radiative_drag = "None" # Particle emission policy for the species # @type: string # @default: "None" @@ -567,7 +609,7 @@ # @note: Only one emission mechanism allowed # @note: Appropriate radiation drag flag will be applied automatically # (unless explicitly set to "None") - emission = "None" + emission = "None" # Timesteps between spatial sorting of particles for given species # @type: uint # @default: 0 @@ -580,7 +622,7 @@ # @default: 100 # @note: Set to 0 to disable re-sorting # @note: Overrides `particles.clear_interval` for the given species - clear_interval = 100 + clear_interval = 100 # Parameters for specific problem generators and setups # @note: Free-form: keys are defined by the problem generator, so nothing here @@ -593,12 +635,12 @@ # @type: string # @default: "bpfile" # @enum: "disabled", "hdf5", "BPFile" - format = "bpfile" + format = "bpfile" # Number of timesteps between all outputs # @type: uint [> 0] # @default: 100 # @note: Value is overriden by output intervals for specific outputs - interval = 100 + interval = 100 # Physical (code) time interval between all outputs # @type: float # @default: -1.0 @@ -612,7 +654,7 @@ # Toggle for the field output # @type: bool # @default: true - enable = true + enable = true # Field quantities to output # @type: array # @default: [] @@ -624,17 +666,17 @@ # "3" # @note: By default, we accumulate moments from all massive species, one # can specify only specific species: `Ttt_1_2`, `Rho_1`, `Rho_3_4` - quantities = [] + quantities = [] # Custom (user-defined) field quantities # @type: array # @default: [] - custom = [] + custom = [] # Number of timesteps between field outputs # @type: uint # @default: 0 # @note: When `!= 0`, overrides `output.interval` # @note: When `== 0`, `output.interval` is used - interval = 0 + interval = 0 # Physical (code) time interval between field outputs # @type: float # @default: -1.0 @@ -646,14 +688,14 @@ # @default: [1, 1, 1] # @note: The output is downsampled by the given factors in each direction # @note: If a scalar is given, it is applied to all directions - downsampling = [1, 1, 1] + downsampling = [1, 1, 1] # Smoothing of the output moments [output.fields.smoothing] # Smoothing order for the output of moments ("Rho", "Charge", "T", ...) # @type: ushort # @default: 0 - order = 0 + order = 0 # Smoothing algorithm # @type: string # @default: "spline" @@ -669,22 +711,22 @@ # Toggle for the particles output # @type: bool # @default: true - enable = true + enable = true # Particle species indices to output # @type: array # @default: [] # @note: If empty, all species are output - species = [] + species = [] # Stride for the output of particles # @type: uint [>= 1] # @default: 100 - stride = 100 + stride = 100 # Number of timesteps between particle outputs # @type: uint # @default: 0 # @note: When `!= 0`, overrides `output.interval` # @note: When `== 0`, `output.interval` is used - interval = 0 + interval = 0 # Physical (code) time interval between particle outputs # @type: float # @default: -1.0 @@ -697,28 +739,28 @@ # Toggle for the spectra output # @type: bool # @default: true - enable = true + enable = true # Minimum energy for the spectra output # @type: float # @default: 1e-3 - e_min = 1e-3 + e_min = 1e-3 # Maximum energy for the spectra output # @type: float # @default: 1e3 - e_max = 1e3 + e_max = 1e3 # Whether to use logarithmic bins for energy # @type: bool # @default: true - log_bins = true + log_bins = true # Number of energy bins for the spectra output # @type: uint [> 0] # @default: 200 # @deprecated: removed in 1.6+, use `num_energy_bins` instead - n_bins = 200 + n_bins = 200 # Number of energy bins for the spectra output # @type: uint [> 0] # @default: 200 - num_energy_bins = 200 + num_energy_bins = 200 # Number of spatial bins for the spectra output # @type: array [size 1 :->: 3] # @default: [1, 1, 1] @@ -728,20 +770,20 @@ # @default: 0 # @note: When `!= 0`, overrides `output.interval` # @note: When `== 0`, `output.interval` is used - interval = 0 + interval = 0 # Physical (code) time interval between spectra outputs # @type: float # @default: -1.0 # @note: When `< 0`, the output is controlled by `interval` # @note: When specified, overrides `output.interval_time` - interval_time = -1.0 + interval_time = -1.0 # Debug output parameters [output.debug] # Output fields "as is" without conversions # @type: bool # @default: false - as_is = false + as_is = false # Output fields with values in ghost cells # @type: bool # @default: false @@ -752,12 +794,12 @@ # Toggle for the stats output # @type: bool # @default: true - enable = true + enable = true # Number of timesteps between stat outputs # @type: uint [> 0] # @default: 100 # @note: Overriden if `output.stats.interval_time != -1` - interval = 100 + interval = 100 # Physical (code) time interval between stat outputs # @type: float # @default: -1.0 @@ -770,18 +812,18 @@ # "Tij" # @note: For particle moments, ... # @note: ... same notation is used as for `output.fields.quantities` - quantities = ["B^2", "E^2", "ExB", "Rho", "T00"] + quantities = ["B^2", "E^2", "ExB", "Rho", "T00"] # Custom (user-defined) stats # @type: array # @default: [] - custom = [] + custom = [] # Checkpointing parameters [checkpoint] # Number of timesteps between checkpoints # @type: uint [> 0] # @default: 1000 - interval = 1000 + interval = 1000 # Physical (code) time interval between checkpoints # @type: float [> 0] # @default: -1.0 @@ -792,23 +834,23 @@ # @default: 2 # @note: 0 = disable checkpointing # @note: -1 = keep all checkpoints - keep = 2 + keep = 2 # Write a checkpoint once after a fixed walltime # @type: string # @default: "00:00:00" # @note: The format is "HH:MM:SS" # @note: Empty string or "00:00:00" disables this functionality # @note: Writing checkpoint at walltime does not stop the simulation - walltime = "00:00:00" + walltime = "00:00:00" # Parent directory to write checkpoints to # @type: string # @default: `.ckpt` # @note: The directory is created if it does not exist - write_path = "" + write_path = "" # Parent directory to use when resuming from a checkpoint # @type: string # @default: inherit `write_path` - read_path = "" + read_path = "" # @inferred: # - is_resuming @@ -836,12 +878,12 @@ # @type: uint # @default: 4294967296 # @note: Lower this on memory-constrained nodes; matches ADIOS2's default - max_shm_size = 4294967296 + max_shm_size = 4294967296 # Internal serialization buffer chunk size, in bytes (BP5 BufferChunkSize) # @type: uint # @default: 16777216 # @note: Scales with per-rank output volume; matches ADIOS2's default - buffer_chunk_size = 16777216 + buffer_chunk_size = 16777216 # In-situ renderer. Renders scalar fields on the GPU and writes PNG images # directly to `/renders/` each cadence -- no field data is written to @@ -861,45 +903,45 @@ # Toggle for the on-the-fly renderer # @type: bool # @default: false - enable = false + enable = false # Number of timesteps between renders # @type: uint # @default: 0 # @note: When `!= 0`, overrides `output.interval` # @note: When `== 0`, `interval_time` (or `output.interval`) is used - interval = 0 + interval = 0 # Physical (code) time interval between renders # @type: float # @default: -1.0 # @note: When `< 0`, the output is controlled by `interval` - interval_time = -1.0 + interval_time = -1.0 # Image width in pixels (the rendered region; the PNG is wider if a colorbar # margin is added, see `colorbar_outside`) # @type: int [> 0] # @default: 1024 - width = 1024 + width = 1024 # Image height in pixels # @type: int [> 0] # @default: 1024 - height = 1024 + height = 1024 # Convenience: force a square frame (sets width == height == resolution), the # natural shape for a dome master. Overrides `width`/`height` when > 0. # @type: int [> 0] # @default: 0 (use width/height) - resolution = 0 + resolution = 0 # Number of entries in the color/opacity lookup table # @type: int [> 1] # @default: 256 - n_lut = 256 + n_lut = 256 # Opaque background RGB (each channel 0..1) shown through # transparent/low-opacity pixels; also fills the colorbar margin # @type: array [size 3] # @default: [0.0, 0.0, 0.0] - background = [0.0, 0.0, 0.0] + background = [0.0, 0.0, 0.0] # Draw a colorbar (gradient + value ticks + label) on each PNG # @type: bool # @default: true - colorbar = true + colorbar = true # Draw the colorbar in an added right margin (the PNG becomes wider by a fixed # strip) instead of overlaying it on the rendered volume # @type: bool @@ -910,13 +952,13 @@ # Cartesian or 3D rendering. # @type: bool # @default: true - mirror = true + mirror = true # Draw the current simulation time as a label ("T = ", fixed to 2 # decimals) in the upper-right corner of the render region, in a contrasting # color, vertically centered between the frame top and the colorbar. # @type: bool # @default: false - time_label = false + time_label = false # Draw a spine (frame) + axis ticks + labels around the rendered region. The # PNG gains left/bottom margins (background-filled) for the tick labels and # axis names, so they never overlap the data. @@ -928,23 +970,23 @@ # spine); # 3D = the global box projected to a wireframe with ticks on the # three silhouette edges (x bottom, y & z on the left) - axes = false + axes = false # Axis names. 3D uses all three; the 2D slice uses the first two. When unset, # the 2D slice defaults to "x","y" (Cartesian) or "X","Z" (spherical). # @type: array [size <= 3] # @default: ["x", "y", "z"] - axis_labels = ["x", "y", "z"] + axis_labels = ["x", "y", "z"] # Target number of ticks per axis (actual count is rounded to nice values) # @type: int [>= 2] # @default: 5 - axis_ticks = 5 + axis_ticks = 5 # 3D only: target width (pixels) of the box wireframe "spine". The spine is # drawn inside the ray-march (opaque, depth-occluded by the volume); its width # is floored by the ray step, so for a crisper thin line raise `samples` as # well. # @type: float [> 0.0] # @default: 2.0 - spine_width = 2.0 + spine_width = 2.0 # Limit the render region to axis-aligned box in physical/world coordinates. # Left unset it spans the full domain. Clamped to the box. @@ -976,13 +1018,13 @@ # @default: 400 # @note: The world-space step is `box_diagonal / samples` unless # `step_size` is set. Higher = better quality, slower. - samples = 400 + samples = 400 # Fixed world-space step between ray samples # @type: float [>= 0.0] # @default: 0.0 # @note: 0 derives the step from `samples`. The step is identical on all # ranks, which is what makes the multi-domain composite seamless. - step_size = 0.0 + step_size = 0.0 # Stop marching a ray once its accumulated opacity reaches this value # @type: float [0.0 -> 1.0] # @default: 0.99 @@ -997,7 +1039,7 @@ # @type: array [size 2 or 3] # @default: [] (static view) # @example: velocity = [0.9, 0.0] # pan along +x1 at 0.9 c - velocity = [] + velocity = [] # Sim time at which the view starts moving (static before it, e.g. to let an # initial ramp-up finish) # @type: float @@ -1020,29 +1062,29 @@ # (unlike ortho/perspective, which need the eye outside the box). # Set a square frame (`resolution`, or width == height). `forward` # is the dome ZENITH (screen-up defaults to +y for a +z zenith). - mode = "orthographic" + mode = "orthographic" # Camera (eye) position in world (physical) coordinates # @type: array [size 3] # @default: box center pushed back ~1.7 box-diagonals along (1, 1, 1); # for `mode = "dome"`, the domain center (interior eye) - position = [0.0, 0.0, 0.0] + position = [0.0, 0.0, 0.0] # Point the camera looks at, in world coordinates (the dome ZENITH target) # @type: array [size 3] # @default: box center; for `mode = "dome"`, the zenith defaults to +z - look_at = [0.0, 0.0, 0.0] + look_at = [0.0, 0.0, 0.0] # Camera up vector (dome: the disk's screen-up) # @type: array [size 3] # @default: [0.0, 0.0, 1.0]; for `mode = "dome"`, [0.0, 1.0, 0.0] - up = [0.0, 0.0, 1.0] + up = [0.0, 0.0, 1.0] # Vertical field of view in degrees (perspective only) # @type: float [> 0.0] # @default: 35.0 - fov = 35.0 + fov = 35.0 # Full dome field of view in degrees (dome mode only): the image rim is at # dome_fov/2 from the zenith (180 = a full hemisphere down to the horizon). # @type: float [> 0.0, <= 360.0] # @default: 180.0 - dome_fov = 180.0 + dome_fov = 180.0 # Dome far-clip radius in world units (dome mode only): each ray stops this # far from the eye, so the sampled region is a half-ball (hemisphere) of # this radius rather than the whole box -> uniform path length and no box @@ -1052,7 +1094,7 @@ # @default: the largest sphere centered in the box (half the shortest # side), so it touches the face centers and never a corner # @note: 0 disables the clip (rays march to the box boundary) - dome_radius = 0.0 + dome_radius = 0.0 # Vertical extent of the view in world units (orthographic only) # @type: float [> 0.0] # @default: the global box diagonal (the whole box fits from any angle) @@ -1077,21 +1119,21 @@ # Build the fisheye dome master instead of the plain slice # @type: bool # @default: false - enable = false + enable = false # (Cartesian only) Full dome field of view in degrees (image radius maps # linearly to the dome zenith angle: the rim is at fov/2) # @type: float [> 0.0, <= 180.0] # @default: 180.0 # a full hemisphere - fov = 180.0 + fov = 180.0 # (Cartesian only) World radius of the circular cutout mapped onto the dome # @type: float [> 0.0] # @default: half the shorter domain side (the largest centered disk that # fits inside the box) - radius = 1.0 + radius = 1.0 # (Cartesian only) World-space center of the cutout # @type: array [size 2] # @default: the domain center - center = [0.0, 0.0] + center = [0.0, 0.0] # (Cartesian only) How the dome zenith angle maps to a world radius on the # flat slice # @type: string @@ -1121,38 +1163,38 @@ # @default: false # @note: implied true if any scene sets `fieldlines = true` or uses `field # = "fieldlines"` - enable = false + enable = false # Vector field to trace # @type: string # @default: "B" # @enum: "B", "E", "J" - field = "B" + field = "B" # Field coarsening factor (simulation cells per coarse cell, per axis) # @type: int [1..16] # @default: 4 # @note: larger = smoother "morphology" lines + cheaper replication (the # coarse field is ~ N_cells / bin^D floats/rank; D = sim dimension) - bin = 4 + bin = 4 # (3D tubes) Seed-lattice spacing in screen pixels (sets line density) # @type: float [> 0] # @default: 8 # @note: capped by `seed_max`; if seed_px asks for more seeds than that, # the spacing grows to fit and seed_px no longer governs - seed_px = 8 + seed_px = 8 # (3D tubes) Hard cap on the seed count (lattice is n^3, 2 lines per seed) # @type: int [> 0] # @default: 4096 # @note: lower this for fewer / more widely spaced lines - seed_max = 4096 + seed_max = 4096 # (2D contours) Number of evenly-spaced flux-function contour levels # @type: int [> 0] # @default: 16 # @note: evenly-spaced psi levels => line density tracks |B| automatically - levels = 16 + levels = 16 # Tube radius (3D) / contour line width (2D), in screen pixels # @type: float [> 0] # @default: 2 - tube_px = 2 + tube_px = 2 # Colormap for the field lines (mapped by |B| along each line) # @type: string # @default: "inferno" @@ -1161,38 +1203,38 @@ # "dusk", "cosmic", "freeze", "apple", "gothic", "sunburst", # "voltage", "ocean", "fusion", "prinsenvlag" (an optional "cmr." # prefix is ok) - colormap = "inferno" + colormap = "inferno" # Monochrome override: draw the lines in a single [r,g,b] color (each 0..1) # instead of the |B| colormap -- reads well as an overlay on another volume # @type: array [size 3] # @default: [] (empty => color by |B|) # @example: [1.0, 1.0, 1.0] # white field lines - color = [] + color = [] # Map the tube color range logarithmically # @type: bool # @default: false # @note: requires min > 0 - log = false + log = false # Tube color range: lower bound on |field| # @type: float # @default: 0.0 # @note: when min >= max, the range is auto-set from |field| along the # lines - min = 0.0 + min = 0.0 # Tube color range: upper bound on |field| # @type: float # @default: 0.0 # @note: when min >= max, the range is auto-set from |field| along the # lines - max = 0.0 + max = 0.0 # (3D tubes) RK4 integration step as a fraction of one coarse cell # @type: float [> 0] # @default: 0.5 - step_frac = 0.5 + step_frac = 0.5 # (3D tubes) Per-direction integration-step cap # @type: int [> 0] # @default: 4000 - max_steps = 4000 + max_steps = 4000 # (3D tubes) Maximum line length, in global box diagonals (per direction) # @type: float [> 0] # @default: 3.0 @@ -1226,28 +1268,28 @@ # diverging colormap ("cool2warm") to center zero # @note: "fieldlines" renders the magnetic field-line tubes on their own # (no scalar volume sampled); see [render.fieldlines] below - field = "" + field = "" # PNG filename prefix; files are `.png` # @type: string # @default: "_" - prefix = "_" + prefix = "_" # Colorbar title # @type: string # @default: `field` - label = "" + label = "" # Lower bound of the value range mapped onto the colormap/opacity # @type: float # @default: 0.0 - min = 0.0 + min = 0.0 # Upper bound of the value range # @type: float # @default: 1.0 - max = 1.0 + max = 1.0 # Map the value range logarithmically # @type: bool # @default: false # @note: Requires min > 0 and max > 0 - log = false + log = false # Colormap name # @type: string # @default: "viridis" @@ -1256,14 +1298,14 @@ # "dusk", "cosmic", "freeze", "apple", "gothic", "sunburst", # "voltage", "ocean", "fusion", "prinsenvlag" (an optional "cmr." # prefix is ok) - colormap = "viridis" + colormap = "viridis" # Opacity transfer function: [position, opacity] control points, both in [0, # 1], piecewise-linear in the normalized value # @type: array> # @default: linear ramp (opacity = normalized value) # @note: Keep the low end near 0 so empty regions stay transparent # @example: [[0.0, 0.0], [0.3, 0.1], [1.0, 0.7]] - alpha = [] + alpha = [] # Explicit value(s) to label on the colorbar # @type: array # @default: 5 evenly-spaced ticks between min and max @@ -1275,14 +1317,14 @@ # @default: false # @note: requires the [render.fieldlines] enabled (3D only). A scene with # field = "fieldlines" instead renders them alone. - fieldlines = false + fieldlines = false # Diagnostic logging parameters [diagnostics] # Number of timesteps between diagnostic logs # @type: int [> 0] # @default: 1 - interval = 1 + interval = 1 # Blocking timers between successive algorithms # @type: bool # @default: false @@ -1290,11 +1332,11 @@ # Enable colored stdout # @type: bool # @default: true - colored_stdout = true + colored_stdout = true # Specify the log level # @type: string # @default: "VERBOSE" # @enum: "VERBOSE", "WARNING", "ERROR" # @note: "VERBOSE" prints all messages, "WARNING" prints only warnings and # errors, "ERROR" prints only errors - log_level = "VERBOSE" + log_level = "VERBOSE" diff --git a/pgens/rms/cfg.rms.toml b/pgens/rms/cfg.rms.toml new file mode 100644 index 000000000..81b319d03 --- /dev/null +++ b/pgens/rms/cfg.rms.toml @@ -0,0 +1,98 @@ +[simulation] + name = "rms" + engine = "srpic" + runtime = 20000.0 + +[grid] + resolution = [32768] + extent = [[0.0, 4000.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["CONDUCTOR", "MATCH"]] + particles = [["REFLECT", "ABSORB"]] + +[scales] + larmor0 = 3.33333 + skindepth0 = 1.0 + +[algorithms] + current_filters = 8 + + [algorithms.timestep] + CFL = 0.7 + +[two_body] + thomson_optical_depth = 1e-3 + + [[two_body.interaction]] + type = "compton" + group1 = [1] + group2 = [3] + interval = 1 + tile_size = 8 + +[particles] + ppc0 = 8.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 2e7 + + [[particles.species]] + label = "i" + mass = 25.0 + charge = 1.0 + maxnpart = 2e7 + + [[particles.species]] + label = "ph" + mass = 0.0 + charge = 0.0 + maxnpart = 2e7 + +[setup] + drift_ux = 0.1 # speed towards the wall [c] + temperature = 0.0001 # temperature of maxwell distribution [T / (m0 c^2)] + temperature_ratio = 1.0 # temperature ratio of electrons to protons + Bmag = 1.0 # magnetic field strength as fraction of magnetisation + Btheta = 60.0 # magnetic field angle in the plane + Bphi = 0.0 # magnetic field angle out of plane + filling_fraction = 1.0 # fraction of the shock piston filled with plasma + injector_velocity = 0.0 # speed of injector [c] + injection_start = 0.0 # start time of moving injector + injection_frequency = 100 + photon_injection_rate = 0.01 + +[output] + interval_time = 100.0 + + [output.fields] + quantities = ["N_1", "N_2", "N_3", "B", "E"] + custom = [ + "Temperature_e", + "Temperature_i", + "Temperature_ph", + "Vbx", + "Vby", + "Vbz", + ] + mom_smooth = 4 + + [output.particles] + enable = true + stride = 10 + + [output.spectra] + enable = true + +[checkpoint] + interval = 1000 + keep = 2 + +[diagnostics] + colored_stdout = true diff --git a/pgens/rms/pgen.hpp b/pgens/rms/pgen.hpp new file mode 100644 index 000000000..980a1b708 --- /dev/null +++ b/pgens/rms/pgen.hpp @@ -0,0 +1,963 @@ +#ifndef PROBLEM_GENERATOR_H +#define PROBLEM_GENERATOR_H + +#include "enums.h" +#include "global.h" + +#include "traits/pgen.h" +#include "utils/error.h" +#include "utils/numeric.h" + +#include "archetypes/utils.h" +#include "framework/containers/particles.h" +#include "framework/domain/metadomain.h" + +#include +#include + +namespace user { + using namespace ntt; + + enum Component : idx_t { + comp_n = 0u, + comp_rho = 1u, + comp_rho_vx = 2u, + comp_rho_vy = 3u, + comp_rho_vz = 4u, + comp_n_t = 5u, + }; + + struct ComputeN { + static constexpr uint8_t N = 1; + + Inline void operator()(const ParticleArrays& /* prtls */, + float /* mass */, + float /* charge */, + prtlidx_t /* p */, + list_t& contribs) const { + contribs[0] = ONE; + } + }; + + struct ComputeRhoV { + static constexpr uint8_t N = 5; + + Inline void operator()(const ParticleArrays& prtls, + float mass, + float /* charge */, + prtlidx_t p, + list_t& contribs) const { + contribs[0] = ONE; + contribs[1] = mass; + const auto gamma = U2GAMMA(prtls.ux1(p), prtls.ux2(p), prtls.ux3(p)); + contribs[2] = mass * prtls.ux1(p) / gamma; + contribs[3] = mass * prtls.ux2(p) / gamma; + contribs[4] = mass * prtls.ux3(p) / gamma; + } + }; + + struct ComputePhotonTmunu { + static constexpr uint8_t N = 5; + + Inline void operator()(const ParticleArrays& prtls, + float /* mass */, + float /* charge */, + prtlidx_t p, + list_t& contribs) const { + const auto energy = math::sqrt( + SQR(prtls.ux1(p)) + SQR(prtls.ux2(p)) + SQR(prtls.ux3(p))); + contribs[0] = ONE; + contribs[1] = energy; + contribs[2] = prtls.ux1(p); + contribs[3] = prtls.ux2(p); + contribs[4] = prtls.ux3(p); + } + }; + + template + struct ComputeN_Rho_Vi { + static constexpr uint8_t N = 1; + + Inline void operator()(const ParticleArrays& prtls, + float mass, + float /* charge */, + prtlidx_t p, + list_t& contribs) const { + const auto gamma = U2GAMMA(prtls.ux1(p), prtls.ux2(p), prtls.ux3(p)); + if constexpr (C == -1) { + contribs[0] = ONE; + } else if constexpr (C == 0) { + contribs[0] = mass; + } else if constexpr (C == 1) { + contribs[0] = mass * prtls.ux1(p) / gamma; + } else if constexpr (C == 2) { + contribs[0] = mass * prtls.ux2(p) / gamma; + } else if constexpr (C == 3) { + contribs[0] = mass * prtls.ux3(p) / gamma; + } + } + }; + + template + struct Normalize { + ndfield_t buffer; + const idx_t comp, comp_norm; + + Normalize(const ndfield_t& buff, idx_t comp, idx_t comp_norm) + : buffer { buff } + , comp { comp } + , comp_norm { comp_norm } {} + + Inline void operator()(cellidx_t i1) const { + if constexpr (D == Dim::_1D) { + if (buffer(i1, comp_norm) < 1e-6) { + buffer(i1, comp) = ZERO; + } else { + buffer(i1, comp) /= buffer(i1, comp_norm); + } + } else { + raise::KernelError(HERE, "Normalize is only implemented for 1D"); + } + } + + Inline void operator()(cellidx_t i1, cellidx_t i2) const { + if constexpr (D == Dim::_2D) { + if (buffer(i1, i2, comp_norm) < 1e-6) { + buffer(i1, i2, comp) = ZERO; + } else { + buffer(i1, i2, comp) /= buffer(i1, i2, comp_norm); + } + } else { + raise::KernelError(HERE, "Normalize is only implemented for 2D"); + } + } + + Inline void operator()(cellidx_t i1, cellidx_t i2, cellidx_t i3) const { + if constexpr (D == Dim::_3D) { + if (buffer(i1, i2, i3, comp_norm) < 1e-6) { + buffer(i1, i2, i3, comp) = ZERO; + } else { + buffer(i1, i2, i3, comp) /= buffer(i1, i2, i3, comp_norm); + } + } else { + raise::KernelError(HERE, "Normalize is only implemented for 3D"); + } + } + }; + + template + struct ComputePressure { + static constexpr uint8_t N = 1; + ndfield_t buffer; + + const idx_t rho, rho_vx, rho_vy, rho_vz; + + ComputePressure(const ndfield_t& buff, + idx_t rho = comp_rho, + idx_t rho_vx = comp_rho_vx, + idx_t rho_vy = comp_rho_vy, + idx_t rho_vz = comp_rho_vz) + : buffer { buff } + , rho { rho } + , rho_vx { rho_vx } + , rho_vy { rho_vy } + , rho_vz { rho_vz } {} + + Inline void operator()(const ParticleArrays& prtls, + float mass, + float, + prtlidx_t p, + list_t& contribs) const { + real_t Vx1 { ZERO }, Vx2 { ZERO }, Vx3 { ZERO }; + if constexpr (D == Dim::_1D) { + Vx1 = buffer(prtls.i1(p) + N_GHOSTS, rho_vx) / + buffer(prtls.i1(p) + N_GHOSTS, rho); + Vx2 = buffer(prtls.i1(p) + N_GHOSTS, rho_vy) / + buffer(prtls.i1(p) + N_GHOSTS, rho); + Vx3 = buffer(prtls.i1(p) + N_GHOSTS, rho_vz) / + buffer(prtls.i1(p) + N_GHOSTS, rho); + } else if constexpr (D == Dim::_2D) { + Vx1 = buffer(prtls.i1(p) + N_GHOSTS, prtls.i2(p) + N_GHOSTS, rho_vx) / + buffer(prtls.i1(p) + N_GHOSTS, prtls.i2(p) + N_GHOSTS, rho); + Vx2 = buffer(prtls.i1(p) + N_GHOSTS, prtls.i2(p) + N_GHOSTS, rho_vy) / + buffer(prtls.i1(p) + N_GHOSTS, prtls.i2(p) + N_GHOSTS, rho); + Vx3 = buffer(prtls.i1(p) + N_GHOSTS, prtls.i2(p) + N_GHOSTS, rho_vz) / + buffer(prtls.i1(p) + N_GHOSTS, prtls.i2(p) + N_GHOSTS, rho); + } else if constexpr (D == Dim::_3D) { + Vx1 = buffer(prtls.i1(p) + N_GHOSTS, + prtls.i2(p) + N_GHOSTS, + prtls.i3(p) + N_GHOSTS, + rho_vx) / + buffer(prtls.i1(p) + N_GHOSTS, + prtls.i2(p) + N_GHOSTS, + prtls.i3(p) + N_GHOSTS, + rho); + Vx2 = buffer(prtls.i1(p) + N_GHOSTS, + prtls.i2(p) + N_GHOSTS, + prtls.i3(p) + N_GHOSTS, + rho_vy) / + buffer(prtls.i1(p) + N_GHOSTS, + prtls.i2(p) + N_GHOSTS, + prtls.i3(p) + N_GHOSTS, + rho); + Vx3 = buffer(prtls.i1(p) + N_GHOSTS, + prtls.i2(p) + N_GHOSTS, + prtls.i3(p) + N_GHOSTS, + rho_vz) / + buffer(prtls.i1(p) + N_GHOSTS, + prtls.i2(p) + N_GHOSTS, + prtls.i3(p) + N_GHOSTS, + rho); + } + const auto Gamma = ONE / math::sqrt(ONE - SQR(Vx1) - SQR(Vx2) - SQR(Vx3)); + const auto gamma = U2GAMMA(prtls.ux1(p), prtls.ux2(p), prtls.ux3(p)); + contribs[0] = mass * + (SQR(Gamma) * SQR(gamma - Vx1 * prtls.ux1(p) - + Vx2 * prtls.ux2(p) - Vx3 * prtls.ux3(p)) - + ONE) / + (THREE * gamma); + } + }; + + template + struct ComputePhotonTemperature { + const ndfield_t t_array; + const idx_t comp_t, comp_n, comp_t00, comp_t01, comp_t02, comp_t03; + + ComputePhotonTemperature(ndfield_t& t_array, + idx_t comp_t, + idx_t comp_n, + idx_t comp_t00, + idx_t comp_t01, + idx_t comp_t02, + idx_t comp_t03) + : t_array { t_array } + , comp_t { comp_t } + , comp_n { comp_n } + , comp_t00 { comp_t00 } + , comp_t01 { comp_t01 } + , comp_t02 { comp_t02 } + , comp_t03 { comp_t03 } {} + + Inline auto Gamma(real_t T00_Sqr, real_t T0i_Sqr) const -> real_t { + return HALF * math::sqrt(T0i_Sqr) * + math::sqrt( + ONE / (T0i_Sqr + math::sqrt(T00_Sqr) * + (math::sqrt(FOUR * T00_Sqr - THREE * T0i_Sqr) - + TWO * math::sqrt(T00_Sqr)))); + } + + Inline void operator()(cellidx_t i1) const { + if constexpr (D == Dim::_1D) { + const auto T0i_Sqr = (SQR(t_array(i1, comp_t01)) + + SQR(t_array(i1, comp_t02)) + + SQR(t_array(i1, comp_t03))); + const auto T00_Sqr = SQR(t_array(i1, comp_t00)); + t_array(i1, comp_t) = (math::sqrt(FOUR * T00_Sqr - THREE * T0i_Sqr) - + math::sqrt(T00_Sqr)) / + (t_array(i1, comp_n) * Gamma(T00_Sqr, T0i_Sqr) * + static_cast(2.7)); + } else { + raise::KernelError(HERE, "ComputePhotonTemperature is only implemented for 1D"); + } + } + + Inline void operator()(cellidx_t i1, cellidx_t i2) const { + if constexpr (D == Dim::_2D) { + const auto T0i_Sqr = (SQR(t_array(i1, i2, comp_t01)) + + SQR(t_array(i1, i2, comp_t02)) + + SQR(t_array(i1, i2, comp_t03))); + const auto T00_Sqr = SQR(t_array(i1, i2, comp_t00)); + t_array(i1, i2, comp_t) = (math::sqrt(FOUR * T00_Sqr - THREE * T0i_Sqr) - + math::sqrt(T00_Sqr)) / + (t_array(i1, i2, comp_n) * + Gamma(T00_Sqr, T0i_Sqr) * + static_cast(2.7)); + } else { + raise::KernelError(HERE, "ComputePhotonTemperature is only implemented for 2D"); + } + } + + Inline void operator()(cellidx_t i1, cellidx_t i2, cellidx_t i3) const { + if constexpr (D == Dim::_3D) { + const auto T0i_Sqr = (SQR(t_array(i1, i2, i3, comp_t01)) + + SQR(t_array(i1, i2, i3, comp_t02)) + + SQR(t_array(i1, i2, i3, comp_t03))); + const auto T00_Sqr = SQR(t_array(i1, i2, i3, comp_t00)); + t_array(i1, i2, i3, comp_t) = (math::sqrt(FOUR * T00_Sqr - THREE * T0i_Sqr) - + math::sqrt(T00_Sqr)) / + (t_array(i1, i2, i3, comp_n) * + Gamma(T00_Sqr, T0i_Sqr) * + static_cast(2.7)); + } else { + raise::KernelError(HERE, "ComputePhotonTemperature is only implemented for 3D"); + } + } + }; + + template + struct PhotonSpatialDistribution { + const ndfield_t n_array; + + const M metric; + const real_t nmax { 4.0 }; + + PhotonSpatialDistribution(const ndfield_t& n_array, const M& metric) + : n_array { n_array } + , metric { metric } {} + + Inline auto operator()(const coord_t& x_Ph) const -> real_t { + coord_t x_Cd { ZERO }; + metric.template convert(x_Ph, x_Cd); + if constexpr (M::Dim == Dim::_1D) { + return n_array(static_cast(x_Cd[0]) + N_GHOSTS, comp_n) / nmax; + } else if constexpr (M::Dim == Dim::_2D) { + return n_array(static_cast(x_Cd[0]) + N_GHOSTS, + static_cast(x_Cd[1]) + N_GHOSTS, + comp_n) / + nmax; + } else if constexpr (M::Dim == Dim::_3D) { + return n_array(static_cast(x_Cd[0]) + N_GHOSTS, + static_cast(x_Cd[1]) + N_GHOSTS, + static_cast(x_Cd[2]) + N_GHOSTS, + comp_n) / + nmax; + } + } + }; + + template + struct PlanckDistribution { + const real_t T_ph_inj; + const real_t boost_beta; + + random_number_pool_t random_pool; + + PlanckDistribution(real_t T_ph_inj, real_t boost_beta, random_number_pool_t& pool) + : T_ph_inj { T_ph_inj } + , boost_beta { boost_beta } + , random_pool { pool } {} + + Inline void operator()(const coord_t&, vec_t& k) const { + real_t prob { ZERO }, n { ZERO }; + auto gen = random_pool.get_state(); + const auto rnd = Random(gen); + const auto rnd1 = Random(gen); + const auto rnd2 = Random(gen); + const auto rnd3 = Random(gen); + const auto rndth = Random(gen); + const auto rndph = Random(gen); + random_pool.free_state(gen); + + while ((prob < rnd) and (n < 40)) { + n += ONE; + prob += ONE / (static_cast(1.20206) * CUBE(n)); + } + const auto energy = -static_cast(2.7) * T_ph_inj * + math::log( + rnd1 * rnd2 * rnd3 + static_cast(1e-16)) / + n; + const auto costh = TWO * rndth - ONE; + const auto phi = static_cast(constant::TWO_PI) * rndph; + + k[0] = energy * math::sqrt(ONE - SQR(costh)) * math::cos(phi); + k[1] = energy * math::sqrt(ONE - SQR(costh)) * math::sin(phi); + k[2] = energy * costh; + + // boost the photon momentum in the -x direction + const auto gamma = ONE / math::sqrt(ONE - SQR(boost_beta)); + const auto kx = k[0]; + const auto ky = k[1]; + const auto kz = k[2]; + k[0] = gamma * (kx - boost_beta * energy); + k[1] = ky; + k[2] = kz; + } + }; + + template + struct InitFields { + /* + Sets up magnetic and electric field components for the simulation. + Must satisfy E = -v x B for Lorentz Force to be zero. + + @param bmag: magnetic field scaling + @param thetaB: Bx = bmag * cos(thetaB) + @param beta_upstream: drift three-velocity in the x direction + */ + InitFields(real_t bmag, real_t thetaB, real_t beta_upstream) + : Bmag { bmag } + , thetaB { thetaB * static_cast(convert::deg2rad) } + , beta_upstream { beta_upstream } {} + + // magnetic field components + Inline auto bx1(const coord_t&) const -> real_t { + return Bmag * math::cos(thetaB); + } + + Inline auto bx2(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto bx3(const coord_t&) const -> real_t { + return Bmag * math::sin(thetaB); + } + + // electric field components + Inline auto ex1(const coord_t&) const -> real_t { + return ZERO; + } + + Inline auto ex2(const coord_t&) const -> real_t { + return -beta_upstream * Bmag * math::sin(thetaB); + } + + Inline auto ex3(const coord_t&) const -> real_t { + return ZERO; + } + + private: + const real_t Bmag, thetaB, beta_upstream; + }; + + template + struct PGen { + static constexpr auto D { M::Dim }; + // compatibility traits for the problem generator + static constexpr auto engines { + ::traits::pgen::compatible_with {} + }; + static constexpr auto metrics { + ::traits::pgen::compatible_with {} + }; + static constexpr auto dimensions { + ::traits::pgen::compatible_with {} + }; + const SimulationParams& params; + Metadomain& metadomain; + + // domain properties + const real_t global_xmin, global_xmax; + // gas properties + const real_t beta_upstream, Te, Te_ovr_Ti; + // magnetic field properties + const real_t Bmag, thetaB; + // photon properties + // const real_t photon_inj_rate; // units of n0 / time + const real_t photon_density; // units of n0 + const real_t T_ph_inj; // photon injection temperature + // plasma injector properties + const real_t filling_fraction, beta_injector; + const int injection_interval; + + InitFields init_flds; + + PGen(const SimulationParams& p, Metadomain& m) + : params { p } + , metadomain { m } + , global_xmin { metadomain.mesh().extent(in::x1).first } + , global_xmax { metadomain.mesh().extent(in::x1).second } + , beta_upstream { params.template get("setup.beta_upstream") } + , Te { params.template get("setup.Te") } + , Te_ovr_Ti { params.template get("setup.Te_ovr_Ti", ONE) } + , Bmag { params.template get("setup.Bmag", ZERO) } + , thetaB { params.template get("setup.thetaB", ZERO) } + // , photon_inj_rate { params.template get("setup.photon_inj_rate", ZERO) } + , photon_density { params.template get("setup.photon_density", ZERO) } + , T_ph_inj { params.template get("setup.T_ph_inj") } + , filling_fraction { params.template get("setup.filling_fraction", + 1.0) } + , beta_injector { params.template get("setup.beta_injector", 1.0) } + , injection_interval { params.template get( + "setup.injection_interval", + 100) } + , init_flds { Bmag, thetaB, beta_upstream } {} + + auto MatchFields(simtime_t) const -> InitFields { + return init_flds; + } + + auto FixFieldsConst(const bc_in&, const em& comp) const + -> std::pair { + if (comp == em::ex1) { + return { init_flds.ex1({ ZERO }), true }; + } else if ((comp == em::ex2) or (comp == em::ex3)) { + return { ZERO, true }; + } else if (comp == em::bx1) { + return { init_flds.bx1({ ZERO }), true }; + } else if (comp == em::bx2) { + return { init_flds.bx2({ ZERO }), true }; + } else if (comp == em::bx3) { + return { init_flds.bx3({ ZERO }), true }; + } else { + raise::Error("Invalid component", HERE); + return { ZERO, false }; + } + } + + void InitPrtls(Domain& domain) { + /* + * Plasma setup as partially filled box + * + * Plasma setup: + * + * global_xmin global_xmax + * | | + * V V + * |:::::::::::|..........................| + * ^ + * | + * filling_fraction + */ + + // minimum and maximum position of particles + real_t xg_min = global_xmin; + real_t xg_max = global_xmin + filling_fraction * (global_xmax - global_xmin); + + // define box to inject into + boundaries_t box; + // loop over all dimensions + for (auto d { 0u }; d < (unsigned int)M::Dim; ++d) { + // compute the range for the x-direction + if (d == static_cast(in::x1)) { + box.emplace_back(xg_min, xg_max); + } else { + // inject into full range in other directions + box.push_back(Range::All); + } + } + + const auto gamma_upstream = ONE / math::sqrt(ONE - SQR(beta_upstream)); + + // inject particles + arch::InjectUniformMaxwellians( + params, + domain, + TWO, + std::make_pair(Te, Te / Te_ovr_Ti), + { 1, 2 }, + std::make_pair( + std::vector { -gamma_upstream * beta_upstream, ZERO, ZERO }, + std::vector { -gamma_upstream * beta_upstream, ZERO, ZERO }), + false, + box); + + const auto planck_dist = PlanckDistribution(T_ph_inj, + beta_upstream, + domain.random_pool()); + arch::InjectUniform(params, domain, 3, planck_dist, photon_density, false, box); + } + + void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { + const auto dt = params.template get("algorithms.timestep.dt"); + + if (step % injection_interval == 0) { + /* + * Replenish plasma in a moving injector + * + * Injector setup: + * + * global_xmin purge/replenish global_xmax + * | x_init | | + * V v V V + * |:::::::::::;::::::::::|\\\\\\\\|......| + * xmin xmax + * ^ + * | + * moving injector + */ + + // initial position of injector + const auto x_init = global_xmin + + filling_fraction * (global_xmax - global_xmin); + + // compute the position of the injector after the current timestep + const auto xmax = std::min(x_init + beta_injector * (step + 1) * dt, + global_xmax); + + // compute the beginning of the injected region + const auto xmin = (step == 0) + ? std::max(x_init - beta_upstream * dt, + global_xmin) + : xmax - injection_interval * dt * beta_injector - + (injection_interval + 1) * dt * beta_upstream; + + // define indice range to reset fields + boundaries_t incl_ghosts; + for (auto d = 0; d < M::Dim; ++d) { + incl_ghosts.emplace_back(false, false); + } + + // define box to reset fields + boundaries_t purge_box; + // loop over all dimension + for (auto d = 0u; d < M::Dim; ++d) { + if (d == 0) { + purge_box.emplace_back(xmin, global_xmax); + } else { + purge_box.push_back(Range::All); + } + } + + const auto extent = domain.mesh.ExtentToRange(purge_box, incl_ghosts); + tuple_t x_min { 0 }, x_max { 0 }; + for (auto d = 0; d < M::Dim; ++d) { + x_min[d] = extent[d].first; + x_max[d] = extent[d].second; + } + + Kokkos::parallel_for("ResetFields", + CreateRangePolicy(x_min, x_max), + arch::SetEMFields_kernel { + domain.fields.em, + init_flds, + domain.mesh.metric }); + metadomain.CommunicateFields(domain, Comm::E | Comm::B); + + /* + tag particles inside the injection zone as dead + */ + // const auto& mesh = domain.mesh; + + // loop over particle species + // for (auto& species : domain.species) { + // // get particle properties + // auto i1 = species.i1; + // auto dx1 = species.dx1; + // auto tag = species.tag; + + // Kokkos::parallel_for( + // "RemoveParticles", + // species.rangeActiveParticles(), + // Lambda(prtlidx_t p) { + // // check if the particle is already dead + // if (tag(p) == ParticleTag::dead) { + // return; + // } + // const auto x_Cd = static_cast(i1(p)) + + // static_cast(dx1(p)); + // const auto x_Ph = mesh.metric.template convert<1, Crd::Cd, Crd::XYZ>( + // x_Cd); + + // if (x_Ph > xmin) { + // tag(p) = ParticleTag::dead; + // } + // }); + // } + + // define box to inject into + boundaries_t inj_box; + // loop over all dimension + for (auto d = 0u; d < M::Dim; ++d) { + if (d == 0) { + inj_box.emplace_back(xmin, xmax); + } else { + inj_box.push_back(Range::All); + } + } + + const auto gamma_upstream = ONE / math::sqrt(ONE - SQR(beta_upstream)); + + // same maxwell distribution as above + arch::InjectUniformMaxwellians( + params, + domain, + TWO, + std::make_pair(Te, Te / Te_ovr_Ti), + { 1, 2 }, + std::make_pair( + std::vector { -gamma_upstream * beta_upstream, ZERO, ZERO }, + std::vector { -gamma_upstream * beta_upstream, ZERO, ZERO }), + false, + inj_box); + + // replenish photons + const auto planck_dist = PlanckDistribution(T_ph_inj, + beta_upstream, + domain.random_pool()); + arch::InjectUniform(params, domain, 3, planck_dist, photon_density, false, inj_box); + } + + // { + // /* + // * Inject photons + // */ + // auto compute_n = ComputeN {}; + // arch::ComputeMomentWithSpeciesNew( + // params, + // domain, + // { 1, 2 }, + // domain.fields.bckp, + // { comp_n }, + // compute_n); + // + // // inject photons with a Planck distribution in energy and spatial distribution following the plasma density + // const auto energy_dist = PlanckDistribution(T_ph_inj, + // domain.random_pool()); + // const auto spatial_dist = PhotonSpatialDistribution(domain.fields.bckp, + // domain.mesh.metric); + // arch::InjectNonUniform( + // params, + // domain, + // 3, + // energy_dist, + // spatial_dist, + // static_cast(photon_inj_rate * dt)); + // } + } + + void CustomFieldOutput(const std::string& label, + ndfield_t& buff, + cellidx_t buff_idx, + timestep_t, + simtime_t, + const Domain& domain) { + const uint8_t smoothing_order = 2u * N_GHOSTS; + if (label == "Vbx") { + /** + * buff_idx + 1 -> rho_e + rho_i + */ + auto compute_rho = ComputeN_Rho_Vi<0> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)((buff_idx + 1) % 6) }, + compute_rho, + smoothing_order); + /** + * buff_idx -> rho_e * vxe + rho_i * vxi + */ + auto compute_rho_vx = ComputeN_Rho_Vi<1> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)buff_idx }, + compute_rho_vx, + smoothing_order); + /** + * buff_idx -> vx = (rho_e * vxe + rho_i * vxi) / (rho_e + rho_i) + */ + Kokkos::parallel_for( + "ComputeVx", + domain.mesh.rangeActiveCells(), + Normalize { buff, (idx_t)(buff_idx), (idx_t)((buff_idx + 1) % 6) }); + } else if (label == "Vby") { + /** + * buff_idx + 1 -> rho_e + rho_i + */ + auto compute_rho = ComputeN_Rho_Vi<0> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)((buff_idx + 1) % 6) }, + compute_rho, + smoothing_order); + /** + * buff_idx -> rho_e * vye + rho_i * vyi + */ + auto compute_rho_vy = ComputeN_Rho_Vi<2> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)(buff_idx) }, + compute_rho_vy, + smoothing_order); + /** + * buff_idx -> vy = (rho_e * vye + rho_i * vyi) / (rho_e + rho_i) + */ + Kokkos::parallel_for( + "ComputeVy", + domain.mesh.rangeActiveCells(), + Normalize { buff, (idx_t)(buff_idx), (idx_t)((buff_idx + 1) % 6) }); + } else if (label == "Vbz") { + /** + * buff_idx + 1 -> rho_e + rho_i + */ + auto compute_rho = ComputeN_Rho_Vi<0> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)((buff_idx + 1) % 6) }, + compute_rho, + smoothing_order); + /** + * buff_idx -> rho_e * vye + rho_i * vyi + */ + auto compute_rho_vz = ComputeN_Rho_Vi<3> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)buff_idx }, + compute_rho_vz, + smoothing_order); + /** + * buff_idx -> vy = (rho_e * vye + rho_i * vyi) / (rho_e + rho_i) + */ + Kokkos::parallel_for( + "ComputeVy", + domain.mesh.rangeActiveCells(), + Normalize { buff, (idx_t)(buff_idx), (idx_t)((buff_idx + 1) % 6) }); + } else if (label == "Temperature_e") { + /** + * buff_idx + 1 -> n_e + n_i + * buff_idx + 2 -> rho_e + rho_i + * buff_idx + 3 -> rho_e * vxe + rho_i * vxi + * buff_idx + 4 -> rho_e * vye + rho_i * vyi + * buff_idx + 5 -> rho_e * vze + rho_i * vzi + */ + auto compute_rho_v = ComputeRhoV {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)((buff_idx + 1) % 6), + (idx_t)((buff_idx + 2) % 6), + (idx_t)((buff_idx + 3) % 6), + (idx_t)((buff_idx + 4) % 6), + (idx_t)((buff_idx + 5) % 6) }, + compute_rho_v, + smoothing_order); + /** + * buff_idx -> n_e * T_e + */ + auto compute_pressure = ComputePressure { buff, + (idx_t)((buff_idx + 2) % 6), + (idx_t)((buff_idx + 3) % 6), + (idx_t)((buff_idx + 4) % 6), + (idx_t)((buff_idx + 5) % 6) }; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1 }, + buff, + { (idx_t)buff_idx }, + compute_pressure, + smoothing_order); + /** + * buff_idx + 1 -> n_e + */ + auto compute_n = ComputeN_Rho_Vi<-1> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1 }, + buff, + { (idx_t)((buff_idx + 1) % 6) }, + compute_n, + smoothing_order); + /** + * buff_idx -> T_e + */ + Kokkos::parallel_for( + "ComputeTe", + domain.mesh.rangeActiveCells(), + Normalize { buff, (idx_t)(buff_idx), (idx_t)((buff_idx + 1) % 6) }); + } else if (label == "Temperature_i") { + /** + * buff_idx + 1 -> n_e + n_i + * buff_idx + 2 -> rho_e + rho_i + * buff_idx + 3 -> rho_e * vxe + rho_i * vxi + * buff_idx + 4 -> rho_e * vye + rho_i * vyi + * buff_idx + 5 -> rho_e * vze + rho_i * vzi + */ + auto compute_rho_v = ComputeRhoV {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 1, 2 }, + buff, + { (idx_t)((buff_idx + 1) % 6), + (idx_t)((buff_idx + 2) % 6), + (idx_t)((buff_idx + 3) % 6), + (idx_t)((buff_idx + 4) % 6), + (idx_t)((buff_idx + 5) % 6) }, + compute_rho_v, + smoothing_order); + /** + * buff_idx -> n_i * T_i + */ + auto compute_pressure = ComputePressure { buff, + (idx_t)((buff_idx + 2) % 6), + (idx_t)((buff_idx + 3) % 6), + (idx_t)((buff_idx + 4) % 6), + (idx_t)((buff_idx + 5) % 6) }; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 2 }, + buff, + { (idx_t)buff_idx }, + compute_pressure, + smoothing_order); + /** + * buff_idx + 1 -> n_i + */ + auto compute_n = ComputeN_Rho_Vi<-1> {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 2 }, + buff, + { (idx_t)((buff_idx + 1) % 6) }, + compute_n, + smoothing_order); + /** + * buff_idx -> T_i + */ + Kokkos::parallel_for( + "ComputeTi", + domain.mesh.rangeActiveCells(), + Normalize { buff, (idx_t)buff_idx, (idx_t)((buff_idx + 1) % 6) }); + } else if (label == "Temperature_ph") { + /** + * buff_idx + 1 -> n_ph + * buff_idx + 2 -> T^00_ph + * buff_idx + 3 -> T^0x_ph + * buff_idx + 4 -> T^0y_ph + * buff_idx + 5 -> T^0z_ph + */ + auto compute_t_munu = ComputePhotonTmunu {}; + arch::ComputeMomentWithSpeciesNew( + params, + domain, + { 3 }, + buff, + { (idx_t)((buff_idx + 1) % 6), + (idx_t)((buff_idx + 2) % 6), + (idx_t)((buff_idx + 3) % 6), + (idx_t)((buff_idx + 4) % 6), + (idx_t)((buff_idx + 5) % 6) }, + compute_t_munu, + smoothing_order); + /** + * buff_idx -> T_ph + */ + Kokkos::parallel_for( + "ComputeT_ph", + domain.mesh.rangeActiveCells(), + ComputePhotonTemperature { buff, + (idx_t)buff_idx, + (idx_t)((buff_idx + 1) % 6), + (idx_t)((buff_idx + 2) % 6), + (idx_t)((buff_idx + 3) % 6), + (idx_t)((buff_idx + 4) % 6), + (idx_t)((buff_idx + 5) % 6) }); + } + } + }; +} // namespace user + +#endif // PROBLEM_GENERATOR_H diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h index 415da851e..b12d4106b 100644 --- a/src/archetypes/particle_injector.h +++ b/src/archetypes/particle_injector.h @@ -326,12 +326,8 @@ namespace arch { "Weights must be used for non-Cartesian coordinates", HERE); raise::ErrorIf( - params.template get("particles.use_weights") and not use_weights, - "Weights are enabled in the input but not enabled in the injector", - HERE); - raise::ErrorIf( - not params.template get("particles.use_weights") and use_weights, - "Weights are not enabled in the input but enabled in the injector", + params.template get("particles.use_weights") != use_weights, + "Mismatch between use_weights in the input file and the injector", HERE); if (domain.species[species.first - 1].charge() + domain.species[species.second - 1].charge() != @@ -386,6 +382,150 @@ namespace arch { } } + /** + * @brief Injects particles based on spatial distribution function + * @param params Simulation parameters + * @param domain Local domain object + * @param species Species index + * @param energy_dist Energy distribution class + * @param spatial_dist Spatial distribution class + * @param number_density Number density (in units of n0) + * @param use_weights Use weights + * @param box Region to inject the particles in + * @tparam S Simulation engine type + * @tparam M Metric type + * @tparam ED Energy distribution type + * @tparam SD Spatial distribution type + */ + template ED, SpatialDistClass SD> + inline void InjectNonUniform(const SimulationParams& params, + Domain& domain, + spidx_t species, + const ED& energy_dist, + const SD& spatial_dist, + real_t number_density, + bool use_weights = (M::CoordType != Coord::Cartesian), + const boundaries_t& box = {}) { + raise::ErrorIf((M::CoordType != Coord::Cartesian) && (not use_weights), + "Weights must be used for non-Cartesian coordinates", + HERE); + raise::ErrorIf( + params.template get("particles.use_weights") != use_weights, + "Mismatch between use_weights in the input file and the injector", + HERE); + { + range_t cell_range; + if (box.empty()) { + cell_range = domain.mesh.rangeActiveCells(); + } else { + boundaries_t reduced_box(box); + if (reduced_box.size() > M::Dim) { + reduced_box.resize(M::Dim); + } + raise::ErrorIf(reduced_box.size() != M::Dim, + "Box must have the same dimension as the mesh", + HERE); + boundaries_t incl_ghosts; + for (auto d = 0; d < M::Dim; ++d) { + incl_ghosts.emplace_back(false, false); + } + const auto extent = domain.mesh.ExtentToRange(reduced_box, incl_ghosts); + tuple_t x_min { 0 }, x_max { 0 }; + for (auto d = 0; d < M::Dim; ++d) { + x_min[d] = extent[d].first; + x_max[d] = extent[d].second; + } + cell_range = CreateRangePolicy(x_min, x_max); + } + const auto ppc = number_density * + params.template get("particles.ppc0") * HALF; + auto injector_kernel = kernel::SingleSpeciesNonUniformInjector_kernel( + ppc, + domain.species[species - 1], + domain.index(), + domain.mesh.metric, + energy_dist, + spatial_dist, + ONE / params.template get("scales.V0"), + domain.random_pool()); + Kokkos::parallel_for("InjectSingleSpeciesNonUniformNumberDensity", + cell_range, + injector_kernel); + const auto n_inj = injector_kernel.number_injected(); + domain.species[species - 1].set_npart( + domain.species[species - 1].npart() + n_inj); + domain.species[species - 1].set_counter( + domain.species[species - 1].counter() + n_inj); + } + } + + /** + * @brief Injects uniform number density of a single species everywhere in the domain + * @param domain Domain object + * @param species Species index + * @param energy_dist Energy distribution objects + * @param number_density Number density (in units of n0) + * @param use_weights Use weights + * @param box Region to inject the particles in global coords + * @tparam S Simulation engine type + * @tparam M Metric type + * @tparam ED Energy distribution type + */ + template ED> + inline void InjectUniform(const SimulationParams& params, + Domain& domain, + spidx_t species, + const ED& energy_dist, + real_t number_density, + bool use_weights = false, + const boundaries_t& box = {}) { + raise::ErrorIf((M::CoordType != Coord::Cartesian) && (not use_weights), + "Weights must be used for non-Cartesian coordinates", + HERE); + raise::ErrorIf((M::CoordType == Coord::Cartesian) && use_weights, + "Weights should not be used for Cartesian coordinates", + HERE); + raise::ErrorIf(params.template get("particles.use_weights") != use_weights, + "Weights must be enabled from the input file to use them in " + "the injector", + HERE); + if (domain.species[species - 1].charge() != 0.0f) { + raise::Warning("Charge of the injected species is non-zero", HERE); + } + + { + boundaries_t nonempty_box; + for (auto d { 0u }; d < M::Dim; ++d) { + if (d < box.size()) { + nonempty_box.emplace_back(box[d].first, box[d].second); + } else { + nonempty_box.push_back(Range::All); + } + } + const auto result = ComputeNumInject(params, domain, number_density, nonempty_box); + if (not std::get<0>(result)) { + return; + } + const auto nparticles = std::get<1>(result); + const auto xi_min = std::get<2>(result); + const auto xi_max = std::get<3>(result); + + Kokkos::parallel_for("InjectUniform", + nparticles, + kernel::SingleSpeciesUniformInjector_kernel( + domain.species[species - 1], + domain.index(), + domain.mesh.metric, + xi_min, + xi_max, + energy_dist, + ONE / params.template get("scales.V0"), + domain.random_pool())); + domain.species[species - 1].set_npart( + domain.species[species - 1].npart() + nparticles); + } + } + } // namespace arch #endif // ARCHETYPES_PARTICLE_INJECTOR_H diff --git a/src/archetypes/qed/compton.h b/src/archetypes/qed/compton.h new file mode 100644 index 000000000..8776d4a12 --- /dev/null +++ b/src/archetypes/qed/compton.h @@ -0,0 +1,305 @@ +/** + * @file archetypes/qed/compton.h + * @brief Two-body collision policy of Compton scattering between leptons and photons + * @implements + * - arch::qed::ComptonScattering<> + * @namespaces: + * - arch::qed:: + */ +#ifndef ARCHETYPES_QED_COMPTON_H +#define ARCHETYPES_QED_COMPTON_H + +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/comparators.h" +#include "utils/error.h" +#include "utils/numeric.h" +#include "utils/param_container.h" + +#include "framework/containers/particles.h" + +#include + +namespace arch::qed { + using namespace ntt; + + template + struct ComptonScattering { + static constexpr spidx_t MAXSP = 16u; + static constexpr int MAX_ITER = 10; + + ParticleArrays species[MAXSP]; + static constexpr real_t low_energy_limit = static_cast(2e-3); + static constexpr real_t Thomson_limit = static_cast(1e-3); + + const real_t nominal_probability_density; + random_number_pool_t random_pool; + + ComptonScattering(const prm::Parameters& params, + random_number_pool_t& random_pool) + : nominal_probability_density { params.template get( + "compton_scattering.nominal_probability_density") } + , random_pool { random_pool } { + if (nominal_probability_density <= ZERO) { + raise::Error("nominal_probability must be > 0", HERE); + } + } + + /* + * Lorentz boost a 4-momentum p of the photon to the frame moving with u + * @param u: 4-velocity of the boost frame + * @param gamma: Lorentz factor of the boost frame + * @param p: 4-momentum of the photon in the lab frame + * @param e: energy of the photon in the lab frame + * @return: 4-momentum of the photon in the boost frame + * @return: energy of the photon in the boost frame + */ + Inline void LorentzBoost(const vec_t& u, + real_t gamma, + const vec_t& p, + real_t e, + vec_t& p_, + real_t& e_) const { + const auto u_dot_p = DOT(u[0], u[1], u[2], p[0], p[1], p[2]); + + e_ = gamma * e - u_dot_p; + p_[0] = p[0] + (u_dot_p / (ONE + gamma) - e) * u[0]; + p_[1] = p[1] + (u_dot_p / (ONE + gamma) - e) * u[1]; + p_[2] = p[2] + (u_dot_p / (ONE + gamma) - e) * u[2]; + } + + /* + * Calculate the Klein-Nishina cross section for a photon with energy e_ in the lepton rest frame + * @param e_: photon energy in the lepton rest frame + * @return f_KN: the Klein-Nishina cross section normalized to the Thomson cross section + * @note for e_ > low_energy_limit, full Klein-Nishina formula + * @note for Thomson_limit < e_ <= low_energy_limit, 2nd order expansion of the Klein-Nishina formula + * @note for e_ <= Thomson_limit, return 1 (Thomson limit) + */ + Inline auto KNCrossSection(double e_) const -> real_t { + if (e_ > Thomson_limit) { + if (e_ < low_energy_limit) { + // correctly handle the e_ << 1 limit using 2nd order expansion of f_KN + return static_cast(1.0 - 2.0 * e_ + 5.2 * SQR(e_)); + } else { + return static_cast( + 0.375 * + ((1.0 - 2.0 / e_ - 2.0 / SQR(e_)) * math::log(1.0 + 2.0 * e_) + + 0.5 + 4.0 / e_ - 0.5 / SQR(1.0 + 2.0 * e_)) / + e_); + } + } else { + return ONE; + } + } + + /* + * Sample a cosine theta value from the Thomson scattering cross section + */ + Inline auto RandomCosTheta_Th() const -> real_t { + auto gen_ = random_pool.get_state(); + const auto rnd_ = Random(gen_); + random_pool.free_state(gen_); + const auto u = math::pow( + FOUR * rnd_ - TWO + + math::sqrt(FIVE + static_cast(16) * rnd_ * (rnd_ - ONE)), + THIRD); + return u - ONE / u; + } + + /* + * Sample a cosine theta from the Klein-Nishina scattering cross section + */ + Inline auto RandomCosTheta_KN(double e_) const -> real_t { + auto gen_ = random_pool.get_state(); + const auto rnd_ = Random(gen_); + random_pool.free_state(gen_); + + auto u = 2.0 * rnd_ - 1.0; + bool converged = false; + for (int iter = 0; iter < MAX_ITER; ++iter) { + const auto CDF = (-((2.0 + e_ * (4.0 + e_ - 4.0 * (-1.0 + u) * u * e_ + + 2.0 * CUBE(-1.0 + u) * SQR(e_))) / + SQR(1.0 + e_ - u * e_)) + + (2.0 + e_ * (4.0 - e_ * (7.0 + 16.0 * e_))) / + SQR(1.0 + 2.0 * e_) + + 2.0 * (-2.0 + (-2.0 + e_) * e_) * + math::log((1.0 + e_ - u * e_) / (1.0 + 2.0 * e_))) / + ((-4.0 * e_ * (2.0 + e_ * (1.0 + e_) * (8.0 + e_))) / + SQR(1.0 + 2.0 * e_) + + (4.0 - 2.0 * (-2.0 + e_) * e_) * + math::log(1.0 + 2.0 * e_)); + const auto dCDF_du = -((CUBE(e_) * SQR(1.0 + 2.0 * e_) * + (1.0 + SQR(u) - (-1.0 + u) * (1.0 + SQR(u)) * e_ + + SQR(-1.0 + u) * SQR(e_))) / + (CUBE(-1.0 + (-1.0 + u) * e_) * + (2.0 * e_ * (2.0 + e_ * (1.0 + e_) * (8.0 + e_)) + + SQR(1.0 + 2.0 * e_) * (-2.0 + (-2.0 + e_) * e_) * + math::log(1.0 + 2.0 * e_)))); + + const auto du = (rnd_ - CDF) / dCDF_du; + + u += du; + if (u > 1.0) { + u = 1.0; + } else if (u < -1.0) { + u = -1.0; + } + if (math::abs(du) < 1e-3) { + converged = true; + break; + } + } // iterative loop for u + return static_cast(u); + } + + /* + * Scatter a photon with initial momentum p_ and energy e_ in the lepton + * rest frame to a new momentum pnew_ and energy enew_ + * @param KN_regime: whether the photon energy is in the Klein-Nishina regime + * @param p_: initial photon momentum in the lepton rest frame + * @param e_: initial photon energy in the lepton rest frame + * @return pnew_: output photon momentum after scattering in the lepton rest frame + * @return enew_: output photon energy after scattering in the lepton rest frame + * @note the scattering angle is sampled from the Klein-Nishina differential + * cross section if KN_regime is true, otherwise it is sampled from the Thomson limit + */ + Inline void ScatterPhoton(bool KN_regime, + const vec_t& p_, + real_t e_, + vec_t& pnew_, + real_t& enew_) const { + auto rand_costheta_ { ZERO }; + if (not KN_regime) { + rand_costheta_ = RandomCosTheta_Th(); + } else { + rand_costheta_ = RandomCosTheta_KN(e_); + } + const auto rand_sintheta_ = math::sqrt(ONE - SQR(rand_costheta_)); + + auto gen_ = random_pool.get_state(); + const auto rand_phi_ = static_cast(constant::TWO_PI) * + Random(gen_); + random_pool.free_state(gen_); + const auto rand_cosphi_ = math::cos(rand_phi_); + const auto rand_sinphi_ = math::sin(rand_phi_); + + // Define an orthonormal basis: {a_, b_, c_} in the lepton frame + const vec_t a_ { p_[0] / e_, p_[1] / e_, p_[2] / e_ }; + vec_t b_ { ONE, ZERO, ZERO }; + if (not cmp::AlmostZero(a_[0])) { + b_[0] = -a_[1] / a_[0]; + b_[1] = ONE / math::sqrt(ONE + SQR(b_[0])); + b_[0] /= math::sqrt(ONE + SQR(b_[0])); + } + const vec_t c_ { + CROSS_x1(a_[0], a_[1], a_[2], b_[0], b_[1], b_[2]), + CROSS_x2(a_[0], a_[1], a_[2], b_[0], b_[1], b_[2]), + CROSS_x3(a_[0], a_[1], a_[2], b_[0], b_[1], b_[2]) + }; + + enew_ = e_ / (ONE + e_ * (ONE - rand_costheta_)); + + pnew_[0] = enew_ * (rand_costheta_ * a_[0] + + rand_sintheta_ * rand_cosphi_ * b_[0] + + rand_sintheta_ * rand_sinphi_ * c_[0]); + pnew_[1] = enew_ * (rand_costheta_ * a_[1] + + rand_sintheta_ * rand_cosphi_ * b_[1] + + rand_sintheta_ * rand_sinphi_ * c_[1]); + pnew_[2] = enew_ * (rand_costheta_ * a_[2] + + rand_sintheta_ * rand_cosphi_ * b_[2] + + rand_sintheta_ * rand_sinphi_ * c_[2]); + } + + Inline auto should_interact(spidx_t sp1, + npart_t p1, + spidx_t sp2, + npart_t p2, + real_t tile_weight) const -> bool { + // values with "_" are in the lepton rest-frame + const vec_t lepton_u { species[sp1 - 1].ux1(p1), + species[sp1 - 1].ux2(p1), + species[sp1 - 1].ux3(p1) }; + const auto lepton_gamma = U2GAMMA(lepton_u[0], lepton_u[1], lepton_u[2]); + const auto lepton_weight = species[sp1 - 1].weight(p1); + + const vec_t photon_p { species[sp2 - 1].ux1(p2), + species[sp2 - 1].ux2(p2), + species[sp2 - 1].ux3(p2) }; + const auto photon_energy = NORM(photon_p[0], photon_p[1], photon_p[2]); + const auto photon_weight = species[sp2 - 1].weight(p2); + + // boost photon momentum to lepton rest frame + vec_t photon_p_ { ZERO, ZERO, ZERO }; + real_t photon_energy_ { ZERO }; + + LorentzBoost(lepton_u, lepton_gamma, photon_p, photon_energy, photon_p_, photon_energy_); + + const auto f_KN = KNCrossSection(photon_energy_); + + auto gen = random_pool.get_state(); + const auto rnd = Random(gen); + random_pool.free_state(gen); + + return rnd < + (tile_weight * nominal_probability_density * f_KN * photon_energy_ * + lepton_weight * photon_weight / (photon_energy * lepton_gamma)); + } + + Inline void operator()(spidx_t sp1, npart_t p1, spidx_t sp2, npart_t p2) const { + // @TODO coord/vec conversion + // values with "_" are in the lepton rest-frame + const vec_t lepton_u { species[sp1 - 1].ux1(p1), + species[sp1 - 1].ux2(p1), + species[sp1 - 1].ux3(p1) }; + const auto lepton_gamma = U2GAMMA(lepton_u[0], lepton_u[1], lepton_u[2]); + const auto lepton_weight = species[sp1 - 1].weight(p1); + + const vec_t photon_p { species[sp2 - 1].ux1(p2), + species[sp2 - 1].ux2(p2), + species[sp2 - 1].ux3(p2) }; + const auto photon_energy = NORM(photon_p[0], photon_p[1], photon_p[2]); + const auto photon_weight = species[sp2 - 1].weight(p2); + + // boost photon momentum to lepton rest frame + vec_t photon_p_ { ZERO, ZERO, ZERO }; + real_t photon_energy_ { ZERO }; + + LorentzBoost(lepton_u, lepton_gamma, photon_p, photon_energy, photon_p_, photon_energy_); + + vec_t photon_pnew_ { ZERO, ZERO, ZERO }, + photon_pnew { ZERO, ZERO, ZERO }; + real_t photon_energy_new_ { ZERO }, photon_energy_new { ZERO }; + + ScatterPhoton(photon_energy_ > Thomson_limit, + photon_p_, + photon_energy_, + photon_pnew_, + photon_energy_new_); + LorentzBoost({ -lepton_u[0], -lepton_u[1], -lepton_u[2] }, + lepton_gamma, + photon_pnew_, + photon_energy_new_, + photon_pnew, + photon_energy_new); + if constexpr (R1) { + species[sp1 - 1].ux1(p1) += (photon_p[0] - photon_pnew[0]) * + photon_weight / lepton_weight; + species[sp1 - 1].ux2(p1) += (photon_p[1] - photon_pnew[1]) * + photon_weight / lepton_weight; + species[sp1 - 1].ux3(p1) += (photon_p[2] - photon_pnew[2]) * + photon_weight / lepton_weight; + } + + if constexpr (R2) { + species[sp2 - 1].ux1(p2) = photon_pnew[0]; + species[sp2 - 1].ux2(p2) = photon_pnew[1]; + species[sp2 - 1].ux3(p2) = photon_pnew[2]; + } + } + }; + +} // namespace arch::qed + +#endif // ARCHETYPES_QED_COMPTON_H diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h index 1ae775309..405459091 100644 --- a/src/archetypes/utils.h +++ b/src/archetypes/utils.h @@ -146,7 +146,17 @@ namespace arch { const auto inv_n0 = ONE / params.template get("scales.n0"); const auto use_weights = params.template get("particles.use_weights"); - Kokkos::deep_copy(buffer, ZERO); + if constexpr (M::Dim == Dim::_1D) { + Kokkos::deep_copy(Kokkos::subview(buffer, Kokkos::ALL(), buffer_idx), ZERO); + } else if constexpr (M::Dim == Dim::_2D) { + Kokkos::deep_copy( + Kokkos::subview(buffer, Kokkos::ALL, Kokkos::ALL, buffer_idx), + ZERO); + } else if constexpr (M::Dim == Dim::_3D) { + Kokkos::deep_copy( + Kokkos::subview(buffer, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL, buffer_idx), + ZERO); + } auto scatter_buff = Kokkos::Experimental::create_scatter_view(buffer); for (const auto sp : species) { const auto& prtl_spec = domain.species[sp - 1]; @@ -168,6 +178,53 @@ namespace arch { Kokkos::Experimental::contribute(buffer, scatter_buff); } + template + inline void ComputeMomentWithSpeciesNew( + const SimulationParams& params, + const Domain& domain, + const std::vector& species, + ndfield_t& buffer, + const Kokkos::Array& buff_indices, + const MF& func, + uint8_t smoothing_order = 0u, + OutputSmoothingTypeFlag smoothing_method = OutputSmoothingType::SPLINE) { + const auto ni2 = domain.mesh.n_active(in::x2); + const auto inv_n0 = ONE / params.template get("scales.n0"); + const auto use_weights = params.template get("particles.use_weights"); + + for (const auto idx : buff_indices) { + if constexpr (M::Dim == Dim::_1D) { + Kokkos::deep_copy(Kokkos::subview(buffer, Kokkos::ALL(), idx), ZERO); + } else if constexpr (M::Dim == Dim::_2D) { + Kokkos::deep_copy(Kokkos::subview(buffer, Kokkos::ALL, Kokkos::ALL, idx), + ZERO); + } else if constexpr (M::Dim == Dim::_3D) { + Kokkos::deep_copy( + Kokkos::subview(buffer, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL, idx), + ZERO); + } + } + auto scatter_buff = Kokkos::Experimental::create_scatter_view(buffer); + for (const auto sp : species) { + const auto& prtl_spec = domain.species[sp - 1]; + Kokkos::parallel_for( + "ComputeMoment", + prtl_spec.rangeActiveParticles(), + kernel::ParticleMomentsNew_kernel(func, + scatter_buff, + buff_indices, + prtl_spec, + use_weights, + domain.mesh.metric, + domain.mesh.flds_bc(), + ni2, + inv_n0, + smoothing_order, + smoothing_method)); + } + Kokkos::Experimental::contribute(buffer, scatter_buff); + } + template inline void UpdateEMFields(Domain& domain, const F& fieldsetter) { if constexpr (S == SimEngine::SRPIC) { diff --git a/src/engines/engine.hpp b/src/engines/engine.hpp index 13b5e8e28..642998139 100644 --- a/src/engines/engine.hpp +++ b/src/engines/engine.hpp @@ -130,6 +130,7 @@ namespace ntt { auto parameters = prm::Parameters {}; parameters.set("dt", static_cast(dt)); parameters.set("time", static_cast(time)); + parameters.set("step", static_cast(step)); parameters.set("tiled_deposit_team_size", static_cast(tiled_deposit_team_size)); return parameters; @@ -254,6 +255,7 @@ namespace ntt { "ParticlePusher", "FieldBoundaries", "ParticleBoundaries", "Communications", "Injector", "Custom", + "TwoBodyInteractions", "ParticleSort", "LoadBalance", "ParticleSort", "Output", "Render", "Checkpoint" }, diff --git a/src/engines/reporter.cpp b/src/engines/reporter.cpp index 7752840f1..6cef73924 100644 --- a/src/engines/reporter.cpp +++ b/src/engines/reporter.cpp @@ -6,8 +6,10 @@ #include "utils/formatting.h" #include "utils/reporter.h" +#include "framework/parameters/extra.h" #include "framework/parameters/parameters.h" +#include #include #include @@ -166,6 +168,48 @@ namespace ntt { params.template get( "radiation.emission.synchrotron.nominal_photon_energy")); } + + report += "\n"; + const auto two_body_interactions = + params.template get>( + "two_body.interaction"); + if (not two_body_interactions.empty()) { + reporter::AddCategory(report, 4, "Two-body interactions"); + reporter::AddParam( + report, + 6, + "Thomson optical depth", + "%.3e", + params.template get("two_body.thomson_optical_depth")); + for (const auto& interaction : two_body_interactions) { + reporter::AddSubcategory( + report, + 6, + TwoBodyInteraction::to_string(interaction.type).c_str()); + std::vector group1(interaction.group1.size()); + std::vector group2(interaction.group2.size()); + for (size_t g1 = 0; g1 < interaction.group1.size(); ++g1) { + group1[g1] = interaction.group1[g1]; + } + for (size_t g2 = 0; g2 < interaction.group2.size(); ++g2) { + group2[g2] = interaction.group2[g2]; + } + reporter::AddParam(report, + 8, + "group #1 species", + "%s (recoil: %s)", + fmt::formatVector(group1).c_str(), + interaction.recoil1 ? "ON" : "OFF"); + reporter::AddParam(report, + 8, + "group #2 species", + "%s (recoil: %s)", + fmt::formatVector(group2).c_str(), + interaction.recoil2 ? "ON" : "OFF"); + reporter::AddParam(report, 8, "tile size [cells]", "%u", interaction.tile_size); + reporter::AddParam(report, 8, "interval [steps]", "%u", interaction.interval); + } + } return report; } diff --git a/src/engines/srpic/srpic.hpp b/src/engines/srpic/srpic.hpp index d0352318e..42c186c33 100644 --- a/src/engines/srpic/srpic.hpp +++ b/src/engines/srpic/srpic.hpp @@ -24,6 +24,7 @@ #include "engines/srpic/fieldsolvers.h" #include "engines/srpic/particle_pusher.h" #include "engines/srpic/particles_bcs.h" +#include "engines/srpic/twobody.h" #include "framework/domain/domain.h" #include "framework/parameters/parameters.h" @@ -182,6 +183,15 @@ namespace ntt { timers.stop("Injector"); } + if constexpr (CartesianMetricClass) { + timers.start("TwoBodyInteractions"); + srpic::TwoBodyInteractions(dom, this->engineParams(), m_params); + timers.stop("TwoBodyInteractions"); + } + + timers.start("ParticleSort"); + m_metadomain.SortParticles(time, step, m_params, dom); + timers.stop("ParticleSort"); // NOTE: particle sorting is intentionally NOT done here. It runs once per // step in the engine loop (Engine::run) after CustomPostStep and // LoadBalance, so the layout the next deposit uses reflects window diff --git a/src/engines/srpic/twobody.h b/src/engines/srpic/twobody.h new file mode 100644 index 000000000..8019cd20b --- /dev/null +++ b/src/engines/srpic/twobody.h @@ -0,0 +1,99 @@ +#ifndef ENGINES_SRPIC_TWOBODY_H +#define ENGINES_SRPIC_TWOBODY_H + +#include "enums.h" + +#include "traits/metric.h" +#include "utils/error.h" +#include "utils/formatting.h" +#include "utils/log.h" +#include "utils/param_container.h" + +#include "archetypes/qed/compton.h" +#include "framework/domain/domain.h" +#include "framework/parameters/extra.h" +#include "framework/parameters/parameters.h" +#include "kernels/twobody_interactions.hpp" + +namespace ntt { + namespace srpic { + + template + void TwoBodyInteractions(Domain& domain, + const prm::Parameters& engine_params, + const SimulationParams& params) { + logger::Checkpoint("Launching TwoBodyInteractions routines", HERE); + const auto dt = engine_params.get("dt"); + const auto step = engine_params.get("step"); + for (const auto& interaction : + params.template get>( + "two_body.interaction")) { + if (step % interaction.interval == 0u) { + const auto thomson_optical_depth = params.template get( + "two_body.thomson_optical_depth"); + const auto nominal_thomson_probability_density = thomson_optical_depth * + dt * + static_cast( + interaction.interval); + if (interaction.type == TwoBodyInteraction::COMPTON) { + prm::Parameters compton_params; + compton_params.set("compton_scattering.nominal_probability_density", + nominal_thomson_probability_density); + auto launch = [&]() { + auto policy = arch::qed::ComptonScattering( + compton_params, + domain.random_pool()); + + std::vector*> group1_species; + std::vector*> group2_species; + + for (const auto& sp_lepton : interaction.group1) { + raise::ErrorIf( + domain.species[sp_lepton - 1].mass() == ZERO, + fmt::format( + "Species %u is massless but is in the lepton group " + "of a Compton interaction", + sp_lepton), + HERE); + group1_species.push_back(&domain.species[sp_lepton - 1]); + } + for (const auto& sp_photon : interaction.group2) { + raise::ErrorIf( + domain.species[sp_photon - 1].mass() != ZERO, + fmt::format( + "Species %u is massive but is in the photon group " + "of a Compton interaction", + sp_photon), + HERE); + group2_species.push_back(&domain.species[sp_photon - 1]); + } + + kernel::mink::TwoBodyInteraction( + group1_species, + group2_species, + domain.mesh.n_active(), + domain.mesh.extent(), + interaction.tile_size, + params.template get("particles.ppc0"), + domain.random_pool(), + policy); + }; + if (interaction.recoil1 and interaction.recoil2) { + launch.template operator()(); + } else if (interaction.recoil1 and not interaction.recoil2) { + launch.template operator()(); + } else if (not interaction.recoil1 and interaction.recoil2) { + launch.template operator()(); + } else { + launch.template operator()(); + } + } else if (interaction.type == TwoBodyInteraction::CUSTOM) { + raise::Error("Custom two-body interactions not implemented yet", HERE); + } + } + } + } + } // namespace srpic +} // namespace ntt + +#endif // ENGINES_SRPIC_TWOBODY_H diff --git a/src/framework/parameters/extra.cpp b/src/framework/parameters/extra.cpp index dbee39c7c..f0d9d127f 100644 --- a/src/framework/parameters/extra.cpp +++ b/src/framework/parameters/extra.cpp @@ -5,12 +5,14 @@ #include "utils/numeric.h" +#include "framework/parameters/extra.h" #include "framework/parameters/parameters.h" #include #include #include +#include namespace ntt { namespace params { @@ -114,6 +116,29 @@ namespace ntt { compton_photon_weight.value(); compton_nominal_photon_energy = ONE / SQR(compton_gamma_qed.value()); } + + twobody_thomson_optical_depth = toml::find_or( + toml_data, + "two_body", + "thomson_optical_depth", + defaults::twobody::thomson_optical_depth); + + // find two-body interactions + const auto twobody_tab = toml::find_or(toml_data, + "two_body", + "interaction", + toml::array {}); + for (const auto& tbint : twobody_tab) { + twobody_interactions.push_back(TwoBodyInteractionParams { + .type = TwoBodyInteraction::from_string( + toml::find(tbint, "type")), + .group1 = toml::find>(tbint, "group1"), + .group2 = toml::find_or>(tbint, "group2", {}), + .interval = toml::find_or(tbint, "interval", 1), + .tile_size = toml::find_or(tbint, "tile_size", 4u), + .recoil1 = toml::find_or(tbint, "recoil1", true), + .recoil2 = toml::find_or(tbint, "recoil2", true) }); + } } void Extra::setParams(const std::map& extra, @@ -159,6 +184,10 @@ namespace ntt { params->set("radiation.emission.compton.nominal_photon_energy", compton_nominal_photon_energy.value()); } + + params->set("two_body.thomson_optical_depth", + twobody_thomson_optical_depth.value()); + params->set("two_body.interaction", twobody_interactions); } } // namespace params } // namespace ntt diff --git a/src/framework/parameters/extra.h b/src/framework/parameters/extra.h index 29de4527a..c3891429f 100644 --- a/src/framework/parameters/extra.h +++ b/src/framework/parameters/extra.h @@ -12,6 +12,7 @@ #ifndef FRAMEWORK_PARAMETERS_EXTRA_H #define FRAMEWORK_PARAMETERS_EXTRA_H +#include "enums.h" #include "global.h" #include "framework/parameters/parameters.h" @@ -25,6 +26,16 @@ namespace ntt { namespace params { + struct TwoBodyInteractionParams { + TwoBodyInteractionFlag type; + std::vector group1; + std::vector group2; + timestep_t interval; + ncells_t tile_size; + bool recoil1; + bool recoil2; + }; + struct Extra { // radiative drag parameters std::optional synchrotron_gamma_rad; @@ -45,6 +56,10 @@ namespace ntt { std::optional compton_nominal_probability; std::optional compton_nominal_photon_energy; + // two-body interaction parameters + std::optional twobody_thomson_optical_depth; + std::vector twobody_interactions; + void read(const std::map&, const toml::value&, const SimulationParams* const); diff --git a/src/global/defaults.h b/src/global/defaults.h index 101e54659..a0d64af4c 100644 --- a/src/global/defaults.h +++ b/src/global/defaults.h @@ -119,6 +119,10 @@ namespace ntt::defaults { const real_t gamma_rad = 1.0; const real_t gamma_qed = 10.0; } // namespace compton + + namespace twobody { + const real_t thomson_optical_depth = 1.0; + } // namespace twobody } // namespace ntt::defaults #endif // GLOBAL_DEFAULTS_H diff --git a/src/global/enums.h b/src/global/enums.h index 5700a3a28..b88f2da95 100644 --- a/src/global/enums.h +++ b/src/global/enums.h @@ -402,6 +402,43 @@ namespace ntt { using EmissionTypeFlag = uint8_t; + namespace TwoBodyInteraction { + enum TwoBodyInteractionFlag_ : uint8_t { + NONE = 0, + COMPTON = 1, + CUSTOM = 2, + }; + + inline auto to_string(uint8_t flags) -> std::string { + switch (flags) { + case NONE: + return "none"; + case COMPTON: + return "compton"; + case CUSTOM: + return "custom"; + default: + return "unknown"; + } + } + + inline auto from_string(const std::string& s) -> uint8_t { + if (fmt::toLower(s) == "none") { + return NONE; + } else if (fmt::toLower(s) == "compton") { + return COMPTON; + } else if (fmt::toLower(s) == "custom") { + return CUSTOM; + } else { + raise::Error(fmt::format("Invalid TwoBodyInteraction type: %s", s.c_str()), + HERE); + return NONE; + } + } + } // namespace TwoBodyInteraction + + using TwoBodyInteractionFlag = uint8_t; + namespace OutputSmoothingType { enum OutputSmoothingTypeFlag_ : uint8_t { NONE = 0, diff --git a/src/global/traits/pgen.h b/src/global/traits/pgen.h index 179c2dbc9..0f84699dd 100644 --- a/src/global/traits/pgen.h +++ b/src/global/traits/pgen.h @@ -134,4 +134,4 @@ namespace traits::pgen { } // namespace traits::pgen -#endif // TRAITS_PGEN_H \ No newline at end of file +#endif // TRAITS_PGEN_H diff --git a/src/global/traits/policies.h b/src/global/traits/policies.h index a317fd687..8be16642c 100644 --- a/src/global/traits/policies.h +++ b/src/global/traits/policies.h @@ -141,4 +141,40 @@ concept CustomParticleUpdatePolicyClass = } or traits::custom_prtl_update::IsNoPolicy; +namespace traits::twobodyinteractions { + + template + concept HasSpecies = requires(I& interaction_policy) { + { interaction_policy.species } -> std::convertible_to; + }; + + template + concept HasShouldInteract = requires(const I& interaction_policy, + spidx_t sp1, + npart_t p1, + spidx_t sp2, + npart_t p2, + real_t tile_weight) { + { + interaction_policy.should_interact(sp1, p1, sp2, p2, tile_weight) + } -> std::same_as; + }; + + template + concept HasInteraction = requires(const I& interaction_policy, + spidx_t sp1, + npart_t p1, + spidx_t sp2, + npart_t p2) { + { interaction_policy(sp1, p1, sp2, p2) } -> std::same_as; + }; + +} // namespace traits::twobodyinteractions + +template +concept TwoBodyInteractionPolicyClass = + traits::twobodyinteractions::HasSpecies and + traits::twobodyinteractions::HasShouldInteract and + traits::twobodyinteractions::HasInteraction; + #endif // TRAITS_POLICIES_H diff --git a/src/global/utils/sorting.h b/src/global/utils/sorting.h index 006c0bf27..845797c8b 100644 --- a/src/global/utils/sorting.h +++ b/src/global/utils/sorting.h @@ -147,6 +147,9 @@ namespace sort { total_tiles *= ntx3; } if constexpr (Count) { + if (num_ppt.extent(0) == 0u) { + num_ppt = array_t { "num_ppt", total_tiles }; + } raise::ErrorIf(num_ppt.extent(0) != total_tiles, "num_ppt must have extent equal to total tiles", HERE); diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp index 5d7904cd9..49ef50a21 100644 --- a/src/kernels/injectors.hpp +++ b/src/kernels/injectors.hpp @@ -5,6 +5,8 @@ * - kernel::UniformInjector_kernel<> * - kernel::GlobalInjector_kernel<> * - kernel::NonUniformInjector_kernel<> + * - kernel::SingleSpeciesNonUniformInjector_kernel<> + * - kernel::SingleSpeciesUniformInjector_kernel<> * @namespaces: * - kernel:: */ @@ -858,6 +860,363 @@ namespace kernel { } }; // struct NonUniformInjector_kernel + template ED, SpatialDistClass SD> + struct SingleSpeciesNonUniformInjector_kernel { + + const real_t ppc0; + + ParticleArrays particles; + + array_t idx { "idx" }; + + const npart_t offset; + const npart_t domain_idx, cntr; + const bool use_tracking; + const M metric; + const ED energy_dist; + const SD spatial_dist; + const real_t inv_V0; + random_number_pool_t random_pool; + + SingleSpeciesNonUniformInjector_kernel(real_t ppc0, + Particles& particles, + npart_t domain_idx, + const M& metric, + const ED& energy_dist, + const SD& spatial_dist, + real_t inv_V0, + random_number_pool_t& random_pool) + : ppc0 { ppc0 } + , particles { static_cast(particles) } + , offset { particles.npart() } + , domain_idx { domain_idx } + , cntr { particles.counter() } + , use_tracking { particles.use_tracking() } + , metric { metric } + , energy_dist { energy_dist } + , spatial_dist { spatial_dist } + , inv_V0 { inv_V0 } + , random_pool { random_pool } {} + + auto number_injected() const -> npart_t { + auto idx_h = Kokkos::create_mirror_view(idx); + Kokkos::deep_copy(idx_h, idx); + return idx_h(); + } + + Inline auto injected_ppc(const coord_t& x_Ph) const + -> Kokkos::pair { + real_t ppc_real = ppc0, weight = ONE; + if constexpr (SimpleSpatialDistClass) { + ppc_real *= spatial_dist(x_Ph); + } else { + const auto sp_dist = spatial_dist(x_Ph); + ppc_real *= sp_dist.first; + weight = sp_dist.second; + } + auto ppc = static_cast(ppc_real); + auto rand_gen = random_pool.get_state(); + if (Random(rand_gen) < (ppc_real - static_cast(ppc))) { + ppc += 1; + } + random_pool.free_state(rand_gen); + return { ppc, weight }; + } + + Inline void inject(const prtlidx_t index, + const tuple_t& xi_Cd, + const tuple_t& dxi_Cd, + const vec_t& v_Cd, + const real_t weight) const { + // clang-format off + if (not use_tracking) { + InjectParticle(index + offset, + particles.i1, particles.i2, particles.i3, + particles.dx1, particles.dx2, particles.dx3, + particles.ux1, particles.ux2, particles.ux3, + particles.phi, particles.weight, particles.tag, particles.pld_i, + xi_Cd, dxi_Cd, v_Cd, weight, ZERO); + } else { + InjectParticle(index + offset, + particles.i1, particles.i2, particles.i3, + particles.dx1, particles.dx2, particles.dx3, + particles.ux1, particles.ux2, particles.ux3, + particles.phi, particles.weight, particles.tag, particles.pld_i, + xi_Cd, dxi_Cd, v_Cd, weight, ZERO, + domain_idx, index + cntr); + } + // clang-format on + } + + Inline void operator()(cellidx_t i1) const { + if constexpr (M::Dim == Dim::_1D) { + const auto i1_ = COORD(i1); + const coord_t x_Cd { i1_ + HALF }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + + auto [ppc, weight] = injected_ppc(x_Ph); + if (ppc == 0) { + return; + } + + if constexpr (M::CoordType != Coord::Cartesian) { + weight *= metric.sqrt_det_h({ i1_ + HALF }) * inv_V0; + } + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + random_pool.free_state(rand_gen); + + vec_t v_XYZ { ZERO }; + { + vec_t v_T { ZERO }; + energy_dist(x_Ph, v_T); + metric.template transform_xyz(x_Cd, v_T, v_XYZ); + } + inject(index, { static_cast(i1_) }, { dx1 }, v_XYZ, weight); + } + } else { + raise::KernelError( + HERE, + "SingleSpeciesNonUniformInjector_kernel 1D called for 2D/3D"); + } + } + + Inline void operator()(cellidx_t i1, cellidx_t i2) const { + if constexpr (M::Dim == Dim::_2D) { + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + const coord_t x_Cd { i1_ + HALF, i2_ + HALF }; + coord_t x_Ph { ZERO }; + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + if constexpr (S == SimEngine::SRPIC and M::CoordType != Coord::Cartesian) { + x_Cd_[2] = ZERO; + } + metric.template convert(x_Cd, x_Ph); + + auto [ppc, weight] = injected_ppc(x_Ph); + if (ppc == 0) { + return; + } + + if constexpr (M::CoordType != Coord::Cartesian) { + weight *= metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF }) * inv_V0; + } + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + random_pool.free_state(rand_gen); + + vec_t v_Cd { ZERO }; + { + vec_t v_T { ZERO }; + energy_dist(x_Ph, v_T); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd_, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd_, v_T, v_Cd); + } + } + inject(index, + { static_cast(i1_), static_cast(i2_) }, + { dx1, dx2 }, + v_Cd, + weight); + } + } + + else { + raise::KernelError( + HERE, + "SingleSpeciesNonUniformInjector_kernel 2D called for 1D/3D"); + } + } + + Inline void operator()(cellidx_t i1, cellidx_t i2, cellidx_t i3) const { + if constexpr (M::Dim == Dim::_3D) { + const auto i1_ = COORD(i1); + const auto i2_ = COORD(i2); + const auto i3_ = COORD(i3); + const coord_t x_Cd { i1_ + HALF, i2_ + HALF, i3_ + HALF }; + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + + auto [ppc, weight] = injected_ppc(x_Ph); + if (ppc == 0) { + return; + } + + if constexpr (M::CoordType != Coord::Cartesian) { + weight *= metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF }) * + inv_V0; + } + for (auto p { 0u }; p < ppc; ++p) { + const auto index = Kokkos::atomic_fetch_add(&idx(), 1); + + auto rand_gen = random_pool.get_state(); + const auto dx1 = Random(rand_gen); + const auto dx2 = Random(rand_gen); + const auto dx3 = Random(rand_gen); + random_pool.free_state(rand_gen); + + vec_t v_Cd { ZERO }; + { + vec_t v_T { ZERO }; + energy_dist(x_Ph, v_T); + if constexpr (S == SimEngine::SRPIC) { + metric.template transform_xyz(x_Cd, v_T, v_Cd); + } else if constexpr (S == SimEngine::GRPIC) { + metric.template transform(x_Cd, v_T, v_Cd); + } + } + inject( + index, + { static_cast(i1_), static_cast(i2_), static_cast(i3_) }, + { dx1, dx2, dx3 }, + v_Cd, + weight); + } + } else { + raise::KernelError( + HERE, + "SingleSpeciesNonUniformInjector_kernel 3D called for 1D/2D"); + } + } + }; // struct SingleSpeciesNonUniformInjector_kernel + + template ED> + struct SingleSpeciesUniformInjector_kernel { + + ParticleArrays particles; + + const npart_t offset; + const npart_t domain_idx, cntr; + const bool use_tracking; + const M metric; + const array_t xi_min, xi_max; + const ED energy_dist; + const real_t inv_V0; + random_number_pool_t random_pool; + + SingleSpeciesUniformInjector_kernel(Particles& particles, + npart_t domain_idx, + const M& metric, + const array_t& xi_min, + const array_t& xi_max, + const ED& energy_dist, + real_t inv_V0, + random_number_pool_t& random_pool) + : particles { particles } + , offset { particles.npart() } + , domain_idx { domain_idx } + , cntr { particles.counter() } + , use_tracking { particles.use_tracking() } + , metric { metric } + , xi_min { xi_min } + , xi_max { xi_max } + , energy_dist { energy_dist } + , inv_V0 { inv_V0 } + , random_pool { random_pool } { + if (use_tracking) { +#if !defined(MPI_ENABLED) + raise::ErrorIf(particles.pld_i.extent(1) < 1, + "Particle tracking is enabled but the " + "particle integer payload size is less " + "than 1", + HERE); +#else + raise::ErrorIf(particles.pld_i.extent(1) < 2, + "Particle tracking is enabled but the " + "particle integer payload size is less " + "than 2", + HERE); +#endif + } + } + + Inline void operator()(prtlidx_t p) const { + coord_t x_Cd { ZERO }; + tuple_t xi_Cd { 0 }; + tuple_t dxi_Cd { static_cast(0) }; + vec_t v { ZERO, ZERO, ZERO }; + { // generate a random coordinate + auto rand_gen = random_pool.get_state(); + if constexpr (M::Dim == Dim::_1D or M::Dim == Dim::_2D or + M::Dim == Dim::_3D) { + x_Cd[0] = xi_min(0) + Random(rand_gen) * (xi_max(0) - xi_min(0)); + xi_Cd[0] = static_cast(x_Cd[0]); + dxi_Cd[0] = static_cast(x_Cd[0] - xi_Cd[0]); + } + if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) { + x_Cd[1] = xi_min(1) + Random(rand_gen) * (xi_max(1) - xi_min(1)); + xi_Cd[1] = static_cast(x_Cd[1]); + xi_Cd[1] = static_cast(x_Cd[1]); + dxi_Cd[1] = static_cast(x_Cd[1] - xi_Cd[1]); + } + if constexpr (M::Dim == Dim::_3D) { + x_Cd[2] = xi_min(2) + Random(rand_gen) * (xi_max(2) - xi_min(2)); + xi_Cd[2] = static_cast(x_Cd[2]); + dxi_Cd[2] = static_cast(x_Cd[2] - xi_Cd[2]); + } + random_pool.free_state(rand_gen); + } + { // generate the velocity + coord_t x_Ph { ZERO }; + metric.template convert(x_Cd, x_Ph); + if constexpr (M::CoordType == Coord::Cartesian) { + energy_dist(x_Ph, v); + } else if constexpr (S == SimEngine::SRPIC) { + coord_t x_Cd_ { ZERO }; + x_Cd_[0] = x_Cd[0]; + x_Cd_[1] = x_Cd[1]; + x_Cd_[2] = ZERO; // phi = 0 + vec_t v_Ph { ZERO }; + energy_dist(x_Ph, v_Ph); + metric.template transform_xyz(x_Cd_, v_Ph, v); + } else if constexpr (S == SimEngine::GRPIC) { + vec_t v_Ph { ZERO, ZERO, ZERO }; + energy_dist(x_Ph, v_Ph); + metric.template transform(x_Cd, v_Ph, v); + } else { + raise::KernelError(HERE, "Unknown simulation engine"); + } + } + real_t weight = ONE; + if constexpr (M::CoordType != Coord::Cartesian) { + const auto sqrt_det_h = metric.sqrt_det_h(x_Cd); + weight = sqrt_det_h * inv_V0; + } + // clang-format off + if (not use_tracking) { + InjectParticle( + p + offset, + particles.i1, particles.i2, particles.i3, + particles.dx1, particles.dx2, particles.dx3, + particles.ux1, particles.ux2, particles.ux3, + particles.phi, particles.weight, particles.tag, particles.pld_i, + xi_Cd, dxi_Cd, v, weight, ZERO); + } else { + InjectParticle( + p + offset, + particles.i1, particles.i2, particles.i3, + particles.dx1, particles.dx2, particles.dx3, + particles.ux1, particles.ux2, particles.ux3, + particles.phi, particles.weight, particles.tag, particles.pld_i, + xi_Cd, dxi_Cd, v, weight, ZERO, + domain_idx, cntr + p); + } + // clang-format on + } + }; // struct SingleSpeciesUniformInjector_kernel + } // namespace kernel #endif // KERNELS_INJECTORS_HPP diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp index bb5e837b0..af9ecf3f6 100644 --- a/src/kernels/particle_moments.hpp +++ b/src/kernels/particle_moments.hpp @@ -3,6 +3,7 @@ * @brief Algorithm for computing different moments from particle distribution * @implements * - kernel::ParticleMoments_kernel<> + * - kernel::ParticleMomentsNew_kernel<> * - kernel::NormalizeVectorByRho_kernel<> * - kernel::Normalize4VelocityByNorm_kernel<> * - kernel::Transform4VelocitySpatialToPhysical_kernel<> @@ -31,6 +32,29 @@ namespace kernel { using namespace ntt; + namespace particle_moment_functor { + + template + concept IsValid = requires(const MF& fm, + const ParticleArrays& prtls, + float mass, + float charge, + prtlidx_t p, + list_t& buffer) { + { fm(prtls, mass, charge, p, buffer) } -> std::same_as; + }; + + template + concept HasN = requires { + { MF::N } -> std::convertible_to; + }; + + } // namespace particle_moment_functor + + template + concept ParticleMomentsFunctor = particle_moment_functor::IsValid && + particle_moment_functor::HasN; + template auto get_contrib(float mass, float charge) -> real_t { if constexpr (F == FldsID::Rho) { @@ -252,7 +276,11 @@ namespace kernel { } else { u0 = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2])); } - return (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / u0; + if (c1 > 0u) { + return (mass == ZERO ? ONE : mass) * u_Phys[c1 - 1] / u0; + } else { + return (mass == ZERO ? ONE : (mass * math::sqrt(ONE - SQR(ONE / u0)))); + } } Inline auto computeEckartVelocityFluxComponent(prtlidx_t p) const -> real_t { @@ -407,6 +435,244 @@ namespace kernel { } }; + /** + * @brief Generic moment-deposition kernel. + * + * Computes a per-particle scalar via the supplied @p Func and deposits it onto + * a scatter buffer, automatically applying volume normalization, weighting, + * smoothing (shape function) and axis reflection. @p Func is any device- + * callable object with signature + * `Inline auto operator()(const ParticleArrays&, prtlidx_t, const M&) const + * -> real_t` + * i.e. it receives the particle container, the particle index, and the metric, + * and returns the raw (un-normalized, un-weighted, un-smoothed) contribution of + * that particle to the moment. + * + * @tparam M Metric type + * @tparam N Last dimension of the buffer + * @tparam Func Per-particle contribution functor (deduced) + */ + template + class ParticleMomentsNew_kernel { + static_assert( + MF::N <= N, + "Buffer size N must be >= number of components N deposited by Func"); + static constexpr auto D = M::Dim; + + const MF func; + scatter_ndfield_t Buff; + const Kokkos::Array buff_indices; + const ParticleArrays particles; + const float mass, charge; + const bool use_weights; + const bool apply_norm; + const M metric; + const int ni2; + const real_t inv_n0; + + const uint8_t order; + const uint8_t window; + const OutputSmoothingTypeFlag smoothing; + + bool is_axis_i2min { false }, is_axis_i2max { false }; + + public: + ParticleMomentsNew_kernel( + const MF& func, + const scatter_ndfield_t& scatter_buff, + const Kokkos::Array& buff_indices, + const Particles& particles, + bool use_weights, + const M& metric, + const boundaries_t& boundaries, + ncells_t ni2, + real_t inv_n0, + uint8_t order = 0u, + OutputSmoothingTypeFlag smoothing = OutputSmoothingType::SPLINE, + bool apply_norm = true) + : func { func } + , Buff { scatter_buff } + , buff_indices { buff_indices } + , particles { static_cast(particles) } + , mass { particles.mass() } + , charge { particles.charge() } + , use_weights { use_weights } + , apply_norm { apply_norm } + , metric { metric } + , ni2 { static_cast(ni2) } + , inv_n0 { inv_n0 } + , order { order } + , smoothing { smoothing } + , window { static_cast( + math::ceil(static_cast(order) / 2.0f)) } { + raise::ErrorIf(window > N_GHOSTS, "Window size too large", HERE); + if constexpr ((M::CoordType != Coord::Cartesian) && + ((D == Dim::_2D) || (D == Dim::_3D))) { + raise::ErrorIf(boundaries.size() < 2, "boundaries defined incorrectly", HERE); + is_axis_i2min = (boundaries[1].first == FldsBC::AXIS); + is_axis_i2max = (boundaries[1].second == FldsBC::AXIS); + } + } + + Inline auto shapeFunction(real_t delta_x) const -> real_t { + if (smoothing == OutputSmoothingType::SPLINE) { + if (order == 0) { + return ONE; + } else if (order == 1) { + return prtl_shape::S1(delta_x); + } else if (order == 2) { + return prtl_shape::S2(delta_x); + } else if (order == 3) { + return prtl_shape::S3(delta_x); + } else if (order == 4) { + return prtl_shape::S4(delta_x); + } else if (order == 5) { + return prtl_shape::S5(delta_x); + } else if (order == 6) { + return prtl_shape::S6(delta_x); + } else if (order == 7) { + return prtl_shape::S7(delta_x); + } else if (order == 8) { + return prtl_shape::S8(delta_x); + } else if (order == 9) { + return prtl_shape::S9(delta_x); + } else if (order == 10) { + return prtl_shape::S10(delta_x); + } else if (order == 11) { + return prtl_shape::S11(delta_x); + } else { + raise::KernelError(HERE, "Unsupported shape function order"); + return ZERO; + } + } else if (smoothing == OutputSmoothingType::CONST) { + return ONE / (TWO * static_cast(window) + ONE); + } else { + raise::KernelError(HERE, "Unsupported smoothing method"); + return ZERO; + } + } + + Inline void operator()(prtlidx_t p) const { + if (particles.tag(p) == ParticleTag::dead) { + return; + } + list_t contributions { ZERO }; + func(particles, mass, charge, p, contributions); + for (uint8_t i = 0; i < MF::N; ++i) { + // apply volume normalization and (optionally) particle weights; + // skipped e.g. for nppc, which counts raw particles per cell + if constexpr (D == Dim::_1D) { + contributions[i] *= inv_n0 / + metric.sqrt_det_h( + { static_cast(particles.i1(p)) + HALF }); + } else if constexpr (D == Dim::_2D) { + contributions[i] *= inv_n0 / + metric.sqrt_det_h( + { static_cast(particles.i1(p)) + HALF, + static_cast(particles.i2(p)) + HALF }); + } else if constexpr (D == Dim::_3D) { + contributions[i] *= inv_n0 / + metric.sqrt_det_h( + { static_cast(particles.i1(p)) + HALF, + static_cast(particles.i2(p)) + HALF, + static_cast(particles.i3(p)) + HALF }); + } + if (use_weights) { + contributions[i] *= particles.weight(p); + } + } + auto buff_access = Buff.access(); + for (uint8_t i = 0; i < MF::N; ++i) { + const auto coeff = contributions[i]; + const auto buff_idx = buff_indices[i]; + if constexpr (D == Dim::_1D) { + for (auto di1 { -window }; di1 <= window; ++di1) { + const real_t delta_x1 = math::abs(static_cast(particles.dx1(p)) - + (static_cast(di1) + HALF)); + buff_access(particles.i1(p) + di1 + N_GHOSTS, + buff_idx) += coeff * shapeFunction(delta_x1); + } + } else if constexpr (D == Dim::_2D) { + for (auto di2 { -window }; di2 <= window; ++di2) { + for (auto di1 { -window }; di1 <= window; ++di1) { + const real_t delta_x1 = math::abs( + static_cast(particles.dx1(p)) - + (static_cast(di1) + HALF)); + const real_t delta_x2 = math::abs( + static_cast(particles.dx2(p)) - + (static_cast(di2) + HALF)); + const auto shape_coeff = shapeFunction(delta_x1) * + shapeFunction(delta_x2); + if constexpr (M::CoordType == Coord::Cartesian) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } else { + // reflect contribution at axes + if (is_axis_i2min && (particles.i2(p) + di2 < 0)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + N_GHOSTS - (particles.i2(p) + di2), + buff_idx) += coeff * shape_coeff; + } else if (is_axis_i2max && (particles.i2(p) + di2 >= ni2)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + 2 * ni2 - (particles.i2(p) + di2) + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } else { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } + } + } + } + } else if constexpr (D == Dim::_3D) { + for (auto di3 { -window }; di3 <= window; ++di3) { + for (auto di2 { -window }; di2 <= window; ++di2) { + for (auto di1 { -window }; di1 <= window; ++di1) { + const auto delta_x1 = math::abs( + static_cast(particles.dx1(p)) - + (static_cast(di1) + HALF)); + const auto delta_x2 = math::abs( + static_cast(particles.dx2(p)) - + (static_cast(di2) + HALF)); + const auto delta_x3 = math::abs( + static_cast(particles.dx3(p)) - + (static_cast(di3) + HALF)); + const auto shape_coeff = shapeFunction(delta_x1) * + shapeFunction(delta_x2) * + shapeFunction(delta_x3); + if constexpr (M::CoordType == Coord::Cartesian) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } else { + // reflect contribution at axes + if (is_axis_i2min && (particles.i2(p) + di2 < 0)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + N_GHOSTS - (particles.i2(p) + di2), + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } else if (is_axis_i2max && (particles.i2(p) + di2 >= ni2)) { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + 2 * ni2 - (particles.i2(p) + di2) + N_GHOSTS, + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } else { + buff_access(particles.i1(p) + di1 + N_GHOSTS, + particles.i2(p) + di2 + N_GHOSTS, + particles.i3(p) + di3 + N_GHOSTS, + buff_idx) += coeff * shape_coeff; + } + } + } + } + } + } + } + } + }; + template class NormalizeVectorByRho_kernel { const ndfield_t Rho; diff --git a/src/kernels/twobody_interactions.hpp b/src/kernels/twobody_interactions.hpp new file mode 100644 index 000000000..9a7df409a --- /dev/null +++ b/src/kernels/twobody_interactions.hpp @@ -0,0 +1,443 @@ +/** + * @file kernels/twobody_interactions.hpp + * @brief Generic two-body interaction kernel that can be used to implement various + * types of collisions between species, e.g. Compton scattering, Breit-Wheeler pair production, etc. + * @implements + * - kernel::mink::TwoBodyInteraction<> + * @namespaces: + * - arch::mink:: + */ +#ifndef KERNELS_TWOBODY_INTERACTIONS_HPP +#define KERNELS_TWOBODY_INTERACTIONS_HPP + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "traits/policies.h" +#include "utils/error.h" +#include "utils/sorting.h" + +#include "framework/containers/particles.h" + +#include +#include + +#include +#include + +namespace kernel::mink { + using namespace ntt; + + namespace { + struct CollisionSpecies { + const spidx_t sp; + const npart_t npart; + ncells_t num_tiles { 0u }; + + array_t tileidx; + array_t num_ppt; + + CollisionSpecies(spidx_t sp, + npart_t npart, + const array_t& tileidx, + const array_t& num_ppt, + ncells_t num_tiles) + : sp { sp } + , npart { npart } + , tileidx { tileidx } + , num_ppt { num_ppt } + , num_tiles { num_tiles } {} + }; + + template + Inline void UnravelTileIdx(ncells_t tile_idx, + ncells_t ntx2, + ncells_t ntx3, + ncells_t& ti, + ncells_t& tj, + ncells_t& tk) { + if constexpr (D == Dim::_1D) { + ti = tile_idx; + } else if constexpr (D == Dim::_2D) { + ti = tile_idx / ntx2; + tj = tile_idx % ntx2; + } else if constexpr (D == Dim::_3D) { + ti = tile_idx / (ntx2 * ntx3); + const auto rem = tile_idx % (ntx2 * ntx3); + tj = rem / ntx3; + tk = rem % ntx3; + } else { + raise::KernelError(HERE, "Wrong D in TileIdxUnravel"); + } + } + + template + Inline auto NCellsOnTile(ncells_t tile_idx, + ncells_t tile_size, + ncells_t ntx2, + ncells_t ntx3, + ncells_t nx1, + ncells_t nx2, + ncells_t nx3) -> ncells_t { + ncells_t ncells_on_tile { 1 }; + ncells_t ti { 0u }, tj { 0u }, tk { 0u }; + UnravelTileIdx(tile_idx, ntx2, ntx3, ti, tj, tk); + if constexpr ((D == Dim::_1D) or (D == Dim::_2D) or (D == Dim::_3D)) { + const auto i1_min = ti * tile_size; + const auto i1_max = math::min(i1_min + tile_size, nx1); + ncells_on_tile *= (i1_max - i1_min); + } + if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) { + const auto i2_min = tj * tile_size; + const auto i2_max = math::min(i2_min + tile_size, nx2); + ncells_on_tile *= (i2_max - i2_min); + } + if constexpr (D == Dim::_3D) { + const auto i3_min = tk * tile_size; + const auto i3_max = math::min(i3_min + tile_size, nx3); + ncells_on_tile *= (i3_max - i3_min); + } + return ncells_on_tile; + } + + namespace { + struct PackIndex { + const npart_t offset; + const CollisionSpecies species; + array_t combined_idx; + array_t combined_tileidx; + + PackIndex(npart_t offset, + const CollisionSpecies& species, + array_t& combined_idx, + array_t& combined_tileidx) + : offset { offset } + , species { species } + , combined_idx { combined_idx } + , combined_tileidx { combined_tileidx } {} + + Inline void operator()(prtlidx_t p) const { + // pack species idx into top 8 bits + prtl index into the remaining 56 bits + combined_idx(offset + p) = (static_cast(species.sp) << 56) | + static_cast(p); + combined_tileidx(offset + p) = species.tileidx(p); + } + }; + + struct CombineNumPpt { + const CollisionSpecies species; + array_t combined_num_ppt; + + CombineNumPpt(const CollisionSpecies& species, + array_t& combined_num_ppt) + : species { species } + , combined_num_ppt { combined_num_ppt } {} + + Inline void operator()(cellidx_t t) const { + combined_num_ppt(t) += species.num_ppt(t); + } + }; + + struct PackRandom { + const array_t combined_tileidx; + array_t shuffle_key; + random_number_pool_t random_pool; + + PackRandom(const array_t& combined_tileidx, + array_t& shuffle_key, + random_number_pool_t& random_pool) + : combined_tileidx { combined_tileidx } + , shuffle_key { shuffle_key } + , random_pool { random_pool } {} + + Inline void operator()(prtlidx_t p) const { + auto gen = random_pool.get_state(); + const auto rnd = static_cast(gen.urand()); + random_pool.free_state(gen); + const auto tile_idx = static_cast(combined_tileidx(p)); + // packing top 32 bits with tile index, and the rest -- random + shuffle_key(p) = (tile_idx << 32) | rnd; + } + }; + + struct TileOffsets { + array_t tile_offsets; + array_t combined_num_ppt; + + TileOffsets(array_t& tile_offsets, + const array_t& combined_num_ppt) + : tile_offsets { tile_offsets } + , combined_num_ppt { combined_num_ppt } {} + + Inline void operator()(cellidx_t t, npart_t& acc, const bool is_final) const { + if (is_final) { + tile_offsets(t) = acc; + } + acc += combined_num_ppt(t); + } + }; + } // namespace + + template + struct CollisionGroup { + std::vector group; + + array_t combined_idx; + array_t combined_tileidx; + array_t combined_num_ppt; + array_t tile_offsets; + + ncells_t num_tiles { 0u }; + + CollisionGroup(const std::vector*>& particles, + const std::vector& ncells, + ncells_t tile_size, + random_number_pool_t& random_pool) { + for (const auto* species : particles) { + const auto npart_s = species->npart(); + array_t tileidx { "tile_idx", npart_s }; + auto tile_indexing_kernel = sort::PositionToTileIndex( + species->i1, + species->i2, + species->i3, + species->tag, + tileidx, + ncells, + tile_size); + Kokkos::parallel_for("TileIndexing", species->npart(), tile_indexing_kernel); + group.emplace_back(species->sp, + npart_s, + tileidx, + tile_indexing_kernel.num_ppt, + tile_indexing_kernel.total_tiles); + if (num_tiles == 0u) { + num_tiles = group.back().num_tiles; + } else if (num_tiles != group.back().num_tiles) { + raise::Error("unequal num_tiles across species in group", HERE); + } + raise::ErrorIf(group.back().tileidx.extent(0) != species->npart(), + "tileidx must have the same extent as npart for all " + "species in group", + HERE); + } + + npart_t tot_npart = 0u; + for (const auto& species : group) { + tot_npart += species.npart; + } + + combined_idx = array_t { "combined_idx", tot_npart }; + combined_tileidx = array_t { "combined_tileidx", tot_npart }; + combined_num_ppt = array_t { "combined_num_ppt", num_tiles }; + tile_offsets = array_t { "tile_offsets", num_tiles }; + + { + // combine particle indices in the group & compute total number in each tile + npart_t offset = 0u; + for (const auto& species : group) { + Kokkos::parallel_for( + "CombineInGroup", + species.npart, + PackIndex { offset, species, combined_idx, combined_tileidx }); + offset += species.npart; + Kokkos::parallel_for("CombineNumPpt", + species.num_tiles, + CombineNumPpt { species, combined_num_ppt }); + Kokkos::fence(); + } + } + { + // randomly shuffle particles within each tile and sort by tiles + array_t shuffle_key { "shuffle_key", tot_npart }; + Kokkos::parallel_for( + "PackRandom", + tot_npart, + PackRandom { combined_tileidx, shuffle_key, random_pool }); + Kokkos::Experimental::sort_by_key(Kokkos::DefaultExecutionSpace {}, + shuffle_key, + combined_idx); + } + { + // compute index offsets for each tile + Kokkos::parallel_scan("TileOffsets", + num_tiles, + TileOffsets { tile_offsets, combined_num_ppt }); + } + } + }; + } // namespace + + template + void TwoBodyInteraction( + const std::vector*>& species1, + const std::vector*>& species2, + const std::vector& ncells, + const boundaries_t& domain_extent, + ncells_t tile_size, + real_t ppc0, + random_number_pool_t& random_pool, + I& interaction_policy) { + raise::ErrorIf(species1.empty() or species2.empty(), + "species groups must be non-empty", + HERE); + raise::ErrorIf(ncells.size() != static_cast(D), + "ncells size must match D", + HERE); + raise::ErrorIf(domain_extent.size() != static_cast(D), + "domain_extent size must match D", + HERE); + // compute base cell volume in physical units + real_t cell_volume { ONE }; + for (int d = 0; d < static_cast(D); ++d) { + cell_volume *= static_cast( + domain_extent[d].second - domain_extent[d].first) / + static_cast(ncells[d]); + } + + const auto group1 = CollisionGroup(species1, ncells, tile_size, random_pool); + const auto group2 = CollisionGroup(species2, ncells, tile_size, random_pool); + raise::ErrorIf(group1.num_tiles != group2.num_tiles, + "number of tiles differ in group1 vs group2", + HERE); + const auto num_tiles = group1.num_tiles; + + const auto& combined_idx1 = group1.combined_idx; + const auto& combined_idx2 = group2.combined_idx; + const auto& combined_num_ppt1 = group1.combined_num_ppt; + const auto& combined_num_ppt2 = group2.combined_num_ppt; + const auto& tile_offsets1 = group1.tile_offsets; + const auto& tile_offsets2 = group2.tile_offsets; + + // fill species in the interaction policy + for (auto& sp1 : species1) { + interaction_policy.species[sp1->sp - 1] = static_cast( + *sp1); + } + for (auto& sp2 : species2) { + interaction_policy.species[sp2->sp - 1] = static_cast( + *sp2); + } + + // total particle weight on each tile + auto weights_on_tile1 = array_t { "weights_on_tile1", num_tiles }; + auto weights_on_tile2 = array_t { "weights_on_tile2", num_tiles }; + Kokkos::parallel_for( + "ComputeWeightsOnTiles", + Kokkos::TeamPolicy<>(num_tiles, Kokkos::AUTO), + Lambda(const Kokkos::TeamPolicy<>::member_type& team) { + const ncells_t t = team.league_rank(); + + const auto o1 = tile_offsets1(t); + const auto num_ppt1 = combined_num_ppt1(t); + real_t weight_on_tile1 = ZERO; + Kokkos::parallel_reduce( + Kokkos::TeamThreadRange(team, num_ppt1), + [&](prtlidx_t i, real_t& lsum) { + const auto sp1 = static_cast(combined_idx1(o1 + i) >> 56); + const auto p1 = static_cast(combined_idx1(o1 + i) & + ((1ull << 56) - 1)); + + lsum += (interaction_policy.species[sp1 - 1].tag(p1) == + ParticleTag::alive) + ? interaction_policy.species[sp1 - 1].weight(p1) + : ZERO; + }, + weight_on_tile1); + weights_on_tile1(t) = weight_on_tile1; + + const auto o2 = tile_offsets2(t); + const auto num_ppt2 = combined_num_ppt2(t); + real_t weight_on_tile2 = ZERO; + Kokkos::parallel_reduce( + Kokkos::TeamThreadRange(team, num_ppt2), + [&](prtlidx_t i, real_t& lsum) { + const auto sp2 = static_cast(combined_idx2(o2 + i) >> 56); + const auto p2 = static_cast(combined_idx2(o2 + i) & + ((1ull << 56) - 1)); + + lsum += (interaction_policy.species[sp2 - 1].tag(p2) == + ParticleTag::alive) + ? interaction_policy.species[sp2 - 1].weight(p2) + : ZERO; + }, + weight_on_tile2); + weights_on_tile2(t) = weight_on_tile2; + }); + + // number of cells in each direction + ncells_t nx1 { 1u }, nx2 { 1u }, nx3 { 1u }; + if constexpr ((D == Dim::_1D) or (D == Dim::_2D) or (D == Dim::_3D)) { + nx1 = ncells[0]; + } + if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) { + nx2 = ncells[1]; + } + if constexpr (D == Dim::_3D) { + nx3 = ncells[2]; + } + ncells_t ntx2 { 1u }, ntx3 { 1u }; + if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) { + ntx2 = static_cast( + math::ceil(static_cast(nx2) / static_cast(tile_size))); + } + if constexpr (D == Dim::_3D) { + ntx3 = static_cast( + math::ceil(static_cast(nx3) / static_cast(tile_size))); + } + + array_t interaction_pairs { "interaction_pairs", + combined_idx1.extent(0) }; + array_t counter { "counter" }; + + Kokkos::parallel_for( + "PopulateInteractionPairs", + Kokkos::TeamPolicy<>(num_tiles, Kokkos::AUTO), + Lambda(const Kokkos::TeamPolicy<>::member_type& team) { + const ncells_t t = team.league_rank(); + const auto tile_weight = + math::max(weights_on_tile1(t), weights_on_tile2(t)) / + (NCellsOnTile(t, tile_size, ntx2, ntx3, nx1, nx2, nx3) * ppc0); + + const auto k = math::min(combined_num_ppt1(t), combined_num_ppt2(t)); + const auto o1 = tile_offsets1(t); + const auto o2 = tile_offsets2(t); + Kokkos::parallel_for(Kokkos::TeamThreadRange(team, k), [&](prtlidx_t i) { + // unpack the higher 8 bits + const auto sp1 = static_cast(combined_idx1(o1 + i) >> 56); + const auto sp2 = static_cast(combined_idx2(o2 + i) >> 56); + + // unpack the lower 56 bits + const auto p1 = static_cast(combined_idx1(o1 + i) & + ((1ull << 56) - 1)); + const auto p2 = static_cast(combined_idx2(o2 + i) & + ((1ull << 56) - 1)); + if ((interaction_policy.species[sp1 - 1].tag(p1) == ParticleTag::alive) and + (interaction_policy.species[sp2 - 1].tag(p2) == ParticleTag::alive)) { + if (interaction_policy.should_interact(sp1, p1, sp2, p2, tile_weight)) { + const auto idx = Kokkos::atomic_fetch_add(&counter(), 1); + interaction_pairs(idx, 0) = combined_idx1(o1 + i); + interaction_pairs(idx, 1) = combined_idx2(o2 + i); + } + } + }); + }); + auto counter_h = Kokkos::create_mirror_view(counter); + Kokkos::deep_copy(counter_h, counter); + Kokkos::parallel_for( + "ProcessInteractions", + counter_h(), + Lambda(prtlidx_t idx) { + const auto sp1 = static_cast(interaction_pairs(idx, 0) >> 56); + const auto sp2 = static_cast(interaction_pairs(idx, 1) >> 56); + const auto p1 = static_cast(interaction_pairs(idx, 0) & + ((1ull << 56) - 1)); + const auto p2 = static_cast(interaction_pairs(idx, 1) & + ((1ull << 56) - 1)); + interaction_policy(sp1, p1, sp2, p2); + }); + } + +} // namespace kernel::mink + +#endif // KERNELS_TWOBODY_INTERACTIONS_HPP diff --git a/tests/archetypes/CMakeLists.txt b/tests/archetypes/CMakeLists.txt index 4a5b501e1..4c1aad515 100644 --- a/tests/archetypes/CMakeLists.txt +++ b/tests/archetypes/CMakeLists.txt @@ -16,7 +16,7 @@ function(gen_test title) set(src ${title}.cpp) add_executable(${exec} ${src}) - set(libs ntt_archetypes ntt_global ntt_metrics) + set(libs ntt_archetypes ntt_framework ntt_global ntt_metrics) add_dependencies(${exec} ${libs}) target_link_libraries(${exec} PRIVATE ${libs}) @@ -28,3 +28,4 @@ gen_test(spatial_dist) gen_test(field_setter) gen_test(powerlaw) gen_test(pgen) +gen_test(qed_compton) diff --git a/tests/archetypes/qed_compton.cpp b/tests/archetypes/qed_compton.cpp new file mode 100644 index 000000000..14350cbf4 --- /dev/null +++ b/tests/archetypes/qed_compton.cpp @@ -0,0 +1,280 @@ +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" + +#include "archetypes/qed/compton.h" +#include "framework/containers/particles.h" +#include "kernels/twobody_interactions.hpp" + +#include + +#include +#include +#include + +using namespace ntt; + +void fill_random(array_t& i1, + array_t& i2, + array_t& ux1, + array_t& ux2, + array_t& ux3, + array_t& weight, + array_t& tag, + npart_t npart, + ncells_t nx1, + ncells_t nx2, + random_number_pool_t& rpool) { + Kokkos::parallel_for( + "FillRandom", + npart, + KOKKOS_LAMBDA(const npart_t p) { + auto gen = rpool.get_state(); + i1(p) = static_cast(gen.urand() % static_cast(nx1)); + i2(p) = static_cast(gen.urand() % static_cast(nx2)); + ux1(p) = Random(gen) * TWO - ONE; + ux2(p) = Random(gen) * TWO - ONE; + ux3(p) = Random(gen) * TWO - ONE; + weight(p) = ONE; + tag(p) = ParticleTag::alive; + rpool.free_state(gen); + }); + Kokkos::fence(); +} + +auto get_total_energy(bool is_massive, + array_t& ux1, + array_t& ux2, + array_t& ux3, + npart_t npart) -> real_t { + real_t total_energy = ZERO; + Kokkos::parallel_reduce( + "TotalEnergy", + npart, + Lambda(const npart_t p, real_t& local_sum) { + if (is_massive) { + local_sum += U2GAMMA(ux1(p), ux2(p), ux3(p)); + + } else { + local_sum += NORM(ux1(p), ux2(p), ux3(p)); + } + }, + total_energy); + return total_energy; +} + +auto get_total_momentum_in(in dir, + array_t& ux1, + array_t& ux2, + array_t& ux3, + npart_t npart) -> real_t { + real_t total_momentum_in = ZERO; + Kokkos::parallel_reduce( + "TotalMomentumIn", + npart, + Lambda(const npart_t p, real_t& local_sum) { + if (dir == in::x1) { + local_sum += ux1(p); + } else if (dir == in::x2) { + local_sum += ux2(p); + } else if (dir == in::x3) { + local_sum += ux3(p); + } + }, + total_momentum_in); + return total_momentum_in; +} + +auto main(int argc, char* argv[]) -> int { + ntt::GlobalInitialize(argc, argv); + + try { + const ncells_t nx1 = 32u; + const ncells_t nx2 = 64u; + const ncells_t tile_size = 3u; + const std::vector ncells = { nx1, nx2 }; + const ncells_t ntx1 = static_cast( + math::ceil(static_cast(nx1) / static_cast(tile_size))); + const ncells_t ntx2 = static_cast( + math::ceil(static_cast(nx2) / static_cast(tile_size))); + const npart_t npart = 1000u; + random_number_pool_t random_pool { 12345u }; + + Particles sp1 { 1u, + "sp1", + 1.0f, + 1.0f, + npart, + 0u, + 0u, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + Particles sp2 { 2u, + "sp2", + 1.0f, + -1.0f, + npart, + 0u, + 0u, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + Particles sp3 { 3u, + "sp3", + 0.0f, + 0.0f, + npart, + 0u, + 0u, + ParticlePusher::PHOTON, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + + for (auto* sp : { &sp1, &sp2, &sp3 }) { + sp->set_npart(npart); + fill_random(sp->i1, + sp->i2, + sp->ux1, + sp->ux2, + sp->ux3, + sp->weight, + sp->tag, + npart, + nx1, + nx2, + random_pool); + } + + boundaries_t extent { + { ZERO, ONE }, + { -ONE, ONE } + }; + const auto ppc0 = static_cast(npart) / (nx1 * nx2); + + prm::Parameters params; + params.set("compton_scattering.nominal_probability_density", + static_cast(1e-3)); + + auto policy = arch::qed::ComptonScattering(params, + random_pool); + policy.species[0] = static_cast(sp1); + policy.species[1] = static_cast(sp2); + policy.species[2] = static_cast(sp3); + + std::array init_energies { + get_total_energy(true, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_energy(true, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_energy(false, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + std::array init_moms_x1 { + get_total_momentum_in(in::x1, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_momentum_in(in::x1, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_momentum_in(in::x1, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + std::array init_moms_x2 { + get_total_momentum_in(in::x2, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_momentum_in(in::x2, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_momentum_in(in::x2, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + std::array init_moms_x3 { + get_total_momentum_in(in::x3, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_momentum_in(in::x3, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_momentum_in(in::x3, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + + for (int i = 0; i < 1000; ++i) { + kernel::mink::TwoBodyInteraction({ &sp1, &sp2 }, + { &sp3 }, + ncells, + extent, + tile_size, + ppc0, + random_pool, + policy); + } + + { + std::array fin_energies { + get_total_energy(true, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_energy(true, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_energy(false, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + std::array fin_moms_x1 { + get_total_momentum_in(in::x1, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_momentum_in(in::x1, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_momentum_in(in::x1, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + std::array fin_moms_x2 { + get_total_momentum_in(in::x2, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_momentum_in(in::x2, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_momentum_in(in::x2, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + std::array fin_moms_x3 { + get_total_momentum_in(in::x3, sp1.ux1, sp1.ux2, sp1.ux3, sp1.npart()), + get_total_momentum_in(in::x3, sp2.ux1, sp2.ux2, sp2.ux3, sp2.npart()), + get_total_momentum_in(in::x3, sp3.ux1, sp3.ux2, sp3.ux3, sp3.npart()) + }; + + const auto fin_energy = fin_energies[0] + fin_energies[1] + fin_energies[2]; + const auto init_energy = init_energies[0] + init_energies[1] + + init_energies[2]; + const auto fin_mom_x1 = fin_moms_x1[0] + fin_moms_x1[1] + fin_moms_x1[2]; + const auto init_mom_x1 = init_moms_x1[0] + init_moms_x1[1] + init_moms_x1[2]; + const auto fin_mom_x2 = fin_moms_x2[0] + fin_moms_x2[1] + fin_moms_x2[2]; + const auto init_mom_x2 = init_moms_x2[0] + init_moms_x2[1] + init_moms_x2[2]; + const auto fin_mom_x3 = fin_moms_x3[0] + fin_moms_x3[1] + fin_moms_x3[2]; + const auto init_mom_x3 = init_moms_x3[0] + init_moms_x3[1] + init_moms_x3[2]; + + const auto err_energy = (fin_energy - init_energy) / init_energy; + const auto err_mom_x1 = (fin_mom_x1 - init_mom_x1) / + (std::abs(init_mom_x1) + 1e-10); + const auto err_mom_x2 = (fin_mom_x2 - init_mom_x2) / + (std::abs(init_mom_x2) + 1e-10); + const auto err_mom_x3 = (fin_mom_x3 - init_mom_x3) / + (std::abs(init_mom_x3) + 1e-10); + + raise::ErrorIf(err_energy > 1e-5, + fmt::format("energy is not conserved %e -> %e [%e]", + init_energy, + fin_energy, + err_energy), + HERE); + raise::ErrorIf(err_mom_x1 > 1e-5, + fmt::format("x1 momentum is not conserved %e -> %e [%e]", + init_mom_x1, + fin_mom_x1, + err_mom_x1), + HERE); + raise::ErrorIf(err_mom_x2 > 1e-5, + fmt::format("x2 momentum is not conserved %e -> %e [%e]", + init_mom_x2, + fin_mom_x2, + err_mom_x2), + HERE); + raise::ErrorIf(err_mom_x3 > 1e-5, + fmt::format("x3 momentum is not conserved %e -> %e [%e]", + init_mom_x3, + fin_mom_x3, + err_mom_x3), + HERE); + } + + } catch (std::exception& e) { + std::cerr << e.what() << '\n'; + ntt::GlobalFinalize(); + return 1; + } + ntt::GlobalFinalize(); + return 0; +} diff --git a/tests/global/tiling.cpp b/tests/global/tiling.cpp index 2f926ce28..993d3806c 100644 --- a/tests/global/tiling.cpp +++ b/tests/global/tiling.cpp @@ -108,8 +108,10 @@ void test_tiling(const array_t& i1, for (auto ts { 1u }; ts <= 11u; ++ts) { ncells_t nt1 = 1u, nt2 = 1u, nt3 = 1u; - nt1 = static_cast( - math::ceil(static_cast(ncells[0]) / static_cast(ts))); + if constexpr ((D == Dim::_1D) or (D == Dim::_2D) or (D == Dim::_3D)) { + nt1 = static_cast( + math::ceil(static_cast(ncells[0]) / static_cast(ts))); + } if constexpr ((D == Dim::_2D) or (D == Dim::_3D)) { nt2 = static_cast( math::ceil(static_cast(ncells[1]) / static_cast(ts))); diff --git a/tests/kernels/twobody_interactions.cpp b/tests/kernels/twobody_interactions.cpp new file mode 100644 index 000000000..8051c5805 --- /dev/null +++ b/tests/kernels/twobody_interactions.cpp @@ -0,0 +1,187 @@ +#include "kernels/twobody_interactions.hpp" + +#include "enums.h" +#include "global.h" + +#include "arch/kokkos_aliases.h" +#include "utils/error.h" + +#include "framework/containers/particles.h" +#include "kernels/twobody_interactions.hpp" + +#include + +#include +#include + +using namespace ntt; + +// Verifies that each paired particle from group1 and group2 lies in the same tile +struct SameTilePolicy { + ParticleArrays species[4]; + const ncells_t tile_size; + const ncells_t ncx1, ncx2; // number of cells in each direction + const ncells_t ntx1, ntx2; // numbers of tiles + array_t diff_tile_errors { "diff_tile_errors" }; + + SameTilePolicy(ncells_t tile_size, + ncells_t ncx1, + ncells_t ncx2, + ncells_t ntx1, + ncells_t ntx2) + : tile_size { tile_size } + , ncx1 { ncx1 } + , ncx2 { ncx2 } + , ntx1 { ntx1 } + , ntx2 { ntx2 } {} + + Inline auto should_interact(spidx_t, npart_t, spidx_t, npart_t, real_t) const + -> bool { + return true; + } + + Inline void operator()(spidx_t sp1, npart_t p1, spidx_t sp2, npart_t p2) const { + const auto x1_1 = species[sp1 - 1].i1(p1); + const auto x2_1 = species[sp1 - 1].i2(p1); + const auto x1_2 = species[sp2 - 1].i1(p2); + const auto x2_2 = species[sp2 - 1].i2(p2); + const auto t1 = static_cast(x1_1 / tile_size) * ntx2 + + static_cast(x2_1 / tile_size); + const auto t2 = static_cast(x1_2 / tile_size) * ntx2 + + static_cast(x2_2 / tile_size); + if (t1 != t2) { + Kokkos::atomic_add(&diff_tile_errors(), 1); + } + } +}; + +void fill_random(array_t& i1, + array_t& i2, + array_t& tag, + npart_t npart, + ncells_t nx1, + ncells_t nx2, + random_number_pool_t& rpool) { + Kokkos::parallel_for( + "FillRandom", + npart, + KOKKOS_LAMBDA(const npart_t p) { + auto gen = rpool.get_state(); + i1(p) = static_cast(gen.urand() % static_cast(nx1)); + i2(p) = static_cast(gen.urand() % static_cast(nx2)); + tag(p) = ParticleTag::alive; + rpool.free_state(gen); + }); + Kokkos::fence(); +} + +auto main(int argc, char* argv[]) -> int { + ntt::GlobalInitialize(argc, argv); + + try { + const ncells_t nx1 = 32u; + const ncells_t nx2 = 64u; + const ncells_t tile_size = 3u; + const std::vector ncells = { nx1, nx2 }; + const boundaries_t extent = { + { ZERO, ONE }, + { ZERO, TWO } + }; + const ncells_t ntx1 = static_cast( + math::ceil(static_cast(nx1) / static_cast(tile_size))); + const ncells_t ntx2 = static_cast( + math::ceil(static_cast(nx2) / static_cast(tile_size))); + const npart_t npart = 1000u; + random_number_pool_t random_pool { 12345u }; + + Particles sp1 { 1u, + "sp1", + 1.0f, + 1.0f, + npart, + 0u, + 0u, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + Particles sp2 { 2u, + "sp2", + 1.0f, + 1.0f, + npart, + 0u, + 0u, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + Particles sp3 { 3u, + "sp3", + 1.0f, + 1.0f, + npart, + 0u, + 0u, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + Particles sp4 { 4u, + "sp4", + 1.0f, + 1.0f, + npart, + 0u, + 0u, + ParticlePusher::BORIS, + false, + RadiativeDrag::NONE, + EmissionType::NONE, + 0u, + 0u }; + + for (auto* sp : { &sp1, &sp2, &sp3, &sp4 }) { + sp->set_npart(npart); + fill_random(sp->i1, sp->i2, sp->tag, npart, nx1, nx2, random_pool); + } + + const std::vector*> group1 = { &sp1, + &sp2 }; + const std::vector*> group2 = { &sp3, + &sp4 }; + + auto policy = SameTilePolicy { tile_size, nx1, nx2, ntx1, ntx2 }; + + kernel::mink::TwoBodyInteraction(group1, + group2, + ncells, + extent, + tile_size, + ONE, + random_pool, + policy); + Kokkos::fence(); + + { + auto errors_h = Kokkos::create_mirror_view(policy.diff_tile_errors); + Kokkos::deep_copy(errors_h, policy.diff_tile_errors); + raise::ErrorIf(errors_h() != 0, + "paired particles from different tiles detected", + HERE); + } + + } catch (std::exception& e) { + std::cerr << e.what() << '\n'; + ntt::GlobalFinalize(); + return 1; + } + ntt::GlobalFinalize(); + return 0; +} diff --git a/tutorials/cartesian_sr/cartesian_sr.py b/tutorials/cartesian_sr/cartesian_sr.py new file mode 100644 index 000000000..289d8c4e1 --- /dev/null +++ b/tutorials/cartesian_sr/cartesian_sr.py @@ -0,0 +1,70 @@ +import matplotlib.pyplot as plt +import nt2 +import numpy as np + + +def get_dipole(xs, ys): + xx, yy = np.meshgrid(xs, ys) + rr = np.sqrt(xx**2 + yy**2) + bx = 3 * xx * yy / rr**5 + by = (3 * yy**2 - rr**2) / rr**5 + return bx, by + + +def plot(t, data): + fig = plt.figure(figsize=(6, 5), dpi=150) + gs = fig.add_gridspec(1, 2, width_ratios=[1, 0.05], wspace=0.05) + ax = fig.add_subplot(gs[0, 0]) + ax_cbar = fig.add_subplot(gs[0, 1]) + + d = data.fields.sel(t=t, method="nearest") + (d.N_1 + d.N_2).plot(ax=ax, vmin=0, vmax=20, cmap="inferno", add_colorbar=False) + cbar_ticks = np.linspace(0, 20, 50) + ax_cbar.pcolormesh( + [0, 1], + cbar_ticks, + np.array([cbar_ticks] * 2).T, + cmap="inferno", + vmin=0, + vmax=20, + rasterized=True, + ) + ax_cbar.yaxis.tick_right() + ax_cbar.yaxis.set_label_position("right") + ax_cbar.set(xticks=[], ylabel=r"$n_\pm$") + for spine in ax_cbar.spines.values(): + spine.set_visible(False) + + ys = np.linspace(-2.9, 2.9, 50) + xs = -0.5 * np.ones_like(ys) + + bx, by = get_dipole(data.fields.x, data.fields.y) + ax.streamplot( + data.fields.x.values, + data.fields.y.values, + d.Bx.values + bx, + d.By.values + by, + color="#ffffff90", + linewidth=0.25, + zorder=10, + arrowstyle="->", + arrowsize=0.5, + density=30, + start_points=np.array([xs, ys]).T, + ) + ax.add_artist( + plt.Circle((0, 0), data.attrs["setup.r_plummet"], color="C0", zorder=20) + ) + + ax.set( + xlim=(-4, 2), + ylim=(-3, 3), + aspect=1, + xlabel="$x$", + ylabel="$y$", + title=f"$t = {t:.2f}$", + ) + + +data = nt2.Data("cartesian_sr") +data.makeMovie() diff --git a/tutorials/cartesian_sr/tutorial_cartesian_sr.toml b/tutorials/cartesian_sr/cartesian_sr.toml similarity index 77% rename from tutorials/cartesian_sr/tutorial_cartesian_sr.toml rename to tutorials/cartesian_sr/cartesian_sr.toml index 8a85767b2..d7fdb90bc 100644 --- a/tutorials/cartesian_sr/tutorial_cartesian_sr.toml +++ b/tutorials/cartesian_sr/cartesian_sr.toml @@ -1,5 +1,5 @@ [simulation] - name = "tutorial_cartesian_sr" + name = "cartesian_sr" engine = "srpic" runtime = 20.0 @@ -60,5 +60,18 @@ enable = false [checkpoint] - interval = 1.0 - keep = 1 + interval_time = 1.0 + keep = 1 + +[render] + enable = true + interval_time = 1.0 + resolution = 2048 + axes = true + + [[render.scene]] + field = "N_1_2" + label = "N / n0" + min = 0.0 + max = 20.0 + colormap = "inferno"