From f88d158b35ab006d2d8b3cdb32ccaa218a91e23d Mon Sep 17 00:00:00 2001 From: lisyarus Date: Fri, 12 Dec 2025 15:35:31 +0300 Subject: [PATCH] Weather simulation wip --- examples/biomes.png | Bin 4367 -> 4595 bytes examples/weather.cpp | 141 +++++++++++++++++++++++++++++++++++-------- 2 files changed, 117 insertions(+), 24 deletions(-) diff --git a/examples/biomes.png b/examples/biomes.png index 00437a7e4e552cf3b489c429dff201810656167a..36b3c9bb5a810424dcca536a930e8be1e05cf8a0 100644 GIT binary patch delta 725 zcmeBI`m8*mpOJau1ZUnqnhXrl3=9lldg22ab|Yge10yTb$;^ya^%M6$OJZPPY|eCc z4)An#RwyXSPs_|nWnie7Q#;Yt~0;9*=5O3eni=w3<(C zdY$9sHkG4C0}FI5er${=EhxNp^SGe_!~WPe4j^ zS!ArC!{uD#xgFj&@2nE3s&|#B6;Mp@NSkIj-AnbM^}ct44_?PkN&h+jck!3imEDYt zofgMFKAIG|TrD{zzDaKNob49bQ>C)+PTp@Y=X&_Goqok4oH|LInjQfG&XS__H##n@ znw+HQ`@2y7f5&V8iQbCgl07oqE1a6&{5tsmLHYKu#H>r+T*n*Vi68$d!R)l#Ve>nS z`^Wn(+b}krKQQ&)w{YnLZ`5CY-6m$q#I~)m`}MYEZy2XrANVj?arP!{)=aKerVn*) z_opzkjC_2JAC;t1*Z+=_2;+z`pxKmzi8@vtC(~K z2EORcSDCt)RLsmxl8uwiQ*=#D4J>p`(oBpQxc8M4b6=W5+^sauH>>X zHa0giFts!?n{38b0*~X3Z0a!f=KE|yjFb7;V^oY&lTA|0ElqS&EDch0O$<_!b(2hu z(sV77%q$EI3=$2DERrVovHQYxK4Q=35bHAA+XzfXY)RhkE9@k83-6wNC;h2zrV8tWRIt-pUXO@geCxZ C!X6<2 delta 671 zcmV;Q0$}~~Bab4Gj{yRakV6eD0RR920000i0RR92k>DH$F*G_eI65?w0|7dJQ1dyX z0004mX+uL$Nkc;*aB^>EX>4Tx0C=2zkv&MmKpe$iTct%S4t5Z6$WWau_=PyuDionY zs1;guFuC*#nlvOSE{=k0!NHHks)LKOt`4q(Aou~|E;uQ=NQwVT3N2zhIPS;0dyl(! zfWJ{;s@W9>RLwHd$%K&2tqQ?^SA@}v5XLbfF;kxt#S}cp*FAiEy^HZI?{j~SJ|%B5 zz$XyLnQmCb8^kl4md<&fIKoPjLVQjC;;dF`taVTR!f;+&S>`&;VI;ALC5R9pqlPknDzFfvRU^ej zn)VYO{vpRNl1nC68H^kYs6vJ0_`(0+ceiF?YSK*#CxGr3+x{2@ctJ^&f&Ds=-K90H>oMao|Hcz39?w|~#H`uhRj)N;~uq(cz^01bGv)B=hEBsXMY zVl*;iH7zhXIAkp}VK*=>IXPirEj2hbHe)z4FgGz z*ON&RAqNX4F8nj=;FEC?C|E;DL_t&t*Hyu>6+kc)0Kns=(Padi(AvTXY{mv`@JGl9 zF+mh6cU@}@QGL>V?nTCWTu9THqXgz0==~)uLx=6DcmXf05F?$hW8DA%002ovPDHLk FV1n$N5$XT{ diff --git a/examples/weather.cpp b/examples/weather.cpp index cb16ac71..8f19d7c7 100644 --- a/examples/weather.cpp +++ b/examples/weather.cpp @@ -43,8 +43,8 @@ auto make_perlin(random::generator & rng, int min_octave, int max_octave, float void make_force_field(random::generator & rng, util::ndarray, 2> & result, float scale) { - auto noise_1 = make_perlin(rng, 2, 6, 2.f); - auto noise_2 = make_perlin(rng, 2, 6, 2.f); + auto noise_1 = make_perlin(rng, 3, 6, 2.f); + auto noise_2 = make_perlin(rng, 3, 6, 2.f); for (int y = 0; y < result.height(); ++y) { @@ -61,7 +61,7 @@ struct weather_app { static constexpr int N = 128; - const float dt = 4.f; + const float dt = 20.f; const float viscosity = 0.f; const float advection_magnification = 1.f; const float temperature_diffusion = 0.001f; @@ -71,13 +71,13 @@ struct weather_app const float coriolis = 0.f; const float coriolis_bands = 2.f; const float friction = 0.f; - const float slope_force = 0.04f; + const float slope_force = 4.f; const float vorticity_confinement = 0.f; - const float evaporation = 0.1f; - const float max_humidity_factor = 0.001f; + const float evaporation = 0.01f; + const float max_humidity_factor = 100.f; const float precipitation_factor = 0.0003f; const float force_field_amplitude = 0.00005f; - const float random_forces = 0.f; + const float random_forces = 1.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); @@ -89,7 +89,7 @@ struct weather_app gfx::pixmap_rgba biomes_map; - float expected_temperature_at(int y) + float expected_temperature_at(int y) const { // float latitude = (y - N * 0.5f) * 2.f / N; // return std::cos(latitude * float(math::pi)); @@ -97,7 +97,7 @@ struct weather_app return temperature_income_at(y) * dt / (std::exp(cooling * dt) - 1.f); } - float temperature_income_at(int y) + float temperature_income_at(int y) const { // float latitude = (y - N * 0.5f) * 2.f / N; float latitude = y * 1.f / N; @@ -105,6 +105,11 @@ struct weather_app return heating * math::lerp(0.8f, 1.f, 1.f - std::abs(latitude)); } + int wrap(int i) const + { + return (i + N) % N; + } + weather_app(options const &, context const &) { simulation_box_ = {{{0.f, N}, {0.f, N}}}; @@ -141,6 +146,9 @@ struct weather_app float latitude = (N * 0.5f - y) * 2.f / N; velocity_(x, y) = random_velocity(rng) * 0.f + 0.f * math::vector{-std::cos(0.5f * float(math::pi) * latitude * coriolis_bands), 0.f}; temperature_(x, y) = expected_temperature_at(y); + + float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor; + humidity_(x, y) = max_humidity + dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f) * (1.f - precipitation_factor * dt) / precipitation_factor / dt; } } @@ -214,6 +222,9 @@ struct weather_app } float const force_field_t = ((frame_ % force_field_switch_frames) + 0.5f) / force_field_switch_frames; + int xmin = periodic_x ? 0 : 1; + int xmax = periodic_x ? N : N - 1; // exclusive + // Temperature source for (int y = 0; y < N; ++y) { @@ -234,22 +245,14 @@ struct weather_app 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; + // 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; - }; - - int xmin = periodic_x ? 0 : 1; - int xmax = periodic_x ? N : N - 1; // exclusive - // Velocity & temperature advection for (int i = 0; i < N; ++i) { @@ -366,10 +369,14 @@ struct weather_app (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)); - velocity_(x, y) *= std::min(1.f, slope_factor); + // [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * math::dot(math::normalized(velocity_(x, y)), terrain_gradient)); + // velocity_(x, y) *= std::min(1.f, slope_factor); // velocity_(x, y) -= terrain_gradient * slope_force * dt; + [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(math::length(terrain_gradient), 4.f)); + // [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(std::max(0.f, terrain_(x, y)), 1.f)); + velocity_(x, y) *= slope_factor; + // Directional external force // velocity_(x, y)[1] += 0.001f * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands); // velocity_(x, y) += math::direction(frame_ * dt * 2.f * float(math::pi) / 10080.f) * 0.001f * dt; @@ -528,6 +535,68 @@ struct weather_app mouseover_cell = {x, y}; } + if (mouseover_cell && (state().mouse_button_down.contains(app::mouse_button::left) || state().mouse_button_down.contains(app::mouse_button::right))) + { + float delta = state().mouse_button_down.contains(app::mouse_button::left) ? 1.f : -1.f; + + int R = 4; + + for (int iy = -R; iy <= R; ++iy) + { + for (int ix = -R; ix <= R; ++ix) + { + int x = (*mouseover_cell)[0] + ix; + int y = (*mouseover_cell)[1] + iy; + + if (x >= 0 && x < N && y >= 0 && y < N) + { + auto d = math::vector{ix, iy} / 2.f; + terrain_(x, y) += delta * 0.05f * std::exp(- math::dot(d, d)); + } + } + } + } + + if (mouseover_cell && state().mouse_button_down.contains(app::mouse_button::middle)) + { + int R = 4; + + float average = 0.f; + int count = 0; + + for (int iy = -R; iy <= R; ++iy) + { + for (int ix = -R; ix <= R; ++ix) + { + int x = (*mouseover_cell)[0] + ix; + int y = (*mouseover_cell)[1] + iy; + + if (x >= 0 && x < N && y >= 0 && y < N) + { + average += terrain_(x, y); + count += 1; + } + } + } + + average /= count; + + for (int iy = -R; iy <= R; ++iy) + { + for (int ix = -R; ix <= R; ++ix) + { + int x = (*mouseover_cell)[0] + ix; + int y = (*mouseover_cell)[1] + iy; + + if (x >= 0 && x < N && y >= 0 && y < N) + { + auto d = math::vector{ix, iy} / 2.f; + terrain_(x, y) += (average - terrain_(x, y)) * 0.05f * std::exp(- math::dot(d, d)); + } + } + } + } + [[maybe_unused]] float const pixel_size = view_box[0].length() / state().size[0]; auto map_color = [](float value, gfx::color_4f const & negative, gfx::color_4f const & positive){ @@ -604,10 +673,34 @@ struct weather_app 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})); + color = gfx::blend(color, map_color(humidity_(x, y) * 0.001f, {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})); + { + // color = gfx::blend(color, map_color(precipitation_(x, y), {1.f, 1.f, 1.f, -1.f}, {1.f, 1.f, 1.f, 1.f})); + + float alpha = 2.f / (1.f + std::exp(- precipitation_(x, y))) - 1.f; + + gfx::color_4f cloud_color{1.f, 1.f, 1.f, alpha * 0.875f}; + + if (y > 0 && y + 1 < N) + { + if (periodic_x || (x > 0 && x + 1 < N)) + { + math::vector gradient + { + (precipitation_(wrap(x + 1), y) - precipitation_(wrap(x - 1), y)) / 2.f, + (precipitation_(x, y + 1) - precipitation_(x, y - 1)) / 2.f, + }; + + auto normal = math::normalized(math::vector{-gradient[0], -gradient[1], 2.f}); + auto lightness = 0.5f + 0.5f * math::dot(normal, math::normalized(math::vector{1.f, 2.f, 3.f})); + cloud_color = gfx::dark(cloud_color, 1.f - lightness); + } + } + + color = gfx::blend(color, cloud_color); + } 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}));