From 253eda65cecdc75686538519cab63d4b83eaab7a Mon Sep 17 00:00:00 2001 From: lisyarus Date: Sat, 6 Dec 2025 14:07:19 +0300 Subject: [PATCH] Weather simulation test: periodic X-boundary --- examples/weather.cpp | 181 +++++++++++++++++++++++-------------------- 1 file changed, 98 insertions(+), 83 deletions(-) diff --git a/examples/weather.cpp b/examples/weather.cpp index b4c0824e..53aa2078 100644 --- a/examples/weather.cpp +++ b/examples/weather.cpp @@ -39,7 +39,7 @@ struct weather_app 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; + temperature_(x, y) = 0.5f - std::abs(y*2.f/N - 1.f); } } } @@ -47,7 +47,16 @@ struct weather_app void on_event(app::key_event const & event) override { if (event.down && event.key == app::keycode::SPACE) - paused_ = !paused_; + paused_ ^= true; + + if (event.down && event.key == app::keycode::V) + show_velocity_ ^= true; + + if (event.down && event.key == app::keycode::T) + show_temperature_ ^= true; + + if (event.down && event.key == app::keycode::P) + show_pressure_ ^= true; } void update() override @@ -55,62 +64,54 @@ struct weather_app if (paused_) return; - float const dt = 0.1f; - float const viscosity = 0.01f; + float const dt = 0.2f; + float const viscosity = 0.f; float const temperature_diffusion = 0.01f; - float const coriolis = 0.5f; - float const friction = 0.001f; - float const friction_factor = std::exp(- friction * dt); + float const heating = 0.001f; + float const coriolis = 0.1f; + float const coriolis_bands = 5.f; + float const friction = 0.1f; + float const friction_factor = 1.f - 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; - // } + for (int x = 0; x < N; ++x) + { + temperature_(x, y) += dt * heating * (0.5f - std::abs(y*2.f/N - 1.f)); + } } + auto wrap = [](int i) + { + return (i + N) % N; + }; + // Velocity & temperature advection for (int y = 1; y < N - 1; ++y) { - for (int x = 1; x < N - 1; ++x) + for (int x = 0; x < N; ++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 ix = 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), + math::lerp(velocity_(wrap(ix + 0), iy + 0), velocity_(wrap(ix + 1), iy + 0), tx), + math::lerp(velocity_(wrap(ix + 0), iy + 1), velocity_(wrap(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), + math::lerp(temperature_(wrap(ix + 0), iy + 0), temperature_(wrap(ix + 1), iy + 0), tx), + math::lerp(temperature_(wrap(ix + 0), iy + 1), temperature_(wrap(ix + 1), iy + 1), tx), ty ); } @@ -118,18 +119,17 @@ struct weather_app std::swap(velocity_, new_velocity_); std::swap(temperature_, new_temperature_); - // Apply velocity diffusion & friction + // Apply velocity diffusion 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) + for (int x = 0; x < N; ++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); + auto laplacian = velocity_(wrap(x + 1), y) + velocity_(wrap(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_); @@ -140,23 +140,27 @@ struct weather_app new_temperature_(x, y) = temperature_(x, y); for (int y = 1; y < N - 1; ++y) { - for (int x = 1; x < N - 1; ++x) + for (int x = 0; x < N; ++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); + // Temperature Laplacian + auto laplacian = temperature_(wrap(x + 1), y) + temperature_(wrap(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 + // Apply forces & friction for (int y = 1; y < N - 1; ++y) { - for (int x = 1; x < N - 1; ++x) + for (int x = 0; x < N; ++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)); + velocity_(x, y) += math::ort(velocity_(x, y)) * (coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands)); + + bool is_land = math::length(math::vector{x - N / 2.f, y - N / 4.f}) < N / 8.f; + if (is_land && false) + velocity_(x, y) -= friction_factor * velocity_(x, y);// * math::length(velocity_(x, y)); } } @@ -165,29 +169,29 @@ struct weather_app { for (int y = 1; y < N - 1; ++y) { - for (int x = 1; x < N - 1; ++x) + for (int x = 0; x < N; ++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; + float divergence = (velocity_(wrap(x + 1), y)[0] - velocity_(wrap(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; + pressure_(x, y) = (pressure_(wrap(x - 1), y) + pressure_(wrap(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) + for (int i = 0; i < N; ++i) { - pressure_(0, i) = pressure_(1, i); - pressure_(N - 1, i) = pressure_(N - 2, 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; + // 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; @@ -201,20 +205,19 @@ struct weather_app // by subtracting pressure gradient for (int y = 1; y < N - 1; ++y) { - for (int x = 1; x < N - 1; ++x) + for (int x = 0; x < N; ++x) { // Pressure gradient - math::vector gradient{pressure_(x + 1, y) - pressure_(x - 1, y), pressure_(x, y + 1) - pressure_(x, y - 1)}; - + math::vector gradient{pressure_(wrap(x + 1), y) - pressure_(wrap(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) + for (int i = 0; i < N; ++i) { - velocity_(0, i)[0] = -velocity_(1, i)[0]; - velocity_(N-1, i)[0] = -velocity_(N-2, i)[0]; + // 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]; } @@ -248,36 +251,43 @@ struct weather_app { 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) + if (show_temperature_ || show_pressure_) { - 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 = 0.f; + if (show_temperature_) + value = 2.f * temperature_(x, y); + else if (show_pressure_) + value = 1000.f * pressure_(x, y); - 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}; + } - 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_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color)); + if (show_velocity_) + { + 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 = 1.5f; + v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l; + } + auto n = math::ort(v) * 0.3f; - 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); + painter_.triangle(center - v - n, center - v + n, center + v, {255, 255, 255, 255}); + } } } @@ -294,7 +304,9 @@ struct weather_app 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])); + push_text(std::format("V = {} {}", velocity_(x, y)[0], velocity_(x, y)[1])); + push_text(std::format("P = {}", pressure_(x, y))); + push_text(std::format("T = {}", temperature_(x, y))); } painter_.render(math::orthographic_camera{view_box}.transform()); @@ -306,6 +318,9 @@ private: math::box simulation_box_; bool paused_ = true; + bool show_velocity_ = true; + bool show_temperature_ = true; + bool show_pressure_ = false; util::ndarray obstacle_; util::ndarray, 2> velocity_;