r.sim: Give each walker random numbers of its own by wenzeslaus · Pull Request #7980 · OSGeo/grass

wenzeslaus · 2026-09-30T20:53:23Z

Builds on the random number layouts (#7979), now merged.

Each walker of r.sim.water and r.sim.sediment now draws from a random number state of its own, placed on the walker's stream of a bounded layout of the seed, instead of from the generator shared by the whole program. A walker's random numbers depend only on the seed and the walker's number, not on the thread which moves it or on the order of the walkers. Without the contention on the shared generator, more threads now make the tools faster instead of slower.

Without a seed, the tools now generate one and record it in the history of the output maps and in the run summary, as r.mapcalc does, instead of using 12345; the -s flag is deprecated. Renaming random_seed to seed is left for a separate PR.

Results for a seed change, also with one thread, and no design avoids that. The old values came from one sequence which all walkers consumed in program order, with a varying number of draws per walker and step, so only replaying that order reproduces them, which is what the shared generator does and what serializes the threads. Per-thread states would keep the single-thread values but give a different result for every thread count. Paying once for values which depend only on the seed and the walker is the better investment, and the outputs of r.sim are sums over millions of walkers, so a different realization of the same seed moves them within their sampling noise, unlike rand() in r.mapcalc, whose values are the output and were kept.

Details (AI generated)

This description and the code were written with Claude Code and reviewed by the author of the PR.

Design

Each walker's record holds its struct G_random_state as a member: the walkers get a type of their own, struct walker (position, weight, state), and struct point3D stays for the output walkers. With the states in an array of their own beside the walkers, as in the first commit, r.sim.water was 7 to 9 % slower at every thread count than with per-thread states; with the state in the record it runs at the same speed, and the random numbers are the same either way. The states are placed in parallel before the walkers, with G_random_init_layout_bounded() over the allocated walkers and a bound of 2 + 8 * steps + 48 values per walker: two values at placement and, per time step, a normal pair by the polar method (8 / pi values on average) plus one in a trap, and a margin of 48 for one step's run of rejected attempts (see below). In the counted runs on the nc_spm elevation (145 and 158 steps, with and without traps) the most a walker drew was 628 values, at most 51 % of its stream.
gasdev() takes the walker's state; gasdev_for_paralel() is gone. The walker loop is unchanged.
Without random_seed, the seed is generated, as in r.mapcalc, instead of being 12345, and recorded in the history of the output maps and in the -p summary as seed, the name the option will get, so that a run can be repeated. The -s flag, which generated a seed, is deprecated with r.mapcalc's message, and random_seed with -s is an error. The option is the standard seed option with its old key. The seed is parsed with strtoll() in simwe_seed(), refusing what is not an integer and a seed outside -2^31 to 2^32 - 1, which was silently reduced to 32 bits before; GRASS_RANDOM_SEED sets the generated seed.
Walker paths are the same for any nprocs as long as the water depth stays below hmax and there is no infiltration. The depth and the sediment concentration still differ slightly at nprocs > 1, because the order of the atomic additions to the depth and the order in which walkers reach a cell with infiltration capacity depend on the threads; this PR does not change that, and the tool docs say so.
The build without OpenMP gives the same results as nprocs=1.

The margin of the bound

The slack of the bound, 8 draws per step against 8 / pi drawn on average, adds up with the steps, while the random excess grows only with their square root, so after a few steps the bound is safe and the exposure is in the first steps: with one step, four rejected attempts in the polar method exceed the stream, which happens with a probability of 0.215^4, about 2e-3 per walker. A walker which exceeds its stream continues into the next walker's, so the two are correlated.

A constant margin on top of the 8 per step is slack for one step's run of rejections, the same at every run length; its protection is 0.215^(margin / 2) per walker. A floor on the steps, max(steps, N), protects short runs by treating them as N steps long, so N must sit where the per-step rule is already safe. The table gives the probability per walker and the probability that a run of N walkers has any overrun, for a one-step run with a trap. The walker counts come from the region: twice the cells are requested by default and one more per cell is allocated, so 10^7 is an everyday run, 10^8 a large one and 2^31 the most the int counts hold.

margin	stride	per walker	10^6 walkers	10^7	10^8	2 * 10^9
none (8 per step)	11	2e-3	100 %	100 %	100 %	100 %
+16	27	1e-8	0.95 %	9.1 %	61 %	100 %
+32	43	4e-14	4e-6 %	4e-5 %	4e-4 %	9e-3 %
+48	59	2e-19	2e-11 %	2e-10 %	2e-9 %	4e-8 %
+64	75	9e-25	0	0	1e-14 %	2e-13 %
max(steps, 30)	243	6e-81	0	0	0	0

The cost is in the batches of the layout, the number of walker sets which fit into the span of 2^46 draws, which r.sim does not use beyond batch 0 today. At the int ceiling of walkers:

steps	none	+32	+48	+64	max(steps, 30)
1	3,198	818	596	469	144
10	423	305	268	239	144
100	43	42	41	40	43
600	7	7	7	7	7
3600	1	1	1	1	1

The margin costs batches only below about ten steps, where hundreds remain even at the ceiling, and nothing measurable elsewhere; the floor costs most of what a short run has. 48 is the first value at which the probability stays negligible at the ceiling, below 1e-6 for a run; 64 buys nothing observable. The margin changes every stream, so the references were regenerated once more.

Results and speed

The gunittest references and the pinned landform arrays were regenerated for the new values; the commands used reproduce the old references exactly with the previous code, so nothing but the random numbers changed.

r.sim.water on the whole nc_spm elevation (6,075,000 walkers, 145 iterations), wall seconds on a shared 36-core machine, medians of 3 runs (of 8 runs pinned to one CPU for "after" at nprocs 1):

	nprocs 1	nprocs 4	nprocs 16
before (shared generator with atomics)	67.9	143.4	158.7
after (a state in each walker's record)	50.0	13.3	4.8

Tests

New r_sim_water_random_test.py (15 tests) and r_sim_sediment_random_test.py (11 tests): the same walkers for nprocs 1, 2 and 4 with a large hmax and no infiltration, different seeds give different results, the seed range and non-integer seeds are refused, GRASS_RANDOM_SEED seeds the run, a run without a seed records the seed which repeats it, the deprecated -s still generates one, random_seed with -s is refused, and the summary reports the seed.
The flow control test of the existing suite passed at seed 42 by luck: with the ten walkers five cells get by default, a trapping probability of 0.3 raises the total depth in 23 of 30 seeds. It now runs with a thousand walkers, where the three cases separate for every seed tried.
The whole r.sim pytest suite (94 tests) runs in about 33 s; both gunittest suites pass.
The run summary examples in the tool docs still show numbers from the previous generator; the region they were made with is not recorded in the docs. The new seed line was added to them with the seed of the command.

Suggested squash commit message

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

Each walker draws from a random number state of its own, placed on the
walker's stream of a bounded layout of the seed, instead of from the
shared generator. A walker's random numbers now depend only on the
seed and the walker's number, not on the thread which moves it or on
the order of the walkers, so walker paths are the same for any nprocs
as long as the depth-dependent diffusion and infiltration do not come
into play. The water depth and the sediment concentration still differ
slightly at nprocs > 1 because the order of the atomic additions to the
gama grid and the order in which walkers reach a cell with infiltration
capacity depend on the threads, which this does not change.

The streams hold 2 + 8 * steps + 48 values, over twice what a walker
draws on average: two values at placement and, per step, a normal pair
by the polar method, 8 / pi values on average, plus one in a cell of
the flow_control map (a trap in the code). The 8 per step hold the sum
over many steps, and the 48 allow a run of 24 rejected attempts within
one step, so that a walker exceeds its stream with a probability below
1e-16 and a run of 2^31 walkers, the most the int counts hold, with a
probability below 1e-6. In counted runs on the nc_spm elevation, the
most a walker drew over 145 or 158 steps was 628 values, at most 51% of
its stream. A run whose walkers and steps need more than one seed's
span provides gets a warning. The states are placed in parallel before
the walkers are, each a member of its walker's record, for which the
walkers get a type of their own, struct walker. A first version kept
the states in an array of their own beside the walkers, which cost 7 to
9% of the wall time at every thread count, likely for the extra stream
of memory reads and writes in each walker's step; in the record, the
state shares the cache line with the walker's position and weight and
the cost is gone. Without the contention on the shared generator,
nprocs > 1 is now faster than nprocs = 1: on the whole nc_spm elevation
(6 million walkers), 50 s at nprocs=1, 13 s at 4 and 4.8 s at 16,
against 68 s, 143 s and 159 s before. The build without OpenMP now
gives the same results as nprocs=1.

Without random_seed, the seed is generated, as in r.mapcalc, instead
of being 12345, and recorded in the history of the output maps and in
the -p summary, so that a run can be repeated. The results therefore
differ from run to run unless random_seed is set. The -s flag is
deprecated, and random_seed with -s is an error. The seed is parsed
with strtoll() and checked right after parsing, so a seed outside
-2^31 to 2^32 - 1 is an error instead of being silently reduced to 32
bits. The tool docs get a section on the random numbers, the seed and
reproducibility with more than one thread.

This change has the price of changed outputs also at nprocs=1, but
invests into results independent of the thread count where the model
allows it. The alternative, one state per thread, would keep the
values at nprocs=1 and give a different result for each thread count.
The outputs of r.sim.water and r.sim.sediment are sums over millions
of walkers, so a different realization of the same seed moves them
within their sampling noise, and from the user's perspective they do
not change much anyway. The gunittest references and the pinned
landform arrays are regenerated; the commands used reproduce the old
references exactly with the previous code.

Written with the help of Claude Code (Fable and Opus).

Tools which draw random numbers in parallel, or run ensembles, or keep several model components reproducible on their own, cannot do so with the one generator shared by the whole program: whichever thread asks first gets the next value, and different seeds are related on this generator (seeds 2^31 or 2^30 apart are shifted by a constant, and consecutive seeds form a lattice). Add generators the program owns, struct G_random_state, and layouts, struct G_random_layout, which place one state per unit of work on the generator's cycle. A layout is a seed plus the way the span is cut into streams of stride draws, one stream per unit and run. The span is the first quarter of the cycle after the seed, 2^46 draws: the multiplier has order 2^46, so positions 2^46, 2^47 or 3 * 2^46 apart give values differing by a constant, and the other three quarters hold the same values plus 0.25, 0.5 and 0.75. Within the span no two positions are that far apart. An exact layout takes the draws of a unit as given, so the streams are the positions of a serial run and a tool keeps its earlier values; a bounded layout takes a bound and rounds it up to odd, so that units a power of two apart relate only a few low bits; a spread layout gives every unit the largest stream the span allows, for one run, dividing the span into an odd number of parts: with an even number, the unit halfway along would start 2^45 draws after unit 0 and the unit a quarter of the way along near 2^44, distances at which the values differ by two or four constants in turn. A spread layout holds at most 2^20 units, because the stride rounded down to odd loses up to two draws per unit, which accumulate along the span; up to that count the units nearest to 2^45 and 2^44 start at least 0.48 and 0.24 of a stride away from them. Runs are appended blocks of streams, never reserved, so an ensemble needs only the run number. The queries return the runs that fit and the stride; the two state calls place a state on a unit's stream, or a unit's stream of a run; G_random_advance() moves a state by any number of draws in one jump; G_random_seed() gives the shared generator's sequence. Add G_random_generate_seed(), which reads GRASS_RANDOM_SEED, or else SOURCE_DATE_EPOCH, as an integer: an empty value counts as not set, a value which is not an integer or does not fit is a fatal error, a value from -2^31 to 2^32 - 1 is used as it is, and a value outside that range is reduced to its low 32 bits with a warning. Without either variable the seed is a weak hash of the time and the process ID. G_srand48_auto() seeds with that value and returns it; before, it returned an environment value unreduced although it had used only its low 32 bits. Add GRASS_RANDOM_SEED to the list of environment variables. The shared generator is unchanged. Tests cover the identity of G_random_seed() with G_srand48(), the three layouts, the runs, the refusals, the absence of constant-shift twins within the span, the distance of spread units from the two- and four-constant levels for every unit count a spread layout accepts, the low-bits rule, the advance, the seed from each environment source, and independence from the shared generator and between threads. Written with the help of Claude Code.

The tests which set GRASS_RANDOM_SEED and SOURCE_DATE_EPOCH with pytest's monkeypatch in the test process failed on Windows, where the Python interpreter and the GRASS library use separate C runtimes with environments of their own, so a variable set from Python is not what the library's getenv() reads. Run these tests in a subprocess which starts with the variables set, as the other environment tests in the file do, and check the shared generator's sequence after G_srand48_auto() through the values the subprocess prints. Written with the help of Claude Code.

Add Vaclav Petras to the authors of lrand48.c and drop the remark in parentheses after Maris Nartiss, since the header now describes the thread safety itself. Add the SPDX lines to the random number test file, whose header changes anyway. Written with the help of Claude Code.

The file header said the states a program owns are never shared, which the library cannot promise; whether a state is used by one thread at a time is up to the program. Say that instead, as the condition under which the functions need neither atomics nor locks and a unit's values do not depend on the thread which draws them. Written with the help of Claude Code.

The environment reader and the check of the seeds given to the layout and state functions each spelled out the range from -2^31 to 2^32 - 1. Test it in one place, with the reason for the range, and let both use it. Written with the help of Claude Code.

The comment said the caller has checked that units * stride does not overflow, which reads as a description rather than as the requirement it is. Say that callers must check it beforehand. Written with the help of Claude Code.

wenzeslaus · 2026-09-30T21:41:56Z

One question here is, since this is changing results, should be change the initial seed (value = 12345;) to let's say 1? Or, should we even move to the r.mapcalc behavior when the seed is different for every run? I would say yes to the first new seed and no to the randomization of seed (unlike r.mapcalc rand(), the stochasticity is not expected for SIMWE and the results are generally the same for different seeds). The argument for making it randomized would be consistency with other random-number driven tools. Either way, I don't think it is that important for SIMWE, so that's why I think it is okay to change it between minor versions (for years, there was only one seed value hardcoded which is the current default).

Each walker draws from a random number state of its own, placed on the walker's stream of a bounded layout of the seed, instead of from the shared generator. A walker's random numbers now depend only on the seed and the walker's number, not on the thread which moves it or on the order of the walkers, so walker paths are the same for any nprocs as long as the depth-dependent diffusion and infiltration do not come into play. The water depth and the sediment concentration still differ slightly at nprocs > 1 because the order of the atomic additions to the gama grid depends on the threads and the updates of the inf grid are not synchronized, which this does not change. The streams hold 2 + 8 * steps values, over twice what a walker draws on average: two values at placement and, per step, a normal pair by the polar method, 8 / pi values on average, plus one in a trap. A single step can draw more than eight values, but over a whole simulation a walker stays well within its stream; on the nc_spm elevation, no walker drew more than 53% of it. The states are placed in parallel before the walkers are. Without the contention on the shared generator, nprocs > 1 is now faster than nprocs = 1: on the whole nc_spm elevation (4 million walkers), 50 s at nprocs=1, 13.5 s at 4 and 5.6 s at 12, against 63 s, 125 s and 140 s before. The build without OpenMP now gives the same results as nprocs=1. The seed is parsed with strtoll(), refusing what is not an integer, and checked right after parsing, so a seed outside -2^31 to 2^32 - 1 is an error instead of being silently reduced to 32 bits. The -s flag uses G_random_generate_seed(), which honors GRASS_RANDOM_SEED. Results for a seed change also at nprocs=1, and no design keeps them. The old values came from one sequence which all walkers consumed in program order, with a varying number of draws per walker and step, so a walker's values depended on wherever the other walkers had left the sequence, and only replaying that order reproduces them, which is what the shared generator does and what serializes the threads. The other design, one state per thread, keeps the values at nprocs=1 and gives a different result for every other thread count. Paying once for values which depend only on the seed and the walker's number is the better investment. The outputs of r.sim.water and r.sim.sediment are sums over millions of walkers, so a different realization of the same seed moves them within their sampling noise, unlike rand() in r.mapcalc, whose values are the output and therefore kept their old values. The gunittest references and the pinned landform arrays are regenerated; the commands used reproduce the old references exactly with the previous code. Written with the help of Claude Code.

The two new test files carried the license identifier only; add the copyright line the style guide asks for in every file. Written with the help of Claude Code.

Make the struct G_random_state of a walker a member of the walker's record instead of an element of an array of its own beside the walkers. The walkers get a type of their own, struct walker, because struct point3D also serves the output walkers, where m holds the elevation and no state is needed. The random numbers are the same; only where the state lives changes. At nprocs=1 the depth, discharge and error maps of r.sim.water, the five output maps of r.sim.sediment and the output walkers are bit-identical to the previous commit, so the tests and their pinned values are unchanged. The array cost time. On the whole nc_spm elevation (6 million walkers, 145 time steps) at nprocs=1, pinned to one CPU, median of 8 runs, r.sim.water took 53.8 s with the array and 50.0 s with the state in the record, close to the 49.4 s of the per-thread states of OSGeo#7884; at nprocs 4 and 16, 13.3 s and 4.8 s against 14.1 s and 5.1 s with the array (median of 3 runs). Why the array costs this much was not measured. Likely it adds one more stream of memory reads and writes to each walker's step, where the state in the record shares the cache line with x, y and m, and the loop passes the address of an array element to a library function, which the compiler cannot see through. Written with the help of Claude Code.

G_random_seed() becomes G_random_state_from_seed(): it places a state as G_random_state_for_unit() does, and the old name read as returning a seed next to G_random_generate_seed(). G_random_init_layout_spread() becomes G_random_init_layout(): the layout which gives the units the whole span is the basic one, of which a single seeded state is the case of one unit, while the exact and bounded layouts are those with a known stride. "Spread" also names a process in tools which would call the function, such as r.spread. A run becomes a batch, so G_random_layout_runs() and G_random_state_for_run() become G_random_layout_batches() and G_random_state_for_batch(). A batch of units says what the second index is whether a program uses it for the members of an ensemble, for repeated passes or for the processes of a model, and "run" keeps its ordinary meaning. The internal init_layout() becomes fill_layout(), the fields of the layout and the limit on the units of a layout of the whole span are renamed to match, and the messages say batch. The longer name of the limit realigns the neighboring definitions. Nothing else changes; the tests are renamed with the functions. Written with the help of Claude Code.

Call G_random_state_from_seed() instead of G_random_seed() to check the seed and G_random_layout_batches() instead of G_random_layout_runs() to check that the walkers fit one seed. Written with the help of Claude Code.

wenzeslaus · 2026-10-01T20:51:04Z

In my benchmarks, I noticed that this PR was slightly slower than the reentrant random number streams (#7884). Moving the RNG state to a new walker structure improved the speed to that of the per-thread RNG states of #7884.

With one thread, this PR takes 26 % less time than current main (49.7 s against 67.2 s), and with more threads it keeps getting faster, down to 4.8 s at 16 threads. Current main, where all threads share one generator advanced with atomics (#6480), gets slower with more threads, up to 148 s at 16 threads, 31 times the time of this PR. The old code before #6480, where the threads shared the generator without any synchronization and so had a data race, did not scale either and went from 61 s with one thread to 74 s with 16.

r.sim.water on the whole nc_spm elevation:

Details (AI generated)

The measurements, tables and this text were made with Claude Code and reviewed by the author of the PR.

Builds

All builds are Autotools with the same configure line, gcc, CFLAGS -g -O2, OpenMP on, and run with random_seed=1.

label	commit	random numbers in the walker loop
this PR	`7c0d776`	one `G_random_state` per walker, a member of the walker record (`struct walker`), placed on a bounded layout
current main, atomic in RNG	`1f83c2d`	shared `G_drand48()` advanced with C11 atomics (the thread-safe generator, #6480)
before #6480, data race on RNG	`55de52a`	shared `G_drand48()` without synchronization; results are racy, only the timing is meaningful; also lacks the walker-loop fixes #7964, #7965 and #7974
this PR, separate array for states	`c75d1c5`	the previous head of this PR: the same states in an array of their own
per-thread states (#7884 on main)	`6db2f46` merged onto main locally	one generator state per thread, seeded by chunk index

The C code of main is the same in all builds except the one before #6480. None of them has the latest r.sim changes on main, the lat-lon refusal (#7978) and the null handling of Manning's n and infiltration (#7981), which do not touch the walker loop.

Protocol

Shared Intel Xeon W-2295 (18 cores, 36 logical CPUs), Linux 6.8. One grass --tmp-mapset session on nc_spm_08_grass7; each run uses the binary and libraries of its build, checked through /proc/<pid>/maps. Builds and thread counts are interleaved within each of 3 repeats; the tables give the median and, in brackets, the minimum and maximum of the 3 runs. The 1-minute load average was 2.6 to 17.9 (median 7.6) during the runs.

Cases:

1500 x 1350 cells: g.region raster=elevation, default walkers (the tool reports 4,050,000 requested; 6,075,000 are simulated), 145 iterations.
100 x 100 cells: n=224000 s=223000 e=637000 w=636000 res=10, nwalkers=1000000, 158 iterations for r.sim.water and 68 for r.sim.sediment (water depth 0.05, detachment and transport coefficients 0.001, shear stress 0).

Results

r.sim.water, 1500 x 1350 cells, 6,075,000 walkers; wall time in seconds, median (min-max) of 3 runs:

	nprocs 1	nprocs 2	nprocs 4	nprocs 8	nprocs 16
this PR	49.7 (49.4-50.3)	25.5 (25.4-25.9)	13.2 (13.2-13.3)	7.4 (7.1-7.4)	4.8 (4.7-5.0)
current main, atomic in RNG	67.2 (66.3-67.2)	106.9 (104.5-117.5)	144.1 (143.9-144.2)	142.0 (140.4-142.9)	147.9 (145.8-148.0)
before #6480, data race on RNG	61.3 (60.9-62.1)	62.2 (61.4-62.7)	64.4 (63.3-64.7)	68.6 (67.7-70.0)	74.3 (73.6-75.9)
this PR, separate array for states	52.7 (52.4-52.7)	27.0 (27.0-27.8)	14.0 (13.9-14.2)	7.8 (7.8-7.8)	5.3 (5.3-5.3)
per-thread states (#7884 on main)	49.2 (49.1-49.7)	25.2 (25.2-25.6)	13.3 (13.0-13.6)	7.3 (7.2-7.4)	4.8 (4.8-5.0)

r.sim.water, 100 x 100 cells, 1,000,000 walkers; wall time in seconds, median (min-max) of 3 runs:

	nprocs 1	nprocs 2	nprocs 4	nprocs 8	nprocs 16
this PR	7.49 (7.46-7.72)	4.09 (4.06-4.14)	2.15 (2.14-2.24)	1.20 (1.18-1.21)	0.83 (0.80-0.84)
current main, atomic in RNG	10.25 (10.19-10.41)	18.86 (18.82-19.30)	23.37 (23.35-23.98)	23.79 (23.52-23.90)	24.82 (24.65-24.90)
before #6480, data race on RNG	9.46 (9.40-9.47)	10.08 (9.96-10.35)	10.86 (10.81-11.05)	11.39 (11.22-11.43)	12.52 (12.19-12.59)
this PR, separate array for states	7.76 (7.75-7.90)	4.27 (4.20-4.34)	2.23 (2.19-2.26)	1.24 (1.22-1.26)	0.85 (0.84-0.87)
per-thread states (#7884 on main)	7.37 (7.37-7.53)	4.06 (4.04-4.10)	2.14 (2.14-2.16)	1.20 (1.18-1.23)	0.88 (0.82-0.88)

r.sim.sediment, 100 x 100 cells, 1,000,000 walkers; wall time in seconds, median (min-max) of 3 runs:

	nprocs 1	nprocs 2	nprocs 4	nprocs 8	nprocs 16
this PR	3.92 (3.90-4.03)	2.01 (2.01-2.05)	1.10 (1.09-1.12)	0.59 (0.58-0.59)	0.40 (0.39-0.41)
current main, atomic in RNG	5.20 (5.14-5.21)	9.05 (8.41-9.20)	11.17 (10.37-11.46)	11.39 (11.34-11.49)	11.52 (11.50-11.67)
before #6480, data race on RNG	4.75 (4.67-4.77)	4.87 (4.77-5.02)	5.10 (5.01-5.12)	5.45 (5.42-5.51)	5.88 (5.79-5.88)
this PR, separate array for states	3.98 (3.93-4.06)	2.09 (2.06-2.15)	1.10 (1.08-1.16)	0.63 (0.60-0.64)	0.42 (0.41-0.42)
per-thread states (#7884 on main)	3.88 (3.84-3.95)	1.98 (1.96-2.04)	1.09 (1.05-1.10)	0.60 (0.56-0.62)	0.39 (0.39-0.40)

State in the walker record

With one thread on the 1500 x 1350 case, pinned to one CPU and with 8 repeats in rotating order, the median wall time was 50.58 s for this PR, 50.60 s for per-thread states and 55.07 s for the previous head with the separate array. Both heads of this PR write the same depth map (maximum absolute difference 0). Why the separate array was slower was not measured (no hardware counters were available); likely causes are the extra memory stream of 8 bytes per walker and step, and that a pointer into the array passed to the library keeps the compiler from holding values in registers across the call.

Peak memory (maximum resident set size) on the 1500 x 1350 case is 444 MB for this PR and 396 MB for per-thread states: the walker record grows from 24 to 32 bytes, 48.6 MB for the 6,075,000 allocated walkers, as much as the separate array took.

Batches started units * stride draws apart. When that number is even, the values which different batches give at the same draw are related, the more simply the more often 2 divides the distance: with a distance divisible by 2^20, the value of every batch is given by a quadratic through those of the three batches before it, much as consecutive seeds give values on a straight line. Make the distance units * stride, or one more when that is even, as the strides of bounded layouts and of layouts of the whole span are already odd for the same reason between units. Batch 0 does not change, a batch costs at most one draw more, and no tool uses further batches yet. This does not remove the relation. It moves it from the draws which carry the same number in every batch to draws whose numbers differ by one from batch to batch. A layout gives streams which do not overlap, not independent ones. The batches that fit are now those whose last stream ends within the span, so a batch of exactly 2^46 draws still fits once. Written with the help of Claude Code.

The infiltration capacity of a cell is now reduced under a lock, so the note on results with more threads names what still depends on them: the order in which walkers reach the cell. Written with the help of Claude Code.

The comment on the trap draw was carried over from the old call it described, and the one on the walker array repeats its type. Written with the help of Claude Code.

Without random_seed, the seed is now generated, as in r.mapcalc, instead of being 12345. The seed used is recorded in the history of the output maps and in the -p summary as random_seed, so that a run can be repeated. The -s flag, which generated a seed, is deprecated with the message r.mapcalc gives, and random_seed with -s is an error, as there. The option is now the standard seed option with its old key. Written with the help of Claude Code.

The type comes with grass/gis.h, but code which uses it includes the header itself rather than relying on that. Written with the help of Claude Code.

The bound of 2 + 8 draws per time step holds the sum over many steps with room to spare, but the first steps are exposed: with one step, a run of four rejected attempts in the polar method exceeds the stream, which happens with a probability of 2e-3 per walker, so a short run on a large region has walkers which continue into the next walker's stream. A constant margin of 48 draws, 24 further attempts, bounds the probability per walker by 0.215^24, below 1e-16, at any number of steps, so a run of 2^31 walkers exceeds the bound with a probability below 1e-6. The streams move, so the pinned landform depths and the gunittest references change; the infiltration reference differs only by rounding and is kept. Written with the help of Claude Code.

With the ten walkers the five cells get by default, the total depth with a trapping probability of 0.3 is within noise of the total without trapping: over 30 seeds it was higher in 23. With a thousand walkers the three cases separate for every seed tried, so the test checks the effect rather than the seed. Written with the help of Claude Code.

The key in the history and the run summary is new, so it takes the name the option will have, seed, which the other tools use. Written with the help of Claude Code.

Written with the help of Claude Code.

simlib.h, included above them, defines both structs. The file no longer mentions point3D at all now that the walkers have a type of their own, and it never used point2D. Written with the help of Claude Code.

The section described the generator; the reader needs what the seed does and when the results depend on the number of threads. Results of r.sim.water are the same for any nprocs except where the water depth exceeds hmax or infiltration is used, since a walker there reacts to what the walkers before it in the same step left in the cell, and that order depends on the threads. Results of r.sim.sediment do not depend on nprocs at all, as its water depth is an input. Both statements are confirmed on the nc_spm elevation at 1, 2, 4 and 8 threads and tested on the synthetic data: the depth is identical below hmax, the sediment flux always. The seed range, the verbose message and the environment variable are not repeated here: the tool refuses a bad seed, and the variable is documented on the environment variables page. The note on seeds below 2^30 goes: it concerns runs with many seeds and belongs to a general page on random numbers, should one be needed. Written with the help of Claude Code.

The condition under which they are the same for any thread count, a water depth below hmax and no infiltration, is too narrow to promise; the rule for the user is nprocs=1. r.sim.sediment refers to r.sim.water for the same rule rather than claiming more. Written with the help of Claude Code.

The figures were made with the head of OSGeo#7980 before it was merged. Rerun the benchmark script on this branch merged with main, which now contains OSGeo#7980 (per-walker random numbers) and OSGeo#8013 (the 48-bit multiplication in a helper), with the same sizes, threads and repeats. The times agree with the earlier run within about 2 % (1 m: 113.1 s with 1 thread, 12.2 s with 16, speedup 9.3). Written with the help of Claude Code.

Measure r.sim.water with nprocs from 1 to 16 on the elevation of the SECREF North Carolina dataset at 4 m, 2 m and 1 m resolution, using grass.benchmark as the benchmarks of other tools do. The script plots time, speedup and efficiency and saves the results as JSON, from which it can redo the plots. In grass.benchmark, nprocs_plot() shades the range of the runs in each line's color for time, speedup and efficiency instead of a gray band for time only (for speedup over the ratios of every serial run to every parallel run), takes figsize and dpi, plots results without all_times, and sets the default title for time. The range is computed by nprocs_ranges(), which has pytest tests. The documentation of the experimental subpackage says that it is changed together with the benchmarks which need more from it. Add the figures to the documentation in a Performance section, as the other parallel tools have, at twice the 600 pixel width, shown at 600 pixels in the HTML page. They were made with the per-walker random number states of #7980. Written with the help of Claude Code.

wenzeslaus added 6 commits September 30, 2026 12:36

This was referenced Sep 30, 2026

r.sim: Give each walker random numbers of its own wenzeslaus/grass#61

Closed

libgis: Add a helper which reads the seed from a tool's options wenzeslaus/grass#62

Closed

github-actions Bot added raster Related to raster data processing Python Related code is in Python C Related code is in C HTML Related code is in HTML libraries module docs markdown Related to markdown, markdown files tests Related to Test Suite labels Sep 30, 2026

wenzeslaus added 3 commits September 30, 2026 22:02

r.sim: Add copyright lines to the random number tests

1cf8a00

The two new test files carried the license identifier only; add the copyright line the style guide asks for in every file. Written with the help of Claude Code.

Merge remote-tracking branch 'upstream/main' into rsim-layouts

efa774c

wenzeslaus force-pushed the rsim-layouts branch from acbdf3e to efa774c Compare October 1, 2026 02:03

wenzeslaus added 4 commits October 1, 2026 02:17

Merge branch 'random-layouts' into rsim-layouts

bc46c1e

r.sim: Follow the renamed random number functions

83a4c95

Call G_random_state_from_seed() instead of G_random_seed() to check the seed and G_random_layout_batches() instead of G_random_layout_runs() to check that the walkers fit one seed. Written with the help of Claude Code.

wenzeslaus mentioned this pull request Oct 2, 2026

r.sim.water: Add a benchmark for the number of threads wenzeslaus/grass#69

Closed

wenzeslaus added 2 commits October 2, 2026 13:56

Merge branch 'random-layouts' into rsim-layouts

f30708d

Merge remote-tracking branch 'upstream/main' into rsim-layouts

5554480

wenzeslaus marked this pull request as ready for review October 4, 2026 03:54

r.sim.water: Describe the locked infiltration in the nprocs note

96cd2a0

The infiltration capacity of a cell is now reduced under a lock, so the note on results with more threads names what still depends on them: the order in which walkers reach the cell. Written with the help of Claude Code.

petrasovaa reviewed Oct 5, 2026

View reviewed changes

Comment thread raster/r.sim/simlib/random.c Outdated

petrasovaa reviewed Oct 5, 2026

View reviewed changes

Comment thread raster/r.sim/simlib/hydro.c Outdated

petrasovaa reviewed Oct 5, 2026

View reviewed changes

Comment thread raster/r.sim/simlib/random.c

petrasovaa reviewed Oct 5, 2026

View reviewed changes

Comment thread raster/r.sim/simlib/hydro.c Outdated

wenzeslaus added 7 commits October 5, 2026 13:05

r.sim: Drop comments which restate the code

9599115

The comment on the trap draw was carried over from the old call it described, and the one on the walker array repeats its type. Written with the help of Claude Code.

r.sim: Include stdint.h where the seed is held as int64_t

e64d51b

The type comes with grass/gis.h, but code which uses it includes the header itself rather than relying on that. Written with the help of Claude Code.

r.sim: Record the seed under the key seed

73d86a0

The key in the history and the run summary is new, so it takes the name the option will have, seed, which the other tools use. Written with the help of Claude Code.

r.sim.water: Name the cases the flow control test separates

1bb9d0f

Written with the help of Claude Code.

wenzeslaus mentioned this pull request Oct 5, 2026

r.sim: Rename the option random_seed to seed wenzeslaus/grass#72

Closed

wenzeslaus added 3 commits October 5, 2026 17:07

r.sim: Drop two forward declarations hydro.c does not need

d831e04

simlib.h, included above them, defines both structs. The file no longer mentions point3D at all now that the walkers have a type of their own, and it never used point2D. Written with the help of Claude Code.

This was referenced Oct 6, 2026

r.sim: Accumulate water depth independently of walker order wenzeslaus/grass#74

Open

r.sim: Accumulate water depth independently of walker order wenzeslaus/grass#58

Closed

petrasovaa approved these changes Oct 7, 2026

View reviewed changes

wenzeslaus merged commit ad9ce6b into OSGeo:main Oct 7, 2026
27 checks passed

github-actions Bot added this to the 8.6.0 milestone Oct 7, 2026

This was referenced Oct 7, 2026

lib/gis: clarify intentional overflow in 48-bit LCG arithmetic #8013

Merged

r.sim: Rename the option random_seed to seed #8032

Merged

libgis: Add a function which reads the seed option of a tool wenzeslaus/grass#75

Closed

wenzeslaus mentioned this pull request Oct 7, 2026

libgis: Add a seed option parsing function #8033

Merged

wenzeslaus mentioned this pull request Oct 12, 2026

libgis: Add reentrant drand48 for independent parallel streams #7884

Closed

Navigation Menu

Repository navigation

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

r.sim: Give each walker random numbers of its own#7980
wenzeslaus merged 36 commits into
OSGeo:mainfrom
wenzeslaus:rsim-layouts

wenzeslaus commented Sep 30, 2026 •

edited

Loading

Uh oh!

wenzeslaus commented Sep 30, 2026

Uh oh!

wenzeslaus commented Oct 1, 2026

Builds

Protocol

Results

State in the walker record

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

wenzeslaus commented Sep 30, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Design

The margin of the bound

Results and speed

Tests

Suggested squash commit message

Uh oh!

wenzeslaus commented Sep 30, 2026

Uh oh!

wenzeslaus commented Oct 1, 2026

Builds

Protocol

Results

State in the walker record

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

wenzeslaus commented Sep 30, 2026 •

edited

Loading