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

More weather experiments

Browse source
This commit is contained in:
Nikita Lisitsa 2026-04-06 23:03:22 +03:00
parent d1a3bf15d4
commit 2a486c6d8f
1 changed files with 90 additions and 53 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

143
examples/weather.cpp
View file

@ -61,23 +61,29 @@ struct weather_app
{ {
static constexpr int N = 128; static constexpr int N = 128;
const bool static_mode = true; const bool static_mode = false;
const float dt = 20.f; const float dt = 20.f;
const float viscosity = static_mode ? 0.005f : 0.f; const float viscosity = static_mode ? 0.01f : 0.f;
const bool temperature_advection = true;
const float advection_magnification = 1.f; const float advection_magnification = 1.f;
const float temperature_diffusion = 0.001f; const float temperature_diffusion = 0.0004f;
const float cooling = 0.01f / 300.f; const float cooling = 0.01f / 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.f; const float water_heating_factor = 0.9f;
const float coriolis_bands = 2.f; const float coriolis = 0.001f;
const float coriolis_bands = 6.f;
const float band_force = 0.00001f;
const float friction = 0.f; const float friction = 0.f;
const float slope_force = 0.05f; const float slope_friction = 1.f;
const float slope_force = 0.001f;
const float land_force = 0.05f;
const float buoyancy_factor = 0.0002f * 0;
const float vorticity_confinement = 0.f; const float vorticity_confinement = 0.f;
const float elevation_temperature_drop = 30.f; const float elevation_temperature_drop = 30.f;
const float evaporation = 1.0f; const float evaporation = 1.0f;
const float max_humidity_factor = 100.f; const float max_humidity_factor = 1.f;
const float precipitation_factor = 0.0001f; const float precipitation_factor = 0.0001f;
const float force_field_amplitude = 0.00005f; const float force_field_amplitude = 0.00005f;
const float random_forces = 0.25f * (static_mode ? 0.f : 1.f); const float random_forces = 0.25f * (static_mode ? 0.f : 1.f);
@ -92,20 +98,20 @@ struct weather_app
gfx::pixmap_rgba biomes_map; gfx::pixmap_rgba biomes_map;
float expected_temperature_at(int y) const float expected_temperature_at(int y, bool water) const
{ {
// float latitude = (y - N * 0.5f) * 2.f / N; // float latitude = (y - N * 0.5f) * 2.f / N;
// return std::cos(latitude * float(math::pi)); // return std::cos(latitude * float(math::pi));
return temperature_income_at(y) * dt / (std::exp(cooling * dt) - 1.f); return temperature_income_at(y) * (water ? water_heating_factor : 1.f) * dt / (std::exp(cooling * dt) - 1.f);
} }
float temperature_income_at(int y) const float temperature_income_at(int y) const
{ {
// float latitude = (y - N * 0.5f) * 2.f / N; // float latitude = (y - N * 0.5f) * 2.f / N;
float latitude = y * 1.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.8f, 1.f, 1.f - std::abs(latitude)); // return heating * math::lerp(0.8f, 1.f, 1.f - std::abs(latitude));
} }
int wrap(int i) const int wrap(int i) const
@ -133,24 +139,6 @@ struct weather_app
precipitation_.resize({N, N}); precipitation_.resize({N, N});
average_precipitation_.resize({N, N}); average_precipitation_.resize({N, N});
// random::generator rng{0, 0};
random::uniform_ball_vector_distribution<float, 2> random_velocity{};
// for (auto & v : velocity_)
// {
// v = random_velocity(rng) * 0.01f;
// v += std::cos(0.5f * float(math::pi) * latitude * coriolis_bands
// }
for (int y = 0; y < N; ++y)
{
for (int x = 0; x < N; ++x)
{
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);
}
}
auto terrain_noise = make_perlin(rng, 2, 10, 1.6f); auto terrain_noise = make_perlin(rng, 2, 10, 1.6f);
@ -167,7 +155,7 @@ struct weather_app
} }
} }
auto heightmap = gfx::read_image<std::uint8_t>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "heightmap_seed_3.png"}); auto heightmap = gfx::read_image<std::uint8_t>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "heightmap_seed_1.png"});
for (int y = 0; y < N; ++y) for (int y = 0; y < N; ++y)
{ {
@ -177,6 +165,24 @@ struct weather_app
} }
} }
random::uniform_ball_vector_distribution<float, 2> random_velocity{};
// for (auto & v : velocity_)
// {
// v = random_velocity(rng) * 0.01f;
// v += std::cos(0.5f * float(math::pi) * latitude * coriolis_bands
// }
for (int y = 0; y < N; ++y)
{
for (int x = 0; x < N; ++x)
{
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, terrain_(x, y) <= 0.f);
}
}
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)
@ -186,6 +192,7 @@ struct weather_app
float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor; 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; humidity_(x, y) = max_humidity + dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f) * (1.f - precipitation_factor * dt) / precipitation_factor / dt;
humidity_(x, y) = 0.f;
} }
} }
@ -242,7 +249,7 @@ struct weather_app
std::swap(force_field_current_, force_field_next_); std::swap(force_field_current_, force_field_next_);
make_force_field(rng, force_field_next_, 0.5f); make_force_field(rng, force_field_next_, 0.5f);
} }
float const force_field_t = ((frame_ % force_field_switch_frames) + 0.5f) / force_field_switch_frames; [[maybe_unused]] float const force_field_t = ((frame_ % force_field_switch_frames) + 0.5f) / force_field_switch_frames;
int xmin = periodic_x ? 0 : 1; int xmin = periodic_x ? 0 : 1;
int xmax = periodic_x ? N : N - 1; // exclusive int xmax = periodic_x ? N : N - 1; // exclusive
@ -252,8 +259,9 @@ struct weather_app
{ {
for (int x = 0; x < N; ++x) for (int x = 0; x < N; ++x)
{ {
bool const is_water = terrain_(x, y) <= 0.f;
// temperature_(x, y) = math::lerp(temperature_(x, y), expected_temperature_at(y), 1.f - std::exp(- heating * dt)); // temperature_(x, y) = math::lerp(temperature_(x, y), expected_temperature_at(y), 1.f - std::exp(- heating * dt));
temperature_(x, y) += dt * temperature_income_at(y); temperature_(x, y) += dt * temperature_income_at(y) * (is_water ? water_heating_factor : 1.f);
temperature_(x, y) *= cooling_factor; temperature_(x, y) *= cooling_factor;
} }
} }
@ -268,7 +276,8 @@ struct weather_app
// float discharge = std::min(humidity_(x, y), precipitation_factor * dt); // float discharge = std::min(humidity_(x, y), precipitation_factor * dt);
// float discharge = 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 max_humidity = temperature_(x, y) * max_humidity_factor;
float max_humidity = temperature_(x, y) * temperature_(x, y) * max_humidity_factor;
float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt; float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt;
humidity_(x, y) -= discharge; humidity_(x, y) -= discharge;
precipitation_(x, y) = discharge / dt; precipitation_(x, y) = discharge / dt;
@ -293,6 +302,7 @@ struct weather_app
new_humidity_(N - 1, i) = humidity_(N - 1, 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)
{ {
for (int x = xmin; x < xmax; ++x) for (int x = xmin; x < xmax; ++x)
@ -334,7 +344,7 @@ struct weather_app
} }
} }
std::swap(velocity_, new_velocity_); std::swap(velocity_, new_velocity_);
std::swap(temperature_, new_temperature_); if (temperature_advection) std::swap(temperature_, new_temperature_);
std::swap(humidity_, new_humidity_); std::swap(humidity_, new_humidity_);
// Apply velocity diffusion // Apply velocity diffusion
@ -377,27 +387,40 @@ struct weather_app
{ {
for (int x = xmin; x < xmax; ++x) for (int x = xmin; x < xmax; ++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; [[maybe_unused]] 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 * 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) + random_forces * 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 // velocity_(x, y)[0] += dt * 0.000001f * std::cos(0.5f * float(math::pi) * latitude * 4.f);
// velocity_(x, y)[0] += dt * 0.000001f;
velocity_(x, y)[0] += dt * band_force * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands);
[[maybe_unused]] math::vector terrain_gradient
{ {
(std::max(0.f, terrain_(x + 1, y)) - std::max(0.f, terrain_(x - 1, y))) / 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, (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]] math::vector temperature_gradient
// velocity_(x, y) *= std::min(1.f, slope_factor); {
(temperature_(x + 1, y) - temperature_(x - 1, y)) / 2.f,
(temperature_(x, y + 1) - temperature_(x, y - 1)) / 2.f,
};
velocity_(x, y) += temperature_gradient * buoyancy_factor * dt;
[[maybe_unused]] float slope_factor = std::exp(- dt * slope_friction * 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; // 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(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)); [[maybe_unused]] float land_factor = std::exp(- dt * land_force * std::pow(std::max(0.f, terrain_(x, y)), 1.f));
velocity_(x, y) *= slope_factor; velocity_(x, y) *= land_factor;
// Directional external force // Directional external force
// velocity_(x, y)[1] += 0.001f * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands); // velocity_(x, y)[1] += 0.001f * dt * std::sin(0.5f * float(math::pi) * latitude * coriolis_bands);
@ -415,8 +438,8 @@ struct weather_app
float local_friction = friction * terrain_(x, y); float local_friction = friction * terrain_(x, y);
// velocity_(x, y) -= local_friction * velocity_(x, y) * math::length(velocity_(x, y)); // velocity_(x, y) -= local_friction * velocity_(x, y) * math::length(velocity_(x, y));
float local_friction_factor = std::exp(- local_friction * dt); [[maybe_unused]] float local_friction_factor = std::exp(- local_friction * dt);
velocity_(x, y) *= local_friction_factor; // velocity_(x, y) *= local_friction_factor;
} }
} }
@ -487,7 +510,7 @@ struct weather_app
{ {
if (!periodic_x) if (!periodic_x)
{ {
float left_boundary_flow = 0.01f;//0.01f * std::sin((i * 1.f / N) * float(math::pi) * 4.f); float left_boundary_flow = 0.f;//0.01f * std::sin((i * 1.f / N) * float(math::pi) * 4.f);
float right_boundary_flow = -left_boundary_flow; float right_boundary_flow = -left_boundary_flow;
velocity_(1, i)[0] = left_boundary_flow; velocity_(1, i)[0] = left_boundary_flow;
@ -499,8 +522,8 @@ struct weather_app
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];
velocity_(i, 1)[0] = 0.01f; // velocity_(i, 1)[0] = 0.01f;
velocity_(i, N-2)[0] = 0.01f; // velocity_(i, N-2)[0] = 0.01f;
} }
// Uncomment to visualize the force field // Uncomment to visualize the force field
@ -523,6 +546,18 @@ struct weather_app
// } // }
// } // }
// Apply boundary conditions for humidity
for (int i = 0; i < N; ++i)
{
if (!periodic_x)
{
humidity_(0, i) = 0.f;
humidity_(N - 1, i) = 0.f;
}
humidity_(i, 0) = 0.f;
humidity_(i, N - 1) = 0.f;
}
++frame_; ++frame_;
// Update all-time average temperature & precipitation // Update all-time average temperature & precipitation
@ -639,7 +674,7 @@ struct weather_app
auto map_biome = [this](float temperature, float precipitation) 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 x = math::clamp<int>(math::unlerp({ -3.f, 5.f}, precipitation) * biomes_map.width() , {0, biomes_map.width() - 1});
auto y = math::clamp<int>(math::unlerp({-10.f, 30.f}, temperature ) * biomes_map.height(), {0, biomes_map.height() - 1}); auto y = math::clamp<int>(math::unlerp({-10.f, 30.f}, temperature ) * biomes_map.height(), {0, biomes_map.height() - 1});
return gfx::to_colorf(biomes_map(x, y)); return gfx::to_colorf(biomes_map(x, y));
@ -690,14 +725,16 @@ struct weather_app
} }
} }
bool const is_water = terrain_(x, y) <= 0.f;
if (show_temperature_) if (show_temperature_)
color = map_temperature(temperature_(x, y) - 273.f); 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, is_water)), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f}));
if (show_average_temperature_delta_) if (show_average_temperature_delta_)
color = gfx::blend(color, map_color((average_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((average_temperature_(x, y) - expected_temperature_at(y, is_water)), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f}));
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}));
@ -750,7 +787,7 @@ struct weather_app
auto color = gfx::color_4f::zero(); auto color = gfx::color_4f::zero();
if (auto l = math::length(v); l > 0.f) if (auto l = math::length(v); l > 0.f)
{ {
float const magnification = 1000.f; float const magnification = 200.f;
float const max_length = 1.5f; float const max_length = 1.5f;
v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l; v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l;
@ -784,7 +821,7 @@ struct weather_app
push_text(std::format("P = {:.3f}", pressure_(x, y) * 1000.f)); 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, terrain_(x, y) <= 0.f) - 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("W = {:.3f}", humidity_(x, y)));
push_text(std::format("R = {:.3f}", precipitation_(x, y))); push_text(std::format("R = {:.3f}", precipitation_(x, y)));