Weather experiments v2
This commit is contained in:
parent
005008b720
commit
dd8781a1b2
1 changed files with 880 additions and 0 deletions
880
examples/weather_v2.cpp
Normal file
880
examples/weather_v2.cpp
Normal file
|
|
@ -0,0 +1,880 @@
|
||||||
|
#include <psemek/app/application_base.hpp>
|
||||||
|
#include <psemek/app/default_application_factory.hpp>
|
||||||
|
#include <psemek/gfx/painter.hpp>
|
||||||
|
#include <psemek/gfx/gl.hpp>
|
||||||
|
#include <psemek/math/camera.hpp>
|
||||||
|
#include <psemek/math/gauss.hpp>
|
||||||
|
#include <psemek/random/generator.hpp>
|
||||||
|
#include <psemek/random/device.hpp>
|
||||||
|
#include <psemek/random/uniform_ball.hpp>
|
||||||
|
#include <psemek/pcg/perlin.hpp>
|
||||||
|
#include <psemek/pcg/fractal.hpp>
|
||||||
|
#include <psemek/util/ndarray.hpp>
|
||||||
|
#include <psemek/util/clock.hpp>
|
||||||
|
#include <psemek/log/log.hpp>
|
||||||
|
#include <psemek/io/file_stream.hpp>
|
||||||
|
|
||||||
|
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_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_2 = 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::normalized(math::vector{2.f * noise_1(p) - 1.f, 2.f * noise_2(p) - 1.f}) * scale;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
void add(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src)
|
||||||
|
{
|
||||||
|
auto src_begin = src.begin();
|
||||||
|
auto src_end = src.end();
|
||||||
|
auto dst_begin = dst.begin();
|
||||||
|
for (; src_begin != src_end;)
|
||||||
|
*dst_begin++ += *src_begin++;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
void scale(util::ndarray<T, 2> & array, float factor)
|
||||||
|
{
|
||||||
|
for (auto & v : array)
|
||||||
|
v *= factor;
|
||||||
|
}
|
||||||
|
|
||||||
|
float do_lerp(float x, float y, float t)
|
||||||
|
{
|
||||||
|
return math::lerp(x, y, 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);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
void lerp(util::ndarray<T, 2> & dst, util::ndarray<T, 2> const & src0, util::ndarray<T, 2> const & src1, float t)
|
||||||
|
{
|
||||||
|
auto src0_begin = src0.begin();
|
||||||
|
auto src0_end = src0.end();
|
||||||
|
auto src1_begin = src1.begin();
|
||||||
|
auto dst_begin = dst.begin();
|
||||||
|
|
||||||
|
for (; src0_begin != src0_end;)
|
||||||
|
*dst_begin++ = do_lerp(*src0_begin++, *src1_begin++, t);
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
T sample_bilinear(util::ndarray<T, 2> const & array, math::point<float, 2> p)
|
||||||
|
{
|
||||||
|
p[0] = math::clamp(p[0] - 0.5f, {0.f, array.width() - 1.f});
|
||||||
|
p[1] = math::clamp(p[1] - 0.5f, {0.f, array.height() - 1.f});
|
||||||
|
|
||||||
|
int ix = std::min<int>(array.width() - 2, std::floor(p[0]));
|
||||||
|
int iy = std::min<int>(array.height() - 2, std::floor(p[1]));
|
||||||
|
|
||||||
|
float tx = p[0] - ix;
|
||||||
|
float ty = p[1] - iy;
|
||||||
|
|
||||||
|
return math::lerp(
|
||||||
|
math::lerp(array(ix + 0, iy + 0), array(ix + 1, iy + 0), tx),
|
||||||
|
math::lerp(array(ix + 0, iy + 1), array(ix + 1, iy + 1), tx),
|
||||||
|
ty
|
||||||
|
);
|
||||||
|
}
|
||||||
|
|
||||||
|
std::string season_str(int season)
|
||||||
|
{
|
||||||
|
switch (season)
|
||||||
|
{
|
||||||
|
case 0: return "Winter";
|
||||||
|
case 1: return "Spring";
|
||||||
|
case 2: return "Summer";
|
||||||
|
case 3: return "Autumn";
|
||||||
|
}
|
||||||
|
|
||||||
|
return "???";
|
||||||
|
}
|
||||||
|
|
||||||
|
std::string day_night_str(int day_night)
|
||||||
|
{
|
||||||
|
switch (day_night)
|
||||||
|
{
|
||||||
|
case 0: return "Night";
|
||||||
|
case 1: return "Day";
|
||||||
|
}
|
||||||
|
|
||||||
|
return "???";
|
||||||
|
}
|
||||||
|
|
||||||
|
struct climate_snapshot
|
||||||
|
{
|
||||||
|
util::ndarray<math::vector<float, 2>, 2> wind_velocity;
|
||||||
|
util::ndarray<float, 2> temperature;
|
||||||
|
util::ndarray<float, 2> humidity;
|
||||||
|
};
|
||||||
|
|
||||||
|
static constexpr int N = 128;
|
||||||
|
|
||||||
|
struct solver
|
||||||
|
{
|
||||||
|
util::ndarray<float, 2> const & terrain;
|
||||||
|
random::generator & rng;
|
||||||
|
int season;
|
||||||
|
int day_night;
|
||||||
|
|
||||||
|
climate_snapshot & snapshot;
|
||||||
|
|
||||||
|
int solve_velocity_iterations = 0;
|
||||||
|
int max_solve_velocity_iterations = N / 2;
|
||||||
|
|
||||||
|
int solve_temperature_iterations = 0;
|
||||||
|
int max_solve_temperature_iterations = N / 4;
|
||||||
|
|
||||||
|
int solve_humidity_iterations = 0;
|
||||||
|
int max_solve_humidity_iterations = 3 * N;
|
||||||
|
|
||||||
|
util::ndarray<math::vector<float, 2>, 2> initial_random_wind_velocity;
|
||||||
|
util::ndarray<float, 2> pressure;
|
||||||
|
util::ndarray<float, 2> new_temperature;
|
||||||
|
util::ndarray<float, 2> new_humidity;
|
||||||
|
|
||||||
|
int stage = 0;
|
||||||
|
bool finished = false;
|
||||||
|
|
||||||
|
solver(util::ndarray<float, 2> const & terrain, random::generator & rng, int season, int day_night, climate_snapshot & snapshot)
|
||||||
|
: terrain(terrain)
|
||||||
|
, rng(rng)
|
||||||
|
, season(season)
|
||||||
|
, day_night(day_night)
|
||||||
|
, snapshot(snapshot)
|
||||||
|
{
|
||||||
|
initial_random_wind_velocity.resize({N, N});
|
||||||
|
snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
|
||||||
|
pressure.resize({N, N}, 0.f);
|
||||||
|
snapshot.temperature.resize({N, N}, 0.f);
|
||||||
|
new_temperature.resize({N, N}, 0.f);
|
||||||
|
snapshot.humidity.resize({N, N}, 0.f);
|
||||||
|
new_humidity.resize({N, N}, 0.f);
|
||||||
|
|
||||||
|
make_random_vector_field(rng, initial_random_wind_velocity, 3, 6, 0.005f);
|
||||||
|
|
||||||
|
for (int y = 0; y < N; ++y)
|
||||||
|
for (int x = 0; x < N; ++x)
|
||||||
|
if (terrain(x, y) < 0.f)
|
||||||
|
initial_random_wind_velocity(x, y) *= 5.f;
|
||||||
|
else
|
||||||
|
initial_random_wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
|
||||||
|
|
||||||
|
snapshot.wind_velocity = initial_random_wind_velocity.copy();
|
||||||
|
|
||||||
|
for (int y = 0; y < N; ++y)
|
||||||
|
for (int x = 0; x < N; ++x)
|
||||||
|
snapshot.temperature(x, y) = expected_temperature(y);
|
||||||
|
}
|
||||||
|
|
||||||
|
float expected_temperature(int y)
|
||||||
|
{
|
||||||
|
return math::lerp(40.f, -20.f, (y + 0.5f) / N) // Base latitude temperature in spring/autumn
|
||||||
|
- 15.f * std::cos(season * math::pi * 0.5f) // Seasonal offset of +/-15C
|
||||||
|
- 5.f * std::cos(day_night * math::pi) // Day-night fluctuation of +/-5C
|
||||||
|
;
|
||||||
|
}
|
||||||
|
|
||||||
|
float expected_humidity(int x, int y)
|
||||||
|
{
|
||||||
|
if (terrain(x, y) < 0.f)
|
||||||
|
return 1.f;
|
||||||
|
else
|
||||||
|
return 0.f;
|
||||||
|
}
|
||||||
|
|
||||||
|
void step()
|
||||||
|
{
|
||||||
|
if (stage == 0)
|
||||||
|
{
|
||||||
|
for (int y = 1; y + 1 < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 1; x + 1 < N; ++x)
|
||||||
|
{
|
||||||
|
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, y + 1)[1] - initial_random_wind_velocity(x, y + 1)[1]) * 0.5f
|
||||||
|
;
|
||||||
|
|
||||||
|
// (sum - 4.f * p(x,y)) = divergence
|
||||||
|
// sum - 4.f * p = divergence
|
||||||
|
// 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;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int y = 1; y + 1 < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 1; x + 1 < N; ++x)
|
||||||
|
{
|
||||||
|
math::vector pressure_gradient
|
||||||
|
{
|
||||||
|
(pressure(x + 1, y) - pressure(x - 1, y)) * 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;
|
||||||
|
|
||||||
|
if (terrain(x, y) > 0.f)
|
||||||
|
snapshot.wind_velocity(x, y) *= std::exp(- 2.f * terrain(x, y));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
++solve_velocity_iterations;
|
||||||
|
if (solve_velocity_iterations == max_solve_velocity_iterations)
|
||||||
|
++stage;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (stage == 1)
|
||||||
|
{
|
||||||
|
for (int y = 0; y < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 0; x < N; ++x)
|
||||||
|
{
|
||||||
|
float step_time = 100.f;
|
||||||
|
|
||||||
|
math::point p{x + 0.5f, y + 0.5f};
|
||||||
|
float advection = sample_bilinear(snapshot.temperature, p - snapshot.wind_velocity(x, y) * step_time);
|
||||||
|
|
||||||
|
float heating_speed = 0.001f;
|
||||||
|
float heating_delta = (expected_temperature(y) - snapshot.temperature(x, y)) * (- std::expm1(- heating_speed * step_time));
|
||||||
|
|
||||||
|
new_temperature(x, y) = heating_delta + advection;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
std::swap(snapshot.temperature, new_temperature);
|
||||||
|
|
||||||
|
for (int y = 1; y + 1 < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 1; x + 1 < N; ++x)
|
||||||
|
{
|
||||||
|
float average = (0.f
|
||||||
|
+ snapshot.temperature(x - 1, y)
|
||||||
|
+ snapshot.temperature(x + 1, y)
|
||||||
|
+ snapshot.temperature(x, y - 1)
|
||||||
|
+ snapshot.temperature(x, y + 1)
|
||||||
|
) * 0.25f;
|
||||||
|
|
||||||
|
new_temperature(x, y) = math::lerp(snapshot.temperature(x, y), average, 0.5f);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
std::swap(snapshot.temperature, new_temperature);
|
||||||
|
|
||||||
|
++solve_temperature_iterations;
|
||||||
|
if (solve_temperature_iterations >= max_solve_temperature_iterations)
|
||||||
|
++stage;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (stage == 2)
|
||||||
|
{
|
||||||
|
for (int y = 0; y < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 0; x < N; ++x)
|
||||||
|
{
|
||||||
|
float step_time = 100.f;
|
||||||
|
|
||||||
|
math::point p{x + 0.5f, y + 0.5f};
|
||||||
|
float advection = sample_bilinear(snapshot.humidity, p - snapshot.wind_velocity(x, y) * step_time);
|
||||||
|
|
||||||
|
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) = evaporation_delta + advection;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
std::swap(snapshot.humidity, new_humidity);
|
||||||
|
|
||||||
|
for (int y = 1; y + 1 < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 1; x + 1 < N; ++x)
|
||||||
|
{
|
||||||
|
float average = (0.f
|
||||||
|
+ snapshot.humidity(x - 1, y)
|
||||||
|
+ snapshot.humidity(x + 1, y)
|
||||||
|
+ snapshot.humidity(x, y - 1)
|
||||||
|
+ snapshot.humidity(x, y + 1)
|
||||||
|
) * 0.25f;
|
||||||
|
|
||||||
|
new_humidity(x, y) = math::lerp(snapshot.humidity(x, y), average, 0.5f);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
std::swap(snapshot.humidity, new_humidity);
|
||||||
|
|
||||||
|
++solve_humidity_iterations;
|
||||||
|
if (solve_humidity_iterations >= max_solve_humidity_iterations)
|
||||||
|
{
|
||||||
|
++stage;
|
||||||
|
finished = true;
|
||||||
|
}
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
struct weather_app
|
||||||
|
: app::application_base
|
||||||
|
{
|
||||||
|
random::generator rng{random::device{}};
|
||||||
|
// random::generator rng{0, 0};
|
||||||
|
|
||||||
|
gfx::painter painter;
|
||||||
|
|
||||||
|
math::box<float, 2> simulation_box;
|
||||||
|
|
||||||
|
bool paused = true;
|
||||||
|
bool show_velocity = false;
|
||||||
|
bool show_temperature = false;
|
||||||
|
bool show_humidity = false;
|
||||||
|
bool show_biomes = true;
|
||||||
|
|
||||||
|
gfx::pixmap_rgba biomes_map;
|
||||||
|
util::ndarray<float, 2> terrain;
|
||||||
|
|
||||||
|
// 0 - Winter
|
||||||
|
// 1 - Spring
|
||||||
|
// 2 - Summer
|
||||||
|
// 3 - Autumn
|
||||||
|
int season = 0;
|
||||||
|
|
||||||
|
// 0 - Night
|
||||||
|
// 1 - Day
|
||||||
|
int day_night = 0;
|
||||||
|
|
||||||
|
climate_snapshot snapshots[4][2];
|
||||||
|
climate_snapshot average;
|
||||||
|
climate_snapshot display_snapshot;
|
||||||
|
|
||||||
|
bool need_update_display_snapshot = false;
|
||||||
|
|
||||||
|
std::optional<struct solver> solver;
|
||||||
|
|
||||||
|
float display_season = 0.f;
|
||||||
|
|
||||||
|
util::clock<> frame_clock;
|
||||||
|
|
||||||
|
weather_app(options const &, context const & context)
|
||||||
|
{
|
||||||
|
simulation_box = {{{0.f, N}, {0.f, N}}};
|
||||||
|
|
||||||
|
terrain.resize({N, N}, 0.f);
|
||||||
|
auto heightmap = gfx::read_image<std::uint8_t>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "heightmap_seed_3.png"});
|
||||||
|
for (int y = 0; y < N; ++y)
|
||||||
|
for (int x = 0; x < N; ++x)
|
||||||
|
terrain(x, y) = ((heightmap(x, y) / 255.f) * 2048.f - 512.f) / 1024.f;
|
||||||
|
|
||||||
|
biomes_map = gfx::read_image<gfx::color_rgba>(io::file_istream{std::filesystem::path{PSEMEK_EXAMPLES_DIR} / "biomes.png"});
|
||||||
|
|
||||||
|
context.vsync(false);
|
||||||
|
}
|
||||||
|
|
||||||
|
void on_event(app::key_event const & event) override
|
||||||
|
{
|
||||||
|
application_base::on_event(event);
|
||||||
|
|
||||||
|
if (event.down && event.key == app::keycode::SPACE)
|
||||||
|
paused ^= true;
|
||||||
|
|
||||||
|
if (event.down && event.key == app::keycode::V)
|
||||||
|
show_velocity ^= true;
|
||||||
|
|
||||||
|
if (event.down && event.key == app::keycode::T)
|
||||||
|
show_temperature ^= true;
|
||||||
|
|
||||||
|
if (event.down && event.key == app::keycode::W)
|
||||||
|
show_humidity ^= true;
|
||||||
|
|
||||||
|
if (event.down && event.key == app::keycode::B)
|
||||||
|
show_biomes ^= true;
|
||||||
|
}
|
||||||
|
|
||||||
|
void update() override
|
||||||
|
{
|
||||||
|
float const dt = frame_clock.restart().count();
|
||||||
|
|
||||||
|
if (solver && solver->finished)
|
||||||
|
{
|
||||||
|
solver.reset();
|
||||||
|
++day_night;
|
||||||
|
if (day_night == 2)
|
||||||
|
{
|
||||||
|
day_night = 0;
|
||||||
|
++season;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (season == 4)
|
||||||
|
{
|
||||||
|
compute_average();
|
||||||
|
need_update_display_snapshot = true;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (!solver && season < 4)
|
||||||
|
solver.emplace(terrain, rng, season, day_night, snapshots[season][day_night]);
|
||||||
|
|
||||||
|
if (!paused && solver)
|
||||||
|
solver->step();
|
||||||
|
|
||||||
|
if (state().key_down.contains(app::keycode::LEFT))
|
||||||
|
{
|
||||||
|
display_season -= 0.5f * dt;
|
||||||
|
need_update_display_snapshot = true;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (state().key_down.contains(app::keycode::RIGHT))
|
||||||
|
{
|
||||||
|
display_season += 0.5f * dt;
|
||||||
|
need_update_display_snapshot = true;
|
||||||
|
}
|
||||||
|
|
||||||
|
while (display_season < 0.f)
|
||||||
|
display_season += 1.f;
|
||||||
|
while (display_season >= 1.f)
|
||||||
|
display_season -= 1.f;
|
||||||
|
}
|
||||||
|
|
||||||
|
void compute_average()
|
||||||
|
{
|
||||||
|
average.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
|
||||||
|
average.temperature.resize({N, N}, 0.f);
|
||||||
|
average.humidity.resize({N, N}, 0.f);
|
||||||
|
|
||||||
|
for (auto const & season : snapshots)
|
||||||
|
for (auto const & s : season)
|
||||||
|
{
|
||||||
|
add(average.wind_velocity, s.wind_velocity);
|
||||||
|
add(average.temperature, s.temperature);
|
||||||
|
add(average.humidity, s.humidity);
|
||||||
|
}
|
||||||
|
|
||||||
|
scale(average.wind_velocity, 1.f / 8.f);
|
||||||
|
scale(average.temperature, 1.f / 8.f);
|
||||||
|
scale(average.humidity, 1.f / 8.f);
|
||||||
|
}
|
||||||
|
|
||||||
|
void update_display_snapshot()
|
||||||
|
{
|
||||||
|
display_snapshot.wind_velocity.resize({N, N}, math::vector{0.f, 0.f});
|
||||||
|
display_snapshot.temperature.resize({N, N}, 0.f);
|
||||||
|
display_snapshot.humidity.resize({N, N}, 0.f);
|
||||||
|
|
||||||
|
int i = std::min<int>(3, std::floor(4.f * display_season));
|
||||||
|
int j = (i + 1) % 4;
|
||||||
|
float t = (4.f * display_season) - i;
|
||||||
|
|
||||||
|
auto & s0 = snapshots[i][0];
|
||||||
|
auto & s1 = snapshots[j][0];
|
||||||
|
|
||||||
|
lerp(display_snapshot.wind_velocity, s0.wind_velocity, s1.wind_velocity, t);
|
||||||
|
lerp(display_snapshot.temperature, s0.temperature, s1.temperature, t);
|
||||||
|
lerp(display_snapshot.humidity, s0.humidity, s1.humidity, t);
|
||||||
|
}
|
||||||
|
|
||||||
|
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::Clear(gl::COLOR_BUFFER_BIT);
|
||||||
|
|
||||||
|
float const aspect_ratio = state().size[0] * 1.f / state().size[1];
|
||||||
|
math::box<float, 2> view_box = math::expand(simulation_box, 1.f);
|
||||||
|
if (view_box[0].length() / view_box[1].length() > aspect_ratio)
|
||||||
|
view_box[1] = math::expand(view_box[1], (view_box[0].length() / aspect_ratio - view_box[1].length()) / 2.f);
|
||||||
|
else
|
||||||
|
view_box[0] = math::expand(view_box[0], (view_box[1].length() * aspect_ratio - view_box[0].length()) / 2.f);
|
||||||
|
|
||||||
|
std::optional<math::vector<int, 2>> mouseover_cell;
|
||||||
|
{
|
||||||
|
auto mouse = math::lerp(view_box, math::vector{state().mouse[0] * 1.f / state().size[0], 1.f - state().mouse[1] * 1.f / state().size[1]});
|
||||||
|
int x = std::floor(mouse[0]);
|
||||||
|
int y = std::floor(mouse[1]);
|
||||||
|
|
||||||
|
if (x >= 0 && x < N && y >= 0 && y < N)
|
||||||
|
mouseover_cell = {x, y};
|
||||||
|
}
|
||||||
|
|
||||||
|
[[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)));
|
||||||
|
};
|
||||||
|
|
||||||
|
auto map_color_positive = [](float value, gfx::color_4f const & zero, gfx::color_4f const & positive){
|
||||||
|
return math::lerp(zero, positive, 2.f / (1.f + std::exp(- value)) - 1.f);
|
||||||
|
};
|
||||||
|
|
||||||
|
auto map_temperature = [&](float temperature)
|
||||||
|
{
|
||||||
|
return map_color(0.5f * std::floor(temperature / 10.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f});
|
||||||
|
};
|
||||||
|
|
||||||
|
auto map_humidity = [&](float humidity)
|
||||||
|
{
|
||||||
|
return map_color_positive(std::round(10.f * humidity), {0.f, 1.f, 1.f, 0.f}, {0.f, 1.f, 1.f, 0.75f});
|
||||||
|
};
|
||||||
|
|
||||||
|
auto map_biome = [this](float temperature, float precipitation)
|
||||||
|
{
|
||||||
|
auto x = math::clamp<int>(math::unlerp({0.f, 0.5f}, 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});
|
||||||
|
|
||||||
|
return gfx::to_colorf(biomes_map(x, y));
|
||||||
|
};
|
||||||
|
|
||||||
|
auto & avg_snapshot = (season < 4) ? snapshots[season][day_night] : average;
|
||||||
|
auto & snapshot = (season < 4) ? snapshots[season][day_night] : display_snapshot;
|
||||||
|
|
||||||
|
for (int y = 0; y < N; ++y)
|
||||||
|
{
|
||||||
|
for (int x = 0; x < N; ++x)
|
||||||
|
{
|
||||||
|
gfx::color_4f color = gfx::color_4f::zero();
|
||||||
|
|
||||||
|
bool const is_water = terrain(x, y) < 0.f;
|
||||||
|
bool const is_land = !is_water;
|
||||||
|
|
||||||
|
if (is_water)
|
||||||
|
color = map_color(8.f * terrain(x, y), {0.f, 0.f, 0.125f, 1.f}, {0.f, 1.f, 1.5f, 1.f});
|
||||||
|
else if (show_biomes)
|
||||||
|
{
|
||||||
|
float temperature = avg_snapshot.temperature(x, y) - std::max(0.f, terrain(x, y)) * 30.f;
|
||||||
|
color = map_biome(temperature, avg_snapshot.humidity(x, y));
|
||||||
|
}
|
||||||
|
else
|
||||||
|
color = gfx::color_4f{0.3f, 0.5f, 0.1f, 1.f};
|
||||||
|
|
||||||
|
if (is_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.125f});
|
||||||
|
|
||||||
|
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_temperature(snapshot.temperature(x, y)));
|
||||||
|
}
|
||||||
|
|
||||||
|
if (show_humidity)
|
||||||
|
{
|
||||||
|
color = gfx::blend(color, map_humidity(snapshot.humidity(x, y)));
|
||||||
|
}
|
||||||
|
|
||||||
|
painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
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};
|
||||||
|
auto v = snapshot.wind_velocity(x, y);
|
||||||
|
if (auto l = math::length(v); l > 0.f)
|
||||||
|
{
|
||||||
|
float const magnification = 200.f;
|
||||||
|
float const max_length = 1.5f;
|
||||||
|
v *= 0.5f * max_length * (1.f - std::exp(- magnification * l)) / l;
|
||||||
|
}
|
||||||
|
auto n = math::ort(v) * 0.3f;
|
||||||
|
|
||||||
|
auto color = gfx::color_4f{1.f, 1.f, 1.f, 1.f};
|
||||||
|
if (!show_temperature)
|
||||||
|
{
|
||||||
|
float t = 0.2f * snapshot.temperature(x, y);
|
||||||
|
if (t >= 0.f)
|
||||||
|
color = map_color_positive( t, color, {1.f, 0.5f, 0.125f, 1.f});
|
||||||
|
else
|
||||||
|
color = map_color_positive(-t, color, {0.125f, 0.5f, 1.f, 1.f});
|
||||||
|
}
|
||||||
|
|
||||||
|
painter.triangle(center - v - n, center - v + n, center + v, gfx::to_coloru8(color));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
int row = 0;
|
||||||
|
auto push_text = [&](std::string const & text) mutable
|
||||||
|
{
|
||||||
|
painter.text(view_box.corner(0.f, 1.f) - math::vector{0.f, row * pixel_size * 2.f * 12.f}, text, {.scale = {2.f * pixel_size, - 2.f * pixel_size}, .x = gfx::painter::x_align::left, .y = gfx::painter::y_align::top, .c = {255, 255, 255, 255}});
|
||||||
|
++row;
|
||||||
|
};
|
||||||
|
|
||||||
|
if (solver)
|
||||||
|
{
|
||||||
|
push_text(std::format("{} {}", season_str(season), day_night_str(day_night)));
|
||||||
|
|
||||||
|
if (solver->stage >= 0)
|
||||||
|
push_text(std::format("Velocity {}/{}", solver->solve_velocity_iterations, solver->max_solve_velocity_iterations));
|
||||||
|
else
|
||||||
|
push_text("");
|
||||||
|
|
||||||
|
if (solver->stage >= 1)
|
||||||
|
push_text(std::format("Temperature {}/{}", solver->solve_temperature_iterations, solver->max_solve_temperature_iterations));
|
||||||
|
else
|
||||||
|
push_text("");
|
||||||
|
|
||||||
|
if (solver->stage >= 2)
|
||||||
|
push_text(std::format("Humidity {}/{}", solver->solve_humidity_iterations, solver->max_solve_humidity_iterations));
|
||||||
|
else
|
||||||
|
push_text("");
|
||||||
|
|
||||||
|
push_text("");
|
||||||
|
}
|
||||||
|
|
||||||
|
if (mouseover_cell)
|
||||||
|
{
|
||||||
|
int x = (*mouseover_cell)[0];
|
||||||
|
int y = (*mouseover_cell)[1];
|
||||||
|
painter.rect({{{x, x + 1.f}, {y, y + 1.f}}}, {255, 255, 255, 127});
|
||||||
|
|
||||||
|
push_text(std::format("{} {}", x, y));
|
||||||
|
push_text(std::format("V = {:.3f}", length(snapshot.wind_velocity(x, y)) * 1000.f));
|
||||||
|
push_text(std::format("H = {:.3f}", terrain(x, y) * 1024.f));
|
||||||
|
push_text(std::format("T = {:.3f}", snapshot.temperature(x, y)));
|
||||||
|
push_text(std::format("W = {:.3f}", snapshot.humidity(x, y)));
|
||||||
|
push_text("");
|
||||||
|
}
|
||||||
|
|
||||||
|
if (season == 4)
|
||||||
|
{
|
||||||
|
++row;
|
||||||
|
|
||||||
|
float axis_width = 300.f;
|
||||||
|
float axis_height = 24.f;
|
||||||
|
|
||||||
|
float axis_left = view_box[0].min + 30.f * pixel_size;
|
||||||
|
float axis_bottom = view_box[1].max - row * pixel_size * 2.f * 12.f - axis_height * pixel_size;
|
||||||
|
|
||||||
|
auto line = [&](float x0, float y0, float x1, float y1, float w, gfx::color_rgba const & c)
|
||||||
|
{
|
||||||
|
painter.line(
|
||||||
|
{axis_left + x0 * axis_width * pixel_size, axis_bottom + y0 * axis_height * pixel_size},
|
||||||
|
{axis_left + x1 * axis_width * pixel_size, axis_bottom + y1 * axis_height * pixel_size},
|
||||||
|
pixel_size * w, c, false
|
||||||
|
);
|
||||||
|
};
|
||||||
|
|
||||||
|
line(0.f, 0.5f, 1.f, 0.5f, 2.f, {255, 255, 255, 127});
|
||||||
|
line(display_season, 0.f, display_season, 1.f, 4.f, {255, 255, 255, 255});
|
||||||
|
|
||||||
|
++row;
|
||||||
|
++row;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (mouseover_cell && season == 4)
|
||||||
|
{
|
||||||
|
int x = (*mouseover_cell)[0];
|
||||||
|
int y = (*mouseover_cell)[1];
|
||||||
|
|
||||||
|
float plot_width = 300.f;
|
||||||
|
float plot_height = 150.f;
|
||||||
|
|
||||||
|
float plot_left = view_box[0].min + 30.f * pixel_size;
|
||||||
|
float plot_bottom = view_box[1].max - row * pixel_size * 2.f * 12.f - plot_height * 1.5f * pixel_size;
|
||||||
|
|
||||||
|
auto line = [&](float x0, float y0, float x1, float y1, gfx::color_rgba const & c)
|
||||||
|
{
|
||||||
|
painter.line(
|
||||||
|
{plot_left + x0 * plot_width * pixel_size, plot_bottom + y0 * plot_height * pixel_size},
|
||||||
|
{plot_left + x1 * plot_width * pixel_size, plot_bottom + y1 * plot_height * pixel_size},
|
||||||
|
pixel_size * 2.f, c, false
|
||||||
|
);
|
||||||
|
};
|
||||||
|
|
||||||
|
auto begin_plot = [&]{
|
||||||
|
painter.rect({{{plot_left, plot_left + plot_width * pixel_size}, {plot_bottom, plot_bottom + plot_height * pixel_size}}}, {255, 255, 255, 31});
|
||||||
|
|
||||||
|
for (int x = 0; x <= 4; ++x)
|
||||||
|
line(x / 4.f, 0.f, x / 4.f, 1.f, {255, 255, 255, (x == 0 || x == 4) ? 127 : 31});
|
||||||
|
};
|
||||||
|
|
||||||
|
auto draw_plot = [&](std::vector<float> const & ys, gfx::color_rgba const & c)
|
||||||
|
{
|
||||||
|
std::vector<float> xs(5);
|
||||||
|
for (int i = 0; i <= 4; ++i)
|
||||||
|
xs[i] = i / 4.f;
|
||||||
|
|
||||||
|
auto system = math::matrix<float, 16, 16>::zero();
|
||||||
|
auto coeffs = math::vector<float, 16>::zero();
|
||||||
|
|
||||||
|
for (int i = 0; i < 4; ++i)
|
||||||
|
{
|
||||||
|
float x0 = xs[i];
|
||||||
|
system[2 * i + 0][4 * i + 0] = 1.f;
|
||||||
|
system[2 * i + 0][4 * i + 1] = x0;
|
||||||
|
system[2 * i + 0][4 * i + 2] = x0 * x0;
|
||||||
|
system[2 * i + 0][4 * i + 3] = x0 * x0 * x0;
|
||||||
|
coeffs[2 * i + 0] = ys[i];
|
||||||
|
|
||||||
|
float x1 = xs[i + 1];
|
||||||
|
system[2 * i + 1][4 * i + 0] = 1.f;
|
||||||
|
system[2 * i + 1][4 * i + 1] = x1;
|
||||||
|
system[2 * i + 1][4 * i + 2] = x1 * x1;
|
||||||
|
system[2 * i + 1][4 * i + 3] = x1 * x1 * x1;
|
||||||
|
coeffs[2 * i + 1] = ys[i + 1];
|
||||||
|
}
|
||||||
|
|
||||||
|
for (int i = 0; i < 4; ++i)
|
||||||
|
{
|
||||||
|
int j = (i + 1) % 4;
|
||||||
|
float x0 = xs[i + 1];
|
||||||
|
float x1 = xs[j];
|
||||||
|
|
||||||
|
system[8 + 2 * i + 0][4 * i + 1] = 1.f;
|
||||||
|
system[8 + 2 * i + 0][4 * i + 2] = 2.f * x0;
|
||||||
|
system[8 + 2 * i + 0][4 * i + 3] = 3.f * x0 * x0;
|
||||||
|
system[8 + 2 * i + 0][4 * j + 1] = - 1.f;
|
||||||
|
system[8 + 2 * i + 0][4 * j + 2] = - 2.f * x1;
|
||||||
|
system[8 + 2 * i + 0][4 * j + 3] = - 3.f * x1 * x1;
|
||||||
|
|
||||||
|
system[8 + 2 * i + 1][4 * i + 2] = 2.f;
|
||||||
|
system[8 + 2 * i + 1][4 * i + 3] = 6.f * x0;
|
||||||
|
system[8 + 2 * i + 1][4 * j + 2] = - 2.f;
|
||||||
|
system[8 + 2 * i + 1][4 * j + 3] = - 6.f * x1;
|
||||||
|
}
|
||||||
|
|
||||||
|
math::gauss(system, coeffs);
|
||||||
|
|
||||||
|
auto eval = [&](int i, float x)
|
||||||
|
{
|
||||||
|
return 0.f
|
||||||
|
+ coeffs[4 * i + 0]
|
||||||
|
+ coeffs[4 * i + 1] * x
|
||||||
|
+ coeffs[4 * i + 2] * x * x
|
||||||
|
+ coeffs[4 * i + 3] * x * x * x
|
||||||
|
;
|
||||||
|
};
|
||||||
|
|
||||||
|
for (int i = 0; i < 4; ++i)
|
||||||
|
{
|
||||||
|
int D = 8;
|
||||||
|
for (int j = 0; j < D; ++j)
|
||||||
|
{
|
||||||
|
float x0 = (i + (j + 0.f) / D) / 4.f;
|
||||||
|
float x1 = (i + (j + 1.f) / D) / 4.f;
|
||||||
|
line(x0, eval(i, x0), x1, eval(i, x1), c);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (i == std::min<int>(3, std::floor(4.f * display_season)))
|
||||||
|
{
|
||||||
|
float x = display_season;
|
||||||
|
float y = eval(i, x);
|
||||||
|
painter.circle({plot_left + x * plot_width * pixel_size, plot_bottom + y * plot_height * pixel_size}, 4.f * pixel_size, c);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
// Temperature
|
||||||
|
{
|
||||||
|
begin_plot();
|
||||||
|
|
||||||
|
for (int y = 0; y <= 10; ++y)
|
||||||
|
{
|
||||||
|
line(0.f, y / 10.f, 1.f, y / 10.f, {255, 255, 255, (y == 5) ? 127 : 31});
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<float> ys[2];
|
||||||
|
|
||||||
|
for (int d = 0; d < 2; ++d)
|
||||||
|
{
|
||||||
|
for (int s = 0; s <= 4; ++s)
|
||||||
|
{
|
||||||
|
float value = 0.5f + snapshots[s % 4][d].temperature(x, y) / 100.f;
|
||||||
|
ys[d].push_back(value);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
draw_plot(ys[0], gfx::color_rgba{0, 127, 255, 255});
|
||||||
|
draw_plot(ys[1], gfx::color_rgba{255, 127, 0, 255});
|
||||||
|
}
|
||||||
|
|
||||||
|
plot_bottom -= (plot_height + 50.f) * pixel_size;
|
||||||
|
|
||||||
|
// Humidity
|
||||||
|
{
|
||||||
|
begin_plot();
|
||||||
|
|
||||||
|
for (int y = 0; y <= 5; ++y)
|
||||||
|
{
|
||||||
|
line(0.f, y / 5.f, 1.f, y / 5.f, {255, 255, 255, (y == 0) ? 127 : 31});
|
||||||
|
}
|
||||||
|
|
||||||
|
std::vector<float> ys;
|
||||||
|
|
||||||
|
for (int s = 0; s <= 4; ++s)
|
||||||
|
{
|
||||||
|
float value = (snapshots[s % 4][0].humidity(x, y) + snapshots[s % 4][1].humidity(x, y)) / 2.f;
|
||||||
|
ys.push_back(value);
|
||||||
|
}
|
||||||
|
|
||||||
|
draw_plot(ys, {0, 255, 255, 255});
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
painter.render(math::orthographic_camera{view_box}.transform());
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
namespace psemek::app
|
||||||
|
{
|
||||||
|
|
||||||
|
std::unique_ptr<application::factory> make_application_factory()
|
||||||
|
{
|
||||||
|
return default_application_factory<weather_app>({.name = "Weather simulation test"});
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
Loading…
Add table
Reference in a new issue