dev/hybrid

dependencies.py

@@ -367,8 +367,10 @@ def on_install_confirmed(settings: Settings) -> None:
                 "kokkos",
                 settings.kokkos_version,
                 settings.kokkos_backend, 
-                settings.kokkos_arch.strip().lower(),
             )
+            if settings.kokkos_arch.strip().lower():
+                kokkos_modfile_file = os.path.join(kokkos_modfile_file, settings.kokkos_arch.strip().lower())
+
             os.makedirs(os.path.dirname(kokkos_modfile_file), exist_ok=True)
             if os.path.exists(kokkos_modfile_file) and not settings.overwrite:
                 raise FileExistsError(

input.example.toml

@@ -7,7 +7,7 @@
   # Simulation engine to use
   #   @required
   #   @type: string
-  #   @enum: "SRPIC", "GRPIC"
+  #   @enum: "SRPIC", "GRPIC", "HYBRID"
   engine = ""
   # Max runtime in physical (code) units
   #   @required
@@ -224,6 +224,78 @@
   #     @from: `scales.larmor0`
   #     @value: `1 / larmor0`
 
+# Electron-fluid closure for the hybrid Ohm's law
+#   @note: Only read by the `HYBRID` engine (ignored by `SRPIC`/`GRPIC`)
+[hybrid]
+  # Adiabatic index of the massless electron fluid
+  #   @type: float [> 0.0]
+  #   @default: 1.66667 (= 5/3)
+  #   @note: `1.0` corresponds to an isothermal closure
+  gamma_ad = ""
+  # Electron temperature `T_e` in physical (code) units
+  #   @type: float [>= 0.0]
+  #   @default: 0.0
+  #   @note: `0.0` corresponds to cold electrons
+  theta0 = ""
+  # Vacuum threshold for the hybrid Ohm's law, in units of `n0`
+  #   @type: float [> 0.0]
+  #   @default: 1e-3
+  #   @note: Above it `E` follows the usual `[...]/N`; below it `E` ramps continuously
+  #          to zero so plasma-free cells keep `B` frozen instead of amplifying the
+  #          right-hand side by `1/dens_min` (statistically-empty or unfilled ghost
+  #          cells, a fraction `e^-ppc` for finite `ppc`)
+  dens_min = ""
+  # Cap on the local whistler Courant the Ohm's-law Hall term may imply
+  #   @type: float
+  #   @default: 0.5
+  #   @note: The implied whistler speed `v_w ~ |B|/N` can outrun `dx/dt` at shock
+  #          overshoots and magnetic cavities; the Hall term is scaled so its
+  #          whistler Courant never exceeds this. Identity in resolvable cells
+  #   @note: `<= 0.0` disables the limiter
+  hall_lim = ""
+  # Vacuum resistivity, in physical (code) units
+  #   @type: float [>= 0.0]
+  #   @default: 0.0
+  #   @note: `> 0.0` makes the sub-`dens_min` vacuum resistive: `E` crossfades to
+  #          `eta * curl B` below the threshold, so vacuum regions diffuse toward a
+  #          curl-free field instead of freezing. Needed when fast plasma sweeps
+  #          along a vacuum boundary, where a hard `E = 0` mask pumps `B` at `~ u/dx`
+  #   @note: Internally clamped to the explicit-diffusion stability limit; any value is safe
+  #   @note: `0.0` disables it (hard `E = 0` vacuum)
+  resist_vac = ""
+  # Hyper-resistivity coefficient `eta_H` (adds `-eta_H grad^4 B` to the field advance)
+  #   @type: float [>= 0.0]
+  #   @default: 0.0
+  #   @note: Damps the few-cell band the moments filter leaves without kinetic
+  #          feedback, while leaving physical scales (`k dx < 0.3`) essentially
+  #          untouched; size so `eta_H (pi/dx)^4` exceeds the drive
+  #   @note: Internally clamped to the explicit `grad^4` stability limit; any value is safe
+  #   @note: `0.0` disables it
+  resist_hyper = ""
+  # Pegasus-style sub-cycled magnetic-field advance
+  #   @type: bool
+  #   @default: true
+  #   @note: Integrates the whistler-stiff Ohm/Faraday loop with adaptive SSP-RK3
+  #          sub-steps at its own CFL; `false` uses the legacy 3-push scheme
+  subcycle = ""
+  # Target whistler Courant per field sub-step
+  #   @type: float [> 0.0]
+  #   @default: 0.5
+  #   @note: Only read when `subcycle = true`
+  subcycle_courant = ""
+  # Cap on the number of field sub-steps per advance
+  #   @type: int [>= 1]
+  #   @default: 64
+  #   @note: Beyond the cap the per-cell Hall limiter (`hall_lim`) takes over
+  subcycle_max = ""
+  # Characteristic flow speed for the hybrid CFL, in physical (code) units
+  #   @type: float [>= 0.0]
+  #   @default: 0.0
+  #   @note: Added to the Alfven + whistler signal speeds in the Pegasus CFL
+  #          `dt <= CFL * min(dx0 / v_max, 2*pi/Omega)` [Kunz, Stone & Bai (2014), §3.5]
+  #   @note: `0.0` lets the Alfven + whistler speeds set `dt`; raise it for super-whistler flows
+  v_max = ""
+
 [radiation]
   [radiation.drag]
     [radiation.drag.synchrotron]

pgens/hybrid/hybrid.toml

@@ -0,0 +1,50 @@
+[simulation]
+  name    = "hybrid"
+  engine  = "hybrid"
+  runtime = 1.0
+
+[grid]
+  resolution = [512, 256]
+  extent     = [[0.0, 2.0], [0.0, 1.0]]
+
+  [grid.metric]
+    metric = "minkowski"
+
+  [grid.boundaries]
+    fields    = [["PERIODIC"], ["PERIODIC"]]
+    particles = [["PERIODIC"], ["PERIODIC"]]
+
+[scales]
+  larmor0    = 1.0
+  skindepth0 = 1.0
+
+[algorithms]
+  [algorithms.timestep]
+    CFL = 0.5
+
+[particles]
+  ppc0 = 10.0
+
+  [[particles.species]]
+    label    = "ions"
+    mass     = 1.0
+    charge   = 1.0
+    maxnpart = 1e5
+
+[setup]
+
+[output]
+  interval_time = 0.1
+
+  [output.fields]
+    quantities = ["N_1", "B"]
+
+  [output.particles]
+    species = [1]
+    stride  = 10
+
+[checkpoint]
+  keep = 0
+
+[diagnostics]
+  colored_stdout = true

pgens/hybrid/pgen.hpp

@@ -0,0 +1,30 @@
+#ifndef PROBLEM_GENERATOR_H
+#define PROBLEM_GENERATOR_H
+
+#include "enums.h"
+#include "global.h"
+
+#include "traits/pgen.h"
+
+#include "framework/domain/metadomain.h"
+#include "framework/parameters/parameters.h"
+
+namespace user {
+
+  template <SimEngine::type S, class M>
+  struct PGen {
+    static constexpr auto D { M::Dim };
+
+    static constexpr auto engines =
+      ::traits::pgen::compatible_with<SimEngine::HYBRID> {};
+    static constexpr auto metrics =
+      ::traits::pgen::compatible_with<Metric::Minkowski> {};
+    static constexpr auto dimensions =
+      ::traits::pgen::compatible_with<Dim::_1D, Dim::_2D, Dim::_3D> {};
+
+    PGen(const SimulationParams&, const Metadomain<S, M>&) {}
+  };
+
+} // namespace user
+
+#endif

pgens/shock/pgen.hpp

@@ -259,7 +259,7 @@ namespace user {
                              domain.fields.em,
                              init_flds,
                              domain.mesh.metric });
-      metadomain.CommunicateFields(domain, Comm::E | Comm::B);
+      metadomain.CommunicateFields(domain, Comm::EM);
 
       /*
         tag particles inside the injection zone as dead

pgens/streaming/pgen.hpp

@@ -54,7 +54,8 @@ namespace user {
   struct PGen {
     static constexpr auto D { M::Dim };
     // compatibility traits for the problem generator
-    static constexpr auto engines = ::traits::pgen::compatible_with<SimEngine::SRPIC> {};
+    static constexpr auto engines =
+      ::traits::pgen::compatible_with<SimEngine::SRPIC, SimEngine::HYBRID> {};
     static constexpr auto metrics =
       ::traits::pgen::compatible_with<Metric::Minkowski> {};
     static constexpr auto dimensions =
@@ -87,13 +88,13 @@ namespace user {
       raise::ErrorIf(nspec % 2 != 0,
                      "Number of species must be even for this setup",
                      HERE);
-      for (auto n = 0u; n < nspec; n += 2) {
-        raise::ErrorIf(
-          global_domain.species_params()[n].charge() !=
-            -global_domain.species_params()[n + 1].charge(),
-          "Charges of i-th and i+1-th species must be opposite for this setup",
-          HERE);
-      }
+      // for (auto n = 0u; n < nspec; n += 2) {
+      //   raise::ErrorIf(
+      //     global_domain.species_params()[n].charge() !=
+      //       -global_domain.species_params()[n + 1].charge(),
+      //     "Charges of i-th and i+1-th species must be opposite for this
+      //     setup", HERE);
+      // }
       for (auto* specs :
            { &drifts_in_x, &drifts_in_y, &drifts_in_z, &temperatures }) {
         if (specs->empty()) {
@@ -126,13 +127,29 @@ namespace user {
         const auto drift_2 = prmvec_t { drifts_in_x[n + 1],
                                         drifts_in_y[n + 1],
                                         drifts_in_z[n + 1] };
-        arch::InjectUniformMaxwellians<S, M>(
+        const auto edist1  = arch::energy_dist::MaxwellianNonRel<M::Dim> {
+          domain.random_pool(),
+          temperatures[n],
+          drift_1
+        };
+        const auto edist2 = arch::energy_dist::MaxwellianNonRel<M::Dim> {
+          domain.random_pool(),
+          temperatures[n + 1],
+          drift_2
+        };
+        arch::InjectUniform<S, M, decltype(edist1), decltype(edist2)>(
           params,
           domain,
-          densities[n / 2],
-          { temperatures[n], temperatures[n + 1] },
-          { n + 1, n + 2 },
-          { drift_1, drift_2 });
+          { 1u, 2u },
+          { edist1, edist2 },
+          ONE);
+        // arch::InjectUniformMaxwellians<S, M>(
+        //   params,
+        //   domain,
+        //   densities[n / 2],
+        //   { temperatures[n], temperatures[n + 1] },
+        //   { n + 1, n + 2 },
+        //   { drift_1, drift_2 });
       }
     }
   };

src/archetypes/energy_dist.h

@@ -76,6 +76,38 @@ namespace arch::energy_dist {
     random_number_pool_t pool;
   };
 
+  Inline void NonRelMaxwellian(vec_t<Dim::_3D>&            v,
+                               real_t                      temp,
+                               const random_number_pool_t& pool) {
+    auto   rand_gen = pool.get_state();
+    // Juttner-Synge distribution using the Box-Muller method - non-relativistic
+    real_t randX1   = Random<real_t>(rand_gen);
+    while (cmp::AlmostZero(randX1)) {
+      randX1 = Random<real_t>(rand_gen);
+    }
+    randX1        = math::sqrt(-TWO * math::log(randX1));
+    real_t randX2 = static_cast<real_t>(constant::TWO_PI) *
+                    Random<real_t>(rand_gen);
+    v[0] = randX1 * math::cos(randX2) * math::sqrt(temp);
+
+    randX1 = Random<real_t>(rand_gen);
+    while (cmp::AlmostZero(randX1)) {
+      randX1 = Random<real_t>(rand_gen);
+    }
+    randX1 = math::sqrt(-TWO * math::log(randX1));
+    randX2 = static_cast<real_t>(constant::TWO_PI) * Random<real_t>(rand_gen);
+    v[1]   = randX1 * math::cos(randX2) * math::sqrt(temp);
+
+    randX1 = Random<real_t>(rand_gen);
+    while (cmp::AlmostZero(randX1)) {
+      randX1 = Random<real_t>(rand_gen);
+    }
+    randX1 = math::sqrt(-TWO * math::log(randX1));
+    randX2 = static_cast<real_t>(constant::TWO_PI) * Random<real_t>(rand_gen);
+    v[2]   = randX1 * math::cos(randX2) * math::sqrt(temp);
+    pool.free_state(rand_gen);
+  }
+
   Inline void JuttnerSinge(vec_t<Dim::_3D>&            v,
                            real_t                      temp,
                            const random_number_pool_t& pool) {
@@ -291,6 +323,48 @@ namespace arch::energy_dist {
     short drift_dir { 0 };
   };
 
+  template <Dimension D>
+  struct MaxwellianNonRel {
+    MaxwellianNonRel(random_number_pool_t& pool,
+                     real_t                temperature,
+                     const std::vector<real_t>& drift_vel = { ZERO, ZERO, ZERO })
+      : pool { pool }
+      , temperature { temperature } {
+      raise::ErrorIf(drift_vel.size() != 3,
+                     "Maxwellian: Drift velocity must be a 3D vector",
+                     HERE);
+      raise::ErrorIf(temperature < ZERO,
+                     "Maxwellian: Temperature must be non-negative",
+                     HERE);
+      if (not cmp::AlmostZero_host(NORM(drift_vel[0], drift_vel[1], drift_vel[2]))) {
+        drift_vel_x1 = drift_vel[0];
+        drift_vel_x2 = drift_vel[1];
+        drift_vel_x3 = drift_vel[2];
+      }
+    }
+
+    Inline void operator()(const coord_t<D>&, vec_t<Dim::_3D>& v) const {
+
+      if (cmp::AlmostZero(temperature)) {
+        v[0] = ZERO;
+        v[1] = ZERO;
+        v[2] = ZERO;
+      } else {
+        NonRelMaxwellian(v, temperature, pool);
+      }
+      v[0] += drift_vel_x1;
+      v[1] += drift_vel_x2;
+      v[2] += drift_vel_x3;
+    }
+
+  private:
+    random_number_pool_t pool;
+
+    const real_t temperature;
+
+    real_t drift_vel_x1 { ZERO }, drift_vel_x2 { ZERO }, drift_vel_x3 { ZERO };
+  };
+
 } // namespace arch::energy_dist
 
 #endif // ARCHETYPES_ENERGY_DIST_HPP