From d6f7724fdf3d0bb9f6f66531d4ec67682d3de401 Mon Sep 17 00:00:00 2001
From: Vaclav Petras <wenzeslaus@gmail.com>
Date: Sun, 27 Sep 2026 00:02:04 -0400
Subject: [PATCH] r.sim: Keep the diffusion coefficient in a walker-local
 variable

The diffusion coefficient of a walker's move was stored in a raster-sized
grid cell only to be read back by the same walker a few lines later, so
concurrent walkers in the same cell could read each other's value. Keep
it in a local variable instead and drop the grid.

The local variable is a float, the type of the grid, rather than a
double, so that single-threaded results do not change: the depth,
discharge, sediment flux and erosion/deposition maps of single-threaded
runs on the nc_spm elevation are byte-identical to those before this
change. With more threads, results can differ from before, because a
walker no longer moves with a coefficient another walker wrote.

Written with the help of Claude Code.
---
 raster/r.sim/simlib/hydro.c  | 14 +++++++-------
 raster/r.sim/simlib/input.c  |  2 --
 raster/r.sim/simlib/simlib.h |  1 -
 3 files changed, 7 insertions(+), 10 deletions(-)

diff --git a/raster/r.sim/simlib/hydro.c b/raster/r.sim/simlib/hydro.c
index 215c22414ab..0ebc27fb075 100644
--- a/raster/r.sim/simlib/hydro.c
+++ b/raster/r.sim/simlib/hydro.c
@@ -236,16 +236,17 @@ void main_loop(const Setup *setup, const Geometry *geometry,
 #endif
                             double hhc = pow(d1, 3. / 5.);
                             double velx, vely;
+                            /* Diffusion coefficient of this walker's move */
+                            float dif;
                             if (hhc > settings->hhmax &&
                                 inputs->wdepth == NULL) { /* increased diffusion
                                                      if w.depth > hhmax */
-                                grids->dif[k][l] =
-                                    (settings->halpha + 1) * deldif;
+                                dif = (settings->halpha + 1) * deldif;
                                 velx = sim->vavg[lw].x;
                                 vely = sim->vavg[lw].y;
                             }
                             else {
-                                grids->dif[k][l] = deldif;
+                                dif = deldif;
                                 velx = grids->v1[k][l];
                                 vely = grids->v2[k][l];
                             }
@@ -263,10 +264,9 @@ void main_loop(const Setup *setup, const Geometry *geometry,
                                 }
                             }
 
-                            sim->w[lw].x +=
-                                (velx +
-                                 grids->dif[k][l] * gaux); /* move the walker */
-                            sim->w[lw].y += (vely + grids->dif[k][l] * gauy);
+                            /* Move the walker. */
+                            sim->w[lw].x += (velx + dif * gaux);
+                            sim->w[lw].y += (vely + dif * gauy);
 
                             if (hhc > settings->hhmax &&
                                 inputs->wdepth == NULL) {
diff --git a/raster/r.sim/simlib/input.c b/raster/r.sim/simlib/input.c
index e651d6d0c11..643cb36f2db 100644
--- a/raster/r.sim/simlib/input.c
+++ b/raster/r.sim/simlib/input.c
@@ -34,7 +34,6 @@ void alloc_grids_water(const Geometry *geometry, const Outputs *outputs,
     grids->gama = G_alloc_matrix(geometry->my, geometry->mx);
     if (outputs->err != NULL)
         grids->gammas = G_alloc_matrix(geometry->my, geometry->mx);
-    grids->dif = G_alloc_fmatrix(geometry->my, geometry->mx);
 }
 
 void alloc_grids_sediment(const Geometry *geometry, const Outputs *outputs,
@@ -47,7 +46,6 @@ void alloc_grids_sediment(const Geometry *geometry, const Outputs *outputs,
 
     /* memory allocation for output grids */
 
-    grids->dif = G_alloc_fmatrix(geometry->my, geometry->mx);
     if (outputs->erdep != NULL || outputs->et != NULL)
         grids->er = G_alloc_fmatrix(geometry->my, geometry->mx);
 }
diff --git a/raster/r.sim/simlib/simlib.h b/raster/r.sim/simlib/simlib.h
index edd2e7302e8..61215aec31a 100644
--- a/raster/r.sim/simlib/simlib.h
+++ b/raster/r.sim/simlib/simlib.h
@@ -120,7 +120,6 @@ typedef struct {
     float **er;      // Erosion [output]
     float **ct;      // Transport capacity coefficient [input]
     float **trap;    // Traps [input]
-    float **dif;     // Diffusion coefficient [internal]
 } Grids;
 
 struct point2D {