Add bulk velocity as moments to output

input.example.toml

@@ -329,7 +329,7 @@
     # Field quantities to output:
     #   @type: array of strings
     #   @valid: fields: "E", "B", "J", "divE"
-    #   @valid: moments: "Rho", "Charge", "N", "Nppc", "T0i", "Tij"
+    #   @valid: moments: "Rho", "Charge", "N", "Nppc", "T0i", "Tij", "Vi"
     #   @valid: for GR: "D", "H", "divD", "A"
     #   @default: []
     #   @note: For T, you can use unspecified indices, e.g., Tij, T0i, or specific ones, e.g., Ttt, T00, T02, T23

src/framework/domain/output.cpp

@@ -294,6 +294,14 @@ namespace ntt {
                                                  {},
                                                  local_domain->fields.bckp,
                                                  c);
+            } else if (fld.id() == FldsID::V) {
+              ComputeMoments<S, M, FldsID::V>(params,
+                                                 local_domain->mesh,
+                                                 local_domain->species,
+                                                 fld.species,
+                                                 fld.comp[0],
+                                                 local_domain->fields.bckp,
+                                                 c);
             } else {
               raise::Error("Wrong moment requested for output", HERE);
             }

src/global/enums.h

@@ -289,16 +289,17 @@ namespace ntt {
       N       = 12,
       Nppc    = 13,
       Custom  = 14,
+      V       = 15,
     };
 
     constexpr FldsID(uint8_t c) : enums_hidden::BaseEnum<FldsID> { c } {}
 
     static constexpr type variants[]      = { E, divE, D,   divD,   B, H,    J,
-                                              A, T,    Rho, Charge, N, Nppc, Custom };
+                                              A, T,    Rho, Charge, N, Nppc, Custom , V};
     static constexpr const char* lookup[] = { "e",    "dive",  "d",      "divd",
                                               "b",    "h",     "j",      "a",
                                               "t",    "rho",   "charge", "n",
-                                              "nppc", "custom" };
+                                              "nppc", "custom", "v" };
     static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
   };
 

src/global/tests/enums.cpp

@@ -68,7 +68,7 @@ auto main() -> int {
 
   enum_str_t all_out_flds = { "e",      "dive", "d",    "divd",  "b",
                               "h",      "j",    "a",    "t",     "rho",
-                              "charge", "n",    "nppc", "custom" };
+                              "charge", "n",    "nppc", "custom" , "v"};
 
   checkEnum<Coord>(all_coords);
   checkEnum<Metric>(all_metrics);

src/kernels/particle_moments.hpp

@@ -41,7 +41,7 @@ namespace kernel {
     static constexpr auto D = M::Dim;
 
     static_assert((F == FldsID::Rho) || (F == FldsID::Charge) ||
-                    (F == FldsID::N) || (F == FldsID::Nppc) || (F == FldsID::T),
+                    (F == FldsID::N) || (F == FldsID::Nppc) || (F == FldsID::T) || (F == FldsID::V),
                   "Invalid field ID");
 
     const unsigned short     c1, c2;
@@ -89,7 +89,7 @@ namespace kernel {
                            std::size_t                        ni2,
                            real_t                             inv_n0,
                            unsigned short                     window)
-      : c1 { (components.size() == 2) ? components[0]
+      : c1 { (components.size() > 0) ? components[0]
                                       : static_cast<unsigned short>(0) }
       , c2 { (components.size() == 2) ? components[1]
                                       : static_cast<unsigned short>(0) }
@@ -205,6 +205,64 @@ namespace kernel {
         coeff = contrib;
       }
 
+      if constexpr (F == FldsID::V) {
+        real_t          gamma { ZERO };
+        // for stress-energy tensor
+        vec_t<Dim::_3D> u_Phys { ZERO };
+        if constexpr (S == SimEngine::SRPIC) {
+          // SR
+          // stress-energy tensor for SR is computed in the tetrad (hatted) basis
+          if constexpr (M::CoordType == Coord::Cart) {
+            u_Phys[0] = ux1(p);
+            u_Phys[1] = ux2(p);
+            u_Phys[2] = ux3(p);
+          } else {
+            static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D");
+            coord_t<M::PrtlDim> x_Code { ZERO };
+            x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
+            x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
+            if constexpr (D == Dim::_3D) {
+              x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
+            } else {
+              x_Code[2] = phi(p);
+            }
+            metric.template transform_xyz<Idx::XYZ, Idx::T>(
+              x_Code,
+              { ux1(p), ux2(p), ux3(p) },
+              u_Phys);
+          }
+          if (mass == ZERO) {
+            gamma = NORM(u_Phys[0], u_Phys[1], u_Phys[2]);
+          } else {
+            gamma = math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]));
+          }
+        } else {
+          // GR
+          // stress-energy tensor for GR is computed in contravariant basis
+          static_assert(D != Dim::_1D, "GRPIC 1D");
+          coord_t<D> x_Code { ZERO };
+          x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
+          x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
+          if constexpr (D == Dim::_3D) {
+            x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
+          }
+          vec_t<Dim::_3D> u_Cntrv { ZERO };
+          // compute u_i u^i for energy
+          metric.template transform<Idx::D, Idx::U>(x_Code,
+                                                    { ux1(p), ux2(p), ux3(p) },
+                                                    u_Cntrv);
+          gamma = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + u_Cntrv[2] * ux3(p);
+          if (mass == ZERO) {
+            gamma = math::sqrt(gamma);
+          } else {
+            gamma = math::sqrt(ONE + gamma);
+          }
+          metric.template transform<Idx::U, Idx::PU>(x_Code, u_Cntrv, u_Phys);
+        }
+        // compute the corresponding moment
+        coeff = u_Phys[c1 - 1] / gamma;
+      } 
+
       if constexpr (F != FldsID::Nppc) {
         // for nppc calculation ...
         // ... do not take volume, weights or smoothing into account

src/output/fields.cpp

@@ -44,6 +44,9 @@ namespace out {
     } else if (id() == FldsID::T) {
       // energy-momentum tensor
       comp = InterpretComponents({ name.substr(1, 1), name.substr(2, 1) });
+    } else if (id() == FldsID::V) {
+      // energy-momentum tensor
+      comp = InterpretComponents({ name.substr(1, 1) });
     } else {
       // scalar (Rho, divE, Custom, etc.)
       comp = {};

src/output/fields.h

@@ -43,7 +43,7 @@ namespace out {
     [[nodiscard]]
     auto is_moment() const -> bool {
       return (id() == FldsID::T || id() == FldsID::Rho || id() == FldsID::Nppc ||
-              id() == FldsID::N || id() == FldsID::Charge);
+              id() == FldsID::N || id() == FldsID::Charge || id() == FldsID::V);
     }
 
     [[nodiscard]]
@@ -94,6 +94,8 @@ namespace out {
           tmp += m_name.substr(1, 2);
         } else if (id() == FldsID::A) {
           tmp += "3";
+        } else if (id() == FldsID::V) {
+          tmp += m_name.substr(1, 1);
         } else if (is_field()) {
           tmp += "i";
         }