v1.2.3 Release Candidate

.gitignore

@@ -39,6 +39,7 @@ venv/
 *.png
 *.mov
 *.mp4
+!pgens/**/*.png
 
 # Accidental files
 *.xc

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) ")

input.example.toml

@@ -467,8 +467,9 @@
     # Field quantities to output
     #   @type: array<string>
     #   @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<string>

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"

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 <string>
+
 namespace user {
   using namespace ntt;
 
+  enum class FieldGeometry {
+    dipole,
+    monopole
+  };
+
   template <Dimension D>
   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<D>& 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<D>& 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 <Dimension D>
   struct DriveFields : public InitFields<D> {
-    DriveFields(real_t time, real_t bsurf, real_t rstar, real_t omega)
-      : InitFields<D> { bsurf, rstar }
+    DriveFields(real_t             time,
+                real_t             bsurf,
+                real_t             rstar,
+                real_t             omega,
+                const std::string& field_geometry)
+      : InitFields<D> { bsurf, rstar, field_geometry }
       , time { time }
       , Omega { omega } {}
 
@@ -73,25 +98,27 @@ namespace user {
     using arch::ProblemGenerator<S, M>::C;
     using arch::ProblemGenerator<S, M>::params;
 
-    const real_t  Bsurf, Rstar, Omega;
-    InitFields<D> init_flds;
+    const real_t      Bsurf, Rstar, Omega;
+    const std::string field_geom;
+    InitFields<D>     init_flds;
 
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
       : arch::ProblemGenerator<S, M>(p)
       , Bsurf { p.template get<real_t>("setup.Bsurf", ONE) }
       , Rstar { m.mesh().extent(in::x1).first }
       , Omega { static_cast<real_t>(constant::TWO_PI) /
                 p.template get<real_t>("setup.period", ONE) }
-      , init_flds { Bsurf, Rstar } {}
+      , field_geom { p.template get<std::string>("setup.field_geometry", "dipole") }
+      , init_flds { Bsurf, Rstar, field_geom } {}
 
     inline PGen() {}
 
     auto AtmFields(real_t time) const -> DriveFields<D> {
-      return DriveFields<D> { time, Bsurf, Rstar, Omega };
+      return DriveFields<D> { time, Bsurf, Rstar, Omega, field_geom };
     }
 
     auto MatchFields(real_t) const -> InitFields<D> {
-      return InitFields<D> { Bsurf, Rstar };
+      return InitFields<D> { Bsurf, Rstar, field_geom };
     }
   };
 

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")

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<S, M>& local_domain) {
       // background
-      const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
-                                                      local_domain.random_pool,
-                                                      bg_temperature);
-      const auto injector    = arch::UniformInjector<S, M, arch::Maxwellian>(
-        energy_dist,
-        { 1, 2 });
-      arch::InjectUniform<S, M, arch::UniformInjector<S, M, arch::Maxwellian>>(
-        params,
-        local_domain,
-        injector,
-        ONE);
+      arch::InjectUniformMaxwellian<S, M>(params,
+                                          local_domain,
+                                          ONE,
+                                          bg_temperature,
+                                          { 1, 2 });
 
       const auto sigma = params.template get<real_t>("scales.sigma0");
       const auto c_omp = params.template get<real_t>("scales.skindepth0");

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