r.sim: Give each walker random numbers of its own

raster/r.sim/r.sim.sediment/main.c

@@ -57,6 +57,7 @@
 /* INCLUDES                     */
 
 /********************************/
+#include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -102,7 +103,6 @@ int main(int argc, char *argv[])
     struct Cell_head cellhd;
     struct options parm;
     struct flags flag;
-    long seed_value;
 
     G_gisinit(argv[0]);
 
@@ -279,21 +279,16 @@ int main(int argc, char *argv[])
     parm.diffc->description = _("Water diffusion constant");
     parm.diffc->guisection = _("Parameters");
 
-    parm.seed = G_define_option();
+    parm.seed = G_define_standard_option(G_OPT_M_SEED);
     parm.seed->key = "random_seed";
-    parm.seed->type = TYPE_INTEGER;
-    parm.seed->required = NO;
-    parm.seed->label = _("Seed for random number generator");
-    parm.seed->description =
-        _("The same seed can be used to obtain same results"
-          " or random seed can be generated by other means.");
 
     flag.generateSeed = G_define_flag();
     flag.generateSeed->key = 's';
-    flag.generateSeed->label = _("Generate random seed");
+    flag.generateSeed->label =
+        _("Generate random seed (result is non-deterministic) [deprecated]");
     flag.generateSeed->description =
-        _("Automatically generates random seed for random number"
-          " generator (use when you don't want to provide the seed option)");
+        _("This flag is deprecated and will be removed in a future release. "
+          "Seeding is automatic or use parameter random_seed.");
 
     parm.threads = G_define_option();
     parm.threads->key = "nprocs";
@@ -329,20 +324,7 @@ int main(int argc, char *argv[])
                              ? SUMMARY_JSON
                              : SUMMARY_PLAIN;
 
-    if (flag.generateSeed->answer) {
-        seed_value = G_srand48_auto();
-        G_verbose_message(_("Generated random seed (-s): %ld"), seed_value);
-    }
-    else if (parm.seed->answer) {
-        seed_value = atol(parm.seed->answer);
-        G_srand48(seed_value);
-        G_verbose_message(_("Read random seed from %s option: %ld"),
-                          parm.seed->key, seed_value);
-    }
-    else {
-        /* default as it used to be */
-        G_srand48(12345);
-    }
+    int64_t seed = simwe_seed(parm.seed, flag.generateSeed);
 
     G_get_set_window(&cellhd);
 
@@ -353,6 +335,7 @@ int main(int argc, char *argv[])
     ObservationPoints points = {0};
     settings.hhmax = settings.halpha = settings.hbeta = 0;
     settings.ts = false;
+    settings.seed = seed;
     Inputs inputs = {0};
     Outputs outputs = {0};
     Grids grids = {0};

raster/r.sim/r.sim.sediment/r.sim.sediment.html

@@ -44,8 +44,9 @@ <h2>DESCRIPTION</h2>
 <br>
 
 <p>
-Increasing the number of threads with <b>nprocs</b> does not really
-speed up the simulation.
+Increasing the number of threads with <b>nprocs</b> speeds up the
+simulation. The random numbers, the seed and the reproducibility of the
+results with more than one thread are as in <em>r.sim.water</em>.
 
 <h2>NOTES</h2>
 
@@ -93,6 +94,7 @@ <h3>Run summary</h3>
     "walkers_requested": 60,
     "walkers_generated": 78,
     "walkers_remaining": 43,
+    "seed": 1,
     "duration": 60,
     "simulated_time": 54.891343113611583,
     "time_step": 5.4891343113611581,

raster/r.sim/r.sim.sediment/r.sim.sediment.md

@@ -37,8 +37,9 @@ controlled by *duration* \[minutes\] parameter. If the resulting
 erosion/deposition map is noisy, higher number of walkers, given by
 *nwalkers* should be used.  
 
-Increasing the number of threads with **nprocs** does not really speed
-up the simulation.
+Increasing the number of threads with **nprocs** speeds up the
+simulation. The random numbers, the seed and the reproducibility of the
+results with more than one thread are as in *r.sim.water*.
 
 ## NOTES
 
@@ -83,6 +84,7 @@ r.sim.sediment elevation=elevation water_depth=water_depth detachment_coeff=deta
     "walkers_requested": 60,
     "walkers_generated": 78,
     "walkers_remaining": 43,
+    "seed": 1,
     "duration": 60,
     "simulated_time": 54.891343113611583,
     "time_step": 5.4891343113611581,

raster/r.sim/r.sim.sediment/tests/r_sim_sediment_random_test.py

@@ -0,0 +1,117 @@
+# SPDX-FileCopyrightText: 2026 GRASS Development Team
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+"""Tests of the random numbers of r.sim.sediment.
+
+Each walker draws from a random number state of its own, set by the seed and
+the walker's number, so its random numbers do not depend on the number of
+threads. The sediment flux still does at nprocs > 1, because walkers in the
+same cell add to the concentration grid without synchronization, so the test
+which compares numbers of threads compares the walker positions instead. A
+walker's path depends only on the slope and its random numbers.
+"""
+
+import os
+from io import StringIO
+
+import numpy as np
+import pytest
+
+import grass.script as gs
+from grass.exceptions import CalledModuleError
+from grass.tools import Tools
+
+
+@pytest.fixture(scope="module")
+def session(tmp_path_factory):
+    """Session in an XY project with a bowl as elevation and constant inputs."""
+    project = tmp_path_factory.mktemp("r_sim_sediment_random") / "xy_test"
+    gs.create_project(project)
+    with (
+        gs.setup.init(project, env=os.environ.copy()) as session,
+        Tools(session=session) as tools,
+    ):
+        tools.g_region(s=0, n=40, w=0, e=40, res=1)
+        tools.r_mapcalc(
+            expression="elevation = 0.01 * ((col() - 20.5)^2 + (row() - 20.5)^2)"
+        )
+        tools.r_mapcalc(expression="water_depth = 0.1")
+        tools.r_mapcalc(expression="detachment = 0.001")
+        tools.r_mapcalc(expression="transport = 0.001")
+        tools.r_mapcalc(expression="shear_stress = 0")
+        yield session
+
+
+def simulate(session, walkers=None, env=None, **kwargs):
+    """Run the simulation, return the flux and the walker coordinates."""
+    tools = Tools(env=env) if env else Tools(session=session)
+    if walkers:
+        kwargs["walkers_output"] = walkers
+    flux = tools.r_sim_sediment(
+        elevation="elevation",
+        water_depth="water_depth",
+        detachment_coeff="detachment",
+        transport_coeff="transport",
+        shear_stress="shear_stress",
+        sediment_flux=np.array,
+        duration=2,
+        **kwargs,
+    )
+    if not walkers:
+        return flux, None
+    text = tools.v_out_ascii(input=walkers, format="point", precision=15).text
+    return flux, np.loadtxt(StringIO(text), delimiter="|")
+
+
+def test_same_seed_gives_same_flux(session):
+    """Two runs with the same seed give the same sediment flux."""
+    first, _ = simulate(session, random_seed=1)
+    second, _ = simulate(session, random_seed=1)
+    assert np.array_equal(first, second)
+
+
+def test_different_seeds_give_different_results(session):
+    """Runs with different seeds give different fluxes and walker positions."""
+    first_flux, first_walkers = simulate(session, "walkers_seed_1", random_seed=1)
+    second_flux, second_walkers = simulate(session, "walkers_seed_2", random_seed=2)
+    assert not np.array_equal(first_flux, second_flux)
+    assert not np.array_equal(first_walkers, second_walkers)
+
+
+def test_generated_seed_is_the_seed_option(session):
+    """Without a seed, GRASS_RANDOM_SEED gives what the seed option gives."""
+    env = session.env.copy()
+    env["GRASS_RANDOM_SEED"] = "3"
+    generated, _ = simulate(session, env=env)
+    given, _ = simulate(session, random_seed=3)
+    assert np.array_equal(generated, given)
+
+
+@pytest.mark.parametrize("nprocs", [2, 4])
+def test_results_do_not_depend_on_nprocs(session, nprocs):
+    """The flux and the walkers are identical with any number of threads.
+
+    The water depth is an input here, so nothing a walker reacts to is
+    changed by the other walkers of the same step.
+    """
+    serial_flux, serial = simulate(session, f"walkers_{nprocs}_serial", random_seed=5)
+    parallel_flux, parallel = simulate(
+        session, f"walkers_{nprocs}_parallel", random_seed=5, nprocs=nprocs
+    )
+    assert serial.size
+    assert np.array_equal(parallel, serial)
+    assert np.array_equal(parallel_flux, serial_flux)
+
+
+@pytest.mark.parametrize("seed", [-(2**31) - 1, 2**32])
+def test_seed_outside_range_is_an_error(session, seed):
+    """A seed the generator cannot use is refused, not silently wrapped."""
+    with pytest.raises(CalledModuleError, match="outside the range"):
+        simulate(session, random_seed=seed)
+
+
+@pytest.mark.parametrize("seed", ["12abc", "1.5", "-", "99999999999999999999"])
+def test_seed_which_is_not_an_integer_is_an_error(session, seed):
+    """A seed which the parser lets through but which is not an integer is refused."""
+    with pytest.raises(CalledModuleError, match="Invalid random seed"):
+        simulate(session, random_seed=seed)

raster/r.sim/r.sim.sediment/tests/r_sim_sediment_summary_test.py

@@ -15,6 +15,7 @@
     "walkers_requested",
     "walkers_generated",
     "walkers_remaining",
+    "seed",
     "duration",
     "simulated_time",
     "time_step",

raster/r.sim/r.sim.water/main.c

@@ -64,6 +64,7 @@
 /* INCLUDES                     */
 
 /********************************/
+#include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -102,7 +103,6 @@ int main(int argc, char *argv[])
     struct Cell_head cellhd;
     struct options parm;
     struct flags flag;
-    long seed_value;
 
     G_gisinit(argv[0]);
 
@@ -314,21 +314,16 @@ int main(int argc, char *argv[])
     flag.tserie->description = _("Time-series output");
     flag.tserie->guisection = _("Output");
 
-    parm.seed = G_define_option();
+    parm.seed = G_define_standard_option(G_OPT_M_SEED);
     parm.seed->key = "random_seed";
-    parm.seed->type = TYPE_INTEGER;
-    parm.seed->required = NO;
-    parm.seed->label = _("Seed for random number generator");
-    parm.seed->description =
-        _("The same seed can be used to obtain same results"
-          " or random seed can be generated by other means.");
 
     flag.generateSeed = G_define_flag();
     flag.generateSeed->key = 's';
-    flag.generateSeed->label = _("Generate random seed");
+    flag.generateSeed->label =
+        _("Generate random seed (result is non-deterministic) [deprecated]");
     flag.generateSeed->description =
-        _("Automatically generates random seed for random number"
-          " generator (use when you don't want to provide the seed option)");
+        _("This flag is deprecated and will be removed in a future release. "
+          "Seeding is automatic or use parameter random_seed.");
 
     parm.threads = G_define_option();
     parm.threads->key = "nprocs";
@@ -365,20 +360,7 @@ int main(int argc, char *argv[])
                              ? SUMMARY_JSON
                              : SUMMARY_PLAIN;
 
-    if (flag.generateSeed->answer) {
-        seed_value = G_srand48_auto();
-        G_verbose_message(_("Generated random seed (-s): %ld"), seed_value);
-    }
-    else if (parm.seed->answer) {
-        seed_value = atol(parm.seed->answer);
-        G_srand48(seed_value);
-        G_verbose_message(_("Read random seed from %s option: %ld"),
-                          parm.seed->key, seed_value);
-    }
-    else {
-        /* default as it used to be */
-        G_srand48(12345);
-    }
+    int64_t seed = simwe_seed(parm.seed, flag.generateSeed);
 
     G_get_set_window(&cellhd);
 
@@ -417,6 +399,7 @@ int main(int argc, char *argv[])
     G_debug(3, "xmax: %f, ymax: %f", geometry.xmax, geometry.ymax);
 
     settings.ts = flag.tserie->answer;
+    settings.seed = seed;
 
     inputs.elevin = parm.elevin->answer;
     inputs.dxin = parm.dxin->answer;

raster/r.sim/r.sim.water/r.sim.water.html

@@ -206,6 +206,7 @@ <h3>Run summary</h3>
  <tr><td><code>walkers_requested</code></td><td>Number of walkers from <b>nwalkers</b>, by default twice the number of cells</td><td>count</td></tr>
  <tr><td><code>walkers_generated</code></td><td>Walkers created, at least one per cell and more where the source rate is higher</td><td>count</td></tr>
  <tr><td><code>walkers_remaining</code></td><td>Walkers still in the domain at the end of the run</td><td>count</td></tr>
+ <tr><td><code>seed</code></td><td>Seed of the random numbers, given or generated</td><td></td></tr>
  <tr><td><code>duration</code></td><td>Requested simulation length (<b>duration</b>)</td><td>s</td></tr>
  <tr><td><code>simulated_time</code></td><td>Simulated time reached at the end of the run</td><td>s</td></tr>
  <tr><td><code>time_step</code></td><td>Simulated time per iteration</td><td>s</td></tr>
@@ -286,6 +287,7 @@ <h3>Run summary</h3>
     "walkers_requested": 100000,
     "walkers_generated": 120000,
     "walkers_remaining": 112724,
+    "seed": 3,
     "duration": 1200,
     "simulated_time": 1199.2085202681737,
     "time_step": 1.0631281208051186,
@@ -320,6 +322,16 @@ <h3>Run summary</h3>
 }
 </pre></div>
 
+<h3>Random numbers and parallel processing</h3>
+
+The walkers are placed and moved using pseudo-random numbers. The seed
+is given by <b>random_seed</b>; without it, a seed is generated and
+recorded in the history of the output maps and in the run summary as
+<code>seed</code>, so that the run can be repeated. With more than one thread, the
+results differ slightly between thread counts and between repeated runs,
+since the order in which the walkers reach a cell depends on the
+threads. Use <b>nprocs=1</b> when results must be reproducible.
+
 <h2>EXAMPLE</h2>
 
 This example uses the

raster/r.sim/r.sim.water/r.sim.water.md

@@ -196,6 +196,7 @@ stored in the history of the output raster maps under the same keys (see
 | `walkers_requested` | Number of walkers from **nwalkers**, by default twice the number of cells | count |
 | `walkers_generated` | Walkers created, at least one per cell and more where the source rate is higher | count |
 | `walkers_remaining` | Walkers still in the domain at the end of the run | count |
+| `seed` | Seed of the random numbers, given or generated | |
 | `duration` | Requested simulation length (**duration**) | s |
 | `simulated_time` | Simulated time reached at the end of the run | s |
 | `time_step` | Simulated time per iteration | s |
@@ -277,6 +278,7 @@ The printed summary:
     "walkers_requested": 100000,
     "walkers_generated": 120000,
     "walkers_remaining": 112724,
+    "seed": 3,
     "duration": 1200,
     "simulated_time": 1199.2085202681737,
     "time_step": 1.0631281208051186,
@@ -311,6 +313,16 @@ The printed summary:
 }
 ```
 
+### Random numbers and parallel processing
+
+The walkers are placed and moved using pseudo-random numbers. The seed
+is given by **random_seed**; without it, a seed is generated and
+recorded in the history of the output maps and in the run summary as
+`seed`, so that the run can be repeated. With more than one thread, the
+results differ slightly between thread counts and between repeated runs,
+since the order in which the walkers reach a cell depends on the
+threads. Use **nprocs=1** when results must be reproducible.
+
 ## EXAMPLE
 
 This example uses the