diff --git a/examples/weather.cpp b/examples/weather.cpp index 1ee6cc01..b26c487a 100644 --- a/examples/weather.cpp +++ b/examples/weather.cpp @@ -13,23 +13,73 @@ using namespace psemek; +auto make_perlin(random::generator & rng, int min_octave, int max_octave, float power) +{ + std::vector> octaves; + std::vector weights; + + random::uniform_sphere_vector_distribution random_vector{}; + + for (int octave = min_octave; octave < max_octave; ++octave) + { + int size = 1 << 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(power, - octave)); + } + + float weight_sum = 0.f; + for (auto w : weights) + weight_sum += w; + for (auto & w : weights) + w /= weight_sum; + + return pcg::fractal>(std::move(octaves), std::move(weights)); +} + +void make_force_field(random::generator & rng, util::ndarray, 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 : app::application_base { static constexpr int N = 128; - const float dt = 0.2f; + const float dt = 4.f; const float viscosity = 0.f; + const float advection_magnification = 1.f; 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 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 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 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) { @@ -43,7 +93,7 @@ struct weather_app { 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)); + return heating * math::lerp(0.65f, 1.f, 1.f - std::abs(latitude)); } weather_app(options const &, context const &) @@ -57,8 +107,11 @@ struct weather_app temperature_.resize({N, N}, -1.f); new_temperature_.resize({N, N}); 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::uniform_ball_vector_distribution random_velocity{}; @@ -73,33 +126,12 @@ struct weather_app for (int x = 0; x < N; ++x) { 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); } } - 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)); + auto terrain_noise = make_perlin(rng, 2, 10, 1.6f); 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); (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; 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 @@ -143,6 +178,14 @@ struct weather_app if (paused_) 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 for (int y = 0; y < N; ++y) { @@ -179,7 +222,7 @@ struct weather_app for (int x = xmin; x < xmax; ++x) { 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[1] = math::clamp(p[1] - 0.5f, {0.f, N - 1.f}); @@ -241,6 +284,11 @@ struct weather_app } 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 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::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); // velocity_(x, y) -= local_friction * velocity_(x, y) * math::length(velocity_(x, y)); float local_friction_factor = std::exp(- local_friction * dt); @@ -261,7 +336,11 @@ struct weather_app // Solve Poisson equation for pressure 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) { @@ -308,7 +387,10 @@ struct weather_app 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)}; + 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; } } @@ -318,13 +400,39 @@ struct weather_app { if (!periodic_x) { - velocity_(0, i)[0] = -velocity_(1, i)[0]; - velocity_(N-1, i)[0] = -velocity_(N-2, i)[0]; + 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; + + 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, 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_; // Update all-time average temperature @@ -332,7 +440,9 @@ struct weather_app { 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]; + 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 x = 0; x < N; ++x) { 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_) 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})); 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)); + } + } - 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}; auto v = velocity_(x, y); + auto color = gfx::color_4f::zero(); 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; 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; - 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 temperature_; util::ndarray new_temperature_; util::ndarray average_temperature_; + util::ndarray, 2> force_field_main_; + util::ndarray, 2> force_field_current_; + util::ndarray, 2> force_field_next_; + util::ndarray vorticity_; int frame_ = 0; };