diff --git a/.github/workflows/cmake-multi-platform.yml b/.github/workflows/cmake-multi-platform.yml index 7ab1b4d..dffe897 100644 --- a/.github/workflows/cmake-multi-platform.yml +++ b/.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" ] diff --git a/CMakeLists.txt b/CMakeLists.txt index 0b34d78..4a86f10 100644 --- a/CMakeLists.txt +++ b/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}) + diff --git a/src/begrun.cpp b/src/begrun.cpp index bcc006e..baf97ec 100644 --- a/src/begrun.cpp +++ b/src/begrun.cpp @@ -2,6 +2,7 @@ #include #include #include"floatdef.h" +#include"dt/adaptive.h" std::vector 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 \n"; + std::cerr << "Usage: next
\n"; return 1; } const char* filename = argv[1]; + real dt = std::stod(argv[2]); std::vector particles = LoadParticlesFromFile(filename); - real dt = 0.01; while (true) { - Step(particles, dt); + real dtAdaptive = computeAdaptiveDt(particles, dt); + Step(particles, dtAdaptive); } return 0; -} \ No newline at end of file +} + diff --git a/src/dt/adaptive.h b/src/dt/adaptive.h new file mode 100644 index 0000000..5af823c --- /dev/null +++ b/src/dt/adaptive.h @@ -0,0 +1,31 @@ +#pragma once +#include "struct/particle.h" +#include +#include +#include +#include "floatdef.h" + +real computeAdaptiveDt(const std::vector& 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; +}