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

Weather simulation: biome distribution

Browse source
This commit is contained in:
Nikita Lisitsa 2025-12-11 16:44:36 +03:00
parent 69a2a04811
commit 02817be114
2 changed files with 144 additions and 14 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

BIN
examples/biomes.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 4.3 KiB

158
examples/weather.cpp
View file

@ -3,6 +3,7 @@
#include <psemek/gfx/painter.hpp> #include <psemek/gfx/painter.hpp>
#include <psemek/gfx/gl.hpp> #include <psemek/gfx/gl.hpp>
#include <psemek/math/camera.hpp> #include <psemek/math/camera.hpp>
#include <psemek/math/gradient.hpp>
#include <psemek/random/generator.hpp> #include <psemek/random/generator.hpp>
#include <psemek/random/device.hpp> #include <psemek/random/device.hpp>
#include <psemek/random/uniform_ball.hpp> #include <psemek/random/uniform_ball.hpp>
@ -10,6 +11,7 @@
#include <psemek/pcg/fractal.hpp> #include <psemek/pcg/fractal.hpp>
#include <psemek/util/ndarray.hpp> #include <psemek/util/ndarray.hpp>
#include <psemek/log/log.hpp> #include <psemek/log/log.hpp>
#include <psemek/io/file_stream.hpp>
using namespace psemek; using namespace psemek;
@ -63,15 +65,19 @@ struct weather_app
const float viscosity = 0.f; const float viscosity = 0.f;
const float advection_magnification = 1.f; const float advection_magnification = 1.f;
const float temperature_diffusion = 0.001f; const float temperature_diffusion = 0.001f;
const float cooling = 0.01f / 300.f; const float cooling = 0.1f / 300.f;
const float cooling_factor = std::exp(- cooling * dt); const float cooling_factor = std::exp(- cooling * dt);
const float heating = 323.f * (std::exp(cooling * dt) - 1.f) / dt; const float heating = 323.f * (std::exp(cooling * dt) - 1.f) / dt;
const float coriolis = 0.01f; const float coriolis = 0.f;
const float coriolis_bands = 2.f; const float coriolis_bands = 2.f;
const float friction = 0.f; const float friction = 0.f;
const float slope_force = 0.04f; const float slope_force = 0.04f;
const float vorticity_confinement = 0.f; const float vorticity_confinement = 0.f;
const float evaporation = 0.1f;
const float max_humidity_factor = 0.001f;
const float precipitation_factor = 0.0003f;
const float force_field_amplitude = 0.00005f; const float force_field_amplitude = 0.00005f;
const float random_forces = 0.f;
const float force_field_switch_duration = 720.f * 7.5f; // 7.5 days const float force_field_switch_duration = 720.f * 7.5f; // 7.5 days
const int force_field_switch_frames = std::round(force_field_switch_duration / dt); const int force_field_switch_frames = std::round(force_field_switch_duration / dt);
// const float friction_factor = 1.f - std::exp(- friction * dt); // const float friction_factor = 1.f - std::exp(- friction * dt);
@ -81,6 +87,8 @@ struct weather_app
random::generator rng{random::device{}}; random::generator rng{random::device{}};
// random::generator rng{0, 0}; // random::generator rng{0, 0};
gfx::pixmap_rgba biomes_map;
float expected_temperature_at(int y) float expected_temperature_at(int y)
{ {
// float latitude = (y - N * 0.5f) * 2.f / N; // float latitude = (y - N * 0.5f) * 2.f / N;
@ -91,9 +99,10 @@ struct weather_app
float temperature_income_at(int y) float temperature_income_at(int y)
{ {
float latitude = (y - N * 0.5f) * 2.f / N; // float latitude = (y - N * 0.5f) * 2.f / N;
float latitude = y * 1.f / N;
// return heating * math::lerp(0.75f, 1.f, std::cos(latitude * float(math::pi) / 2.f)); // return heating * math::lerp(0.75f, 1.f, std::cos(latitude * float(math::pi) / 2.f));
return heating * math::lerp(0.65f, 1.f, 1.f - std::abs(latitude)); return heating * math::lerp(0.8f, 1.f, 1.f - std::abs(latitude));
} }
weather_app(options const &, context const &) weather_app(options const &, context const &)
@ -111,6 +120,10 @@ struct weather_app
force_field_current_.resize({N, N}); force_field_current_.resize({N, N});
force_field_next_.resize({N, N}); force_field_next_.resize({N, N});
vorticity_.resize({N, N}); vorticity_.resize({N, N});
humidity_.resize({N, N});
new_humidity_.resize({N, N});
precipitation_.resize({N, N});
average_precipitation_.resize({N, N});
// random::generator rng{0, 0}; // random::generator rng{0, 0};
random::uniform_ball_vector_distribution<float, 2> random_velocity{}; random::uniform_ball_vector_distribution<float, 2> random_velocity{};
@ -141,12 +154,15 @@ struct weather_app
(void)d; (void)d;
float value = terrain_noise((x + 0.5f) / N, (y + 0.5f) / N); float value = terrain_noise((x + 0.5f) / N, (y + 0.5f) / N);
value = pow(value, 4.f) - d / 4.f; value = pow(value, 4.f) - d / 4.f;
terrain_(x, y) = std::max(0.f, math::lerp(1.f, 16.f, value)); value = math::lerp(1.f, 16.f, value);
terrain_(x, y) = value;
} }
} }
make_force_field(rng, force_field_main_, 0.5f); make_force_field(rng, force_field_main_, 0.5f);
make_force_field(rng, force_field_next_, 0.5f); make_force_field(rng, force_field_next_, 0.5f);
biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});
} }
void on_event(app::key_event const & event) override void on_event(app::key_event const & event) override
@ -171,6 +187,18 @@ struct weather_app
if (event.down && event.key == app::keycode::H) if (event.down && event.key == app::keycode::H)
show_land_ ^= true; show_land_ ^= true;
if (event.down && event.key == app::keycode::W)
show_water_vapor_ ^= true;
if (event.down && event.key == app::keycode::R)
show_precipitation_ ^= true;
if (event.down && event.key == app::keycode::Q)
show_average_precipitation_ ^= true;
if (event.down && event.key == app::keycode::B)
show_biomes_ ^= true;
} }
void update() override void update() override
@ -197,6 +225,23 @@ struct weather_app
} }
} }
// Evaporation
for (int y = 0; y < N; ++y)
{
for (int x = 0; x < N; ++x)
{
if (terrain_(x, y) <= 0.f)
humidity_(x, y) += dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f);
// float discharge = std::min(humidity_(x, y), precipitation_factor * dt);
float discharge = humidity_(x, y) * precipitation_factor * dt;
// float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor;
// float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt;
humidity_(x, y) -= discharge;
precipitation_(x, y) = discharge;
}
}
auto wrap = [](int i) auto wrap = [](int i)
{ {
return (i + N) % N; return (i + N) % N;
@ -211,10 +256,16 @@ struct weather_app
new_temperature_(i, 0) = temperature_(i, 0); new_temperature_(i, 0) = temperature_(i, 0);
new_temperature_(i, N - 1) = temperature_(i, N - 1); new_temperature_(i, N - 1) = temperature_(i, N - 1);
new_humidity_(i, 0) = humidity_(i, 0);
new_humidity_(i, N - 1) = humidity_(i, N - 1);
if (!periodic_x) if (!periodic_x)
{ {
new_temperature_(0, i) = temperature_(0, i); new_temperature_(0, i) = temperature_(0, i);
new_temperature_(N - 1, i) = temperature_(N - 1, i); new_temperature_(N - 1, i) = temperature_(N - 1, i);
new_humidity_(0, i) = humidity_(0, i);
new_humidity_(N - 1, i) = humidity_(N - 1, i);
} }
} }
for (int y = 1; y < N - 1; ++y) for (int y = 1; y < N - 1; ++y)
@ -249,10 +300,17 @@ struct weather_app
math::lerp(temperature_(wrap(ix + 0), iy + 1), temperature_(wrap(ix + 1), iy + 1), tx), math::lerp(temperature_(wrap(ix + 0), iy + 1), temperature_(wrap(ix + 1), iy + 1), tx),
ty ty
); );
new_humidity_(x, y) = math::lerp(
math::lerp(humidity_(wrap(ix + 0), iy + 0), humidity_(wrap(ix + 1), iy + 0), tx),
math::lerp(humidity_(wrap(ix + 0), iy + 1), humidity_(wrap(ix + 1), iy + 1), tx),
ty
);
} }
} }
std::swap(velocity_, new_velocity_); std::swap(velocity_, new_velocity_);
std::swap(temperature_, new_temperature_); std::swap(temperature_, new_temperature_);
std::swap(humidity_, new_humidity_);
// Apply velocity diffusion // Apply velocity diffusion
for (int y = 0; y < N; ++y) for (int y = 0; y < N; ++y)
@ -299,13 +357,13 @@ struct weather_app
// velocity_(x, y) += math::ort(velocity_(x, y)) * (coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands)); // velocity_(x, y) += math::ort(velocity_(x, y)) * (coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands));
velocity_(x, y) = math::rotate(velocity_(x, y), coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands)); velocity_(x, y) = math::rotate(velocity_(x, y), coriolis * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands));
auto force = force_field_main_(x, y) + math::lerp(force_field_current_(x, y), force_field_next_(x, y), force_field_t); auto force = force_field_main_(x, y) + random_forces * math::lerp(force_field_current_(x, y), force_field_next_(x, y), force_field_t);
velocity_(x, y) += (dt * force_field_amplitude) * force; velocity_(x, y) += (dt * force_field_amplitude) * force;
math::vector terrain_gradient math::vector terrain_gradient
{ {
(terrain_(x + 1, y) - terrain_(x - 1, y)) / 2.f, (std::max(0.f, terrain_(x + 1, y)) - std::max(0.f, terrain_(x - 1, y))) / 2.f,
(terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f, (std::max(0.f, terrain_(x, y + 1)) - std::max(0.f, terrain_(x, y - 1))) / 2.f,
}; };
[[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * math::dot(math::normalized(velocity_(x, y)), terrain_gradient)); [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * math::dot(math::normalized(velocity_(x, y)), terrain_gradient));
@ -435,7 +493,7 @@ struct weather_app
++frame_; ++frame_;
// Update all-time average temperature // Update all-time average temperature & precipitation
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)
@ -443,6 +501,7 @@ struct weather_app
float t = 1.f / frame_; float t = 1.f / frame_;
// float t = 1.f / std::min(8192, frame_); // float t = 1.f / std::min(8192, frame_);
average_temperature_(x, y) = math::lerp(average_temperature_(x, y), temperature_(x, y), t); average_temperature_(x, y) = math::lerp(average_temperature_(x, y), temperature_(x, y), t);
average_precipitation_(x, y) = math::lerp(average_precipitation_(x, y), precipitation_(x, y), t);
} }
} }
} }
@ -475,17 +534,65 @@ struct weather_app
return math::lerp(negative, positive, 1.f/ (1.f + std::exp(- value))); return math::lerp(negative, positive, 1.f/ (1.f + std::exp(- value)));
}; };
auto map_temperature = [&](float value) {
return map_color(2.f * std::round(value / 20.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f});
};
auto map_biome = [this](float temperature, float precipitation)
{
auto x = math::clamp<int>(math::unlerp({ -3.f, 3.f}, precipitation) * biomes_map.width() , {0, biomes_map.width() - 1});
auto y = math::clamp<int>(math::unlerp({-10.f, 40.f}, temperature ) * biomes_map.height(), {0, biomes_map.height() - 1});
return gfx::to_colorf(biomes_map(x, y));
};
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)
{ {
gfx::color_4f color = gfx::color_4f::zero(); gfx::color_4f color = gfx::color_4f::zero();
if (show_land_) if (show_land_ || show_biomes_)
color = gfx::blend(color, map_color(terrain_(x, y), {-1.f, -1.f, -1.f, 1.f}, {1.f, 1.f, 1.f, 1.f})); {
if (!show_biomes_)
{
if (terrain_(x, y) <= 0.f)
color = {0.5f, 0.5f, 1.f, 1.f};
else
color = {1.f, 1.f, 1.f, 1.f};
}
else
{
if (terrain_(x, y) <= 0.f)
color = map_color(4.f * terrain_(x, y), {0.f, 0.f, 0.125f, 1.f}, {0.f, 1.f, 1.5f, 1.f});
else
{
float temperature = average_temperature_(x, y) - 273.f - std::max(0.f, terrain_(x, y)) * 6.5f;
// float precipitation = (std::log10(std::max(1e-9f, average_precipitation_(x, y))) + 3.f) * 2.f;
// float precipitation = std::pow(average_precipitation_(x, y) * 125.f, 2.f) * 8.f;
float precipitation = std::log2(std::max(1e-9f, average_precipitation_(x, y)));
color = map_biome(temperature, precipitation);
}
}
if (show_land_ && x > 0 && x + 1 < N && y > 0 && y + 1 < N)
{
math::vector terrain_gradient
{
(terrain_(x + 1, y) - terrain_(x - 1, y)) / 2.f,
(terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f,
};
auto terrain_normal = math::normalized(math::vector{-terrain_gradient[0], -terrain_gradient[1], 0.5f});
float lightness = 0.5f + 0.5f * math::dot(terrain_normal, math::normalized(math::vector{1.f, 2.f, 3.f}));
color = gfx::dark(color, 1.f - lightness);
}
}
if (show_temperature_) if (show_temperature_)
color = gfx::blend(color, map_color(0.1f * (temperature_(x, y) - 273.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f})); color = map_temperature(temperature_(x, y) - 273.f);
if (show_temperature_delta_) if (show_temperature_delta_)
color = gfx::blend(color, map_color((temperature_(x, y) - expected_temperature_at(y)), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f})); color = gfx::blend(color, map_color((temperature_(x, y) - expected_temperature_at(y)), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f}));
@ -496,6 +603,15 @@ struct weather_app
if (show_pressure_) if (show_pressure_)
color = gfx::blend(color, map_color(10000.f * pressure_(x, y), {0.f, 0.f, 1.f, 0.75f}, {1.f, 0.f, 0.f, 0.75f})); color = gfx::blend(color, map_color(10000.f * pressure_(x, y), {0.f, 0.f, 1.f, 0.75f}, {1.f, 0.f, 0.f, 0.75f}));
if (show_water_vapor_)
color = gfx::blend(color, map_color(humidity_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
if (show_precipitation_)
color = gfx::blend(color, map_color(precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
if (show_average_precipitation_)
color = gfx::blend(color, map_color(average_precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
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));
} }
} }
@ -542,11 +658,14 @@ struct weather_app
push_text(std::format("{} {}", x, y)); push_text(std::format("{} {}", x, y));
push_text(std::format("V = {:.3f} {:.3f}", velocity_(x, y)[0] * 1000.f, velocity_(x, y)[1] * 1000.f)); push_text(std::format("V = {:.3f} {:.3f}", velocity_(x, y)[0] * 1000.f, velocity_(x, y)[1] * 1000.f));
push_text(std::format("P = {:.3f}", pressure_(x, y))); push_text(std::format("P = {:.3f}", pressure_(x, y) * 1000.f));
push_text(std::format("T = {:.3f}", temperature_(x, y) - 273.f)); push_text(std::format("T = {:.3f}", temperature_(x, y) - 273.f));
push_text(std::format("A = {:.3f}", average_temperature_(x, y) - 273.f)); push_text(std::format("A = {:.3f}", average_temperature_(x, y) - 273.f));
push_text(std::format("E = {:.3f}", expected_temperature_at(y) - 273.f)); push_text(std::format("E = {:.3f}", expected_temperature_at(y) - 273.f));
push_text(std::format("H = {:.3f}", terrain_(x, y))); push_text(std::format("H = {:.3f}", terrain_(x, y)));
push_text(std::format("W = {:.3f}", humidity_(x, y)));
push_text(std::format("R = {:.3f}", precipitation_(x, y)));
push_text(std::format("AR= {:.3f}", average_precipitation_(x, y)));
} }
painter_.render(math::orthographic_camera{view_box}.transform()); painter_.render(math::orthographic_camera{view_box}.transform());
@ -564,18 +683,29 @@ private:
bool show_average_temperature_delta_ = false; bool show_average_temperature_delta_ = false;
bool show_pressure_ = false; bool show_pressure_ = false;
bool show_land_ = true; bool show_land_ = true;
bool show_biomes_ = true;
bool show_water_vapor_ = false;
bool show_precipitation_ = false;
bool show_average_precipitation_ = false;
util::ndarray<float, 2> terrain_; util::ndarray<float, 2> terrain_;
util::ndarray<math::vector<float, 2>, 2> velocity_; util::ndarray<math::vector<float, 2>, 2> velocity_;
util::ndarray<math::vector<float, 2>, 2> new_velocity_; util::ndarray<math::vector<float, 2>, 2> new_velocity_;
util::ndarray<float, 2> pressure_; util::ndarray<float, 2> pressure_;
util::ndarray<float, 2> vorticity_;
util::ndarray<float, 2> temperature_; util::ndarray<float, 2> temperature_;
util::ndarray<float, 2> new_temperature_; util::ndarray<float, 2> new_temperature_;
util::ndarray<float, 2> average_temperature_; util::ndarray<float, 2> average_temperature_;
util::ndarray<float, 2> humidity_;
util::ndarray<float, 2> new_humidity_;
util::ndarray<float, 2> precipitation_;
util::ndarray<float, 2> average_precipitation_;
util::ndarray<math::vector<float, 2>, 2> force_field_main_; util::ndarray<math::vector<float, 2>, 2> force_field_main_;
util::ndarray<math::vector<float, 2>, 2> force_field_current_; util::ndarray<math::vector<float, 2>, 2> force_field_current_;
util::ndarray<math::vector<float, 2>, 2> force_field_next_; util::ndarray<math::vector<float, 2>, 2> force_field_next_;
util::ndarray<float, 2> vorticity_;
int frame_ = 0; int frame_ = 0;
}; };