From 428b1d08dee85564dbabed5f8459b29e37d1a797 Mon Sep 17 00:00:00 2001
From: Ryker Fish <rfish@mines.edu>
Date: Fri, 24 Jan 2025 13:33:11 -0700
Subject: [PATCH] fixed bug in damping code that made it not symmetric

---
 solvers/NBody/extra/hydrodynamicKernels.cuh | 16 ++++++++--------
 tests/utils.py                              |  2 --
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/solvers/NBody/extra/hydrodynamicKernels.cuh b/solvers/NBody/extra/hydrodynamicKernels.cuh
index 04c65f6..71f93fb 100644
--- a/solvers/NBody/extra/hydrodynamicKernels.cuh
+++ b/solvers/NBody/extra/hydrodynamicKernels.cuh
@@ -317,21 +317,21 @@ __device__ real3 RPY_WF(real3 rij,real r, real rh){
       // implements damping from Appendix 1 in [2] so the matrix is positive definite when a particle overlaps the wall
       real bi = min(pi.z/rh, real(1.0));
       bi = max(bi, real(0.0));
-      if(bi == real(0.0)){
-        return result;
-      }
-      real bi2 = bi*bi;
+      real bj = min(pj.z/rh, real(1.0));
+      bj = max(bj, real(0.0));
+      real bij = bi*bj;
+
       pi.z = max(pi.z, rh);
       pj.z = max(pj.z, rh);
 
       real3 rij = make_real3(pi)-make_real3(pj);
       const real r = sqrt(dot(rij, rij));
 
-      result.MF += bi2*dotProduct_UF(rij, r, fj, pj.z);
+      result.MF += bij*dotProduct_UF(rij, r, fj, pj.z);
       if (haveTorque){
-        result.MF += bi2*dotProduct_UT(rij, r, tj, pi.z); // this is correct with pi.z: see note on dotProduct_UT
-        result.MT += bi2*dotProduct_WF(rij, r, fj, pj.z);
-        result.MT += bi2*dotProduct_WT(rij, r, tj, pj.z);
+        result.MF += bij*dotProduct_UT(rij, r, tj, pi.z); // this is correct with pi.z: see note on dotProduct_UT
+        result.MT += bij*dotProduct_WF(rij, r, fj, pj.z);
+        result.MT += bij*dotProduct_WT(rij, r, tj, pj.z);
       }
 
       return result;
diff --git a/tests/utils.py b/tests/utils.py
index 3665f9d..74a07cb 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -91,7 +91,5 @@ def generate_positions_in_box(parameters, numberParticles):
     if "algorithm" in parameters:
         positions[:, 2] += 0.5  # [-0.5, 0.5] -> [0, 1]
         positions *= 10  # [0, 1] -> [0, 10]
-        # Move particles to at least one hydrodynamic radius away from the bottom wall
-        # positions[:, 2] += 1  # [0, 10] -> [1, 11]
 
     return positions