Weather simulation experiments 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, 3, 6, 2.f);
+	auto noise_1 = make_perlin(rng, 4, 6, 2.f);
-	auto noise_2 = make_perlin(rng, 3, 6, 2.f);
+	auto noise_2 = make_perlin(rng, 4, 6, 2.f);
 
 	for (int y = 0; y < result.height(); ++y)
 	{
@@ -61,31 +61,34 @@ struct weather_app
 {
 	static constexpr int N = 128;
 
+	const bool static_mode = true;
+
 	const float dt = 20.f;
-	const float viscosity = 0.f;
+	const float viscosity = static_mode ? 0.005f : 0.f;
 	const float advection_magnification = 1.f;
 	const float temperature_diffusion = 0.001f;
-	const float cooling = 0.1f / 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_bands = 2.f;
 	const float friction = 0.f;
-	const float slope_force = 4.f;
+	const float slope_force = 0.05f;
 	const float vorticity_confinement = 0.f;
-	const float evaporation = 0.01f;
+	const float elevation_temperature_drop = 30.f;
+	const float evaporation = 1.0f;
 	const float max_humidity_factor = 100.f;
-	const float precipitation_factor = 0.0003f;
+	const float precipitation_factor = 0.0001f;
 	const float force_field_amplitude = 0.00005f;
-	const float random_forces = 1.f;
+	const float random_forces = 0.25f * (static_mode ? 0.f : 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);
 
 	const bool periodic_x = true;
 
-	random::generator rng{random::device{}};
+	// random::generator rng{random::device{}};
-	// random::generator rng{0, 0};
+	random::generator rng{0, 0};
 
 	gfx::pixmap_rgba biomes_map;
 
@@ -146,9 +149,6 @@ 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;
 			}
 		}
 
@@ -167,6 +167,28 @@ struct weather_app
 			}
 		}
 
+		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;
+			}
+		}
+
+		for (int y = 0; y < N; ++y)
+		{
+			for (int x = 0; x < N; ++x)
+			{
+				if (terrain_(x, y) > 0.f)
+					continue;
+
+				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;
+			}
+		}
+
 		make_force_field(rng, force_field_main_, 0.5f);
 		make_force_field(rng, force_field_next_, 0.5f);
 
@@ -249,7 +271,7 @@ struct weather_app
 				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;
 				humidity_(x, y) -= discharge;
-				precipitation_(x, y) = discharge;
+				precipitation_(x, y) = discharge / dt;
 			}
 		}
 
@@ -373,8 +395,8 @@ struct weather_app
 				// 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(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));
+				[[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
@@ -465,7 +487,7 @@ struct weather_app
 		{
 			if (!periodic_x)
 			{
-				float left_boundary_flow = 0.f;//0.01f * std::sin((i * 1.f / N) * float(math::pi) * 4.f);
+				float left_boundary_flow = 0.01f;//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;
@@ -476,6 +498,9 @@ struct weather_app
 			}
 			velocity_(i, 0)[1] = -velocity_(i, 1)[1];
 			velocity_(i, N-2)[1] = -velocity_(i, N-2)[1];
+
+			velocity_(i, 1)[0] = 0.01f;
+			velocity_(i, N-2)[0] = 0.01f;
 		}
 
 		// Uncomment to visualize the force field
@@ -505,8 +530,13 @@ struct weather_app
 		{
 			for (int x = 0; x < N; ++x)
 			{
-				float t = 1.f / frame_;
+				// float t = 1.f / frame_;
 				// float t = 1.f / std::min(8192, frame_);
+				float t = 1.f / std::min<float>(86400.f / dt, frame_); // year average
+
+				if (static_mode)
+					t = 1.f;
+
 				average_temperature_(x, y) = math::lerp(average_temperature_(x, y), temperature_(x, y), t);
 				average_precipitation_(x, y) = math::lerp(average_precipitation_(x, y), precipitation_(x, y), t);
 			}
@@ -610,7 +640,7 @@ struct weather_app
 		auto map_biome = [this](float temperature, float precipitation)
 		{
 			auto x = math::clamp<int>(math::unlerp({  -3.f,  3.f}, precipitation) * biomes_map.width() , {0, biomes_map.width()  - 1});
-			auto y = math::clamp<int>(math::unlerp({-10.f, 40.f}, temperature  ) * biomes_map.height(), {0, biomes_map.height() - 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));
 		};
@@ -633,10 +663,10 @@ struct weather_app
 					else
 					{
 						if (terrain_(x, y) <= 0.f)
-							color = map_color(4.f * terrain_(x, y), {0.f, 0.f, 0.125f, 1.f}, {0.f, 1.f, 1.5f, 1.f});
+							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
 						{
-							float temperature = average_temperature_(x, y) - 273.f - std::max(0.f, terrain_(x, y)) * 6.5f;
+							float temperature = average_temperature_(x, y) - 273.f - std::max(0.f, terrain_(x, y)) * elevation_temperature_drop;
 							// float precipitation = (std::log10(std::max(1e-9f, average_precipitation_(x, y))) + 3.f) * 2.f;
 							// float precipitation = std::pow(average_precipitation_(x, y) * 125.f, 2.f) * 8.f;
 							float precipitation = std::log2(std::max(1e-9f, average_precipitation_(x, y)));
@@ -652,7 +682,7 @@ struct weather_app
 							(terrain_(x, y + 1) - terrain_(x, y - 1)) / 2.f,
 						};
 
-						auto terrain_normal = math::normalized(math::vector{-terrain_gradient[0], -terrain_gradient[1], 0.5f});
+						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}));
 
@@ -673,13 +703,13 @@ 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.001f, {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.00001f, {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), {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;
+					float alpha = 2.f / (1.f + std::exp(- 0.1f * precipitation_(x, y))) - 1.f;
 
 					gfx::color_4f cloud_color{1.f, 1.f, 1.f, alpha * 0.875f};