Merge remote-tracking branch 'upstream/dev' into picmi_mw

1 files changed, 39 insertions, 0 deletions

diff --git a/Source/Particles/Collision/ComputeTemperature.H b/Source/Particles/Collision/ComputeTemperature.H
new file mode 100644
index 000000000..3cc96fb52
--- /dev/null
+++ b/Source/Particles/Collision/ComputeTemperature.H
@@ -0,0 +1,39 @@
+#ifndef WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
+#define WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_
+
+#include <WarpXConst.H>
+
+template <typename T_index, typename T_R>
+T_R ComputeTemperature (
+    T_index const Is, T_index const Ie, T_index const *I,
+    T_R const *ux, T_R const *uy, T_R const *uz, T_R const m )
+{
+
+    T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c );
+
+    int N = Ie - Is;
+    if ( N == 0 ) { return T_R(0.0); }
+
+    T_R vx = T_R(0.0);    T_R vy = T_R(0.0);
+    T_R vz = T_R(0.0);    T_R vs = T_R(0.0);
+    T_R gm = T_R(0.0);    T_R us = T_R(0.0);
+
+    for (int i = Is; i < Ie; ++i)
+    {
+        us = ( ux[ I[i] ] * ux[ I[i] ] +
+               uy[ I[i] ] * uy[ I[i] ] +
+               uz[ I[i] ] * uz[ I[i] ] );
+        gm = std::sqrt( T_R(1.0) + us*inv_c2 );
+        vx += ux[ I[i] ] / gm;
+        vy += uy[ I[i] ] / gm;
+        vz += uz[ I[i] ] / gm;
+        vs += us / gm / gm;
+    }
+
+    vx = vx / N;    vy = vy / N;
+    vz = vz / N;    vs = vs / N;
+
+    return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz));
+}
+
+#endif // WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_