From e3cbbb5b4703e86ff5e1c714f482c9d4cf44c4e9 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Fri, 5 Dec 2025 17:10:52 +0300 Subject: [PATCH] Weather simulation test (wip) --- examples/weather.cpp | 326 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 326 insertions(+) create mode 100644 examples/weather.cpp diff --git a/examples/weather.cpp b/examples/weather.cpp new file mode 100644 index 00000000..b4c0824e --- /dev/null +++ b/examples/weather.cpp @@ -0,0 +1,326 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace psemek; + +struct weather_app + : app::application_base +{ + static constexpr int N = 128; + + weather_app(options const &, context const &) + { + simulation_box_ = {{{0.f, N}, {0.f, N}}}; + + obstacle_.resize({N, N}, 0.f); + velocity_.resize({N, N}); + new_velocity_.resize({N, N}); + pressure_.resize({N, N}, 0.f); + temperature_.resize({N, N}, -1.f); + new_temperature_.resize({N, N}); + + random::generator rng{random::device{}}; + // random::generator rng{0, 0}; + random::uniform_ball_vector_distribution random_velocity{}; + + for (auto & v : velocity_) + v = random_velocity(rng); + + for (int y = 0; y < N; ++y) + { + for (int x = 0; x < N; ++x) + { + // velocity_(x, y)[0] += (y*1.f/N-0.5f); + temperature_(x, y) = (x < N / 2) ? -1.f : 1.f; + } + } + } + + void on_event(app::key_event const & event) override + { + if (event.down && event.key == app::keycode::SPACE) + paused_ = !paused_; + } + + void update() override + { + if (paused_) + return; + + float const dt = 0.1f; + float const viscosity = 0.01f; + float const temperature_diffusion = 0.01f; + float const coriolis = 0.5f; + float const friction = 0.001f; + float const friction_factor = std::exp(- friction * dt); + + // Temperature source + // temperature_(16, 16) += 10.f * dt; + + // Boundary values + for (int y = 0; y < N; ++y) + { + // velocity_(1, y) = {(y*1.f/N)-0.5f, 0.f}; + // velocity_(N-2, y) = {(y*1.f/N)-0.5f, 0.f}; + + // if (y < N/2) + // { + // velocity_(1, y) = {-1.f, 0.f}; + // velocity_(N-2, y) = {-1.f, 0.f}; + // temperature_(N-1, y) = -1.f; + // } + // else + // { + // velocity_(1, y) = {1.f, 0.f}; + // velocity_(N-2, y) = {1.f, 0.f}; + + // temperature_(0, y) = 1.f; + // } + } + + // Velocity & temperature advection + for (int y = 1; y < N - 1; ++y) + { + for (int x = 1; x < N - 1; ++x) + { + auto v = velocity_(x, y); + auto p = math::point{x + 0.5f, y + 0.5f} - dt * v; + p[0] = math::clamp(p[0] - 0.5f, {0.f, N - 1.f}); + p[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f}); + + int ix = std::min(N - 1, std::floor(p[0])); + int iy = std::min(N - 1, std::floor(p[1])); + + float tx = p[0] - ix; + float ty = p[1] - iy; + + new_velocity_(x, y) = math::lerp( + math::lerp(velocity_(ix + 0, iy + 0), velocity_(ix + 1, iy + 0), tx), + math::lerp(velocity_(ix + 0, iy + 1), velocity_(ix + 1, iy + 1), tx), + ty + ); + + new_temperature_(x, y) = math::lerp( + math::lerp(temperature_(ix + 0, iy + 0), temperature_(ix + 1, iy + 0), tx), + math::lerp(temperature_(ix + 0, iy + 1), temperature_(ix + 1, iy + 1), tx), + ty + ); + } + } + std::swap(velocity_, new_velocity_); + std::swap(temperature_, new_temperature_); + + // Apply velocity diffusion & friction + for (int y = 0; y < N; ++y) + for (int x = 0; x < N; ++x) + new_velocity_(x, y) = velocity_(x, y); + for (int y = 1; y < N - 1; ++y) + { + for (int x = 1; x < N - 1; ++x) + { + // Velocity Laplacian + auto laplacian = velocity_(x + 1, y) + velocity_(x - 1, y) + velocity_(x, y + 1) + velocity_(x, y - 1) - 4.f * velocity_(x, y); + new_velocity_(x, y) = velocity_(x, y) + viscosity * dt * laplacian; + new_velocity_(x, y) *= friction_factor; + } + } + std::swap(velocity_, new_velocity_); + + // Apply temperature diffusion + for (int y = 0; y < N; ++y) + for (int x = 0; x < N; ++x) + new_temperature_(x, y) = temperature_(x, y); + for (int y = 1; y < N - 1; ++y) + { + for (int x = 1; x < N - 1; ++x) + { + // Velocity Laplacian + auto laplacian = temperature_(x + 1, y) + temperature_(x - 1, y) + temperature_(x, y + 1) + temperature_(x, y - 1) - 4.f * temperature_(x, y); + new_temperature_(x, y) = temperature_(x, y) + temperature_diffusion * dt * laplacian; + } + } + std::swap(temperature_, new_temperature_); + + // Apply forces + for (int y = 1; y < N - 1; ++y) + { + for (int x = 1; x < N - 1; ++x) + { + float latitude = (N * 0.5f - y) * 2.f / N; + // float latitude = (N - y) * 1.f / N; + velocity_(x, y) += math::ort(velocity_(x, y)) * (coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * 2.f)); + } + } + + // Solve Poisson equation for pressure + for (int iteration = 0; iteration < 16; ++iteration) + { + for (int y = 1; y < N - 1; ++y) + { + for (int x = 1; x < N - 1; ++x) + { + // Velocity divergence + float divergence = (velocity_(x + 1, y)[0] - velocity_(x - 1, y)[0] + velocity_(x, y + 1)[1] - velocity_(x, y - 1)[1]) / 2.f; + + // Gauss-Seidel iteration step + pressure_(x, y) = (pressure_(x - 1, y) + pressure_(x + 1, y) + pressure_(x, y - 1) + pressure_(x, y + 1) - divergence) / 4.f; + } + } + } + + // Apply boundary conditions for pressure + for (int i = 1; i < N - 1; ++i) + { + pressure_(0, i) = pressure_(1, i); + pressure_(N - 1, i) = pressure_(N - 2, i); + pressure_(i, 0) = pressure_(i, 1); + pressure_(i, N - 1) = pressure_(i, N - 2); + } + pressure_(0, 0) = (pressure_(0, 1) + pressure_(1, 0)) / 2.f; + pressure_(N-1, 0) = (pressure_(N-1, 1) + pressure_(N-2, 0)) / 2.f; + pressure_(0, N-1) = (pressure_(0, N-2) + pressure_(1, N-2)) / 2.f; + pressure_(N-1, N-1) = (pressure_(N-1, N-2) + pressure_(N-2, N-1)) / 2.f; + + // Normalize pressure + float average_pressure = 0.f; + for (auto const & value : pressure_) + average_pressure += value; + average_pressure /= (1.f * N * N); + for (auto & value : pressure_) + value -= average_pressure; + + // Project velocity into divergence-free space + // by subtracting pressure gradient + for (int y = 1; y < N - 1; ++y) + { + for (int x = 1; x < N - 1; ++x) + { + // Pressure gradient + math::vector gradient{pressure_(x + 1, y) - pressure_(x - 1, y), pressure_(x, y + 1) - pressure_(x, y - 1)}; + + velocity_(x, y) -= gradient; + } + } + + // Apply boundary conditions for velocity + for (int i = 1; i < N - 1; ++i) + { + velocity_(0, i)[0] = -velocity_(1, i)[0]; + velocity_(N-1, i)[0] = -velocity_(N-2, i)[0]; + velocity_(i, 0)[1] = -velocity_(i, 1)[1]; + velocity_(i, N-2)[1] = -velocity_(i, N-2)[1]; + } + } + + void present() override + { + gl::ClearColor(0.f, 0.f, 0.f, 0.f); + gl::Clear(gl::COLOR_BUFFER_BIT); + + float const aspect_ratio = state().size[0] * 1.f / state().size[1]; + math::box view_box = math::expand(simulation_box_, 1.f); + if (view_box[0].length() / view_box[1].length() > aspect_ratio) + view_box[1] = math::expand(view_box[1], (view_box[0].length() / aspect_ratio - view_box[1].length()) / 2.f); + else + view_box[0] = math::expand(view_box[0], (view_box[1].length() * aspect_ratio - view_box[0].length()) / 2.f); + + std::optional> mouseover_cell; + { + auto mouse = math::lerp(view_box, math::vector{state().mouse[0] * 1.f / state().size[0], 1.f - state().mouse[1] * 1.f / state().size[1]}); + int x = std::floor(mouse[0]); + int y = std::floor(mouse[1]); + + if (x >= 0 && x < N && y >= 0 && y < N) + mouseover_cell = {x, y}; + } + + [[maybe_unused]] float const pixel_size = view_box[0].length() / state().size[0]; + + for (int y = 0; y < N; ++y) + { + for (int x = 0; x < N; ++x) + { + math::point center{x + 0.5f, y + 0.5f}; + + auto v = velocity_(x, y); + if (auto l = math::length(v); l > 0.f) + { + float const magnification = 40.f; + float const max_length = 2.f; + v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l; + } + auto n = math::ort(v) * 0.25f; + + float value = 1000.f * pressure_(x, y); + // float value = 4.f * temperature_(x, y); + + gfx::color_4f color; + if (value > 0.f) + { + value = 1.f - std::exp(- value); + color = {value, value * 0.5f, value * 0.125f, 1.f}; + } + else + { + value = 1.f - std::exp(value); + color = {value * 0.125f, value * 0.5f, value, 1.f}; + } + + painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color)); + + painter_.triangle(center - v - n, center - v + n, center + v, {255, 255, 255, 255}); + // painter_.line(center - v, center + v, 1.f * pixel_size, {255, 255, 255, 255}, false); + } + } + + auto push_text = [&, row = 0](std::string const & text) mutable + { + painter_.text(view_box.corner(0.f, 1.f) - math::vector{0.f, row * pixel_size * 2.f * 12.f}, text, {.scale = {2.f * pixel_size, - 2.f * pixel_size}, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}}); + ++row; + }; + + if (mouseover_cell) + { + int x = (*mouseover_cell)[0]; + int y = (*mouseover_cell)[1]; + painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127}); + + push_text(std::format("{} {}", x, y)); + push_text(std::format("{} {}", velocity_(x, y)[0], velocity_(x, y)[1])); + } + + painter_.render(math::orthographic_camera{view_box}.transform()); + } + +private: + gfx::painter painter_; + + math::box simulation_box_; + + bool paused_ = true; + + util::ndarray obstacle_; + util::ndarray, 2> velocity_; + util::ndarray, 2> new_velocity_; + util::ndarray pressure_; + util::ndarray temperature_; + util::ndarray new_temperature_; +}; + +namespace psemek::app +{ + + std::unique_ptr make_application_factory() + { + return default_application_factory({.name = "Weather simulation test"}); + } + +}