Output to separate files

.gitignore

@@ -51,6 +51,7 @@ venv/
 # CMake testing files
 Testing/
 
+tags
 .clangd
 .schema.json
 *_old/

CMakeLists.txt

@@ -55,7 +55,9 @@ if(${DEBUG} STREQUAL "OFF")
   set(CMAKE_BUILD_TYPE
       Release
       CACHE STRING "CMake build type")
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG")
+  set(CMAKE_CXX_FLAGS
+      "${CMAKE_CXX_FLAGS} -DNDEBUG -Wno-unused-local-typedefs -Wno-unknown-cuda-version"
+  )
 else()
   set(CMAKE_BUILD_TYPE
       Debug
@@ -64,8 +66,6 @@ else()
       "${CMAKE_CXX_FLAGS} -DDEBUG -Wall -Wextra -Wno-unknown-pragmas")
 endif()
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unused-local-typedefs")
-
 # options
 set(precisions
     "single" "double"

input.example.toml

@@ -90,10 +90,10 @@
     # Boundary conditions for fields:
     #   @required
     #   @type: 1/2/3-size array of string tuples, each of size 1 or 2
-    #   @valid: "PERIODIC", "ABSORB", "ATMOSPHERE", "CUSTOM", "HORIZON"
+    #   @valid: "PERIODIC", "ABSORB", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
     #   @example: [["CUSTOM", "ABSORB"]] (for 2D spherical [[rmin, rmax]])
-    #   @note: When periodic in any of the directions, you should only set one value [..., ["PERIODIC"], ...]
-    #   @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]) [["ATMOSPHERE", "ABSORB"]]
+    #   @note: When periodic in any of the directions, you should only set one value: [..., ["PERIODIC"], ...]
+    #   @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "ABSORB"]]
     #   @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]]
     fields = ""
     # Boundary conditions for fields:
@@ -319,6 +319,10 @@
   #   @default: -1.0 (disabled)
   #   @note: When `interval_time` < 0, the output is controlled by `interval`, otherwise by `interval_time`
   interval_time = ""
+  # Whether to output each timestep into separate files:
+  #   @type: bool
+  #   @default: true
+  separate_files = ""
 
   [output.fields]
     # Toggle for the field output:

setups/srpic/shock/pgen.hpp

@@ -5,6 +5,8 @@
 #include "global.h"
 
 #include "arch/traits.h"
+#include "utils/error.h"
+#include "utils/numeric.h"
 
 #include "archetypes/energy_dist.h"
 #include "archetypes/particle_injector.h"
@@ -14,6 +16,53 @@
 namespace user {
   using namespace ntt;
 
+  template <Dimension D>
+  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 btheta: magnetic field polar angle
+      @param bphi: magnetic field azimuthal angle
+      @param drift_ux: drift velocity in the x direction
+    */
+    InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux)
+      : Bmag { bmag * static_cast<real_t>(convert::deg2rad) }
+      , Btheta { btheta * static_cast<real_t>(convert::deg2rad) }
+      , Bphi { bphi * static_cast<real_t>(convert::deg2rad) }
+      , Vx { drift_ux } {}
+
+    // magnetic field components
+    Inline auto bx1(const coord_t<D>&) const -> real_t {
+      return Bmag * math::cos(Btheta);
+    }
+
+    Inline auto bx2(const coord_t<D>&) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+    Inline auto bx3(const coord_t<D>&) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    // electric field components
+    Inline auto ex1(const coord_t<D>&) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex2(const coord_t<D>&) const -> real_t {
+      return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    Inline auto ex3(const coord_t<D>&) const -> real_t {
+      return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+  private:
+    const real_t Btheta, Bphi, Vx, Bmag;
+  };
+
   template <SimEngine::type S, class M>
   struct PGen : public arch::ProblemGenerator<S, M> {
     // compatibility traits for the problem generator
@@ -30,20 +79,42 @@ namespace user {
 
     const real_t drift_ux, temperature;
 
+    const real_t  Btheta, Bphi, Bmag;
+    InitFields<D> init_flds;
+
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
-      : arch::ProblemGenerator<S, M>(p)
+      : arch::ProblemGenerator<S, M> { p }
       , drift_ux { p.template get<real_t>("setup.drift_ux") }
-      , temperature { p.template get<real_t>("setup.temperature") } {}
+      , temperature { p.template get<real_t>("setup.temperature") }
+      , Bmag { p.template get<real_t>("setup.Bmag", ZERO) }
+      , Btheta { p.template get<real_t>("setup.Btheta", ZERO) }
+      , Bphi { p.template get<real_t>("setup.Bphi", ZERO) }
+      , init_flds { Bmag, Btheta, Bphi, drift_ux } {}
 
     inline PGen() {}
 
+    auto FixField(const em& comp) const -> real_t {
+      if (comp == em::ex2) {
+        return init_flds.ex2({ ZERO });
+      } else if (comp == em::ex3) {
+        return init_flds.ex3({ ZERO });
+      } else if (comp == em::bx1) {
+        return init_flds.bx1({ ZERO });
+      } else {
+        raise::Error("Other components should not be requested when BC is in X",
+                     HERE);
+        return ZERO;
+      }
+    }
+
     inline void InitPrtls(Domain<S, M>& local_domain) {
       const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
                                                       local_domain.random_pool,
                                                       temperature,
                                                       -drift_ux,
                                                       in::x1);
-      const auto injector    = arch::UniformInjector<S, M, arch::Maxwellian>(
+
+      const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
         energy_dist,
         { 1, 2 });
       arch::InjectUniform<S, M, arch::UniformInjector<S, M, arch::Maxwellian>>(

setups/srpic/shock/shock.py

@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 import matplotlib as mpl
 
-data = nt2r.Data("shock-03.h5")
+data = nt2r.Data("shock.h5")
 
 
 def frame(ti, f):
@@ -55,7 +55,7 @@ def frame(ti, f):
     axs = [fig.add_subplot(gs[i]) for i in range(len(quantities))]
 
     for ax, q in zip(axs, quantities):
-        q["compute"](f).coarsen(x=2, y=2).mean().plot(
+        q["compute"](f.isel(t=ti)).plot(
             ax=ax,
             cmap=q["cmap"],
             norm=q["norm"],

setups/srpic/shock/shock.toml

@@ -11,7 +11,7 @@
     metric = "minkowski"
 
   [grid.boundaries]
-    fields    = [["CONDUCTOR", "ABSORB"], ["PERIODIC"]]
+    fields = [["FIXED", "ABSORB"], ["PERIODIC"]]
     particles = [["REFLECT", "ABSORB"], ["PERIODIC"]]
 
 [scales]
@@ -42,6 +42,9 @@
 [setup]
   drift_ux    = 0.1
   temperature = 1e-3
+  Bmag = 1.0
+  Btheta = 0.0
+  Bphi = 0.0
 
 [output]
   interval_time = 0.1