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

More weather v2 tests

Browse source
This commit is contained in:
Nikita Lisitsa 2026-04-12 17:13:23 +03:00
parent 5b98171283
commit 2d9355f829
1 changed files with 269 additions and 77 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

346
examples/weather_v2.cpp
View file

@ -44,6 +44,20 @@ auto make_perlin(random::generator & rng, int min_octave, int max_octave, float
return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights)); return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights));
} }
void add_random_scalar_field(random::generator & rng, util::ndarray<float, 2> & result, int min_octave, int max_octave, float scale)
{
auto noise = make_perlin(rng, min_octave, max_octave, 2.f);
for (int y = 0; y < result.height(); ++y)
{
for (int x = 0; x < result.width(); ++x)
{
math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
result(x, y) += noise(p) * scale;
}
}
}
void add_random_vector_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, int min_octave, int max_octave, float scale) void add_random_vector_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, int min_octave, int max_octave, float scale)
{ {
auto noise_1 = make_perlin(rng, min_octave, max_octave, 2.f); auto noise_1 = make_perlin(rng, min_octave, max_octave, 2.f);
@ -59,6 +73,21 @@ void add_random_vector_field(random::generator & rng, util::ndarray<math::vector
} }
} }
template <typename T>
void mult_random_scalar_field(random::generator & rng, util::ndarray<T, 2> & result, int min_octave, int max_octave, float min, float max)
{
auto noise = make_perlin(rng, min_octave, max_octave, 2.f);
for (int y = 0; y < result.height(); ++y)
{
for (int x = 0; x < result.width(); ++x)
{
math::point p{(x + 0.5f) / result.height(), (y + 0.5f) / result.width()};
result(x, y) *= math::lerp(min, max, noise(p));
}
}
}
template <typename T> template <typename T>
void add(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src) void add(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src)
{ {
@ -83,7 +112,11 @@ float do_lerp(float x, float y, float t)
math::vector<float, 2> do_lerp(math::vector<float, 2> const & x, math::vector<float, 2> const & y, float t) math::vector<float, 2> do_lerp(math::vector<float, 2> const & x, math::vector<float, 2> const & y, float t)
{ {
return math::lerp(length(x), length(y), t) * math::slerp(math::normalized(x), math::normalized(y), t); auto lx = length(x);
auto ly = length(y);
if (lx > 0.f && ly > 0.f)
return math::lerp(lx, ly, t) * math::slerp(x / lx, y / ly, t);
return math::lerp(x, y, t);
} }
template <typename T> template <typename T>
@ -160,7 +193,7 @@ struct solver
climate_snapshot & snapshot; climate_snapshot & snapshot;
int solve_velocity_iterations = 0; int solve_velocity_iterations = 0;
int max_solve_velocity_iterations = N / 2; int max_solve_velocity_iterations = 1;
int solve_temperature_iterations = 0; int solve_temperature_iterations = 0;
int max_solve_temperature_iterations = N / 4; int max_solve_temperature_iterations = N / 4;
@ -169,6 +202,7 @@ struct solver
int max_solve_humidity_iterations = 3 * N; int max_solve_humidity_iterations = 3 * N;
util::ndarray<math::vector<float, 2>, 2> initial_random_wind_velocity; util::ndarray<math::vector<float, 2>, 2> initial_random_wind_velocity;
util::ndarray<float, 2> initial_random_wind_velocity_potential;
util::ndarray<float, 2> pressure; util::ndarray<float, 2> pressure;
util::ndarray<float, 2> new_temperature; util::ndarray<float, 2> new_temperature;
util::ndarray<float, 2> new_humidity; util::ndarray<float, 2> new_humidity;
@ -177,14 +211,14 @@ struct solver
bool finished = false; bool finished = false;
solver(util::ndarray<float, 2> const & terrain, random::generator & rng, int season, int day_night, climate_snapshot & snapshot, solver(util::ndarray<float, 2> const & terrain, random::generator & rng, int season, int day_night, climate_snapshot & snapshot,
util::ndarray<math::vector<float, 2>, 2> const & base_wind_velocity) util::ndarray<float, 2> const & base_wind_velocity_potential)
: terrain(terrain) : terrain(terrain)
, rng(rng) , rng(rng)
, season(season) , season(season)
, day_night(day_night) , day_night(day_night)
, snapshot(snapshot) , snapshot(snapshot)
{ {
initial_random_wind_velocity.resize({N, N}); initial_random_wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f}); snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
pressure.resize({N, N}, 0.f); pressure.resize({N, N}, 0.f);
snapshot.temperature.resize({N, N}, 0.f); snapshot.temperature.resize({N, N}, 0.f);
@ -192,8 +226,23 @@ struct solver
snapshot.humidity.resize({N, N}, 0.f); snapshot.humidity.resize({N, N}, 0.f);
new_humidity.resize({N, N}, 0.f); new_humidity.resize({N, N}, 0.f);
initial_random_wind_velocity = base_wind_velocity.copy(); initial_random_wind_velocity_potential = base_wind_velocity_potential.copy();
add_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.002f); add_random_scalar_field(rng, initial_random_wind_velocity_potential, 2, 5, 0.4f);
// add_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.002f);
for (int y = 1; y + 1 < N; ++y)
for (int x = 1; x + 1 < N; ++x)
{
math::vector potential_gradient
{
(initial_random_wind_velocity_potential(x + 1, y) - initial_random_wind_velocity_potential(x - 1, y)) * 0.5f,
(initial_random_wind_velocity_potential(x, y + 1) - initial_random_wind_velocity_potential(x, y - 1)) * 0.5f,
};
initial_random_wind_velocity(x, y) = math::ort(potential_gradient) * 0.5f;
}
// mult_random_scalar_field(rng, initial_random_wind_velocity, 2, 5, 0.01f, 1.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)
@ -220,7 +269,7 @@ struct solver
float expected_humidity(int x, int y) float expected_humidity(int x, int y)
{ {
if (terrain(x, y) < 0.f) if (terrain(x, y) < 0.f)
return 1.f; return 1.5f;
else else
return 0.f; return 0.f;
} }
@ -229,37 +278,37 @@ struct solver
{ {
if (stage == 0) if (stage == 0)
{ {
for (int y = 1; y + 1 < N; ++y) // for (int y = 1; y + 1 < N; ++y)
{ // {
for (int x = 1; x + 1 < N; ++x) // for (int x = 1; x + 1 < N; ++x)
{ // {
float divergence = 0.f // float divergence = 0.f
+ (initial_random_wind_velocity(x + 1, y)[0] - initial_random_wind_velocity(x - 1, y)[0]) * 0.5f // + (initial_random_wind_velocity(x + 1, y)[0] - initial_random_wind_velocity(x - 1, y)[0]) * 0.5f
+ (initial_random_wind_velocity(x, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f // + (initial_random_wind_velocity(x, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f
; // ;
// (sum - 4.f * p(x,y)) = divergence // // (sum - 4.f * p(x,y)) = divergence
// sum - 4.f * p = divergence // // sum - 4.f * p = divergence
// p = (sum - divergence) / 4.f // // p = (sum - divergence) / 4.f
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(x - 1, y) + pressure(x + 1, y) + pressure(x, y - 1) + pressure(x, y + 1) - divergence) / 4.f;
} // }
} // }
for (int y = 1; y + 1 < N; ++y) // for (int y = 1; y + 1 < N; ++y)
{ // {
for (int x = 1; x + 1 < N; ++x) // for (int x = 1; x + 1 < N; ++x)
{ // {
math::vector pressure_gradient // math::vector pressure_gradient
{ // {
(pressure(x + 1, y) - pressure(x - 1, y)) * 0.5f, // (pressure(x + 1, y) - pressure(x - 1, y)) * 0.5f,
(pressure(x, y + 1) - pressure(x, y - 1)) * 0.5f, // (pressure(x, y + 1) - pressure(x, y - 1)) * 0.5f,
}; // };
snapshot.wind_velocity(x, y) = initial_random_wind_velocity(x, y) - pressure_gradient; // snapshot.wind_velocity(x, y) = initial_random_wind_velocity(x, y) - pressure_gradient;
if (terrain(x, y) > 0.f) // if (terrain(x, y) > 0.f)
snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y)); // snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
} // }
} // }
++solve_velocity_iterations; ++solve_velocity_iterations;
if (solve_velocity_iterations == max_solve_velocity_iterations) if (solve_velocity_iterations == max_solve_velocity_iterations)
@ -319,10 +368,17 @@ struct solver
math::point p{x + 0.5f, y + 0.5f}; math::point p{x + 0.5f, y + 0.5f};
float advection = sample_bilinear(snapshot.humidity, p - snapshot.wind_velocity(x, y) * step_time); float advection = sample_bilinear(snapshot.humidity, p - snapshot.wind_velocity(x, y) * step_time);
float evaporation_speed = 0.0001f; float income = 0.f;
float evaporation_delta = (expected_humidity(x, y) - snapshot.humidity(x, y)) * (- std::expm1(- evaporation_speed * step_time)); if (terrain(x, y) < 0.f)
income = step_time * 0.00012f;
else
// income = snapshot.humidity(x, y) * std::expm1(- 0.00006f * step_time);
income = -std::min(snapshot.humidity(x, y), step_time * 0.00001f);
new_humidity(x, y) = evaporation_delta + advection; // float evaporation_speed = 0.0001f;
// float evaporation_delta = (expected_humidity(x, y) - snapshot.humidity(x, y)) * (- std::expm1(- evaporation_speed * step_time));
new_humidity(x, y) = income + advection;
} }
} }
std::swap(snapshot.humidity, new_humidity); std::swap(snapshot.humidity, new_humidity);
@ -386,6 +442,7 @@ struct weather_app
bool show_humidity = false; bool show_humidity = false;
bool show_biomes = true; bool show_biomes = true;
bool show_rivers = true; bool show_rivers = true;
bool show_clouds = true;
gfx::pixmap_rgba biomes_map; gfx::pixmap_rgba biomes_map;
util::ndarray<float, 2> terrain; util::ndarray<float, 2> terrain;
@ -401,17 +458,27 @@ struct weather_app
int day_night = 0; int day_night = 0;
util::ndarray<math::vector<float, 2>, 2> base_wind_velocity; util::ndarray<math::vector<float, 2>, 2> base_wind_velocity;
util::ndarray<float, 2> base_wind_velocity_potential;
climate_snapshot snapshots[4][2]; climate_snapshot snapshots[4][2];
climate_snapshot average; climate_snapshot average;
climate_snapshot display_snapshot; climate_snapshot display_snapshot;
bool need_update_display_snapshot = false;
std::optional<struct solver> solver; std::optional<struct solver> solver;
river_net rivers; river_net rivers;
util::ndarray<float, 2> cloud_density;
util::ndarray<float, 2> new_cloud_density;
struct cloud_particle
{
math::point<float, 2> position;
int lifetime;
};
std::vector<cloud_particle> cloud_particles;
float display_season = 0.f; float display_season = 0.f;
util::clock<> frame_clock; util::clock<> frame_clock;
@ -428,11 +495,16 @@ struct weather_app
biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"}); biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});
base_wind_velocity.resize({N, N}, math::vector{0.f, 0.f}); // base_wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
// add_random_vector_field(rng, base_wind_velocity, 3, 6, 0.003f);
add_random_vector_field(rng, base_wind_velocity, 3, 6, 0.003f); base_wind_velocity_potential.resize({N, N}, 0.f);
add_random_scalar_field(rng, base_wind_velocity_potential, 1, 3, 0.6f);
context.vsync(false); cloud_density.resize({N, N}, 0.f);
new_cloud_density.resize({N, N}, 0.f);
context.vsync(true);
} }
void on_event(app::key_event const & event) override void on_event(app::key_event const & event) override
@ -456,12 +528,17 @@ struct weather_app
if (event.down && event.key == app::keycode::R) if (event.down && event.key == app::keycode::R)
show_rivers ^= true; show_rivers ^= true;
if (event.down && event.key == app::keycode::C)
show_clouds ^= true;
} }
void update() override void update() override
{ {
float const dt = frame_clock.restart().count(); float const dt = frame_clock.restart().count();
bool need_update_display_snapshot = false;
if (solver && solver->finished) if (solver && solver->finished)
{ {
solver.reset(); solver.reset();
@ -480,34 +557,6 @@ struct weather_app
} }
} }
if (!solver && season < 4)
solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night], base_wind_velocity);
if (!paused)
{
if (solver)
{
for (int i = 0; i < 64; ++i)
{
solver->step();
if (solver->finished)
break;
}
}
else if (season == 4 && !rivers.queue.empty())
{
for (int i = 0; i < 64; ++i)
{
propagate_rivers();
if (rivers.queue.empty())
{
compute_river_flow();
break;
}
}
}
}
if (state().key_down.contains(app::keycode::LEFT)) if (state().key_down.contains(app::keycode::LEFT))
{ {
display_season -= 0.5f * dt; display_season -= 0.5f * dt;
@ -524,6 +573,144 @@ struct weather_app
display_season += 1.f; display_season += 1.f;
while (display_season >= 1.f) while (display_season >= 1.f)
display_season -= 1.f; display_season -= 1.f;
if (season == 4 && need_update_display_snapshot)
{
update_display_snapshot();
need_update_display_snapshot = false;
}
if (!solver && season < 4)
solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night], base_wind_velocity_potential);
if (!paused)
{
if (solver)
{
for (int i = 0; i < 64; ++i)
// for (int i = 0; i < 1; ++i)
{
solver->step();
if (solver->finished)
break;
}
}
else if (season == 4 && !rivers.queue.empty())
{
for (int i = 0; i < 64; ++i)
{
propagate_rivers();
if (rivers.queue.empty())
{
compute_river_flow();
paused = true;
break;
}
}
}
else
{
float const dt = 10.f;
// if (false)
// for (int y = 0; y < N; ++y)
// {
// for (int x = 0; x < N; ++x)
// {
// new_cloud_density(x, y) = sample_bilinear(cloud_density, math::point{x + 0.5f, y + 0.5f} - display_snapshot.wind_velocity(x, y) * dt);
// }
// }
// for (int y = 0; y < N; ++y)
// {
// for (int x = 0; x < N; ++x)
// {
// new_cloud_density(x, y) += (terrain(x, y) < 0.f ? 0.0001f : 0.f) * dt;
// }
// }
// // if (false)
// for (int y = 1; y + 1 < N; ++y)
// {
// for (int x = 1; x + 1 < N; ++x)
// {
// auto gradient = math::vector {
// (new_cloud_density(x + 1, y) - new_cloud_density(x - 1, y)) / 2.f,
// (new_cloud_density(x, y + 1) - new_cloud_density(x, y - 1)) / 2.f,
// };
// if (math::length(gradient) > 0.f)
// cloud_density(x, y) = sample_bilinear(new_cloud_density, math::point{x + 0.5f, y + 0.5f} - math::normalized(gradient) * 1.f);
// else
// cloud_density(x, y) = new_cloud_density(x, y);
// }
// }
// std::swap(new_cloud_density, cloud_density);
int total_cloud_particles = N * N / 8;
auto new_cloud_particle = [&]
{
while (true)
{
cloud_particle p{{random::uniform<float>(rng, 0.f, N), random::uniform<float>(rng, 0.f, N)}, 0};
if (sample_bilinear(terrain, p.position) < 0.f)
return p;
}
};
if (cloud_particles.size() < total_cloud_particles)
for (int i = 0; i < 16; ++i)
cloud_particles.push_back(new_cloud_particle());
util::ndarray<std::vector<cloud_particle *>, 2> cloud_particle_index({N, N});
for (auto & p : cloud_particles)
{
p.position += sample_bilinear(display_snapshot.wind_velocity, p.position) * dt;
p.lifetime += 1;
auto r = math::distance_sqr(p.position, math::point{N/2.f, N/2.f}) / math::sqr(N/2.f);
auto t = std::min(1.f, 2.f * std::exp(- 1.f * r));
if (p.lifetime >= 1024 || p.position[0] < 0.f || p.position[0] > N || p.position[1] < 0.f || p.position[1] > N || (p.lifetime > 16 && random::uniform<float>(rng) > t))
{
p = new_cloud_particle();
}
cloud_particle_index(std::min<int>(N-1, std::floor(p.position[0])), std::min<int>(N-1, std::floor(p.position[1]))).push_back(&p);
}
for (auto & p : cloud_particles)
{
int x = std::min<int>(N-1, std::floor(p.position[0]));
int y = std::min<int>(N-1, std::floor(p.position[1]));
for (int tx = -1; tx <= 1; ++tx)
{
for (int ty = -1; ty <= 1; ++ty)
{
int xx = x + tx;
int yy = y + ty;
if (xx < 0 || xx + 1 >= N || yy < 0 || yy + 1 >= N) continue;
for (auto q : cloud_particle_index(xx, yy))
{
if (q == &p) continue;
auto r = p.position - q->position;
if (auto l = math::length(r); l < 1.f)
{
auto n = r / l;
p.position += n * dt * 0.02f;
}
}
}
}
}
}
}
} }
void compute_average() void compute_average()
@ -679,12 +866,6 @@ struct weather_app
void present() override void present() override
{ {
if (need_update_display_snapshot)
{
update_display_snapshot();
need_update_display_snapshot = false;
}
gl::ClearColor(0.f, 0.f, 0.f, 0.f); gl::ClearColor(0.f, 0.f, 0.f, 0.f);
gl::Clear(gl::COLOR_BUFFER_BIT); gl::Clear(gl::COLOR_BUFFER_BIT);
@ -780,6 +961,11 @@ struct weather_app
color = gfx::blend(color, map_humidity(snapshot.humidity(x, y))); color = gfx::blend(color, map_humidity(snapshot.humidity(x, y)));
} }
// if (show_clouds)
// {
// color = gfx::blend(color, map_color_positive(cloud_density(x, y), {1.f, 1.f, 1.f, 0.f}, {1.f, 1.f, 1.f, 1.f}));
// }
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));
} }
} }
@ -839,6 +1025,12 @@ struct weather_app
} }
} }
if (show_clouds)
{
for (auto const & p : cloud_particles)
painter.circle(p.position, 2.f, {255, 255, 255, std::min(255, p.lifetime)}, {255, 255, 255, 0}, 12);
}
int row = 0; int row = 0;
auto push_text = [&](std::string const & text) mutable auto push_text = [&](std::string const & text) mutable
{ {