Weather simulation wip

examples/weather.cpp

@@ -43,8 +43,8 @@ auto make_perlin(random::generator & rng, int min_octave, int max_octave, float
 
 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_1 = make_perlin(rng, 3, 6, 2.f);
-	auto noise_2 = make_perlin(rng, 2, 6, 2.f);
+	auto noise_2 = make_perlin(rng, 3, 6, 2.f);
 
 	for (int y = 0; y < result.height(); ++y)
 	{
@@ -61,7 +61,7 @@ struct weather_app
 {
 	static constexpr int N = 128;
 
-	const float dt = 4.f;
+	const float dt = 20.f;
 	const float viscosity = 0.f;
 	const float advection_magnification = 1.f;
 	const float temperature_diffusion = 0.001f;
@@ -71,13 +71,13 @@ struct weather_app
 	const float coriolis = 0.f;
 	const float coriolis_bands = 2.f;
 	const float friction = 0.f;
-	const float slope_force = 0.04f;
+	const float slope_force = 4.f;
 	const float vorticity_confinement = 0.f;
-	const float evaporation = 0.1f;
+	const float evaporation = 0.01f;
-	const float max_humidity_factor = 0.001f;
+	const float max_humidity_factor = 100.f;
 	const float precipitation_factor = 0.0003f;
 	const float force_field_amplitude = 0.00005f;
-	const float random_forces = 0.f;
+	const float random_forces = 1.f;
 	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);
@@ -89,7 +89,7 @@ struct weather_app
 
 	gfx::pixmap_rgba biomes_map;
 
-	float expected_temperature_at(int y)
+	float expected_temperature_at(int y) const
 	{
 		// float latitude = (y - N * 0.5f) * 2.f / N;
 		// return std::cos(latitude * float(math::pi));
@@ -97,7 +97,7 @@ struct weather_app
 		return temperature_income_at(y) * dt / (std::exp(cooling * dt) - 1.f);
 	}
 
-	float temperature_income_at(int y)
+	float temperature_income_at(int y) const
 	{
 		// float latitude = (y - N * 0.5f) * 2.f / N;
 		float latitude = y * 1.f / N;
@@ -105,6 +105,11 @@ struct weather_app
 		return heating * math::lerp(0.8f, 1.f, 1.f - std::abs(latitude));
 	}
 
+	int wrap(int i) const
+	{
+		return (i + N) % N;
+	}
+
 	weather_app(options const &, context const &)
 	{
 		simulation_box_ = {{{0.f, N}, {0.f, N}}};
@@ -141,6 +146,9 @@ struct weather_app
 				float latitude = (N * 0.5f - y) * 2.f / N;
 				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);
+
+				float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor;
+				humidity_(x, y) = max_humidity + dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f) * (1.f - precipitation_factor * dt) / precipitation_factor / dt;
 			}
 		}
 
@@ -214,6 +222,9 @@ struct weather_app
 		}
 		float const force_field_t = ((frame_ % force_field_switch_frames) + 0.5f) / force_field_switch_frames;
 
+		int xmin = periodic_x ? 0 : 1;
+		int xmax = periodic_x ? N : N - 1; // exclusive
+
 		// Temperature source
 		for (int y = 0; y < N; ++y)
 		{
@@ -234,22 +245,14 @@ struct weather_app
 					humidity_(x, y) += dt * evaporation * std::max(0.f, temperature_(x, y) - 273.f);
 
 				// float discharge = std::min(humidity_(x, y), precipitation_factor * dt);
-				float discharge = humidity_(x, y) * precipitation_factor * dt;
+				// float discharge = humidity_(x, y) * precipitation_factor * dt;
-				// float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor;
+				float max_humidity = std::max(0.f, temperature_(x, y) - 223.f) * max_humidity_factor;
-				// float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt;
+				float discharge = std::max(0.f, humidity_(x, y) - max_humidity) * precipitation_factor * dt;
 				humidity_(x, y) -= discharge;
 				precipitation_(x, y) = discharge;
 			}
 		}
 
-		auto wrap = [](int i)
-		{
-			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)
 		{
@@ -366,10 +369,14 @@ struct weather_app
 					(std::max(0.f, terrain_(x, y + 1)) - std::max(0.f, 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));
+				// [[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) *= std::min(1.f, slope_factor);
 				// velocity_(x, y) -= terrain_gradient * slope_force * dt;
 
+				[[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(math::length(terrain_gradient), 4.f));
+				// [[maybe_unused]] float slope_factor = std::exp(- dt * slope_force * std::pow(std::max(0.f, terrain_(x, y)), 1.f));
+				velocity_(x, y) *= slope_factor;
+
 				// 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;
@@ -528,6 +535,68 @@ struct weather_app
 				mouseover_cell = {x, y};
 		}
 
+		if (mouseover_cell && (state().mouse_button_down.contains(app::mouse_button::left) || state().mouse_button_down.contains(app::mouse_button::right)))
+		{
+			float delta = state().mouse_button_down.contains(app::mouse_button::left) ? 1.f : -1.f;
+
+			int R = 4;
+
+			for (int iy = -R; iy <= R; ++iy)
+			{
+				for (int ix = -R; ix <= R; ++ix)
+				{
+					int x = (*mouseover_cell)[0] + ix;
+					int y = (*mouseover_cell)[1] + iy;
+
+					if (x >= 0 && x < N && y >= 0 && y < N)
+					{
+						auto d = math::vector<float, 2>{ix, iy} / 2.f;
+						terrain_(x, y) += delta * 0.05f * std::exp(- math::dot(d, d));
+					}
+				}
+			}
+		}
+
+		if (mouseover_cell && state().mouse_button_down.contains(app::mouse_button::middle))
+		{
+			int R = 4;
+
+			float average = 0.f;
+			int count = 0;
+
+			for (int iy = -R; iy <= R; ++iy)
+			{
+				for (int ix = -R; ix <= R; ++ix)
+				{
+					int x = (*mouseover_cell)[0] + ix;
+					int y = (*mouseover_cell)[1] + iy;
+
+					if (x >= 0 && x < N && y >= 0 && y < N)
+					{
+						average += terrain_(x, y);
+						count += 1;
+					}
+				}
+			}
+
+			average /= count;
+
+			for (int iy = -R; iy <= R; ++iy)
+			{
+				for (int ix = -R; ix <= R; ++ix)
+				{
+					int x = (*mouseover_cell)[0] + ix;
+					int y = (*mouseover_cell)[1] + iy;
+
+					if (x >= 0 && x < N && y >= 0 && y < N)
+					{
+						auto d = math::vector<float, 2>{ix, iy} / 2.f;
+						terrain_(x, y) += (average - terrain_(x, y)) * 0.05f * std::exp(- math::dot(d, d));
+					}
+				}
+			}
+		}
+
 		[[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){
@@ -604,10 +673,34 @@ struct weather_app
 					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}));
 
 				if (show_water_vapor_)
-					color = gfx::blend(color, map_color(humidity_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
+					color = gfx::blend(color, map_color(humidity_(x, y) * 0.001f, {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
 
 				if (show_precipitation_)
-					color = gfx::blend(color, map_color(precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));
+				{
+					// color = gfx::blend(color, map_color(precipitation_(x, y), {1.f, 1.f, 1.f, -1.f}, {1.f, 1.f, 1.f, 1.f}));
+
+					float alpha = 2.f / (1.f + std::exp(- precipitation_(x, y))) - 1.f;
+
+					gfx::color_4f cloud_color{1.f, 1.f, 1.f, alpha * 0.875f};
+
+					if (y > 0 && y + 1 < N)
+					{
+						if (periodic_x || (x > 0 && x + 1 < N))
+						{
+							math::vector gradient
+							{
+								(precipitation_(wrap(x + 1), y) - precipitation_(wrap(x - 1), y)) / 2.f,
+								(precipitation_(x, y + 1) - precipitation_(x, y - 1)) / 2.f,
+							};
+
+							auto normal = math::normalized(math::vector{-gradient[0], -gradient[1], 2.f});
+							auto lightness = 0.5f + 0.5f * math::dot(normal, math::normalized(math::vector{1.f, 2.f, 3.f}));
+							cloud_color = gfx::dark(cloud_color, 1.f - lightness);
+						}
+					}
+
+					color = gfx::blend(color, cloud_color);
+				}
 
 				if (show_average_precipitation_)
 					color = gfx::blend(color, map_color(average_precipitation_(x, y), {0.f, 1.f, 1.f, -0.75f}, {0.f, 1.f, 1.f, 0.75f}));