Try fix.

GPGPUCourse · Dec 19, 2024 · 0d27767 · 0d27767
1 parent 3ced50c
commit 0d27767
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 29 deletions.
diff --git a/src/cl/lbvh.cl b/src/cl/lbvh.cl
@@ -464,15 +464,11 @@ void calculateForce(float x0, float y0, float m0, __global const struct Node *no
             if (!contains(&child->bbox, x0, y0) && barnesHutCondition(x0, y0, child)) {
                 float dx = child->cmsx - x0;
                 float dy = child->cmsy - y0;
-                float dr2 = max(100.f, dx*dx + dy*dy);
-                float dr2_inv = 1.f / dr2;
-                float dr_inv = sqrt(dr2_inv);
-                float ex = dx * dr_inv;
-                float ey = dy * dr_inv;
-                float fx = ex * dr2_inv * GRAVITATIONAL_FORCE;
-                float fy = ey * dr2_inv * GRAVITATIONAL_FORCE;
-                *force_x += child->mass * fx;
-                *force_y += child->mass * fy;
+                float dr = 1.f / max(100.f, dx*dx + dy*dy);
+                float ex = dx * sqrt(dr);
+                float ey = dy * sqrt(dr);
+                *force_x += child->mass * ex * dr * GRAVITATIONAL_FORCE;
+                *force_y += child->mass * ey * dr * GRAVITATIONAL_FORCE;
             } else {
                 stack[stack_size++] = i_child;
             }

diff --git a/src/main_lbvh.cpp b/src/main_lbvh.cpp
@@ -375,17 +375,6 @@ void calculateForce(float x0, float y0, float m0, const std::vector<Node> &nodes
         const Node &node = nodes[i_node];
 
         if (node.isLeaf()) {
-            // float dx = node.cmsx - x0;
-            // float dy = node.cmsy - y0;
-            // float dr2 = std::max(100.f, dx*dx + dy*dy);
-            // float dr2_inv = 1.f / dr2;
-            // float dr_inv = std::sqrt(dr2_inv);
-            // float ex = dx * dr_inv;
-            // float ey = dy * dr_inv;
-            // float fx = ex * dr2_inv * GRAVITATIONAL_FORCE;
-            // float fy = ey * dr2_inv * GRAVITATIONAL_FORCE;
-            // *force_x += node.mass * fx;
-            // *force_y += node.mass * fy;
             continue;
         }
 
@@ -412,15 +401,11 @@ void calculateForce(float x0, float y0, float m0, const std::vector<Node> &nodes
             if (!child.bbox.contains(x0, y0) && barnesHutCondition(x0, y0, child)) {
                 float dx = child.cmsx - x0;
                 float dy = child.cmsy - y0;
-                float dr2 = std::max(100.f, dx*dx + dy*dy);
-                float dr2_inv = 1.f / dr2;
-                float dr_inv = std::sqrt(dr2_inv);
-                float ex = dx * dr_inv;
-                float ey = dy * dr_inv;
-                float fx = ex * dr2_inv * GRAVITATIONAL_FORCE;
-                float fy = ey * dr2_inv * GRAVITATIONAL_FORCE;
-                *force_x += child.mass * fx;
-                *force_y += child.mass * fy;
+                float dr = 1.f / std::max(100.f, dx*dx + dy*dy);
+                float ex = dx * std::sqrt(dr);
+                float ey = dy * std::sqrt(dr);
+                *force_x += child.mass * ex * dr * GRAVITATIONAL_FORCE;
+                *force_y += child.mass * ey * dr * GRAVITATIONAL_FORCE;
             } else {
                 stack[stack_size] = i_child;
                 stack_size++;