diff --git a/examples/weather.cpp b/examples/weather.cpp index 53aa2078..1ee6cc01 100644 --- a/examples/weather.cpp +++ b/examples/weather.cpp @@ -6,6 +6,8 @@ #include #include #include +#include +#include #include #include @@ -16,30 +18,98 @@ struct weather_app { static constexpr int N = 128; + const float dt = 0.2f; + const float viscosity = 0.f; + const float temperature_diffusion = 0.001f; + const float cooling = 1.f / 300.f; + const float cooling_factor = std::exp(- cooling * dt); + const float heating = 323.f * (std::exp(cooling * dt) - 1.f) / dt; + const float coriolis = 0.f; + const float coriolis_bands = 2.f; + const float friction = 0.f; + // const float friction_factor = 1.f - std::exp(- friction * dt); + + const bool periodic_x = false; + + float expected_temperature_at(int y) + { + // float latitude = (y - N * 0.5f) * 2.f / N; + // return std::cos(latitude * float(math::pi)); + + return temperature_income_at(y) * dt / (std::exp(cooling * dt) - 1.f); + } + + float temperature_income_at(int y) + { + 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.6f, 1.f, 1.f - std::abs(latitude)); + } + weather_app(options const &, context const &) { simulation_box_ = {{{0.f, N}, {0.f, N}}}; - obstacle_.resize({N, N}, 0.f); + terrain_.resize({N, N}, 0.f); velocity_.resize({N, N}); new_velocity_.resize({N, N}); pressure_.resize({N, N}, 0.f); temperature_.resize({N, N}, -1.f); new_temperature_.resize({N, N}); + average_temperature_.resize({N, N}, 0.f); random::generator rng{random::device{}}; // random::generator rng{0, 0}; random::uniform_ball_vector_distribution random_velocity{}; - for (auto & v : velocity_) - v = random_velocity(rng); + // 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) { - // velocity_(x, y)[0] += (y*1.f/N-0.5f); - temperature_(x, y) = 0.5f - std::abs(y*2.f/N - 1.f); + 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}; + temperature_(x, y) = expected_temperature_at(y); + } + } + + std::vector> octaves; + std::vector weights; + + random::uniform_sphere_vector_distribution random_vector{}; + + for (int octave = 0; octave < 8; ++octave) + { + int size = 4 << octave; + util::ndarray, 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> noise(std::move(octaves), std::move(weights)); + + for (int y = 0; y < N; ++y) + { + for (int x = 0; x < N; ++x) + { + auto d = math::length(math::vector{x - N / 2.f, y - N / 2.f}) / (N / 2.f); + (void)d; + float value = noise((x + 0.5f) / N, (y + 0.5f) / N); + value = pow(value, 4.f) - d / 4.f; + terrain_(x, y) = std::max(0.f, math::lerp(1.f, 16.f, value)); } } } @@ -55,8 +125,17 @@ struct weather_app if (event.down && event.key == app::keycode::T) show_temperature_ ^= true; + if (event.down && event.key == app::keycode::D) + show_temperature_delta_ ^= true; + + if (event.down && event.key == app::keycode::A) + show_average_temperature_delta_ ^= true; + if (event.down && event.key == app::keycode::P) show_pressure_ ^= true; + + if (event.down && event.key == app::keycode::H) + show_land_ ^= true; } void update() override @@ -64,21 +143,14 @@ struct weather_app if (paused_) return; - float const dt = 0.2f; - float const viscosity = 0.f; - float const temperature_diffusion = 0.01f; - float const heating = 0.001f; - float const coriolis = 0.1f; - float const coriolis_bands = 5.f; - float const friction = 0.1f; - float const friction_factor = 1.f - std::exp(- friction * dt); - // Temperature source for (int y = 0; y < N; ++y) { for (int x = 0; x < N; ++x) { - temperature_(x, y) += dt * heating * (0.5f - std::abs(y*2.f/N - 1.f)); + // 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) *= cooling_factor; } } @@ -87,19 +159,39 @@ struct weather_app 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) + { + new_temperature_(i, 0) = temperature_(i, 0); + new_temperature_(i, N - 1) = temperature_(i, N - 1); + + if (!periodic_x) + { + new_temperature_(0, i) = temperature_(0, i); + new_temperature_(N - 1, i) = temperature_(N - 1, i); + } + } for (int y = 1; y < N - 1; ++y) { - for (int x = 0; x < N; ++x) + for (int x = xmin; x < xmax; ++x) { auto v = velocity_(x, y); auto p = math::point{x + 0.5f, y + 0.5f} - dt * v; - p[0] = math::clamp(p[0] - 0.5f, {0.f, N - 1.f}); + p[0] = p[0] - 0.5f; p[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f}); + if (!periodic_x) + p[0] = math::clamp(p[0], {0.f, N - 1.f}); + int ix = std::floor(p[0]); int iy = std::min(N - 1, std::floor(p[1])); + if (!periodic_x) + ix = std::min(N - 1, ix); + float tx = p[0] - ix; float ty = p[1] - iy; @@ -125,7 +217,7 @@ struct weather_app new_velocity_(x, y) = velocity_(x, y); for (int y = 1; y < N - 1; ++y) { - for (int x = 0; x < N; ++x) + for (int x = xmin; x < xmax; ++x) { // Velocity Laplacian auto laplacian = velocity_(wrap(x + 1), y) + velocity_(wrap(x - 1), y) + velocity_(x, y + 1) + velocity_(x, y - 1) - 4.f * velocity_(x, y); @@ -140,7 +232,7 @@ struct weather_app new_temperature_(x, y) = temperature_(x, y); for (int y = 1; y < N - 1; ++y) { - for (int x = 0; x < N; ++x) + for (int x = xmin; x < xmax; ++x) { // Temperature Laplacian auto laplacian = temperature_(wrap(x + 1), y) + temperature_(wrap(x - 1), y) + temperature_(x, y + 1) + temperature_(x, y - 1) - 4.f * temperature_(x, y); @@ -152,15 +244,17 @@ struct weather_app // Apply forces & friction for (int y = 1; y < N - 1; ++y) { - for (int x = 0; x < N; ++x) + for (int x = xmin; x < xmax; ++x) { float latitude = (N * 0.5f - y) * 2.f / N; // 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)); - bool is_land = math::length(math::vector{x - N / 2.f, y - N / 4.f}) < N / 8.f; - if (is_land && false) - velocity_(x, y) -= friction_factor * velocity_(x, y);// * math::length(velocity_(x, y)); + float local_friction = friction * terrain_(x, y); + // velocity_(x, y) -= local_friction * velocity_(x, y) * math::length(velocity_(x, y)); + float local_friction_factor = std::exp(- local_friction * dt); + velocity_(x, y) *= local_friction_factor; } } @@ -169,7 +263,7 @@ struct weather_app { for (int y = 1; y < N - 1; ++y) { - for (int x = 0; x < N; ++x) + for (int x = xmin; x < xmax; ++x) { // Velocity divergence float divergence = (velocity_(wrap(x + 1), y)[0] - velocity_(wrap(x - 1), y)[0] + velocity_(x, y + 1)[1] - velocity_(x, y - 1)[1]) / 2.f; @@ -183,15 +277,21 @@ struct weather_app // Apply boundary conditions for pressure for (int i = 0; i < N; ++i) { - // pressure_(0, i) = pressure_(1, i); - // pressure_(N - 1, i) = pressure_(N - 2, i); + if (!periodic_x) + { + pressure_(0, i) = pressure_(1, i); + pressure_(N - 1, i) = pressure_(N - 2, i); + } pressure_(i, 0) = pressure_(i, 1); pressure_(i, N - 1) = pressure_(i, N - 2); } - // pressure_(0, 0) = (pressure_(0, 1) + pressure_(1, 0)) / 2.f; - // pressure_(N-1, 0) = (pressure_(N-1, 1) + pressure_(N-2, 0)) / 2.f; - // pressure_(0, N-1) = (pressure_(0, N-2) + pressure_(1, N-2)) / 2.f; - // pressure_(N-1, N-1) = (pressure_(N-1, N-2) + pressure_(N-2, N-1)) / 2.f; + if (!periodic_x) + { + pressure_(0, 0) = (pressure_(0, 1) + pressure_(1, 0)) / 2.f; + pressure_(N-1, 0) = (pressure_(N-1, 1) + pressure_(N-2, 0)) / 2.f; + pressure_(0, N-1) = (pressure_(0, N-2) + pressure_(1, N-2)) / 2.f; + pressure_(N-1, N-1) = (pressure_(N-1, N-2) + pressure_(N-2, N-1)) / 2.f; + } // Normalize pressure float average_pressure = 0.f; @@ -205,7 +305,7 @@ struct weather_app // by subtracting pressure gradient for (int y = 1; y < N - 1; ++y) { - for (int x = 0; x < N; ++x) + for (int x = xmin; x < xmax; ++x) { // Pressure gradient math::vector gradient{pressure_(wrap(x + 1), y) - pressure_(wrap(x - 1), y), pressure_(x, y + 1) - pressure_(x, y - 1)}; @@ -216,11 +316,25 @@ struct weather_app // Apply boundary conditions for velocity for (int i = 0; i < N; ++i) { - // velocity_(0, i)[0] = -velocity_(1, i)[0]; - // velocity_(N-1, i)[0] = -velocity_(N-2, i)[0]; + if (!periodic_x) + { + 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, N-2)[1] = -velocity_(i, N-2)[1]; } + + ++frame_; + + // Update all-time average temperature + for (int y = 0; y < N; ++y) + { + for (int x = 0; x < N; ++x) + { + average_temperature_(x, y) = math::lerp(average_temperature_(x, y), temperature_(x, y), 1.f / frame_); + } + } } void present() override @@ -251,28 +365,28 @@ struct weather_app { for (int x = 0; x < N; ++x) { - if (show_temperature_ || show_pressure_) - { - float value = 0.f; - if (show_temperature_) - value = 2.f * temperature_(x, y); - else if (show_pressure_) - value = 1000.f * pressure_(x, y); + gfx::color_4f color = gfx::color_4f::zero(); - gfx::color_4f color; - if (value > 0.f) - { - value = 1.f - std::exp(- value); - color = {value, value * 0.5f, value * 0.125f, 1.f}; - } - else - { - value = 1.f - std::exp(value); - color = {value * 0.125f, value * 0.5f, value, 1.f}; - } + 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))); + }; - painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color)); - } + 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})); + + if (show_temperature_) + color = gfx::blend(color, map_color(0.1f * (temperature_(x, y) - 273.f), {0.125f, 0.5f, 1.f, 0.75f}, {1.f, 0.5f, 0.125f, 0.75f})); + + 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})); + + 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})); + + 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})); + + painter_.rect({{{x, x + 1.f}, {y, y + 1.f}}}, gfx::to_coloru8(color)); if (show_velocity_) { @@ -297,6 +411,9 @@ struct weather_app ++row; }; + push_text(std::format("Frame {}", frame_)); + push_text(std::format("Day {:.2f}", frame_ / (720.f / dt))); + if (mouseover_cell) { int x = (*mouseover_cell)[0]; @@ -304,9 +421,12 @@ struct weather_app 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 = {} {}", velocity_(x, y)[0], velocity_(x, y)[1])); - push_text(std::format("P = {}", pressure_(x, y))); - push_text(std::format("T = {}", temperature_(x, y))); + push_text(std::format("V = {:.3f} {:.3f}", velocity_(x, y)[0] * 1000.f, velocity_(x, y)[1] * 1000.f)); + push_text(std::format("P = {:.3f}", pressure_(x, y))); + 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("E = {:.3f}", expected_temperature_at(y) - 273.f)); + push_text(std::format("H = {:.3f}", terrain_(x, y))); } painter_.render(math::orthographic_camera{view_box}.transform()); @@ -319,15 +439,21 @@ private: bool paused_ = true; bool show_velocity_ = true; - bool show_temperature_ = true; + bool show_temperature_ = false; + bool show_temperature_delta_ = false; + bool show_average_temperature_delta_ = false; bool show_pressure_ = false; + bool show_land_ = true; - util::ndarray obstacle_; + util::ndarray terrain_; util::ndarray, 2> velocity_; util::ndarray, 2> new_velocity_; util::ndarray pressure_; util::ndarray temperature_; util::ndarray new_temperature_; + util::ndarray average_temperature_; + + int frame_ = 0; }; namespace psemek::app