Weather simulation test wip

examples/weather.cpp

@@ -6,6 +6,8 @@
 #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/log/log.hpp>
 
@@ -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<float, 2> random_velocity{};
 
-		for (auto & v : velocity_)
+		// for (auto & v : velocity_)
-			v = random_velocity(rng);
+		// {
+		// 	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);
+				float latitude = (N * 0.5f - y) * 2.f / N;
-				temperature_(x, y) = 0.5f - std::abs(y*2.f/N - 1.f);
+				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<pcg::perlin<float, 2>> octaves;
+		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 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<int>(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;
+				float local_friction = friction * terrain_(x, y);
-				if (is_land && false)
+				// velocity_(x, y) -= local_friction * velocity_(x, y) * math::length(velocity_(x, y));
-					velocity_(x, y) -= friction_factor * 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);
+			if (!periodic_x)
-			// pressure_(N - 1, i) = pressure_(N - 2, i);
+			{
+				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;
+		if (!periodic_x)
-		// 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_(0, 0) = (pressure_(0, 1) + pressure_(1, 0)) / 2.f;
-		// pressure_(N-1, N-1) = (pressure_(N-1, N-2) + pressure_(N-2, N-1)) / 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];
+			if (!periodic_x)
-			// velocity_(N-1, i)[0] = -velocity_(N-2, i)[0];
+			{
+				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_)
+				gfx::color_4f color = gfx::color_4f::zero();
-				{
-					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;
+				auto map_color = [](float value, gfx::color_4f const & negative, gfx::color_4f const & positive){
-					if (value > 0.f)
+					return math::lerp(negative, positive, 1.f/ (1.f + std::exp(- value)));
-					{
+				};
-						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};
-					}
 
-					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("V = {:.3f} {:.3f}", velocity_(x, y)[0] * 1000.f, velocity_(x, y)[1] * 1000.f));
-			push_text(std::format("P = {}", pressure_(x, y)));
+			push_text(std::format("P = {:.3f}", pressure_(x, y)));
-			push_text(std::format("T = {}", temperature_(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<float, 2> obstacle_;
+	util::ndarray<float, 2> terrain_;
 	util::ndarray<math::vector<float, 2>, 2> velocity_;
 	util::ndarray<math::vector<float, 2>, 2> new_velocity_;
 	util::ndarray<float, 2> pressure_;
 	util::ndarray<float, 2> temperature_;
 	util::ndarray<float, 2> new_temperature_;
+	util::ndarray<float, 2> average_temperature_;
+
+	int frame_ = 0;
 };
 
 namespace psemek::app