diff --git a/.gitignore b/.gitignore index 9a167b9d..52d372f4 100644 --- a/.gitignore +++ b/.gitignore @@ -39,6 +39,7 @@ venv/ *.png *.mov *.mp4 +!pgens/**/*.png # Accidental files *.xc diff --git a/CMakeLists.txt b/CMakeLists.txt index b088e15c..29566913 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ set(PROJECT_NAME entity) project( ${PROJECT_NAME} - VERSION 1.2.1 + VERSION 1.2.3 LANGUAGES CXX C) add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"") set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ") diff --git a/input.example.toml b/input.example.toml index 6444a16f..028102b0 100644 --- a/input.example.toml +++ b/input.example.toml @@ -467,8 +467,9 @@ # Field quantities to output # @type: array # @default: ["B^2", "E^2", "ExB", "Rho", "T00"] - # @enum: "B^2", "E^2", "ExB", "N", "Charge", "Rho", "T00", "T0i", "Tij" - # @note: Same notation as for `output.fields.quantities` + # @enum: "B^2", "E^2", "ExB", "N", "Npart", "Charge", "Rho", "T00", "T0i", "Tij" + # @note: For particle moments, ... + # @note: ... same notation is used as for `output.fields.quantities` quantities = "" # Custom (user-defined) stats # @type: array diff --git a/pgens/magnetosphere/magnetosphere.toml b/pgens/magnetosphere/magnetosphere.toml index 1a4af8a0..0eb9a03f 100644 --- a/pgens/magnetosphere/magnetosphere.toml +++ b/pgens/magnetosphere/magnetosphere.toml @@ -58,8 +58,9 @@ pusher = "Boris,GCA" [setup] - Bsurf = 1.0 - period = 60.0 + Bsurf = 1.0 + field_geometry = "dipole" # can be "dipole" or "monopole" (default: dipole) + period = 60.0 [output] format = "hdf5" diff --git a/pgens/magnetosphere/pgen.hpp b/pgens/magnetosphere/pgen.hpp index 64fe13cf..1afc8ccc 100644 --- a/pgens/magnetosphere/pgen.hpp +++ b/pgens/magnetosphere/pgen.hpp @@ -6,33 +6,58 @@ #include "arch/kokkos_aliases.h" #include "arch/traits.h" +#include "utils/numeric.h" #include "archetypes/problem_generator.h" #include "framework/domain/metadomain.h" +#include + namespace user { using namespace ntt; + enum class FieldGeometry { + dipole, + monopole + }; + template struct InitFields { - InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {} + InitFields(real_t bsurf, real_t rstar, const std::string& field_geometry) + : Bsurf { bsurf } + , Rstar { rstar } + , field_geom { field_geometry == "monopole" ? FieldGeometry::monopole + : FieldGeometry::dipole } {} Inline auto bx1(const coord_t& x_Ph) const -> real_t { - return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); + if (field_geom == FieldGeometry::monopole) { + return Bsurf / SQR(x_Ph[0] / Rstar); + } else { + return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); + } } Inline auto bx2(const coord_t& x_Ph) const -> real_t { - return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); + if (field_geom == FieldGeometry::monopole) { + return ZERO; + } else { + return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar); + } } private: - const real_t Bsurf, Rstar; + const real_t Bsurf, Rstar; + const FieldGeometry field_geom; }; template struct DriveFields : public InitFields { - DriveFields(real_t time, real_t bsurf, real_t rstar, real_t omega) - : InitFields { bsurf, rstar } + DriveFields(real_t time, + real_t bsurf, + real_t rstar, + real_t omega, + const std::string& field_geometry) + : InitFields { bsurf, rstar, field_geometry } , time { time } , Omega { omega } {} @@ -73,8 +98,9 @@ namespace user { using arch::ProblemGenerator::C; using arch::ProblemGenerator::params; - const real_t Bsurf, Rstar, Omega; - InitFields init_flds; + const real_t Bsurf, Rstar, Omega; + const std::string field_geom; + InitFields init_flds; inline PGen(const SimulationParams& p, const Metadomain& m) : arch::ProblemGenerator(p) @@ -82,16 +108,17 @@ namespace user { , Rstar { m.mesh().extent(in::x1).first } , Omega { static_cast(constant::TWO_PI) / p.template get("setup.period", ONE) } - , init_flds { Bsurf, Rstar } {} + , field_geom { p.template get("setup.field_geometry", "dipole") } + , init_flds { Bsurf, Rstar, field_geom } {} inline PGen() {} auto AtmFields(real_t time) const -> DriveFields { - return DriveFields { time, Bsurf, Rstar, Omega }; + return DriveFields { time, Bsurf, Rstar, Omega, field_geom }; } auto MatchFields(real_t) const -> InitFields { - return InitFields { Bsurf, Rstar }; + return InitFields { Bsurf, Rstar, field_geom }; } }; diff --git a/pgens/magnetosphere/sketch.png b/pgens/magnetosphere/sketch.png new file mode 100644 index 00000000..fe33228e Binary files /dev/null and b/pgens/magnetosphere/sketch.png differ diff --git a/pgens/magnetosphere/sketch.py b/pgens/magnetosphere/sketch.py new file mode 100644 index 00000000..648c5ce7 --- /dev/null +++ b/pgens/magnetosphere/sketch.py @@ -0,0 +1,163 @@ +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches +import myplotlib +import numpy as np + +plt.style.use("latex") + +fig = plt.figure(dpi=600) +ax = fig.add_subplot(111) + +r0 = 0.15 + +ax.plot([0, 0], [r0, 1], c="k", lw=0.5) +ax.plot([0, 0], [-r0, -1], c="k", lw=0.5) +ax.add_patch( + mpatches.Arc((0, 0), 2, 2, angle=0, theta1=-90, theta2=90, color="k", lw=0.5) +) +ax.add_patch( + mpatches.Arc( + (0, 0), 2 * r0, 2 * r0, angle=0, theta1=-90, theta2=90, color="k", lw=0.5 + ) +) +ax.add_patch( + mpatches.Arc( + (0, 0), + 2 - 0.25, + 2 - 0.25, + angle=0, + theta1=-90, + theta2=90, + color="b", + lw=0.5, + ls="--", + ) +) +ax.add_patch( + mpatches.Arc( + (0, 0), + 2 - 0.1, + 2 - 0.1, + angle=0, + theta1=-90, + theta2=90, + color="r", + lw=0.5, + ls="--", + ) +) +ax.text( + 1.05 / np.sqrt(2), 1.05 / np.sqrt(2), "absorbing particle boundaries", c="r", size=5 +) +ax.text( + 1.05 / np.sqrt(2), + 1.05 / np.sqrt(2) - 0.04, + r"$\mathtt{grid.boundaries.absorb.ds}$", + c="r", + size=5, +) +ax.annotate( + "", + xy=(1.02 / np.sqrt(2), 1.02 / np.sqrt(2)), + xytext=(0.93 / np.sqrt(2), 0.93 / np.sqrt(2)), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="r"), + size=5, +) +ax.annotate( + "", + xy=(1.01 * np.cos(np.pi / 6), 1.01 * np.sin(np.pi / 6)), + xytext=(0.87 * np.cos(np.pi / 6), 0.87 * np.sin(np.pi / 6)), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="b"), + size=5, +) +ax.text( + 1.05 * np.cos(np.pi / 6), + 1.05 * np.sin(np.pi / 6), + "matching field boundaries", + c="b", + size=5, +) +ax.text( + 1.05 * np.cos(np.pi / 6), + 1.05 * np.sin(np.pi / 6) - 0.04, + r"$\mathtt{grid.boundaries.match.ds}$", + c="b", + size=5, +) +ax.add_patch( + mpatches.Arc( + (0, 0), + 0.4, + 0.4, + angle=0, + theta1=-90, + theta2=90, + color="g", + lw=0.5, + ls=":", + ) +) +ax.add_patch( + mpatches.Arc( + (0, 0), + 0.5, + 0.5, + angle=0, + theta1=-90, + theta2=90, + color="g", + lw=0.5, + ls=":", + ) +) +ax.annotate( + "", + xy=(0.27 / np.sqrt(2), 0.27 / np.sqrt(2)), + xytext=(0.18 / np.sqrt(2), 0.18 / np.sqrt(2)), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="g"), + size=5, +) +ax.text( + 0.27 / np.sqrt(2), + 0.27 / np.sqrt(2) + 0.04, + "particle atmosphere injection", + c="g", + size=5, +) +ax.text( + 0.27 / np.sqrt(2), + 0.27 / np.sqrt(2), + r"$\mathtt{grid.boundaries.atmosphere.height}$", + c="g", + size=5, +) +ax.annotate( + "", + xy=(0.13 / np.sqrt(2), -0.13 / np.sqrt(2)), + xytext=(0.22 / np.sqrt(2), -0.22 / np.sqrt(2)), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="b"), + size=5, +) +ax.text( + 0.23 / np.sqrt(2), + -0.23 / np.sqrt(2), + "buffer zone for resetting fields", + c="b", + size=5, +) +ax.text( + 0.23 / np.sqrt(2), + -0.23 / np.sqrt(2)-0.04, + r"size in cells = \# of filters", + c="b", + size=5, +) +ax.set( + xlim=(-0.05, 1.05), + ylim=(-1.05, 1.05), + aspect=1, + xticks=[], + yticks=[], + frame_on=False, +) +plt.savefig("sketch.png", bbox_inches="tight") diff --git a/pgens/reconnection/pgen.hpp b/pgens/reconnection/pgen.hpp index 8e8804b2..91aa4639 100644 --- a/pgens/reconnection/pgen.hpp +++ b/pgens/reconnection/pgen.hpp @@ -13,6 +13,7 @@ #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" #include "archetypes/spatial_dist.h" +#include "archetypes/utils.h" #include "framework/domain/metadomain.h" #include "kernels/particle_moments.hpp" @@ -205,17 +206,11 @@ namespace user { inline void InitPrtls(Domain& local_domain) { // background - const auto energy_dist = arch::Maxwellian(local_domain.mesh.metric, - local_domain.random_pool, - bg_temperature); - const auto injector = arch::UniformInjector( - energy_dist, - { 1, 2 }); - arch::InjectUniform>( - params, - local_domain, - injector, - ONE); + arch::InjectUniformMaxwellian(params, + local_domain, + ONE, + bg_temperature, + { 1, 2 }); const auto sigma = params.template get("scales.sigma0"); const auto c_omp = params.template get("scales.skindepth0"); diff --git a/pgens/reconnection/reconnection.toml b/pgens/reconnection/reconnection.toml index 76acbf08..8205a499 100644 --- a/pgens/reconnection/reconnection.toml +++ b/pgens/reconnection/reconnection.toml @@ -1,14 +1,14 @@ [simulation] name = "reconnection" engine = "srpic" - runtime = 10.0 + runtime = 10000.0 [simulation.domain] decomposition = [-1, 2] [grid] - resolution = [512, 512] - extent = [[-1.0, 1.0], [-1.0, 1.0]] + resolution = [4096, 2048] + extent = [[-500.0, 500.0], [-250.0, 250.0]] [grid.metric] metric = "minkowski" @@ -18,11 +18,11 @@ particles = [["PERIODIC"], ["ABSORB", "ABSORB"]] [grid.boundaries.match] - ds = [[0.04], [0.1]] + ds = [[10.0], [20.0]] [scales] - larmor0 = 2e-4 - skindepth0 = 2e-3 + larmor0 = 0.1 + skindepth0 = 1.0 [algorithms] current_filters = 8 @@ -49,23 +49,19 @@ bg_B = 1.0 bg_Bguide = 0.0 bg_temperature = 1e-4 - inj_ypad = 0.25 - cs_width = 0.05 + inj_ypad = 50.0 + cs_width = 10.0 cs_overdensity = 3.0 [output] - format = "hdf5" - interval_time = 0.1 + format = "bpfile" + interval_time = 100.0 [output.fields] quantities = ["N_1", "N_2", "E", "B", "J"] [output.particles] - enable = false + stride = 25 [output.spectra] enable = false - -[diagnostics] - colored_stdout = true - interval = 10 diff --git a/pgens/reconnection/sketch.png b/pgens/reconnection/sketch.png new file mode 100644 index 00000000..6b25615d Binary files /dev/null and b/pgens/reconnection/sketch.png differ diff --git a/pgens/reconnection/sketch.py b/pgens/reconnection/sketch.py new file mode 100644 index 00000000..6b40a3e0 --- /dev/null +++ b/pgens/reconnection/sketch.py @@ -0,0 +1,167 @@ +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches +import myplotlib + +plt.style.use("latex") + +fig = plt.figure(dpi=600) +ax = fig.add_subplot(111) +ax.plot() +ax.set( + xlim=(-1, 1), + ylim=(-0.5, 0.5), + xlabel="$x$", + ylabel="$y$", + xticks=[], + yticks=[], + aspect=1, +) +ax.axhline(0, color="black", lw=0.5, ls="--") +ax.axvline(0.95, color="red", lw=0.5, ls=":") +ax.axvline(-0.95, color="red", lw=0.5, ls=":") +ax.axhline(-0.45, color="blue", lw=0.5, ls=":") +ax.axhline(0.45, color="blue", lw=0.5, ls=":") +ax.axhline(-0.05, color="g", lw=0.5, ls="--") +ax.axhline(0.05, color="g", lw=0.5, ls="--") +ax.axhline(0.35, color="magenta", lw=0.5, ls="--") +ax.axhline(-0.35, color="magenta", lw=0.5, ls="--") +ax.annotate( + "", + xy=(-1, 0.25), + xytext=(-0.95, 0.25), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="red"), + size=3, +) +ax.annotate( + "", + xy=(-0.5, 0.45), + xytext=(-0.5, 0.5), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="b"), + size=3, +) +ax.text( + -0.94, + 0.25, + r"$\mathtt{grid.boundaries.match.ds[0]}$", + color="red", + size=5, + ha="left", + va="center", +) +ax.text( + -0.5, + 0.44, + r"$\mathtt{grid.boundaries.match.ds[1]}$", + color="blue", + size=5, + ha="center", + va="top", +) +ax.annotate( + "", + xy=(-0.5, -0.05), + xytext=(-0.5, 0.05), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="g"), + size=3, +) +ax.annotate( + "", + xy=(0.5, 0.35), + xytext=(0.5, 0.5), + arrowprops=dict(arrowstyle="<->", lw=0.5, color="magenta"), + size=3, +) +ax.text( + -0.5, + -0.06, + r"$\mathtt{setup.cs\_width}$", + color="g", + size=5, + ha="center", + va="top", +) +ax.text( + 0.5, + 0.33, + r"$\mathtt{setup.inj\_ypad}$", + color="magenta", + size=5, + ha="center", + va="top", +) +ax.text( + 0.0, + 0.0, + r"current layer plasma with peak density $\mathtt{setup.cs\_overdensity} \times n_0$", + color="green", + size=5, + ha="center", + va="center", + bbox=dict(boxstyle="round,pad=0.2", fc="white", ec="none"), +) +ax.text( + 0.0, + -0.25, + "background plasma with density $n_0$", + color="green", + size=5, + ha="center", + va="center", +) +ax.annotate( + "", + xy=(0.7, 0.2), + xytext=(0.9, 0.2), + arrowprops=dict(arrowstyle="<-", lw=0.5, color="red"), + size=5, +) +ax.annotate( + "", + xy=(0.7, -0.2), + xytext=(0.9, -0.2), + arrowprops=dict(arrowstyle="->", lw=0.5, color="red"), + size=5, +) +ax.text( + 0.9, + 0.22, + r"magnetic field of strength $\mathtt{setup.bg\_B}\times B_0$", + color="red", + size=5, + ha="right", + va="bottom", +) +ax.text( + -0.93, + -0.22, + "matching fields", + color="red", + size=5, + ha="left", + va="center", + rotation=90, +) +ax.text( + 0, + -0.43, + "matching fields", + color="blue", + size=5, + ha="center", + va="bottom", +) +ax.add_patch(mpatches.Circle((0.75, 0.15), 0.02, color="red", fill=False, lw=0.5)) +ax.add_patch(mpatches.Circle((0.75, 0.15), 0.005, color="red", fill=True, lw=0.5)) + +ax.add_patch(mpatches.Circle((0.75, -0.15), 0.02, color="red", fill=False, lw=0.5)) +ax.add_patch(mpatches.Circle((0.75, -0.15), 0.005, color="red", fill=True, lw=0.5)) +ax.text( + 0.71, + 0.15, + r"$\mathtt{setup.bg\_Bguide} \times B_0$", + color="r", + size=5, + ha="right", + va="center", +) +plt.savefig("sketch.png", bbox_inches="tight") diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp index 6bd6f21a..f93bbfd0 100644 --- a/pgens/shock/pgen.hpp +++ b/pgens/shock/pgen.hpp @@ -12,6 +12,7 @@ #include "archetypes/field_setter.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" +#include "archetypes/utils.h" #include "framework/domain/metadomain.h" #include @@ -67,7 +68,6 @@ namespace user { const real_t Btheta, Bphi, Vx, Bmag; }; - template struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator @@ -136,7 +136,7 @@ namespace user { } } - inline void InitPrtls(Domain& local_domain) { + inline void InitPrtls(Domain& domain) { /* * Plasma setup as partially filled box @@ -173,28 +173,19 @@ namespace user { // species #2 -> protons // energy distribution of the particles - const auto energy_dist = arch::TwoTemperatureMaxwellian( - local_domain.mesh.metric, - local_domain.random_pool, - { temperature_ratio * temperature * local_domain.species[1].mass() , - temperature }, - { 1, 2 }, - -drift_ux, - in::x1); - - // we want to set up a uniform density distribution - const auto injector = arch::UniformInjector( - energy_dist, - { 1, 2 }); - - // inject uniformly within the defined box - arch::InjectUniform>( - params, - local_domain, - injector, - 1.0, // target density - false, // no weights - box); + const auto temperatures = std::make_pair(temperature, + temperature_ratio * temperature); + const auto drifts = std::make_pair( + std::vector { -drift_ux, ZERO, ZERO }, + std::vector { -drift_ux, ZERO, ZERO }); + arch::InjectUniformMaxwellians(params, + domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + box); } void CustomPostStep(timestep_t step, simtime_t time, Domain& domain) { @@ -315,28 +306,19 @@ namespace user { } // same maxwell distribution as above - const auto energy_dist = arch::TwoTemperatureMaxwellian( - domain.mesh.metric, - domain.random_pool, - { temperature_ratio * temperature * domain.species[1].mass(), - temperature }, - { 1, 2 }, - -drift_ux, - in::x1); - - // we want to set up a uniform density distribution - const auto injector = arch::UniformInjector( - energy_dist, - { 1, 2 }); - - // inject uniformly within the defined box - arch::InjectUniform>( - params, - domain, - injector, - 1.0, // target density - false, // no weights - inj_box); + const auto temperatures = std::make_pair(temperature, + temperature_ratio * temperature); + const auto drifts = std::make_pair( + std::vector { -drift_ux, ZERO, ZERO }, + std::vector { -drift_ux, ZERO, ZERO }); + arch::InjectUniformMaxwellians(params, + domain, + ONE, + temperatures, + { 1, 2 }, + drifts, + false, + inj_box); } }; } // namespace user diff --git a/pgens/shock/sketch.png b/pgens/shock/sketch.png new file mode 100644 index 00000000..9df56272 Binary files /dev/null and b/pgens/shock/sketch.png differ diff --git a/pgens/shock/sketch.py b/pgens/shock/sketch.py new file mode 100644 index 00000000..fa78ae6d --- /dev/null +++ b/pgens/shock/sketch.py @@ -0,0 +1,141 @@ +import matplotlib.pyplot as plt +import myplotlib + +plt.style.use("latex") + +fig = plt.figure(figsize=(14, 6), dpi=600) +ax = fig.add_subplot(111) + +ax.plot() +ax.set( + xlim=(0, 10), + ylim=(-2, 2), + aspect=1, +) +ax.text( + 0.1, + 0.5 + 0.15, + r'$\mathtt{grid.boundaries.fields[0][0]: ``CONDUCTOR"}$', + ha="left", + va="center", + c="C3", +) +ax.text( + 0.1, + 0.5 - 0.1, + r'$\mathtt{grid.boundaries.particles[0][0]: ``REFLECT"}$', + ha="left", + va="center", + c="C3", +) +ax.text( + 0.1, + 1.95, + r'$\mathtt{grid.boundaries.fields[1]: ``PERIODIC"}$', + ha="left", + va="top", + c="C2", +) +ax.text( + 0.1, + 1.75, + r'$\mathtt{grid.boundaries.particles[1]: ``PERIODIC"}$', + ha="left", + va="top", + c="C2", +) +ax.text( + 10 - 0.1, + 0.15 - 0.5, + r'$\mathtt{grid.boundaries.fields[0][1]: ``MATCH"}$', + ha="right", + va="center", + c="C1", +) +ax.text( + 10 - 0.1, + -0.1 - 0.5, + r'$\mathtt{grid.boundaries.particles[0][1]: ``ABSORB"}$', + ha="right", + va="center", + c="C1", +) +ax.spines["left"].set_color("C3") +ax.spines["top"].set_color("C2") +ax.spines["bottom"].set_color("C2") +ax.spines["right"].set_color("C1") + +ax.text( + 2, -1.45, r"$\mathtt{setup.filling\_fraction}$", ha="center", va="bottom", c="C0" +) +ax.annotate("", (0, -1.5), (4, -1.5), arrowprops=dict(arrowstyle="<->", color="C0")) +ax.annotate("", (4, -1.35), (4.5, -1.35), arrowprops=dict(arrowstyle="<|-", color="C0")) +ax.text( + 4.25, -1.3, r"$\mathtt{setup.injector\_velocity}$", ha="left", va="bottom", c="C0" +) +ax.annotate("", (4.5, 0.5), (6, 1.75), arrowprops=dict(arrowstyle="<-", color="C4")) +ax.plot([4.5, 5], [0.5, 0.5], c="C4", ls="--", lw=1) +ax.text( + 5.5, + 1.2, + r"initial (upstream) B-field of strength $\mathtt{setup.Bmag}\times B_0$", + ha="left", + va="center", + c="C4", +) +ax.text( + 5, + 0.7, + r"angle w.r.t. the shock normal: $\mathtt{setup.Btheta}$", + ha="left", + va="center", + c="C4", +) +ax.text( + 5, + 0.3, + r"out-of-plane angle: $\mathtt{setup.Bphi}$", + ha="left", + va="center", + c="C4", +) + +ax.axvline(4, c="C0", ls="--") + +ax.annotate("", (2, -0.75), (0.75, -0.75), arrowprops=dict(arrowstyle="<|-", color="C5")) +ax.text( + 0.75, + -0.1, + r"upstream plasma:", + ha="left", + va="center", + c="C5", +) +ax.text( + 0.75, + -0.1-0.15, + r"$u_x = \mathtt{setup.drift\_ux}$", + ha="left", + va="center", + c="C5", +) +ax.text( + 0.75, + -0.1-2 * 0.15, + r"$T_-/m_0 c^2 = \mathtt{setup.temperature}$", + ha="left", + va="center", + c="C5", +) +ax.text( + 0.75, + -0.1-3 * 0.15, + r"$T_i/m_0 c^2 = \mathtt{setup.temperature} \times \mathtt{setup.temperature\_ratio}$", + ha="left", + va="center", + c="C5", +) +ax.set( + xticks=[], yticks=[] +) +plt.savefig("sketch.png", bbox_inches="tight") diff --git a/pgens/streaming/pgen.hpp b/pgens/streaming/pgen.hpp index ee14712d..1b8311b3 100644 --- a/pgens/streaming/pgen.hpp +++ b/pgens/streaming/pgen.hpp @@ -9,9 +9,8 @@ #include "utils/error.h" #include "utils/numeric.h" -#include "archetypes/energy_dist.h" -#include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" +#include "archetypes/utils.h" #include "framework/domain/domain.h" #include "framework/domain/metadomain.h" @@ -81,28 +80,19 @@ namespace user { inline void InitPrtls(Domain& domain) { const auto nspec = domain.species.size(); for (auto n = 0u; n < nspec; n += 2) { - const auto drift_1 = prmvec_t { drifts_in_x[n], + const auto drift_1 = prmvec_t { drifts_in_x[n], drifts_in_y[n], drifts_in_z[n] }; - const auto drift_2 = prmvec_t { drifts_in_x[n + 1], + const auto drift_2 = prmvec_t { drifts_in_x[n + 1], drifts_in_y[n + 1], drifts_in_z[n + 1] }; - const auto injector = arch::experimental:: - UniformInjector( - arch::experimental::Maxwellian(domain.mesh.metric, - domain.random_pool, - temperatures[n], - drift_1), - arch::experimental::Maxwellian(domain.mesh.metric, - domain.random_pool, - temperatures[n + 1], - drift_2), - { n + 1, n + 2 }); - arch::experimental::InjectUniform( + arch::InjectUniformMaxwellians( params, domain, - injector, - densities[n / 2]); + densities[n / 2], + { temperatures[n], temperatures[n + 1] }, + { n + 1, n + 2 }, + { drift_1, drift_2 }); } } }; diff --git a/pgens/streaming/sketch.png b/pgens/streaming/sketch.png new file mode 100644 index 00000000..a2cc6610 Binary files /dev/null and b/pgens/streaming/sketch.png differ diff --git a/pgens/streaming/sketch.py b/pgens/streaming/sketch.py new file mode 100644 index 00000000..632d2f43 --- /dev/null +++ b/pgens/streaming/sketch.py @@ -0,0 +1,33 @@ +import matplotlib.pyplot as plt +import myplotlib + +plt.style.use("latex") + +fig = plt.figure(dpi=600) +ax = fig.add_subplot(111) + +ax.set( + xlim=(0, 2), + ylim=(0, 1), + aspect=1, +) +ax.text(0.2, 0.5, r"i-th species:", c="C0") +ax.text(0.2, 0.5 - 0.06, r"$u_D^i = \mathtt{setup.drifts\_in\_*[i]}$", c="C0") +ax.text(0.2, 0.5 - 2 * 0.06, r"$n_i=\mathtt{setup.densities[i / 2]}$", c="C0") +ax.text(0.2, 0.5 - 3 * 0.06, r"$T_i/m_0c^2=\mathtt{setup.temperatures[i]}$", c="C0") +ax.annotate("", (0.5, 0.2), (0.1, 0.25), arrowprops=dict(arrowstyle="->", color="C0")) + +ax.text(0.95, 0.8, r"i+1-th species:", c="C2") +ax.text(0.95, 0.8 - 0.06, r"$u_D^{i+1} = \mathtt{setup.drifts\_in\_*[i+1]}$", c="C2") +ax.text(0.95, 0.8 - 2 * 0.06, r"$n_{i+1}=\mathtt{setup.densities[(i + 1) / 2]}$", c="C2") +ax.text(0.95, 0.8 - 3 * 0.06, r"$T_{i+1}/m_0c^2=\mathtt{setup.temperatures[i+1]}$", c="C2") +ax.annotate("", (1.25, 0.4), (1.75, 0.5), arrowprops=dict(arrowstyle="->", color="C2")) + +ax.text(1.0, 0.98, r"periodic boundaries", c="C3", ha="center", va="top") +ax.spines["top"].set_color("C3") +ax.spines["bottom"].set_color("C3") +ax.spines["left"].set_color("C3") +ax.spines["right"].set_color("C3") +ax.set(xticks=[], yticks=[]) + +plt.savefig("sketch.png", bbox_inches="tight") diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h new file mode 100644 index 00000000..d3d6ae47 --- /dev/null +++ b/src/archetypes/utils.h @@ -0,0 +1,115 @@ +/** + * @file archetypes/utils.h + * @brief Utility functions that use prefabricated archetypes for the most common tasks + * @implements + * - arch::InjectUniformMaxwellians<> -> void + * - arch::InjectUniformMaxwellian<> -> void + * @namespaces: + * - arch:: + */ + +#ifndef ARCHETYPES_UTILS_H +#define ARCHETYPES_UTILS_H + +#include "enums.h" +#include "global.h" + +#include "archetypes/energy_dist.h" +#include "archetypes/particle_injector.h" +#include "framework/domain/domain.h" + +#include + +namespace arch { + + /** + * @brief Injects uniform number density of particles everywhere in the domain + * following two Maxwellian distributions + * + * @param domain Domain object + * @param tot_number_density Total number density (in units of n0) + * @param temperatures Temperatures of the two species (in units of m0 c^2) + * @param species Pair of species indices + * @param drift_four_vels Pair of drift four-velocities for the two species + * @param use_weights Use weights + * @param box Region to inject the particles in global coords (or empty for the whole domain) + */ + template + inline void InjectUniformMaxwellians( + const SimulationParams& params, + Domain& domain, + real_t tot_number_density, + const std::pair& temperatures, + const std::pair& species, + const std::pair, std::vector>& drift_four_vels = {{ ZERO, ZERO, ZERO }, { ZERO, ZERO, ZERO }}, + bool use_weights = false, + const boundaries_t& box = {}) { + static_assert(M::is_metric, "M must be a metric class"); + + const auto mass_1 = domain.species[species.first - 1].mass(); + const auto mass_2 = domain.species[species.second - 1].mass(); + const auto temperature_1 = temperatures.first / mass_1; + const auto temperature_2 = temperatures.second / mass_2; + + const auto maxwellian_1 = arch::experimental::Maxwellian( + domain.mesh.metric, + domain.random_pool, + temperature_1, + drift_four_vels.first); + const auto maxwellian_2 = arch::experimental::Maxwellian( + domain.mesh.metric, + domain.random_pool, + temperature_2, + drift_four_vels.second); + + const auto injector = arch::experimental:: + UniformInjector( + maxwellian_1, + maxwellian_2, + species); + + arch::experimental::InjectUniform(params, + domain, + injector, + tot_number_density, + use_weights, + box); + } + + /** + * @brief Injects uniform number density of particles everywhere in the domain + * with the same temperature for both species + * + * @param domain Domain object + * @param tot_number_density Total number density (in units of n0) + * @param temperature Temperature (in units of m0 c^2) + * @param species Pair of species indices + * @param drift_four_vels Pair of drift four-velocities for the two species + * @param use_weights Use weights + * @param box Region to inject the particles in global coords (or empty for the whole domain) + */ + template + inline void InjectUniformMaxwellian( + const SimulationParams& params, + Domain& domain, + real_t tot_number_density, + real_t temperature, + const std::pair& species, + const std::pair, std::vector>& drift_four_vels = {{ ZERO, ZERO, ZERO }, { ZERO, ZERO, ZERO }}, + bool use_weights = false, + const boundaries_t& box = {}) { + static_assert(M::is_metric, "M must be a metric class"); + + InjectUniformMaxwellians(params, + domain, + tot_number_density, + { temperature, temperature }, + species, + drift_four_vels, + use_weights, + box); + } + +} // namespace arch + +#endif // ARCHETYPES_UTILS_H diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 960b2c71..6dad452e 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -202,12 +202,12 @@ namespace ntt { const ncells_t k_max = (i2 - (N_GHOSTS)); real_t A3 = ZERO; for (auto k { k_min }; k <= k_max; ++k) { - real_t k_ = static_cast(k); - real_t sqrt_detH_ij1 { metric.sqrt_det_h({ i1_, k_ - HALF }) }; - real_t sqrt_detH_ij2 { metric.sqrt_det_h({ i1_, k_ + HALF }) }; - auto k1 { k + N_GHOSTS }; - A3 += HALF * (sqrt_detH_ij1 * EM(i1, k + N_GHOSTS - 1, em::bx1) + - sqrt_detH_ij2 * EM(i1, k + N_GHOSTS, em::bx1)); + const real_t k_ = static_cast(k); + const real_t sqrt_detH_ij1 { metric.sqrt_det_h({ i1_, k_ - HALF }) }; + const real_t sqrt_detH_ij2 { metric.sqrt_det_h({ i1_, k_ + HALF }) }; + const auto k1 { k + N_GHOSTS }; + A3 += HALF * (sqrt_detH_ij1 * EM(i1, k1 - 1, em::bx1) + + sqrt_detH_ij2 * EM(i1, k1, em::bx1)); } buffer(i1, i2, buff_idx) = A3; }); diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp index 079ad615..884551e2 100644 --- a/src/framework/parameters.cpp +++ b/src/framework/parameters.cpp @@ -31,10 +31,10 @@ namespace ntt { template - auto get_dx0_V0( - const std::vector& resolution, - const boundaries_t& extent, - const std::map& params) -> std::pair { + auto get_dx0_V0(const std::vector& resolution, + const boundaries_t& extent, + const std::map& params) + -> std::pair { const auto metric = M(resolution, extent, params); const auto dx0 = metric.dxMin(); coord_t x_corner { ZERO }; @@ -445,6 +445,35 @@ namespace ntt { /* [particles] ---------------------------------------------------------- */ set("particles.clear_interval", toml::find_or(toml_data, "particles", "clear_interval", defaults::clear_interval)); + const auto species_tab = toml::find_or(toml_data, + "particles", + "species", + toml::array {}); + std::vector species = get>( + "particles.species"); + raise::ErrorIf(species_tab.size() != species.size(), + "number of species changed after restart", + HERE); + + std::vector new_species; + + spidx_t idxM1 = 0; + for (const auto& sp : species_tab) { + const auto maxnpart_real = toml::find(sp, "maxnpart"); + const auto maxnpart = static_cast(maxnpart_real); + const auto particle_species = species[idxM1]; + new_species.emplace_back(particle_species.index(), + particle_species.label(), + particle_species.mass(), + particle_species.charge(), + maxnpart, + particle_species.pusher(), + particle_species.use_gca(), + particle_species.cooling(), + particle_species.npld()); + idxM1++; + } + set("particles.species", new_species); /* [output] ------------------------------------------------------------- */ // fields