fixed ComputeTemperature routine called from Coulomb collisions to wo…

…rk for particles of varying weight.
JustinRayAngus · Jul 10, 2024 · 26ed88e · 26ed88e
1 parent ce9f4af
commit 26ed88e
Show file tree

Hide file tree

Showing 2 changed files with 11 additions and 8 deletions.
diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/ComputeTemperature.H b/Source/Particles/Collision/BinaryCollision/Coulomb/ComputeTemperature.H
@@ -14,6 +14,7 @@ template <typename T_index, typename T_R>
 AMREX_GPU_HOST_DEVICE
 T_R ComputeTemperature (
     T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I,
+    T_R const * AMREX_RESTRICT w,
     T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT uz,
     T_R const m )
 {
@@ -26,21 +27,23 @@ T_R ComputeTemperature (
     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);
+    T_R wtot = T_R(0.0);
 
     for (int i = Is; i < static_cast<int>(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;
+        wtot += w[ I[i] ];
+        vx += w[ I[i] ] * ux[ I[i] ] / gm;
+        vy += w[ I[i] ] * uy[ I[i] ] / gm;
+        vz += w[ I[i] ] * uz[ I[i] ] / gm;
+        vs += w[ I[i] ] * us / gm / gm;
     }
 
-    vx = vx / N;    vy = vy / N;
-    vz = vz / N;    vs = vs / N;
+    vx = vx / wtot;    vy = vy / wtot;
+    vz = vz / wtot;    vs = vs / wtot;
 
     return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz));
 }

diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H b/Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H
@@ -74,12 +74,12 @@ void ElasticCollisionPerez (
     T_PR T1t; T_PR T2t;
     if ( T1 <= T_PR(0.0) && L <= T_PR(0.0) )
     {
-        T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1);
+        T1t = ComputeTemperature(I1s,I1e,I1,w1,u1x,u1y,u1z,m1);
     }
     else { T1t = T1; }
     if ( T2 <= T_PR(0.0) && L <= T_PR(0.0) )
     {
-        T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2);
+        T2t = ComputeTemperature(I2s,I2e,I2,w2,u2x,u2y,u2z,m2);
     }
     else { T2t = T2; }