Merge pull request #10 from TimGoTheCreator/SoA

Update to SoA

Author: TimGoTheCreator

CMakeLists.txt

@@ -68,7 +68,16 @@ add_executable(next ${SRC_FILES} ${ARGPARSE_FILES})
 find_package(OpenMP QUIET)
 if(OpenMP_CXX_FOUND)
     message(STATUS "OpenMP detected — enabling multithreading.")
-    target_link_libraries(next PRIVATE OpenMP::OpenMP_CXX)
+    
+    if(MSVC)
+        # Force MSVC to use the modern LLVM backend (supports OpenMP 3.0+ and size_t)
+        # We add it as a compiler option because MSVC's default find_package 
+        # often defaults to the legacy /openmp flag.
+        target_compile_options(next PRIVATE /openmp:llvm)
+    else()
+        target_link_libraries(next PRIVATE OpenMP::OpenMP_CXX)
+    endif()
+    
 else()
     message(STATUS "OpenMP not found — building in single-threaded mode.")
 endif()

docs/install_linux.md

@@ -40,7 +40,7 @@ cd ..
 ```bash
 cd examples/TwoBodies
 python two_body.py
-../../next two_body.txt 8 0.01 0.1 vtu
+../../next two_body.txt 8 0.001 0.1 vtu
 ```
 
 5. **View results**

docs/install_mac.md

@@ -40,7 +40,7 @@ cd ..
 ```bash
 cd examples/TwoBodies
 python two_body.py
-../../next two_body.txt 8 0.01 0.1 vtu
+../../next two_body.txt 8 0.001 0.1 vtu
 ```
 
 5. **View results**

docs/install_mingw.md

@@ -42,9 +42,14 @@ cd ..
 ```bash
 cd examples/TwoBodies
 python two_body.py
-../../next.exe two_body.txt 8 0.01 0.1 vtu
 ```
 
+#### If the copy command worked (copies to project root by default in CMake):
+../../next.exe two_body.txt 8 0.001 0.1 vtu
+
+#### If your executable is in build/
+../../build/next.exe two_body.txt 8 0.001 0.1 vtu
+
 5. **View results**
 
 Open the `.vtu` output in ParaView.

docs/install_msvc.md

@@ -8,6 +8,7 @@
 - vcpkg
 - HDF5 (installed via vcpkg)
 - Python 3
+- LLVM OpenMP Component (C++ Clang-cl for v143 build tools (x64/x86)
 
 ## Steps
 
@@ -41,9 +42,14 @@ cd ..
 ```powershell
 cd examples/TwoBodies
 python two_body.py
-..\..\next.exe two_body.txt 8 0.01 0.1 vtu
 ```
 
+#### If the copy command worked (copies to project root by default in CMake):
+../../next.exe two_body.txt 8 0.001 0.1 vtu
+
+#### If your executable is in build/
+../../build/Release/next.exe two_body.txt 8 0.001 0.1 vtu
+
 5. **View results**
 
 Open the `.vtu` output in ParaView.

src/begrun.cpp

@@ -28,76 +28,80 @@ using next::OutputFormat;
 int main(int argc, char **argv) {
   auto args = next::parse_arguments(argc, argv);
 
-  // ASCII banner (raw string literal preserves backslashes and spacing)
   static constexpr const char *BANNER = R"NEXTBANNER(
-_   _ ________   _________ 
+ _  _ ________   _________ 
 | \ | |  ____\ \ / /__   __|
 |  \| | |__   \ V /   | |   
 | . ` |  __|   > <    | |   
 | |\  | |____ / . \   | |   
-|_| \_|______/_/ \_\  |_|  
+|_| \_|______/_/ \_\  |_|   
 )NEXTBANNER";
 
-  // Print banner once at startup
   std::cout << BANNER << '\n';
 
   // Set threads and report
   omp_set_num_threads(args.threads);
-
-  std::cout << " Threads: " << args.threads << "\n";
+  std::cout << " Threads:   " << args.threads << "\n";
 
 #ifdef NEXT_FP64
   std::cout << " Precision: FP64\n";
 #elif defined(NEXT_FP32)
   std::cout << " Precision: FP32\n";
 #endif
 
-  // Load particles
-  std::vector<Particle> particles = LoadParticlesFromFile(args.input_file);
+  // --- SoA UPDATE ---
+  // LoadParticlesFromFile now returns a single 'Particle' struct
+  // which internally contains all the parallel vectors.
+  Particle particles = LoadParticlesFromFile(args.input_file);
+
+  std::cout << " Particles: " << particles.size() << "\n";
 
   real simTime = 0;
   real nextDump = 0;
   int step = 0;
   char command;
 
   while (true) {
+    // Both computeAdaptiveDt and Step now take the SoA object by reference
     real dtAdaptive = computeAdaptiveDt(particles, args.dt);
     Step(particles, dtAdaptive);
+    
     simTime += dtAdaptive;
 
-        if (simTime >= nextDump) {
-            std::string out = "dump_" + std::to_string(step);
+    if (simTime >= nextDump) {
+      std::string out = "dump_" + std::to_string(step);
 
-            switch (args.format) {
-                case OutputFormat::VTK:
-                    out += ".vtk";
-                    SaveVTK(particles, out);
-                    break;
+      switch (args.format) {
+        case OutputFormat::VTK:
+          out += ".vtk";
+          SaveVTK(particles, out);
+          break;
 
-                case OutputFormat::VTU:
-                    out += ".vtu";
-                    SaveVTU(particles, out);
-                    break;
+        case OutputFormat::VTU:
+          out += ".vtu";
+          SaveVTU(particles, out);
+          break;
 
-                case OutputFormat::HDF5:
-                    out += ".hdf5";
-                    SaveHDF5(particles, out);
-                    break;
-            }
+        case OutputFormat::HDF5:
+          out += ".hdf5";
+          SaveHDF5(particles, out);
+          break;
+      }
 
-            std::cout << "[Dump " << step << "] t = " << simTime
-                      << ", file: " << out << "\n";
-
-            nextDump += args.dump_interval;
-            step++;
-        }
+      std::cout << "[Dump " << step << "] t = " << simTime
+                << ", file: " << out << "\n";
 
+      nextDump += args.dump_interval;
+      step++;
+    }
 
+    // Non-blocking exit check
     if (std::cin.rdbuf()->in_avail() > 0) {
       std::cin >> command;
-      if (command == 'q' || command == 'Q')
+      if (command == 'q' || command == 'Q') {
         std::cout << "Exiting...\n";
-      break;
+        break;
+      }
     }
   }
 

src/dt/adaptive.h

@@ -16,27 +16,37 @@
 #include <cmath>
 #include <vector>
 
-real computeAdaptiveDt(const std::vector<Particle> &p, real base_dt) {
-  real maxSpeed = 0;
-
-  for (const auto &a : p) {
-    real speed = std::sqrt(a.vx * a.vx + a.vy * a.vy + a.vz * a.vz);
-    if (speed > maxSpeed)
-      maxSpeed = speed;
-  }
-
-  maxSpeed = std::min(maxSpeed, real(1e4));
-
-  // If everything is basically stationary, use base dt
-  if (maxSpeed < real(1e-8))
-    return base_dt;
-
-  // Smaller dt when speeds are high
-  real dt = base_dt / (1 + maxSpeed);
-
-  // Clamp dt to a reasonable range
-  dt = std::max(dt, base_dt * real(0.01)); // never smaller than 1% of base
-  dt = std::min(dt, base_dt * real(1.0));  // never larger than base
-
-  return dt;
+/**
+ * @brief Computes a global adaptive time-step based on the maximum velocity in the system.
+ * Updated for SoA (Structure of Arrays) for better cache performance.
+ */
+real computeAdaptiveDt(const Particle &p, real base_dt) {
+    real maxSpeedSq = 0;
+    const size_t N = p.size();
+
+    // High-speed linear scan through velocity arrays
+    // The compiler can easily vectorize this with SIMD (AVX/SSE)
+    for (size_t i = 0; i < N; ++i) {
+        real speedSq = p.vx[i] * p.vx[i] + p.vy[i] * p.vy[i] + p.vz[i] * p.vz[i];
+        if (speedSq > maxSpeedSq)
+            maxSpeedSq = speedSq;
+    }
+
+    real maxSpeed = std::sqrt(maxSpeedSq);
+    
+    // Safety clamp to prevent dt from exploding or becoming zero
+    maxSpeed = std::min(maxSpeed, real(1e4));
+
+    // If everything is basically stationary, use base dt
+    if (maxSpeed < real(1e-8))
+        return base_dt;
+
+    // Standard Courant-like condition: smaller dt when speeds are high
+    real dt = base_dt / (1 + maxSpeed);
+
+    // Clamp dt to a reasonable range
+    dt = std::max(dt, base_dt * real(0.01)); // never smaller than 1% of base
+    dt = std::min(dt, base_dt * real(1.0));  // never larger than base
+
+    return dt;
 }

src/dt/softening.h

@@ -10,43 +10,43 @@
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #pragma once
+#include "floatdef.h"
 #include <cmath>
+#include <algorithm>
 
-/* Softening for BarnesHut */
+/**
+ * @brief Softening for Barnes-Hut node interactions.
+ * Note: nodeMass and dist are passed as individual reals from the SoA arrays.
+ */
 inline real nextSoftening(real nodeSize, real nodeMass, real dist) {
-  // Base softening from node size
-  real eps_size = nodeSize * real(0.015);
+    // std::pow(x, 1/3) or std::cbrt is slow. 
+    // In SoA loops, this is often the bottleneck.
+    real eps_size = nodeSize * real(0.015);
+    real eps_mass = std::cbrt(nodeMass) * real(0.002);
 
-  // Mass-based softening (physical radius proxy)
-  real eps_mass = std::cbrt(nodeMass) * real(0.002);
+    // Distance taper: strong at r->0, fades smoothly
+    real eps_taper = real(1.0) / (real(1.0) + dist * real(10.0));
 
-  // Distance taper: strong at r→0, fades smoothly
-  real eps_taper = real(1.0) / (real(1.0) + dist * real(10.0));
+    // Combine
+    real eps = (eps_size + eps_mass) * eps_taper;
 
-  // Combine
-  real eps = (eps_size + eps_mass) * eps_taper;
-
-  // Minimum floor
-  const real eps_min = real(1e-4);
-  if (eps < eps_min)
-    eps = eps_min;
-
-  return eps;
+    // Minimum floor - using std::max is cleaner and easier for the compiler to optimize
+    return std::max(eps, real(1e-4));
 }
 
-/* Softening for Gravity kernel */
+/**
+ * @brief Softening for direct particle-particle gravity kernels.
+ * ma and mb are the masses pulled from the ps.m[i] and ps.m[j] arrays.
+ */
 inline real pairSoftening(real ma, real mb) {
-  // Physical radius proxy
-  real ea = std::cbrt(ma) * real(0.002);
-  real eb = std::cbrt(mb) * real(0.002);
-
-  // Symmetric combination (quadrature)
-  real eps = std::sqrt(ea * ea + eb * eb);
+    // Physical radius proxy
+    real ea = std::cbrt(ma) * real(0.002);
+    real eb = std::cbrt(mb) * real(0.002);
 
-  // Minimum floor
-  const real eps_min = real(1e-4);
-  if (eps < eps_min)
-    eps = eps_min;
+    // Symmetric combination (quadrature)
+    real epsSq = ea * ea + eb * eb;
+    real eps = std::sqrt(epsSq);
 
-  return eps; // return epsilon
+    // Minimum floor
+    return std::max(eps, real(1e-4));
 }