From cfb2234754946ebe6b8b476cfaecf5f34dc648b8 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Tue, 29 Sep 2026 16:23:24 -0400 Subject: [PATCH 1/2] r.sim: Add walker weights to the water depth atomically With nprocs > 1, walkers on different threads added their weights to the same cells of the water depth (or sediment concentration) grid without synchronization. This is a data race: of two concurrent additions to a cell, one can be lost, so parallel runs lose water. On a small bowl with 2000 walkers, the total was 13% to 21% lower at nprocs=4 than at nprocs=1. Add the weight with an atomic capture, so that no addition is lost and the diffusion decision uses the updated sum, as in the serial code. Parallel results still depend on the order in which threads add the weights, as before, but no water is lost. Single-threaded results do not change: the depth, discharge, sediment flux and erosion/deposition maps computed from the nc_spm elevation, with and without infiltration, are byte-identical to those before this change. Written with the help of Claude Code. --- raster/r.sim/simlib/hydro.c | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/raster/r.sim/simlib/hydro.c b/raster/r.sim/simlib/hydro.c index 0ea8db6bfd4..da888e693fd 100644 --- a/raster/r.sim/simlib/hydro.c +++ b/raster/r.sim/simlib/hydro.c @@ -221,11 +221,17 @@ void main_loop(const Setup *setup, const Geometry *geometry, } } - grids->gama[k][l] += - (addac * sim->w[lw].m); /* add walker weigh to - water depth or conc. */ + /* Add walker weight to water depth or + * concentration. The captured sum includes the + * weights added before on any thread. */ + double gama; +#pragma omp atomic capture + { + grids->gama[k][l] += addac * sim->w[lw].m; + gama = grids->gama[k][l]; + } - double d1 = grids->gama[k][l] * conn; + double d1 = gama * conn; double gaux, gauy; #if defined(_OPENMP) gasdev_for_paralel(&gaux, &gauy); From 408dc9b25bb93754d685a76d780fc87298ef6fd2 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Tue, 29 Sep 2026 16:23:25 -0400 Subject: [PATCH 2/2] r.sim.water: Test that no water is lost with nprocs > 1 Compare the total water accumulated with 4 threads against a single-threaded run on a bowl which keeps all walkers in the region. Without diffusion, the random numbers do not affect the walkers, and all of them end up in the few cells at the bottom of the bowl, so unsynchronized updates of the same cell are likely. Also test a walker count which is not a multiple of the thread count, which lost the last walker before the chunk size was rounded up. Without the atomic addition of walker weights, both tests fail with totals 12% to 14% short at nprocs=4 in every run. The tests take about 2 s. Written with the help of Claude Code. --- .../tests/r_sim_water_nprocs_test.py | 78 +++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 raster/r.sim/r.sim.water/tests/r_sim_water_nprocs_test.py diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_nprocs_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_nprocs_test.py new file mode 100644 index 00000000000..8510050fa6c --- /dev/null +++ b/raster/r.sim/r.sim.water/tests/r_sim_water_nprocs_test.py @@ -0,0 +1,78 @@ +"""Tests of r.sim.water results with more than one thread.""" + +import os + +import numpy as np +import pytest + +import grass.script as gs +from grass.tools import Tools + +SIZE = 21 +# Walkers start in the middle of a bowl, far enough from the edges that +# diffusion does not take any of them out of the region. The total of the +# water accumulated over the simulation then depends only on the rainfall +# excess and the time, not on the paths of the walkers or the random numbers. +row, col = np.mgrid[0:SIZE, 0:SIZE] +BOWL = ((row - SIZE // 2) ** 2 + (col - SIZE // 2) ** 2) / 10.0 +RAIN = np.where( + (abs(row - SIZE // 2) <= SIZE // 4) & (abs(col - SIZE // 2) <= SIZE // 4), 50, 0 +) + + +@pytest.fixture(scope="module") +def session(tmp_path_factory): + project = tmp_path_factory.mktemp("nprocs") / "project" + gs.create_project(project) + with gs.setup.init(project, env=os.environ.copy()) as session: + Tools(session=session).g_region(w=0, e=SIZE, s=0, n=SIZE, res=1) + yield session + + +def depth(session, nprocs, rain=RAIN, **kwargs): + return Tools(session=session).r_sim_water( + elevation=BOWL, + depth=np.array, + rain=rain, + man_value=0.05, + nwalkers=400, + duration=2, + random_seed=1, + nprocs=nprocs, + **kwargs, + ) + + +def total_water(depth): + """Return the sum over cells of the accumulated walker weights. + + The depth output is the accumulated weight to the power of 3/5. + """ + return np.sum(depth.astype(np.float64) ** (5 / 3)) + + +def test_total_does_not_depend_on_nprocs(session): + """No water is lost when walkers are processed by several threads. + + Without diffusion, all walkers end up in the few cells at the bottom of + the bowl, so unsynchronized updates of the same cell by several threads + are likely. + """ + options = {"infil_value": 0, "diffusion_coeff": 0} + single = total_water(depth(session, nprocs=1, **options)) + multiple = total_water(depth(session, nprocs=4, **options)) + assert multiple == pytest.approx(single, rel=1e-6) + + +def test_all_walkers_move_with_nprocs(session): + """Every walker is simulated whatever the number of threads. + + With rain in all cells, the walkers requested are created as one per + cell, 441 in total, which is not a multiple of the 4 threads. The last + walker is in a corner of the bowl, so it has to be moved to add its + water. Without diffusion, no walker leaves the region. + """ + options = {"rain": np.full((SIZE, SIZE), 50), "diffusion_coeff": 0} + single = total_water(depth(session, nprocs=1, **options)) + multiple = total_water(depth(session, nprocs=4, **options)) + assert multiple == pytest.approx(single, rel=1e-6)