diff --git a/HighMap/include/highmap/hydrology/hydrology.hpp b/HighMap/include/highmap/hydrology/hydrology.hpp index 9d85fa122..a1c8fb7b6 100644 --- a/HighMap/include/highmap/hydrology/hydrology.hpp +++ b/HighMap/include/highmap/hydrology/hydrology.hpp @@ -1037,6 +1037,9 @@ Array flow_direction_d8(const Array &z); * @brief GPU hydraulic flow simulation using a virtual-pipes model; simulates * shallow-water transport over a height field using iterative flux computation * and water transport passes, with optional flux diffusion and post-simulation + * dry-out. The whole iteration loop runs on the device: the terrain is + * uploaded once, depth and fluxes ping-pong between device images, and only + * the final depth (and velocity, if requested) are read back. * @param z Terrain height field. * @param water_height Global water scaling factor. * @param depth_map Initial relative water distribution. diff --git a/HighMap/src/gpu_opencl/kernels/hydraulic_vpipes.cl b/HighMap/src/gpu_opencl/kernels/hydraulic_vpipes.cl index cd3c5a74a..6a29357e7 100644 --- a/HighMap/src/gpu_opencl/kernels/hydraulic_vpipes.cl +++ b/HighMap/src/gpu_opencl/kernels/hydraulic_vpipes.cl @@ -37,8 +37,10 @@ void kernel hydraulic_vpipes_flow_pass(read_only image2d_t z, float ft_val = TGET(ft, i, j); float fb_val = TGET(fb, i, j); - // fast early-exit: if cell is dry and has no existing flux, outgoing fluxes are strictly zero - if (d0 <= 0.f && fl_val <= 0.f && fr_val <= 0.f && ft_val <= 0.f && fb_val <= 0.f) + // fast early-exit: if cell is dry and has no existing flux, outgoing fluxes + // are strictly zero + if (d0 <= 0.f && fl_val <= 0.f && fr_val <= 0.f && ft_val <= 0.f && + fb_val <= 0.f) { TSET(fl_out, i, j, 0.f); TSET(fr_out, i, j, 0.f); @@ -134,8 +136,7 @@ void kernel hydraulic_vpipes_water_pass(read_only image2d_t z, // fast early-exit: if cell is dry and has no incoming or outgoing fluxes if (d0 <= 0.f && fl_val <= 0.f && fr_val <= 0.f && ft_val <= 0.f && - fb_val <= 0.f && in_l <= 0.f && in_r <= 0.f && in_b <= 0.f && - in_t <= 0.f) + fb_val <= 0.f && in_l <= 0.f && in_r <= 0.f && in_b <= 0.f && in_t <= 0.f) { TSET(d2_out, i, j, 0.f); TSET(u_out, i, j, 0.f); @@ -143,7 +144,8 @@ void kernel hydraulic_vpipes_water_pass(read_only image2d_t z, return; } - float dv = dt * (in_l + in_b + in_r + in_t - fl_val - fr_val - ft_val - fb_val); + float dv = dt * + (in_l + in_b + in_r + in_t - fl_val - fr_val - ft_val - fb_val); float d2_new = max(0.f, d0 + dv / (plength * plength)); if (evap_rate > 0.f) @@ -153,10 +155,14 @@ void kernel hydraulic_vpipes_water_pass(read_only image2d_t z, TSET(d2_out, i, j, d2_new); - float u_new = 0.5f * (in_l - fl_val + fr_val - - ((outflow_boundaries && i == nx - 1) ? 0.f : TGET(fl, i + 1, j))); - float v_new = 0.5f * (in_b - fb_val + ft_val - - ((outflow_boundaries && j == ny - 1) ? 0.f : TGET(fb, i, j + 1))); + float u_new = 0.5f * + (in_l - fl_val + fr_val - + ((outflow_boundaries && i == nx - 1) ? 0.f + : TGET(fl, i + 1, j))); + float v_new = 0.5f * + (in_b - fb_val + ft_val - + ((outflow_boundaries && j == ny - 1) ? 0.f + : TGET(fb, i, j + 1))); float dmean = max(0.001f * water_height, d2_new); @@ -262,4 +268,29 @@ void kernel hydraulic_vpipes_sediment_transport_pass(read_only image2d_t u, TSET(s_out, i, j, s_new); } + +// Adds one rainfall increment. FP contraction is disabled so the result is +// bit-identical to the former host-side `d += rain_map * amount`. +#pragma OPENCL FP_CONTRACT OFF +void kernel hydraulic_vpipes_rain_pass(read_only image2d_t d_in, + read_only image2d_t rain_map, + write_only image2d_t d_out, + const int nx, + const int ny, + const float amount, + const int use_map) +{ + const int2 g = {get_global_id(0), get_global_id(1)}; + + if (g.x >= nx || g.y >= ny) return; + + const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | + CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST; + + float d = TGET(d_in, g.x, g.y); + float r = use_map ? TGET(rain_map, g.x, g.y) * amount : amount; + + TSET(d_out, g.x, g.y, d + r); +} +#pragma OPENCL FP_CONTRACT DEFAULT )"" diff --git a/HighMap/src/gpu_opencl/kernels/shallow_viscous_flow.cl b/HighMap/src/gpu_opencl/kernels/shallow_viscous_flow.cl index 08ab4bad5..39e5ab1cb 100644 --- a/HighMap/src/gpu_opencl/kernels/shallow_viscous_flow.cl +++ b/HighMap/src/gpu_opencl/kernels/shallow_viscous_flow.cl @@ -20,7 +20,7 @@ kernel void shallow_viscous_flow(read_only image2d_t z, const sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST; - const float diag = 0.70710678f; + const float diag = 0.70710678f; float h = TGET(h_in, i, j); float hxp = (outflow_boundaries && i == nx - 1) ? 0.f : TGET(h_in, i + 1, j); @@ -28,7 +28,8 @@ kernel void shallow_viscous_flow(read_only image2d_t z, float hyp = (outflow_boundaries && j == ny - 1) ? 0.f : TGET(h_in, i, j + 1); float hym = (outflow_boundaries && j == 0) ? 0.f : TGET(h_in, i, j - 1); - // fast early-exit: skip all elevation fetches, pow(), and flux math if cell and neighbors are dry + // fast early-exit: skip all elevation fetches, pow(), and flux math if cell + // and neighbors are dry if (h <= 0.f && hxp <= 0.f && hxm <= 0.f && hyp <= 0.f && hym <= 0.f) { TSET(h_out, i, j, 0.f); diff --git a/HighMap/src/hydrology/flow_simulation.cpp b/HighMap/src/hydrology/flow_simulation.cpp index 716514921..59dd34df0 100644 --- a/HighMap/src/hydrology/flow_simulation.cpp +++ b/HighMap/src/hydrology/flow_simulation.cpp @@ -2,7 +2,9 @@ * Public License. The full license is in the file LICENSE, distributed with * this software. */ +#include #include +#include #include #include "cl_wrapper/run.hpp" @@ -46,36 +48,49 @@ Array flow_simulation(const Array &z, const glm::ivec2 shape = z.shape; Array d = water_height * depth_map; - Array fl(shape); // left flux - Array fr(shape); // right - Array ft(shape); // top - Array fb(shape); // bottom - Array u(shape); Array v(shape); - // continuous rainfall increment - Array rain_step; - if (p_rain_map && rain_rate > 0.f) + auto finalize = [&]() { - rain_step = (*p_rain_map) * (rain_rate * dt); - } + if (p_vel_u) *p_vel_u = u; + if (p_vel_v) *p_vel_v = v; - // --- instantiate runners once outside the iteration loop + // remove thin layer of remaining water + if (dry_out_ratio != 0.f) + { + float dmax = d.max(); + water_depth_dry_out(d, dry_out_ratio, nullptr, dmax); + } + return d; // water depth + }; - auto run_fp = clwrapper::Run("hydraulic_vpipes_flow_pass"); + if (iterations <= 0) return finalize(); - run_fp.bind_imagef("z", z.vector, shape.x, shape.y); // inputs - run_fp.bind_imagef("fl", fl.vector, shape.x, shape.y); - run_fp.bind_imagef("fr", fr.vector, shape.x, shape.y); - run_fp.bind_imagef("ft", ft.vector, shape.x, shape.y); - run_fp.bind_imagef("fb", fb.vector, shape.x, shape.y); - run_fp.bind_imagef("d1", d.vector, shape.x, shape.y); + // --- device state: depth and the four fluxes are ping-pong pairs (A/B) + // that stay on the GPU for the whole loop; only the terrain is uploaded + // once and only the results are read back at the end. - run_fp.bind_imagef("fl_out", fl.vector, shape.x, shape.y, true); // outputs - run_fp.bind_imagef("fr_out", fr.vector, shape.x, shape.y, true); - run_fp.bind_imagef("ft_out", ft.vector, shape.x, shape.y, true); - run_fp.bind_imagef("fb_out", fb.vector, shape.x, shape.y, true); + Array zeros(shape); + + const bool use_rain = rain_rate > 0.f; + const bool use_map = use_rain && p_rain_map; + + using clwrapper::Direction; + + // flux pass: (z, fl, fr, ft, fb, d1, fl_out, fr_out, ft_out, fb_out, ...) + auto run_fp = clwrapper::Run("hydraulic_vpipes_flow_pass"); + + run_fp.bind_imagef("z", z.vector, shape.x, shape.y); + run_fp.bind_imagef("fl_a", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("fr_a", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("ft_a", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("fb_a", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("d_a", d.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("fl_b", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("fr_b", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("ft_b", zeros.vector, shape.x, shape.y, Direction::INOUT); + run_fp.bind_imagef("fb_b", zeros.vector, shape.x, shape.y, Direction::INOUT); run_fp.bind_arguments(shape.x, shape.y, @@ -84,18 +99,19 @@ Array flow_simulation(const Array &z, flux_diffusion_strength, outflow_boundaries ? 1 : 0); - auto run_wa = clwrapper::Run("hydraulic_vpipes_water_pass"); - - run_wa.bind_imagef("z", z.vector, shape.x, shape.y); // inputs - run_wa.bind_imagef("fl", fl.vector, shape.x, shape.y); - run_wa.bind_imagef("fr", fr.vector, shape.x, shape.y); - run_wa.bind_imagef("ft", ft.vector, shape.x, shape.y); - run_wa.bind_imagef("fb", fb.vector, shape.x, shape.y); - run_wa.bind_imagef("d1", d.vector, shape.x, shape.y); + // water pass: (z, fl, fr, ft, fb, d1, d2_out, u_out, v_out, ...) + auto run_wa = clwrapper::Run("hydraulic_vpipes_water_pass", + run_fp.get_queue()); - run_wa.bind_imagef("d2_out", d.vector, shape.x, shape.y, true); // outputs - run_wa.bind_imagef("u_out", u.vector, shape.x, shape.y, true); - run_wa.bind_imagef("v_out", v.vector, shape.x, shape.y, true); + run_wa.bind_image2d("z", run_fp.get_image2d("z")); + run_wa.bind_image2d("fl_b", run_fp.get_image2d("fl_b")); + run_wa.bind_image2d("fr_b", run_fp.get_image2d("fr_b")); + run_wa.bind_image2d("ft_b", run_fp.get_image2d("ft_b")); + run_wa.bind_image2d("fb_b", run_fp.get_image2d("fb_b")); + run_wa.bind_image2d("d_a", run_fp.get_image2d("d_a")); + run_wa.bind_imagef("d_b", d.vector, shape.x, shape.y, Direction::INOUT); + run_wa.bind_imagef("u", u.vector, shape.x, shape.y, Direction::OUT); + run_wa.bind_imagef("v", v.vector, shape.x, shape.y, Direction::OUT); run_wa.bind_arguments(shape.x, shape.y, @@ -104,65 +120,81 @@ Array flow_simulation(const Array &z, evap_rate, outflow_boundaries ? 1 : 0); + // rain pass: (d_in, rain_map, d_out, nx, ny, amount, use_map) + std::unique_ptr run_rain; + + if (use_rain) + { + run_rain = std::make_unique("hydraulic_vpipes_rain_pass", + run_fp.get_queue()); + + run_rain->bind_image2d("d_a", run_fp.get_image2d("d_a")); + run_rain->bind_imagef("rain", + use_map ? p_rain_map->vector : zeros.vector, + shape.x, + shape.y); + run_rain->bind_image2d("d_b", run_wa.get_image2d("d_b")); + run_rain->bind_arguments(shape.x, shape.y, rain_rate * dt, use_map ? 1 : 0); + } + + // ping-pong handles + const std::array img_d = {run_fp.get_image2d("d_a").cl_image, + run_wa.get_image2d("d_b").cl_image}; + + const std::array, 2> img_f = { + {{run_fp.get_image2d("fl_a").cl_image, + run_fp.get_image2d("fr_a").cl_image, + run_fp.get_image2d("ft_a").cl_image, + run_fp.get_image2d("fb_a").cl_image}, + {run_fp.get_image2d("fl_b").cl_image, + run_fp.get_image2d("fr_b").cl_image, + run_fp.get_image2d("ft_b").cl_image, + run_fp.get_image2d("fb_b").cl_image}}}; + + int dc = 0; // index of the current depth image + int fc = 0; // index of the current flux images + for (int it = 0; it < iterations; ++it) { - // add continuous rainfall - if (p_rain_map && rain_rate > 0.f) - { - d += rain_step; - } - else if (rain_rate > 0.f) + // continuous rainfall: d[dc] + rain -> d[1 - dc] + if (use_rain) { - d += rain_rate * dt; + run_rain->set_argument(0, img_d[dc]); + run_rain->set_argument(2, img_d[1 - dc]); + run_rain->execute_async({shape.x, shape.y}); + dc = 1 - dc; } - // --- flux update - run_fp.write_imagef("fl"); - run_fp.write_imagef("fr"); - run_fp.write_imagef("ft"); - run_fp.write_imagef("fb"); - run_fp.write_imagef("d1"); - - run_fp.execute({shape.x, shape.y}); + // flux update: reads f[fc], d[dc]; writes f[1 - fc] + for (int k = 0; k < 4; ++k) + run_fp.set_argument(1 + k, img_f[fc][k]); + run_fp.set_argument(5, img_d[dc]); + for (int k = 0; k < 4; ++k) + run_fp.set_argument(6 + k, img_f[1 - fc][k]); - // update flux (from GPU to CPU) - run_fp.read_imagef("fl_out"); - run_fp.read_imagef("fr_out"); - run_fp.read_imagef("ft_out"); - run_fp.read_imagef("fb_out"); + run_fp.execute_async({shape.x, shape.y}); - // --- water transport - run_wa.write_imagef("fl"); - run_wa.write_imagef("fr"); - run_wa.write_imagef("ft"); - run_wa.write_imagef("fb"); - run_wa.write_imagef("d1"); + // water transport: reads f[1 - fc], d[dc]; writes d[1 - dc], u, v + for (int k = 0; k < 4; ++k) + run_wa.set_argument(1 + k, img_f[1 - fc][k]); + run_wa.set_argument(5, img_d[dc]); + run_wa.set_argument(6, img_d[1 - dc]); - run_wa.execute({shape.x, shape.y}); + run_wa.execute_async({shape.x, shape.y}); - run_wa.read_imagef("d2_out"); + fc = 1 - fc; + dc = 1 - dc; } - // retrieve velocity field if requested - if (p_vel_u) - { - run_wa.read_imagef("u_out"); - *p_vel_u = u; - } - if (p_vel_v) - { - run_wa.read_imagef("v_out"); - *p_vel_v = v; - } + run_wa.finish(); - // remove thin layer of remaining water - if (dry_out_ratio != 0.f) - { - float dmax = d.max(); - water_depth_dry_out(d, dry_out_ratio, nullptr, dmax); - } + // retrieve results (both depth images map to the host array `d`) + run_wa.read_imagef(dc == 0 ? "d_a" : "d_b"); + + if (p_vel_u) run_wa.read_imagef("u"); + if (p_vel_v) run_wa.read_imagef("v"); - return d; // water depth + return finalize(); } Array flow_simulation_viscous(const Array &z, diff --git a/external/CLWrapper b/external/CLWrapper index dd3af0e1e..64d740e70 160000 --- a/external/CLWrapper +++ b/external/CLWrapper @@ -1 +1 @@ -Subproject commit dd3af0e1ea6a505cdc4bce77e67a2229d0669e99 +Subproject commit 64d740e706154bd33db157df3c55499e22462043 diff --git a/tests/src/test_export_asset.cpp b/tests/src/test_export_asset.cpp index 19857ec5e..9432fa73f 100644 --- a/tests/src/test_export_asset.cpp +++ b/tests/src/test_export_asset.cpp @@ -1,10 +1,10 @@ #include #include -#include - #include "highmap.hpp" +#include + using namespace hmap; // Regression: the TRI_OPTIMIZED (Delaunay) mesh path wrote texture diff --git a/tests/src/test_flow_simulation.cpp b/tests/src/test_flow_simulation.cpp index 50c869a60..f9a28d20e 100644 --- a/tests/src/test_flow_simulation.cpp +++ b/tests/src/test_flow_simulation.cpp @@ -215,3 +215,106 @@ TEST(FlowSimulation, RainAndEvaporation) EXPECT_LT(d_evap.sum(), (0.1f * d_wet).sum()); } + +TEST(FlowSimulation, IterationParityAndDeterminism) +{ + gpu::init_opencl(); + + const glm::ivec2 shape = {64, 64}; + Array z = noise_fbm(NoiseType::PERLIN, shape, {2.f, 2.f}, 42); + remap(z); + + Array depth_map(shape, 1.f); + const float water_height = 0.05f; + + Array d_even = gpu::flow_simulation(z, water_height, depth_map, 50); + Array d_odd = gpu::flow_simulation(z, water_height, depth_map, 51); + Array d_even_again = gpu::flow_simulation(z, water_height, depth_map, 50); + + // deterministic + EXPECT_EQ(d_even.vector, d_even_again.vector); + + // one more step changes little, but the field must have evolved + const float s_even = d_even.sum(); + const float s_odd = d_odd.sum(); + EXPECT_GT(s_even, 0.f); + EXPECT_NEAR(s_odd, s_even, 0.05f * s_even); + EXPECT_NE(d_even.vector, (water_height * depth_map).vector); + EXPECT_GE(d_odd.min(), 0.f); +} + +TEST(FlowSimulation, ZeroIterationsReturnsInitialState) +{ + gpu::init_opencl(); + + const glm::ivec2 shape = {32, 48}; + Array z = noise_fbm(NoiseType::PERLIN, shape, {2.f, 2.f}, 3); + remap(z); + + Array depth_map = noise_fbm(NoiseType::PERLIN, shape, {4.f, 4.f}, 5); + remap(depth_map); + + Array u(shape, 1.f); + Array v(shape, 1.f); + + Array d = gpu::flow_simulation(z, + 0.1f, + depth_map, + 0, + 0.5f, + true, + 0.01f, + 0.f, + nullptr, + 0.f, + 0.f, + false, + &u, + &v); + + EXPECT_EQ(d.vector, (0.1f * depth_map).vector); + EXPECT_EQ(u.shape, shape); + EXPECT_EQ(v.shape, shape); + EXPECT_FLOAT_EQ(u.max(), 0.f); + EXPECT_FLOAT_EQ(v.max(), 0.f); +} + +TEST(FlowSimulation, RainMapOfOnesEqualsUniformRain) +{ + gpu::init_opencl(); + + const glm::ivec2 shape = {64, 64}; + Array z = noise_fbm(NoiseType::PERLIN, shape, {2.f, 2.f}, 42); + remap(z); + + Array depth_map(shape, 0.f); // start dry + Array ones(shape, 1.f); + + Array d_map = gpu::flow_simulation(z, + 1.f, + depth_map, + 30, + 0.1f, + false, + 0.f, + 0.f, + &ones, + 0.02f, + 0.f, + false); + Array d_uniform = gpu::flow_simulation(z, + 1.f, + depth_map, + 30, + 0.1f, + false, + 0.f, + 0.f, + nullptr, + 0.02f, + 0.f, + false); + + EXPECT_GT(d_map.sum(), 0.f); + EXPECT_EQ(d_map.vector, d_uniform.vector); +}