perf(hydrology): keep flow_simulation state on the GPU across iterations

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.

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,16 +136,16 @@ 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);
     TSET(v_out, i, j, 0.f);
     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
 )""

HighMap/src/gpu_opencl/kernels/shallow_viscous_flow.cl

@@ -20,15 +20,16 @@ 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);
   float hxm = (outflow_boundaries && i == 0) ? 0.f : TGET(h_in, i - 1, j);
   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);

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 <array>
 #include <cmath>
+#include <memory>
 #include <vector>
 
 #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<clwrapper::Run> run_rain;
+
+  if (use_rain)
+  {
+    run_rain = std::make_unique<clwrapper::Run>("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<cl::Image2D, 2> img_d = {run_fp.get_image2d("d_a").cl_image,
+                                            run_wa.get_image2d("d_b").cl_image};
+
+  const std::array<std::array<cl::Image2D, 4>, 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,

tests/src/test_export_asset.cpp

@@ -1,10 +1,10 @@
 #include <cstdio>
 #include <filesystem>
 
-#include <gtest/gtest.h>
-
 #include "highmap.hpp"
 
+#include <gtest/gtest.h>
+
 using namespace hmap;
 
 // Regression: the TRI_OPTIMIZED (Delaunay) mesh path wrote texture