lisyarus/psemek
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Weather simulation test: periodic X-boundary

Browse source
This commit is contained in:
Nikita Lisitsa 2025-12-06 14:07:19 +03:00
parent e3cbbb5b47
commit 253eda65ce
1 changed files with 98 additions and 83 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

159
examples/weather.cpp
View file

@ -39,7 +39,7 @@ struct weather_app
for (int x = 0; x < N; ++x) for (int x = 0; x < N; ++x)
{ {
// velocity_(x, y)[0] += (y*1.f/N-0.5f); // 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 void on_event(app::key_event const & event) override
{ {
if (event.down && event.key == app::keycode::SPACE) 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 void update() override
@ -55,62 +64,54 @@ struct weather_app
if (paused_) if (paused_)
return; return;
float const dt = 0.1f; float const dt = 0.2f;
float const viscosity = 0.01f; float const viscosity = 0.f;
float const temperature_diffusion = 0.01f; float const temperature_diffusion = 0.01f;
float const coriolis = 0.5f; float const heating = 0.001f;
float const friction = 0.001f; float const coriolis = 0.1f;
float const friction_factor = std::exp(- friction * dt); float const coriolis_bands = 5.f;
float const friction = 0.1f;
float const friction_factor = 1.f - std::exp(- friction * dt);
// Temperature source // Temperature source
// temperature_(16, 16) += 10.f * dt;
// Boundary values
for (int y = 0; y < N; ++y) for (int y = 0; y < N; ++y)
{ {
// velocity_(1, y) = {(y*1.f/N)-0.5f, 0.f}; for (int x = 0; x < N; ++x)
// velocity_(N-2, y) = {(y*1.f/N)-0.5f, 0.f}; {
temperature_(x, y) += dt * heating * (0.5f - std::abs(y*2.f/N - 1.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;
// }
} }
}
auto wrap = [](int i)
{
return (i + N) % N;
};
// Velocity & temperature advection // Velocity & temperature advection
for (int y = 1; y < N - 1; ++y) 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 v = velocity_(x, y);
auto p = math::point{x + 0.5f, y + 0.5f} - dt * v; 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[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}); p[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f});
int ix = std::min<int>(N - 1, std::floor(p[0])); int ix = std::floor(p[0]);
int iy = std::min<int>(N - 1, std::floor(p[1])); int iy = std::min<int>(N - 1, std::floor(p[1]));
float tx = p[0] - ix; float tx = p[0] - ix;
float ty = p[1] - iy; float ty = p[1] - iy;
new_velocity_(x, y) = math::lerp( new_velocity_(x, y) = math::lerp(
math::lerp(velocity_(ix + 0, iy + 0), velocity_(ix + 1, iy + 0), tx), math::lerp(velocity_(wrap(ix + 0), iy + 0), velocity_(wrap(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 + 1), velocity_(wrap(ix + 1), iy + 1), tx),
ty ty
); );
new_temperature_(x, y) = math::lerp( new_temperature_(x, y) = math::lerp(
math::lerp(temperature_(ix + 0, iy + 0), temperature_(ix + 1, iy + 0), tx), math::lerp(temperature_(wrap(ix + 0), iy + 0), temperature_(wrap(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 + 1), temperature_(wrap(ix + 1), iy + 1), tx),
ty ty
); );
} }
@ -118,18 +119,17 @@ struct weather_app
std::swap(velocity_, new_velocity_); std::swap(velocity_, new_velocity_);
std::swap(temperature_, new_temperature_); std::swap(temperature_, new_temperature_);
// Apply velocity diffusion & friction // Apply velocity diffusion
for (int y = 0; y < N; ++y) for (int y = 0; y < N; ++y)
for (int x = 0; x < N; ++x) for (int x = 0; x < N; ++x)
new_velocity_(x, y) = velocity_(x, y); new_velocity_(x, y) = velocity_(x, y);
for (int y = 1; y < N - 1; ++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 // 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) = velocity_(x, y) + viscosity * dt * laplacian;
new_velocity_(x, y) *= friction_factor;
} }
} }
std::swap(velocity_, new_velocity_); std::swap(velocity_, new_velocity_);
@ -140,23 +140,27 @@ struct weather_app
new_temperature_(x, y) = temperature_(x, y); new_temperature_(x, y) = temperature_(x, y);
for (int y = 1; y < N - 1; ++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 // Temperature Laplacian
auto laplacian = temperature_(x + 1, y) + temperature_(x - 1, y) + temperature_(x, y + 1) + temperature_(x, y - 1) - 4.f * temperature_(x, y); 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; new_temperature_(x, y) = temperature_(x, y) + temperature_diffusion * dt * laplacian;
} }
} }
std::swap(temperature_, new_temperature_); std::swap(temperature_, new_temperature_);
// Apply forces // Apply forces & friction
for (int y = 1; y < N - 1; ++y) 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 * 0.5f - y) * 2.f / N;
// float latitude = (N - y) * 1.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 y = 1; y < N - 1; ++y)
{ {
for (int x = 1; x < N - 1; ++x) for (int x = 0; x < N; ++x)
{ {
// Velocity divergence // 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 // 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 // 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_(0, i) = pressure_(1, i);
pressure_(N - 1, i) = pressure_(N - 2, i); // pressure_(N - 1, i) = pressure_(N - 2, i);
pressure_(i, 0) = pressure_(i, 1); pressure_(i, 0) = pressure_(i, 1);
pressure_(i, N - 1) = pressure_(i, N - 2); pressure_(i, N - 1) = pressure_(i, N - 2);
} }
pressure_(0, 0) = (pressure_(0, 1) + pressure_(1, 0)) / 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_(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_(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_(N-1, N-1) = (pressure_(N-1, N-2) + pressure_(N-2, N-1)) / 2.f;
// Normalize pressure // Normalize pressure
float average_pressure = 0.f; float average_pressure = 0.f;
@ -201,20 +205,19 @@ struct weather_app
// by subtracting pressure gradient // by subtracting pressure gradient
for (int y = 1; y < N - 1; ++y) 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 // 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; velocity_(x, y) -= gradient;
} }
} }
// Apply boundary conditions for velocity // 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_(0, i)[0] = -velocity_(1, i)[0];
velocity_(N-1, i)[0] = -velocity_(N-2, i)[0]; // velocity_(N-1, i)[0] = -velocity_(N-2, i)[0];
velocity_(i, 0)[1] = -velocity_(i, 1)[1]; velocity_(i, 0)[1] = -velocity_(i, 1)[1];
velocity_(i, N-2)[1] = -velocity_(i, N-2)[1]; velocity_(i, N-2)[1] = -velocity_(i, N-2)[1];
} }
@ -248,19 +251,13 @@ struct weather_app
{ {
for (int x = 0; x < N; ++x) for (int x = 0; x < N; ++x)
{ {
math::point center{x + 0.5f, y + 0.5f}; if (show_temperature_ || show_pressure_)
auto v = velocity_(x, y);
if (auto l = math::length(v); l > 0.f)
{ {
float const magnification = 40.f; float value = 0.f;
float const max_length = 2.f; if (show_temperature_)
v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l; value = 2.f * temperature_(x, y);
} else if (show_pressure_)
auto n = math::ort(v) * 0.25f; value = 1000.f * pressure_(x, y);
float value = 1000.f * pressure_(x, y);
// float value = 4.f * temperature_(x, y);
gfx::color_4f color; gfx::color_4f color;
if (value > 0.f) if (value > 0.f)
@ -275,9 +272,22 @@ struct weather_app
} }
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_.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); }
} }
} }
@ -294,7 +304,9 @@ struct weather_app
painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127}); painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127});
push_text(std::format("{} {}", x, y)); 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()); painter_.render(math::orthographic_camera{view_box}.transform());
@ -306,6 +318,9 @@ private:
math::box<float, 2> simulation_box_; math::box<float, 2> simulation_box_;
bool paused_ = true; bool paused_ = true;
bool show_velocity_ = true;
bool show_temperature_ = true;
bool show_pressure_ = false;
util::ndarray<float, 2> obstacle_; util::ndarray<float, 2> obstacle_;
util::ndarray<math::vector<float, 2>, 2> velocity_; util::ndarray<math::vector<float, 2>, 2> velocity_;