From 02817be114aa6c9a584b10ebbc3e04a5352ecd75 Mon Sep 17 00:00:00 2001 From: lisyarus Date: Thu, 11 Dec 2025 16:44:36 +0300 Subject: [PATCH] Weather simulation: biome distribution --- examples/biomes.png | Bin 0 -> 4367 bytes examples/weather.cpp | 158 +++++++++++++++++++++++++++++++++++++++---- 2 files changed, 144 insertions(+), 14 deletions(-) create mode 100644 examples/biomes.png diff --git a/examples/biomes.png b/examples/biomes.png new file mode 100644 index 0000000000000000000000000000000000000000..00437a7e4e552cf3b489c429dff201810656167a GIT binary patch literal 4367 zcmeHKdu&rx7{4vd&CQXCL!B&o>3GQ8efD`TEwHtPRVWUqY`(qE?XtGF_TH`CgeV9q zfr!qxfT)uM38ILBz@ly_8=4@32I8BbF+M@#iy)%kz3sY9NX!!bZ(UB$x!?JHzweyi z`OdfARdw@b7Z+Yph#*LDs5Y*;l==VUQVO$FR zerU>lnYf?8*!UEffBvMz>5_=16uk+Z9SYYWRS?LC)!!RNl*TC%w|jt*q@659vlN9= z4%Y5r9WKO*+}n|ll!4D_XqION)(elh&$*H6er7rWze4f*>q34%nnDc}6M^HQJn^sZ&AzuP(aq~-Bx-Nzor%Q`CZnobwx9s6_5>r?YD`o!vMj;k*& zT|XTu`?j)WRC&3zGx*M}zo%21ckH_B;H#(quzdbO^04;2|6<7Fv=V918^Uv0L5<;@ zsPYn?jwQ@<^v*~pIH6h6QC@0P;y&BSt~YF`BKmCg&M+BH_@ySLwmm5=Y@Zhq+M5MN zw9S}a=uNXg5R-HcO~;~fjZOP(CNB$P!%Wyv6QVc!Yz>fYel;nf4%~s0STL=$(zfY^ zs5dFf?1I3oECoFIY)!hJU@yXUn~(0=jh{FRt3<5NOK8-!bu_)BL;hDdaxBpvI{z| zhZbo|NW=n3Q(KaP6l|5^dc|N0QONlxT9Q$7JEA~HQ7Hyc4Wd#*LK^on=V3_DsKgSc z7l=IssVnk;SVPuk%$VC5yby5D@eV=v@7)9gDja45s?cIAFBI_EjQFgm3W~^@Lylo6 zmtA&ZGAY|Jhd{Y7MwDfYcaS_!akAaxrUp@k;+oFI1<9ZSa$Er(x5)DjT5w_{;}J0j z=O!_R;W*6cak>N#O}d>7Gl=54qykyVMF+3Spb~+~DKmDLz%v-fIDyK|*fHMj^k9tJ zB~c>H^Da^{sYHRDr6yw>9H$cF8YLnTZ!{+i!P)A%kk3ZrSmQy?(yuou=W zqiQmus!^XUSDSOwVYp5rr*i?Q>Oh*S){*{deGO@6$;|Q%y09uLa@)V5jmv|20S`xB zt7y8KY%@)n8?{inBlDDbj4I|;LQ(Tlu$++Dg2uHXA*?e!Fs> zRN!Tgho(6PX6H#a<{%v$#=GH)JM0V(sxltRo(5uaZ znmQDnYLbjo0KqVd#whBaf)Uvb6GqMG-?5iCM-#6J7*u4yFJpt^1+|dK6~nA%M%nov zKiRYRA0q(uNRVOaJ0jPJT*FdeSm2TD8j)*Q3JeQ8l3o9qT!jO#Q&Jp$1*PC+>9+Z8 z0=#Nnz|WZ-KuE-cRBJuGN5OJoqV@(2LGTKrS<;(kw1P2T4~2vI-&|BwG#Ni*|Lz!= z%Jg7^?pI^RuX5y^ckcRSV=e=bKNP5r6db+2XZyIBWGT}9#im2j=ymUXw6`GIdqMQa z+4~Eq>e9;8EzQOAmmXevBG}u0bM5<2bdLV2{Odm_Rr&Ax;f1t?*jiGt;U(cyphH5z Ld4X*;w=e%2tXsi9 literal 0 HcmV?d00001 diff --git a/examples/weather.cpp b/examples/weather.cpp index b26c487a..cb16ac71 100644 --- a/examples/weather.cpp +++ b/examples/weather.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -10,6 +11,7 @@ #include #include #include +#include using namespace psemek; @@ -63,15 +65,19 @@ struct weather_app const float viscosity = 0.f; const float advection_magnification = 1.f; 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 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 friction = 0.f; const float slope_force = 0.04f; 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 random_forces = 0.f; 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 float friction_factor = 1.f - std::exp(- friction * dt); @@ -81,6 +87,8 @@ struct weather_app random::generator rng{random::device{}}; // random::generator rng{0, 0}; + gfx::pixmap_rgba biomes_map; + float expected_temperature_at(int y) { // float latitude = (y - N * 0.5f) * 2.f / N; @@ -91,9 +99,10 @@ struct weather_app 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.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 &) @@ -111,6 +120,10 @@ struct weather_app force_field_current_.resize({N, N}); force_field_next_.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::uniform_ball_vector_distribution random_velocity{}; @@ -141,12 +154,15 @@ struct weather_app (void)d; float value = terrain_noise((x + 0.5f) / N, (y + 0.5f) / N); 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_next_, 0.5f); + + biomes_map = gfx::read_image(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"}); } void on_event(app::key_event const & event) override @@ -171,6 +187,18 @@ struct weather_app if (event.down && event.key == app::keycode::H) 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 @@ -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) { return (i + N) % N; @@ -211,10 +256,16 @@ struct weather_app new_temperature_(i, 0) = temperature_(i, 0); 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) { new_temperature_(0, i) = temperature_(0, 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) @@ -249,10 +300,17 @@ struct weather_app math::lerp(temperature_(wrap(ix + 0), iy + 1), temperature_(wrap(ix + 1), iy + 1), tx), 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(temperature_, new_temperature_); + std::swap(humidity_, new_humidity_); // Apply velocity diffusion 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::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; math::vector terrain_gradient { - (terrain_(x + 1, y) - terrain_(x - 1, y)) / 2.f, - (terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f, + (std::max(0.f, terrain_(x + 1, y)) - std::max(0.f, terrain_(x - 1, y))) / 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)); @@ -435,7 +493,7 @@ struct weather_app ++frame_; - // Update all-time average temperature + // Update all-time average temperature & precipitation for (int y = 0; y < N; ++y) { for (int x = 0; x < N; ++x) @@ -443,6 +501,7 @@ struct weather_app float t = 1.f / frame_; // float t = 1.f / std::min(8192, frame_); 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))); }; + 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(math::unlerp({ -3.f, 3.f}, precipitation) * biomes_map.width() , {0, biomes_map.width() - 1}); + auto y = math::clamp(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 x = 0; x < N; ++x) { gfx::color_4f color = gfx::color_4f::zero(); - if (show_land_) - 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_land_ || show_biomes_) + { + 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_) - 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_) 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_) 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)); } } @@ -542,11 +658,14 @@ struct weather_app 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("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("A = {:.3f}", average_temperature_(x, 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("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()); @@ -564,18 +683,29 @@ private: bool show_average_temperature_delta_ = false; bool show_pressure_ = false; bool show_land_ = true; + bool show_biomes_ = true; + bool show_water_vapor_ = false; + bool show_precipitation_ = false; + bool show_average_precipitation_ = false; util::ndarray terrain_; util::ndarray, 2> velocity_; util::ndarray, 2> new_velocity_; util::ndarray pressure_; + util::ndarray vorticity_; + util::ndarray temperature_; util::ndarray new_temperature_; util::ndarray average_temperature_; + + util::ndarray humidity_; + util::ndarray new_humidity_; + util::ndarray precipitation_; + util::ndarray average_precipitation_; + util::ndarray, 2> force_field_main_; util::ndarray, 2> force_field_current_; util::ndarray, 2> force_field_next_; - util::ndarray vorticity_; int frame_ = 0; };