Add Adaptive dt

.github/workflows/cmake-multi-platform.yml

@@ -4,7 +4,9 @@ name: CMake on multiple platforms
 
 on:
   push:
-    branches: [ "main" ]
+    branches:
+    - main
+    - adaptiveDt
   pull_request:
     branches: [ "main" ]
 

CMakeLists.txt

@@ -68,7 +68,9 @@ endif()
 
 # Includes and sources
 include_directories(${CMAKE_SOURCE_DIR}/include)
+include_directories(${CMAKE_SOURCE_DIR}/src)
 file(GLOB SRC_FILES ${CMAKE_SOURCE_DIR}/src/*.cpp)
 
 # Build executable
 add_executable(next ${SRC_FILES})
+

src/begrun.cpp

@@ -2,6 +2,7 @@
 #include<iostream>
 #include<fstream>
 #include"floatdef.h"
+#include"dt/adaptive.h"
 
 std::vector<Particle> LoadParticlesFromFile(const std::string& filename)
 {
@@ -26,23 +27,25 @@ int main(int argc, char** argv) {
     << "NN   NN  EEEEEEE  XX   XX    TTT  \n"
     << "Newtonian EXact Trajectories\n";
 
-    if (argc != 2)
+    if (argc != 3)
     {
-        std::cerr << "Usage: next <initial.txt>\n";
+        std::cerr << "Usage: next <initial.txt> <dt>\n";
         return 1;
     }
 
     const char* filename = argv[1];
+    real dt = std::stod(argv[2]);
 
     std::vector<Particle> particles = LoadParticlesFromFile(filename);
 
-    real dt = 0.01;
     while (true)
     {
-        Step(particles, dt);
+        real dtAdaptive = computeAdaptiveDt(particles, dt);
+        Step(particles, dtAdaptive);
     }
 
     return 0;
 
 
-}
+}
+

src/dt/adaptive.h

@@ -0,0 +1,31 @@
+#pragma once
+#include "struct/particle.h"
+#include <vector>
+#include <algorithm>
+#include <cmath>
+#include "floatdef.h"
+
+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;
+    }
+
+    // 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;
+}