diff --git a/CMakeLists.txt b/CMakeLists.txt index 10c5f8d..2d2650e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,36 +1,89 @@ cmake_minimum_required(VERSION 3.28) + project(CMakeSFMLProject LANGUAGES CXX) -set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/bin) +# ============================================================ +# Project settings +# ============================================================ set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY + ${CMAKE_SOURCE_DIR}/bin +) + +# ============================================================ +# OpenMP +# ============================================================ + +find_package(OpenMP REQUIRED) -add_executable(sim src/main.cpp) +# ============================================================ +# SFML +# ============================================================ -# 1. Fetch SFML normally (Without OpenMP polluting it) include(FetchContent) -FetchContent_Declare(SFML + +FetchContent_Declare( + SFML GIT_REPOSITORY https://github.com/SFML/SFML.git GIT_TAG 3.1.0 GIT_SHALLOW ON EXCLUDE_FROM_ALL SYSTEM ) + FetchContent_MakeAvailable(SFML) -# 2. Link SFML to your executable -target_compile_features(sim PRIVATE cxx_std_20) -target_link_libraries(sim PRIVATE sfml-graphics sfml-window sfml-system) +# ============================================================ +# Executable +# ============================================================ + +add_executable( + sim + src/main.cpp +) + +# ============================================================ +# Compiler features +# ============================================================ + +target_compile_features( + sim + PRIVATE + cxx_std_20 +) -# 3. Target-Specific OpenMP for Mac (Applies ONLY to 'sim') -if(APPLE) - target_include_directories(sim PRIVATE /opt/homebrew/opt/libomp/include) - target_link_directories(sim PRIVATE /opt/homebrew/opt/libomp/lib) - target_compile_options(sim PRIVATE -Xpreprocessor -fopenmp) - target_link_options(sim PRIVATE -Xpreprocessor -fopenmp) # Ensure linker knows about threads - target_link_libraries(sim PRIVATE omp) -endif() +# ============================================================ +# Compile options +# ============================================================ + +target_compile_options( + sim + PRIVATE + $<$:-O3> +) -set_property(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY VS_STARTUP_PROJECT main) +# ============================================================ +# Libraries +# ============================================================ -set(CMAKE_CXX_FLAGS_RELEASE "-O3 -Xpreprocessor -fopenmp") +target_link_libraries( + sim + PRIVATE + sfml-graphics + sfml-window + sfml-system + OpenMP::OpenMP_CXX +) + +# ============================================================ +# IDE settings +# ============================================================ + +set_property( + DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + PROPERTY VS_STARTUP_PROJECT sim +) diff --git a/src/main.cpp b/src/main.cpp index 9f309bf..85877c1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,8 +1,14 @@ #include -#include -#include + +#include +#include +#include +#include #include -#include +#include +#include +#include +#include #include @@ -10,141 +16,604 @@ #include "barnes_hut_tree.hpp" #include "diagnostics.hpp" -int main() { - sf::Color p_color = sf::Color::Cyan; - const int screen_size = 1380; - sf::RenderWindow window(sf::VideoMode({screen_size, screen_size}) , "Simulation"); - const int frame_rate = 60; - window.setFramerateLimit(frame_rate); - - const double dt = 0.001; - - sf::Clock fpstimer; - const sf::Font font("/Library/Fonts/JetBrainsMonoNL-Regular.ttf"); - - std::mt19937 gen(42); - std::normal_distribution cluster(1.0 * screen_size / 2, 1.0 * screen_size / 10); - std::uniform_real_distribution uniform(100.0, screen_size * 1.0 - 100.0); - - std::vector sources; - Complex center{1.0 * screen_size / 2, 1.0 * screen_size / 2}; - for (size_t i = 0; i < 40000; i++) { - double x = cluster(gen), y = cluster(gen); - while (x < 100 or x > screen_size - 100) x = cluster(gen); - while (y < 100 or y > screen_size - 100) y = cluster(gen); - sources.emplace_back(x, y, 10.0); - } - for (size_t i = 0; i < 10000; i++) - sources.emplace_back(uniform(gen), uniform(gen), 1.0); - - const double orbital_speed = 200.0; - for (auto &s : sources) { - Complex disp = s.position - center; - double r = std::abs(disp); - if (r > 1.0) { - Complex tangent{-disp.imag(), disp.real()}; - tangent /= r; - s.velocity = tangent * orbital_speed; - } - } - - fmm::FmmTree tree(sources, 60, 10); - // fmm::BhTree tree(sources, 1, 2.5); - tree.buildTree(); - - bool drawBoxes = false; - - int frame_limit = 1800, tot_frames = 0; - float tot_time = 0, max_time = 0, min_time = 100; - float times[frame_limit]; - - double error_l2[frame_limit]; - double error_abs[frame_limit]; - double avg_error_l2 = 0, avg_error_abs = 0, max_error = 0; - - while (window.isOpen()) { - while (const std::optional event = window.pollEvent()) { - if (event->is()) window.close(); - - else if (const auto* keyPressed = event->getIf()) { - if (keyPressed->scancode == sf::Keyboard::Scancode::Escape) window.close(); - } - } - fpstimer.restart(); - window.clear(sf::Color::Black); - - #pragma omp parallel for - for (size_t i = 0; i < sources.size(); i++) { - sources[i].velocity += 0.5 * tree.forces[i] * dt; - sources[i].position += sources[i].velocity * dt; - } - - // forces = fmm::computeExactForces(sources); - tree.buildTree(); - // fmm::ErrorData error = fmm::evaluateSimulationError(sources, tree.forces, 1000); - - if (drawBoxes) { - auto boxes = tree.getBoxGeometries(); - sf::VertexArray box_va(sf::PrimitiveType::Lines, boxes.size() * 8); - sf::Color box_color(66, 191, 245); - - #pragma omp parallel for schedule(static) - for (size_t i = 0; i < (int)boxes.size(); i++) { - int id_b = i * 8; - - float cx = static_cast(boxes[i].first.real()); - float cy = static_cast(boxes[i].first.imag()); - float len = static_cast(boxes[i].second / 2.0); - - sf::Vector2f tl(cx - len, cy - len); - sf::Vector2f tr(cx + len, cy - len); - sf::Vector2f br(cx + len, cy + len); - sf::Vector2f bl(cx - len, cy + len); - - box_va[id_b ] = sf::Vertex{tl, box_color}; - box_va[id_b + 1] = sf::Vertex{tr, box_color}; - - box_va[id_b + 2] = sf::Vertex{tr, box_color}; - box_va[id_b + 3] = sf::Vertex{br, box_color}; - - box_va[id_b + 4] = sf::Vertex{br, box_color}; - box_va[id_b + 5] = sf::Vertex{bl, box_color}; - - box_va[id_b + 6] = sf::Vertex{bl, box_color}; - box_va[id_b + 7] = sf::Vertex{tl, box_color}; - } - window.draw(box_va); - } - - #pragma omp parallel for - for (size_t i = 0; i < sources.size(); i++) - sources[i].velocity += 0.5 * tree.forces[i] * dt; - - sf::VertexArray particle_va(sf::PrimitiveType::Points, sources.size()); - - #pragma omp parallel for schedule(static) - for (int i = 0; i < (int)sources.size(); ++i) { - float cx = static_cast(sources[i].position.real()); - float cy = static_cast(sources[i].position.imag()); - - particle_va[i] = sf::Vertex{sf::Vector2f(cx, cy), p_color}; - } - - window.draw(particle_va); - - float ms = 1.0 * fpstimer.getElapsedTime().asMicroseconds() / 1000; - sf::Text number(font, std::to_string(ms) + "ms", 24); - number.setFillColor(sf::Color::White); - window.draw(number); - window.display(); - - // min_time = std::min(min_time, ms); - // max_time = std::max(max_time, ms); - // tot_time += ms; - // times[tot_frames] = ms; - - tot_frames++; - // if (tot_frames == frame_limit) break; - } - return 0; -} \ No newline at end of file +int main() +{ + // ============================================================ + // Window settings + // ============================================================ + + constexpr unsigned int screen_size = 1380; + constexpr unsigned int frame_rate = 60; + + sf::RenderWindow window( + sf::VideoMode({screen_size, screen_size}), + "Simulation" + ); + + window.setFramerateLimit(frame_rate); + + // ============================================================ + // Simulation settings + // ============================================================ + + constexpr double dt = 0.001; + constexpr double orbital_speed = 200.0; + + constexpr std::size_t cluster_particle_count = 40000; + constexpr std::size_t uniform_particle_count = 10000; + + // ============================================================ + // Font + // ============================================================ + + const std::string font_path = + "/usr/share/fonts/truetype/jetbrains-mono/JetBrainsMono-Regular.ttf"; + + if (!std::filesystem::exists(font_path)) + { + std::cerr + << "ERROR: Font not found:\n" + << " " << font_path << "\n\n" + << "Install JetBrains Mono with:\n" + << " sudo apt install fonts-jetbrains-mono\n"; + + return EXIT_FAILURE; + } + + sf::Font font(font_path); + + // ============================================================ + // Random number generation + // ============================================================ + + std::mt19937 gen(42); + + const double center_coordinate = + static_cast(screen_size) / 2.0; + + const double cluster_sigma = + static_cast(screen_size) / 10.0; + + std::normal_distribution cluster_distribution( + center_coordinate, + cluster_sigma + ); + + std::uniform_real_distribution uniform_distribution( + 100.0, + static_cast(screen_size) - 100.0 + ); + + // ============================================================ + // Create particles + // ============================================================ + + std::vector sources; + + sources.reserve( + cluster_particle_count + + uniform_particle_count + ); + + Complex center{ + center_coordinate, + center_coordinate + }; + + // ------------------------------------------------------------ + // Dense central cluster + // ------------------------------------------------------------ + + for (std::size_t i = 0; + i < cluster_particle_count; + ++i) + { + double x = cluster_distribution(gen); + double y = cluster_distribution(gen); + + while ( + x < 100.0 || + x > static_cast(screen_size) - 100.0 + ) + { + x = cluster_distribution(gen); + } + + while ( + y < 100.0 || + y > static_cast(screen_size) - 100.0 + ) + { + y = cluster_distribution(gen); + } + + sources.emplace_back( + x, + y, + 10.0 + ); + } + + // ------------------------------------------------------------ + // Uniform particles + // ------------------------------------------------------------ + + for (std::size_t i = 0; + i < uniform_particle_count; + ++i) + { + sources.emplace_back( + uniform_distribution(gen), + uniform_distribution(gen), + 1.0 + ); + } + + std::cout + << "Created " + << sources.size() + << " particles.\n"; + + // ============================================================ + // Initial orbital velocities + // ============================================================ + + for (auto& source : sources) + { + Complex displacement = + source.position - center; + + const double radius = + std::abs(displacement); + + if (radius > 1.0) + { + Complex tangent{ + -displacement.imag(), + displacement.real() + }; + + tangent /= radius; + + source.velocity = + tangent * orbital_speed; + } + } + + // ============================================================ + // FMM tree + // ============================================================ + + fmm::FmmTree tree( + sources, + 60, + 10 + ); + + // Alternative Barnes-Hut implementation: + // + // fmm::BhTree tree( + // sources, + // 1, + // 2.5 + // ); + + std::cout << "Building initial FMM tree...\n"; + + tree.buildTree(); + + std::cout << "FMM tree ready.\n"; + + // ============================================================ + // Rendering settings + // ============================================================ + + constexpr bool draw_boxes = false; + + const sf::Color particle_color = + sf::Color::Cyan; + + const sf::Color box_color( + 66, + 191, + 245 + ); + + // ============================================================ + // Diagnostics + // ============================================================ + + constexpr int frame_limit = 1800; + + int total_frames = 0; + + float total_time = 0.0f; + float max_time = 0.0f; + float min_time = 1000000.0f; + + std::vector times( + frame_limit + ); + + std::vector error_l2( + frame_limit + ); + + std::vector error_abs( + frame_limit + ); + + double average_error_l2 = 0.0; + double average_error_abs = 0.0; + double max_error = 0.0; + + // Prevent unused-variable warnings until diagnostics are enabled. + (void)error_l2; + (void)error_abs; + (void)average_error_l2; + (void)average_error_abs; + (void)max_error; + + // ============================================================ + // Frame timer + // ============================================================ + + sf::Clock fpstimer; + + // ============================================================ + // Main simulation loop + // ============================================================ + + while (window.isOpen()) + { + // -------------------------------------------------------- + // Events + // -------------------------------------------------------- + + while (const std::optional event = + window.pollEvent()) + { + if (event->is()) + { + window.close(); + } + else if ( + const auto* key_pressed = + event->getIf() + ) + { + if ( + key_pressed->scancode == + sf::Keyboard::Scancode::Escape + ) + { + window.close(); + } + } + } + + // -------------------------------------------------------- + // Start frame timer + // -------------------------------------------------------- + + fpstimer.restart(); + + // -------------------------------------------------------- + // Clear window + // -------------------------------------------------------- + + window.clear( + sf::Color::Black + ); + + // ======================================================== + // First half of velocity Verlet integration + // ======================================================== + + #pragma omp parallel for schedule(static) + for (int i = 0; + i < static_cast(sources.size()); + ++i) + { + sources[i].velocity += + 0.5 * + tree.forces[i] * + dt; + + sources[i].position += + sources[i].velocity * + dt; + } + + // ======================================================== + // Rebuild FMM tree + // ======================================================== + + tree.buildTree(); + + // ======================================================== + // Optional exact-force diagnostics + // ======================================================== + + /* + auto exact_forces = + fmm::computeExactForces(sources); + + fmm::ErrorData error = + fmm::evaluateSimulationError( + sources, + tree.forces, + 1000 + ); + */ + + // ======================================================== + // Draw FMM boxes + // ======================================================== + + if (draw_boxes) + { + const auto boxes = + tree.getBoxGeometries(); + + sf::VertexArray box_vertices( + sf::PrimitiveType::Lines, + boxes.size() * 8 + ); + + #pragma omp parallel for schedule(static) + for (int i = 0; + i < static_cast(boxes.size()); + ++i) + { + const std::size_t id = + static_cast(i) * 8; + + const float cx = + static_cast( + boxes[i].first.real() + ); + + const float cy = + static_cast( + boxes[i].first.imag() + ); + + const float half_size = + static_cast( + boxes[i].second / 2.0 + ); + + const sf::Vector2f top_left{ + cx - half_size, + cy - half_size + }; + + const sf::Vector2f top_right{ + cx + half_size, + cy - half_size + }; + + const sf::Vector2f bottom_right{ + cx + half_size, + cy + half_size + }; + + const sf::Vector2f bottom_left{ + cx - half_size, + cy + half_size + }; + + box_vertices[id] = + sf::Vertex{ + top_left, + box_color + }; + + box_vertices[id + 1] = + sf::Vertex{ + top_right, + box_color + }; + + box_vertices[id + 2] = + sf::Vertex{ + top_right, + box_color + }; + + box_vertices[id + 3] = + sf::Vertex{ + bottom_right, + box_color + }; + + box_vertices[id + 4] = + sf::Vertex{ + bottom_right, + box_color + }; + + box_vertices[id + 5] = + sf::Vertex{ + bottom_left, + box_color + }; + + box_vertices[id + 6] = + sf::Vertex{ + bottom_left, + box_color + }; + + box_vertices[id + 7] = + sf::Vertex{ + top_left, + box_color + }; + } + + window.draw( + box_vertices + ); + } + + // ======================================================== + // Second half of velocity Verlet integration + // ======================================================== + + #pragma omp parallel for schedule(static) + for (int i = 0; + i < static_cast(sources.size()); + ++i) + { + sources[i].velocity += + 0.5 * + tree.forces[i] * + dt; + } + + // ======================================================== + // Create particle vertex array + // ======================================================== + + sf::VertexArray particle_vertices( + sf::PrimitiveType::Points, + sources.size() + ); + + // ======================================================== + // Fill particle vertex array + // ======================================================== + + #pragma omp parallel for schedule(static) + for (int i = 0; + i < static_cast(sources.size()); + ++i) + { + const float x = + static_cast( + sources[i].position.real() + ); + + const float y = + static_cast( + sources[i].position.imag() + ); + + particle_vertices[i] = + sf::Vertex{ + sf::Vector2f{x, y}, + particle_color + }; + } + + // -------------------------------------------------------- + // Draw particles + // -------------------------------------------------------- + + window.draw( + particle_vertices + ); + + // ======================================================== + // Frame timing + // ======================================================== + + const float milliseconds = + static_cast( + fpstimer + .getElapsedTime() + .asMicroseconds() + ) / 1000.0f; + + // ======================================================== + // Draw timing information + // ======================================================== + + sf::Text timing_text( + font, + std::to_string(milliseconds) + " ms", + 24 + ); + + timing_text.setFillColor( + sf::Color::White + ); + + timing_text.setPosition( + {10.0f, 10.0f} + ); + + window.draw( + timing_text + ); + + // ======================================================== + // Display frame + // ======================================================== + + window.display(); + + // ======================================================== + // Timing statistics + // ======================================================== + + if (total_frames < frame_limit) + { + times[total_frames] = + milliseconds; + } + + total_time += milliseconds; + + max_time = + std::max( + max_time, + milliseconds + ); + + min_time = + std::min( + min_time, + milliseconds + ); + + ++total_frames; + + // Uncomment to stop after frame_limit frames: + // + // if (total_frames >= frame_limit) + // { + // break; + // } + } + + // ============================================================ + // Final statistics + // ============================================================ + + if (total_frames > 0) + { + const float average_time = + total_time / + static_cast(total_frames); + + std::cout + << "\nSimulation finished.\n" + << "Frames: " + << total_frames + << '\n' + << "Average frame time: " + << average_time + << " ms\n" + << "Minimum frame time: " + << min_time + << " ms\n" + << "Maximum frame time: " + << max_time + << " ms\n"; + } + + return EXIT_SUCCESS; +}