Merge branch 'dev/main' of github.com:haykh/entity into dev/output

Author: haykh

src/framework/archetypes.hpp

@@ -5,6 +5,7 @@
 
 #include "fields.h"
 #include "meshblock.h"
+#include "qmath.h"
 #include "sim_params.h"
 
 #ifdef NTTINY_ENABLED
@@ -94,6 +95,88 @@ namespace ntt {
     }
   };
 
+  template <Dimension D, SimulationEngine S>
+  struct Maxwellian {
+    Maxwellian(const Meshblock<D, S>& mblock) : pool { *(mblock.random_pool_ptr) } {}
+    // Juttner-Synge distribution
+    Inline void JS(vec_t<Dim3>& v, const real_t& temp) const {
+      typename RandomNumberPool_t::generator_type rand_gen = pool.get_state();
+      real_t                                      u { ZERO }, eta { ZERO }, theta { ZERO };
+      real_t                                      X1 { ZERO }, X2 { ZERO };
+      if (temp < 0.1) {
+        // Juttner-Synge distribution using the Box-Muller method - non-relativistic
+        while (AlmostEqual(u, ZERO)) {
+          u = rand_gen.frand();
+        }
+        eta = math::sqrt(-TWO * math::log(u));
+        while (AlmostEqual(theta, ZERO)) {
+          theta = constant::TWO_PI * rand_gen.frand();
+        }
+        u = eta * math::cos(theta) * math::sqrt(temp);
+      } else {
+        // Juttner-Synge distribution using the Sobol method - relativistic
+        u = ONE;
+        while (SQR(eta) <= SQR(u)) {
+          while (AlmostEqual(X1, ZERO)) {
+            X1 = rand_gen.frand() * rand_gen.frand() * rand_gen.frand();
+          }
+          u  = -temp * math::log(X1);
+          X1 = rand_gen.frand();
+          while (AlmostEqual(X1, 0)) {
+            X1 = rand_gen.frand();
+          }
+          eta = u - temp * math::log(X1);
+        }
+      }
+      X1   = rand_gen.frand();
+      X2   = rand_gen.frand();
+      v[0] = u * (TWO * X1 - ONE);
+      v[2] = TWO * u * math::sqrt(X1 * (ONE - X1));
+      v[1] = v[2] * math::cos(constant::TWO_PI * X2);
+      v[2] = v[2] * math::sin(constant::TWO_PI * X2);
+      pool.free_state(rand_gen);
+    }
+    // Boost a symmetric distribution to a relativistic speed using flipping method
+    // https://arxiv.org/pdf/1504.03910.pdf
+    Inline void boost(vec_t<Dim3>&  v,
+                      const real_t& boost_vel,
+                      const short&  boost_direction) const {
+      typename RandomNumberPool_t::generator_type rand_gen = pool.get_state();
+      real_t boost_beta { boost_vel / math::sqrt(ONE + SQR(boost_vel)) };
+      real_t boost_gamma { boost_vel / boost_beta };
+      real_t ut { math::sqrt(ONE + SQR(v[0]) + SQR(v[1]) + SQR(v[2])) };
+      if (-boost_beta * v[boost_direction] > ut * rand_gen.frand()) {
+        v[boost_direction] = -v[boost_direction];
+      }
+      pool.free_state(rand_gen);
+      v[boost_direction] = boost_gamma * (v[boost_direction] + boost_beta * ut);
+    }
+
+    Inline void operator()(vec_t<Dim3>&  v,
+                           const real_t& temp,
+                           const real_t& boost_vel       = ZERO,
+                           const short&  boost_direction = 0) const {
+      if (AlmostEqual(temp, ZERO)) {
+        v[0] = ZERO;
+        v[1] = ZERO;
+        v[2] = ZERO;
+      } else {
+        JS(v, temp);
+      }
+      if (!AlmostEqual(boost_vel, ZERO)) {
+        if (boost_direction < 0) {
+          boost(v, -boost_vel, -boost_direction - 1);
+        } else if (boost_direction > 0) {
+          boost(v, boost_vel, boost_direction - 1);
+        }
+        // no boost when boost_direction == 0
+      }
+    }
+
+  private:
+    RandomNumberPool_t pool;
+  };
+
   /* -------------------------------------------------------------------------- */
   /*                            Spatial distribution                            */
   /* -------------------------------------------------------------------------- */

src/framework/meshblock/meshblock.h

@@ -11,6 +11,10 @@
 #include <vector>
 
 namespace ntt {
+  namespace dir {
+    enum Direction { x = 1, y = 2, z = 3, r = 1, theta = 2, phi = 3, x1 = 1, x2 = 2, x3 = 3 };
+  }    // namespace dir
+
   enum class CellLayer {
     allLayer,
     activeLayer,

src/pic/pgen/inputs/weibel.toml

@@ -15,7 +15,7 @@ deposit_ON      = true
 current_filters = 25
 
 [units]
-ppc0       = 32.0
+ppc0       = 16.0
 larmor0    = 1.0
 skindepth0 = 1.0
 
@@ -51,10 +51,10 @@ maxnpart = 1e8
 pusher   = "Boris"
 
 [problem]
-udrift_1 = -10.0
-udrift_2 = -10.0
-udrift_3 = 10.0
-udrift_4 = 10.0
+drift_p  = 10.0
+drift_b  = 10.0
+temperature_p = 0.1
+temperature_b = 0.1
 
 [visualization]
 fields        = ["Bx", "By", "mass_density"]

src/pic/pgen/weibel.hpp

@@ -12,29 +12,31 @@
 namespace ntt {
 
   template <Dimension D, SimulationEngine S>
-  struct Drift : public EnergyDistribution<D, S> {
-    Drift(const SimulationParams& params, const Meshblock<D, S>& mblock)
+  struct WeibelInit : public EnergyDistribution<D, S> {
+    WeibelInit(const SimulationParams& params, const Meshblock<D, S>& mblock)
       : EnergyDistribution<D, S>(params, mblock),
-        udrift1 { params.get<real_t>("problem", "udrift1", -10.0) },
-        udrift2 { params.get<real_t>("problem", "udrift2", -10.0) },
-        udrift3 { params.get<real_t>("problem", "udrift3", 10.0) },
-        udrift4 { params.get<real_t>("problem", "udrift4", 10.0) } {}
+        maxwellian { mblock },
+        drift_p { params.get<real_t>("problem", "drift_p", 10.0) },
+        drift_b { params.get<real_t>("problem", "drift_b", 10.0) },
+        temp_p { params.get<real_t>("problem", "temperature_p", 0.0) },
+        temp_b { params.get<real_t>("problem", "temperature_b", 0.0) } {}
     Inline void operator()(const coord_t<D>&,
                            vec_t<Dim3>& v,
                            const int&   species) const override {
       if (species == 1) {
-        v[2] = udrift1;
+        maxwellian(v, temp_p, drift_p, -dir::z);
       } else if (species == 2) {
-        v[2] = udrift2;
+        maxwellian(v, temp_p, drift_p, -dir::z);
       } else if (species == 3) {
-        v[2] = udrift3;
+        maxwellian(v, temp_b, drift_b, dir::z);
       } else if (species == 4) {
-        v[2] = udrift4;
+        maxwellian(v, temp_b, drift_b, dir::z);
       }
     }
 
   private:
-    const real_t udrift1, udrift2, udrift3, udrift4;
+    const Maxwellian<D, S> maxwellian;
+    const real_t           drift_p, drift_b, temp_p, temp_b;
   };
 
   template <Dimension D, SimulationEngine S>