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

Weather simulation force field experiments

Browse source
This commit is contained in:
Nikita Lisitsa 2025-12-10 17:32:29 +03:00
parent 9e4a96f622
commit 69a2a04811
1 changed files with 168 additions and 44 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

212
examples/weather.cpp
View file

@ -13,23 +13,73 @@
using namespace psemek; using namespace psemek;
auto make_perlin(random::generator & rng, int min_octave, int max_octave, float power)
{
std::vector<pcg::perlin<float, 2>> octaves;
std::vector<float> weights;
random::uniform_sphere_vector_distribution<float, 2> random_vector{};
for (int octave = min_octave; octave < max_octave; ++octave)
{
int size = 1 << octave;
util::ndarray<math::vector<float, 2>, 2> gradients({size + 1, size + 1});
for (auto & g : gradients)
g = random_vector(rng);
octaves.emplace_back(std::move(gradients));
weights.push_back(std::pow(power, - octave));
}
float weight_sum = 0.f;
for (auto w : weights)
weight_sum += w;
for (auto & w : weights)
w /= weight_sum;
return pcg::fractal<pcg::perlin<float, 2>>(std::move(octaves), std::move(weights));
}
void make_force_field(random::generator & rng, util::ndarray<math::vector<float, 2>, 2> & result, float scale)
{
auto noise_1 = make_perlin(rng, 2, 6, 2.f);
auto noise_2 = make_perlin(rng, 2, 6, 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) = scale * (math::vector{noise_1(p), noise_2(p)} * 2.f - math::vector{1.f, 1.f});
}
}
}
struct weather_app struct weather_app
: app::application_base : app::application_base
{ {
static constexpr int N = 128; static constexpr int N = 128;
const float dt = 0.2f; const float dt = 4.f;
const float viscosity = 0.f; const float viscosity = 0.f;
const float advection_magnification = 1.f;
const float temperature_diffusion = 0.001f; const float temperature_diffusion = 0.001f;
const float cooling = 1.f / 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 coriolis = 0.01f;
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 vorticity_confinement = 0.f;
const float force_field_amplitude = 0.00005f;
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); // const float friction_factor = 1.f - std::exp(- friction * dt);
const bool periodic_x = false; const bool periodic_x = true;
random::generator rng{random::device{}};
// random::generator rng{0, 0};
float expected_temperature_at(int y) float expected_temperature_at(int y)
{ {
@ -43,7 +93,7 @@ struct weather_app
{ {
float latitude = (y - N * 0.5f) * 2.f / N; float latitude = (y - N * 0.5f) * 2.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.6f, 1.f, 1.f - std::abs(latitude)); return heating * math::lerp(0.65f, 1.f, 1.f - std::abs(latitude));
} }
weather_app(options const &, context const &) weather_app(options const &, context const &)
@ -57,8 +107,11 @@ struct weather_app
temperature_.resize({N, N}, -1.f); temperature_.resize({N, N}, -1.f);
new_temperature_.resize({N, N}); new_temperature_.resize({N, N});
average_temperature_.resize({N, N}, 0.f); average_temperature_.resize({N, N}, 0.f);
force_field_main_.resize({N, N});
force_field_current_.resize({N, N});
force_field_next_.resize({N, N});
vorticity_.resize({N, N});
random::generator rng{random::device{}};
// 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{};
@ -73,33 +126,12 @@ struct weather_app
for (int x = 0; x < N; ++x) for (int x = 0; x < N; ++x)
{ {
float latitude = (N * 0.5f - y) * 2.f / N; float latitude = (N * 0.5f - y) * 2.f / N;
velocity_(x, y) = random_velocity(rng) * 1.1f + 0.001f * math::vector{-std::cos(0.5f * float(math::pi) * latitude * coriolis_bands), 0.f}; 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); temperature_(x, y) = expected_temperature_at(y);
} }
} }
std::vector<pcg::perlin<float, 2>> octaves; auto terrain_noise = make_perlin(rng, 2, 10, 1.6f);
std::vector<float> weights;
random::uniform_sphere_vector_distribution<float, 2> random_vector{};
for (int octave = 0; octave < 8; ++octave)
{
int size = 4 << octave;
util::ndarray<math::vector<float, 2>, 2> gradients({size + 1, size + 1});
for (auto & g : gradients)
g = random_vector(rng);
octaves.emplace_back(std::move(gradients));
weights.push_back(std::pow(2.f, - 0.75f * octave));
}
float weight_sum = 0.f;
for (auto w : weights)
weight_sum += w;
for (auto & w : weights)
w /= weight_sum;
pcg::fractal<pcg::perlin<float, 2>> noise(std::move(octaves), std::move(weights));
for (int y = 0; y < N; ++y) for (int y = 0; y < N; ++y)
{ {
@ -107,11 +139,14 @@ struct weather_app
{ {
auto d = math::length(math::vector{x - N / 2.f, y - N / 2.f}) / (N / 2.f); auto d = math::length(math::vector{x - N / 2.f, y - N / 2.f}) / (N / 2.f);
(void)d; (void)d;
float value = 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)); terrain_(x, y) = std::max(0.f, math::lerp(1.f, 16.f, value));
} }
} }
make_force_field(rng, force_field_main_, 0.5f);
make_force_field(rng, force_field_next_, 0.5f);
} }
void on_event(app::key_event const & event) override void on_event(app::key_event const & event) override
@ -143,6 +178,14 @@ struct weather_app
if (paused_) if (paused_)
return; return;
// Update force field
if ((frame_ % force_field_switch_frames) == 0)
{
std::swap(force_field_current_, force_field_next_);
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;
// Temperature source // Temperature source
for (int y = 0; y < N; ++y) for (int y = 0; y < N; ++y)
{ {
@ -179,7 +222,7 @@ struct weather_app
for (int x = xmin; x < xmax; ++x) for (int x = xmin; x < xmax; ++x)
{ {
auto v = velocity_(x, y); auto v = velocity_(x, y);
auto p = math::point{x + 0.5f, y + 0.5f} - dt * v; auto p = math::point{x + 0.5f, y + 0.5f} - (advection_magnification * dt) * v;
p[0] = p[0] - 0.5f; p[0] = p[0] - 0.5f;
p[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f}); p[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f});
@ -241,6 +284,11 @@ struct weather_app
} }
std::swap(temperature_, new_temperature_); std::swap(temperature_, new_temperature_);
// Compute vorticity
for (int y = 1; y < N - 1; ++y)
for (int x = xmin; x < xmax; ++x)
vorticity_(x, y) = (velocity_(x, y + 1)[0] - velocity_(x, y - 1)[0]) / 2.f - (velocity_(wrap(x + 1), y)[1] - velocity_(wrap(x - 1), y)[1]) / 2;
// Apply forces & friction // Apply forces & friction
for (int y = 1; y < N - 1; ++y) for (int y = 1; y < N - 1; ++y)
{ {
@ -251,6 +299,33 @@ 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);
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,
};
[[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;
// 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;
math::vector vorticity_gradient
{
(vorticity_(wrap(x + 1), y) - vorticity_(wrap(x - 1), y)) / 2.f,
(vorticity_(x, y + 1) - vorticity_(x, y - 1)) / 2.f,
};
if (auto l = math::length(vorticity_gradient); l > 0.f)
vorticity_gradient /= l;
velocity_(x, y) += vorticity_confinement * dt * vorticity_(x, y) * math::ort(vorticity_gradient);
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); float local_friction_factor = std::exp(- local_friction * dt);
@ -261,7 +336,11 @@ struct weather_app
// Solve Poisson equation for pressure // Solve Poisson equation for pressure
for (int iteration = 0; iteration < 16; ++iteration) for (int iteration = 0; iteration < 16; ++iteration)
{ {
for (int y = 1; y < N - 1; ++y) int ymin = ((iteration % 2) == 0) ? 1 : N - 2;
int ymax = ((iteration % 2) == 0) ? N - 1 : 0;
int ystep = ((iteration % 2) == 0) ? 1 : -1;
for (int y = ymin; y != ymax; y += ystep)
{ {
for (int x = xmin; x < xmax; ++x) for (int x = xmin; x < xmax; ++x)
{ {
@ -308,7 +387,10 @@ struct weather_app
for (int x = xmin; x < xmax; ++x) for (int x = xmin; x < xmax; ++x)
{ {
// Pressure gradient // Pressure gradient
math::vector gradient{pressure_(wrap(x + 1), y) - pressure_(wrap(x - 1), y), pressure_(x, y + 1) - pressure_(x, y - 1)}; math::vector gradient{
(pressure_(wrap(x + 1), y) - pressure_(wrap(x - 1), y)) / 2.f,
(pressure_(x, y + 1) - pressure_(x, y - 1)) / 2.f
};
velocity_(x, y) -= gradient; velocity_(x, y) -= gradient;
} }
} }
@ -318,13 +400,39 @@ struct weather_app
{ {
if (!periodic_x) if (!periodic_x)
{ {
velocity_(0, i)[0] = -velocity_(1, i)[0]; float left_boundary_flow = 0.f;//0.01f * std::sin((i * 1.f / N) * float(math::pi) * 4.f);
velocity_(N-1, i)[0] = -velocity_(N-2, i)[0]; float right_boundary_flow = -left_boundary_flow;
velocity_(1, i)[0] = left_boundary_flow;
velocity_(N-2, i)[0] = right_boundary_flow;
velocity_(0, i)[0] = - velocity_(1, i)[0];
velocity_(N-1, i)[0] = - velocity_(N-2, i)[0];
} }
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];
} }
// Uncomment to visualize the force field
// for (int y = 0; y < N; ++y)
// for (int x = 0; x < N; ++x)
// velocity_(x, y) = 100000.f * force_field_(x, y);
// Uncomment to visualize the terrain gradient field
// for (int y = 1; y < N - 1; ++y)
// {
// for (int x = 1; x < N - 1; ++x)
// {
// math::vector terrain_gradient
// {
// (terrain_(x + 1, y) - terrain_(x - 1, y)) / 2.f,
// (terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f,
// };
// velocity_(x, y) = terrain_gradient;
// }
// }
++frame_; ++frame_;
// Update all-time average temperature // Update all-time average temperature
@ -332,7 +440,9 @@ struct weather_app
{ {
for (int x = 0; x < N; ++x) for (int x = 0; x < N; ++x)
{ {
average_temperature_(x, y) = math::lerp(average_temperature_(x, y), temperature_(x, y), 1.f / frame_); 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);
} }
} }
} }
@ -361,16 +471,16 @@ struct weather_app
[[maybe_unused]] float const pixel_size = view_box[0].length() / state().size[0]; [[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){
return math::lerp(negative, positive, 1.f/ (1.f + std::exp(- value)));
};
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();
auto map_color = [](float value, gfx::color_4f const & negative, gfx::color_4f const & positive){
return math::lerp(negative, positive, 1.f/ (1.f + std::exp(- value)));
};
if (show_land_) 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})); 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}));
@ -384,23 +494,33 @@ struct weather_app
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)), {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(1000.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}));
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));
}
}
if (show_velocity_) if (show_velocity_)
{
for (int y = 0; y < N; ++y)
{
for (int x = 0; x < N; ++x)
{ {
math::point center{x + 0.5f, y + 0.5f}; math::point center{x + 0.5f, y + 0.5f};
auto v = velocity_(x, y); auto v = velocity_(x, y);
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 = 40.f; float const magnification = 1000.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;
// color = gfx::lerp(gfx::color_4f{0.5f, 1.f, 0.f, 1.f}, gfx::color_4f{1.f, 0.f, 0.f, 1.f}, 1.f - std::exp(- 0.25f * magnification * l));
color = map_color(0.1f * (temperature_(x, y) - 273.f), {0.125f, 0.5f, 1.f, 1.f}, {1.f, 0.5f, 0.125f, 1.f});
} }
auto n = math::ort(v) * 0.3f; auto n = math::ort(v) * 0.3f;
painter_.triangle(center - v - n, center - v + n, center + v, {255, 255, 255, 255}); painter_.triangle(center - v - n, center - v + n, center + v, gfx::to_coloru8(color));
} }
} }
} }
@ -452,6 +572,10 @@ private:
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<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_next_;
util::ndarray<float, 2> vorticity_;
int frame_ = 0; int frame_ = 0;
}; };