mean_source_ratemean_infiltrationthreadsoutputsoutputs
-Each entry of outputs contains the simulated_time (s) and timestamp
-of the written maps, the number of walkers_remaining at that time, and
+Each entry of outputs contains the simulated_time (s) when the maps
+were written, their timestamp, the number of walkers_remaining at that time, and
the names of the depth, discharge, error and walkers maps, or
null for maps which were not requested.
diff --git a/raster/r.sim/r.sim.water/r.sim.water.md b/raster/r.sim/r.sim.water/r.sim.water.md
index 692f0d06c90..ec2fded33af 100644
--- a/raster/r.sim/r.sim.water/r.sim.water.md
+++ b/raster/r.sim/r.sim.water/r.sim.water.md
@@ -80,7 +80,17 @@ water depth and discharge maps can be saved during simulation using the
time series flag **-t** and **output_step** parameter defining the time
step in minutes for writing output files. Files are saved with a suffix
representing time since the start of simulation in minutes (e.g.
-wdepth.05, wdepth.10). Monitoring of water depth at specific points is
+wdepth.05, wdepth.10) and are timestamped with that time. The
+simulation advances in time steps which usually do not fall exactly on
+the output times. A map holds the state at the time step closest to the
+time in its name, so at most half a time step earlier or later. When the
+time step is longer than **output_step**, there are fewer time steps
+than output times, and a time step writes only the maps for the output
+time closest to it. The series always ends with maps named by the
+**duration** which hold the state at the end of the run, also when the
+duration is not a multiple of **output_step** or when the simulation
+stopped early.
+Monitoring of water depth at specific points is
supported. A vector map with observation points and a path to a logfile
must be provided. For each point in the vector map which is located in
the computational region the water depth is logged each time step in the
@@ -197,12 +207,12 @@ stored in the history of the output raster maps under the same keys (see
| `mean_source_rate` | Mean rainfall excess | m/s |
| `mean_infiltration` | Mean infiltration rate, 0 without infiltration input | m/s |
| `threads` | Threads used for the computation | count |
-| `outputs` | One entry per set of written maps: one per **output_step** with **-t**, otherwise a single entry | |
+| `outputs` | One entry per set of written maps: with **-t**, one per written output step, the last one named by **duration**, otherwise a single entry | |
-Each entry of `outputs` contains the `simulated_time` (s) and `timestamp`
-of the written maps, the number of `walkers_remaining` at that time, and
-the names of the `depth`, `discharge`, `error` and `walkers` maps, or
-`null` for maps which were not requested.
+Each entry of `outputs` contains the `simulated_time` (s) when the maps
+were written, their `timestamp`, the number of `walkers_remaining` at
+that time, and the names of the `depth`, `discharge`, `error` and
+`walkers` maps, or `null` for maps which were not requested.
Summary of a time series run with two output steps in JSON:
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py
index a7c3f0fd74f..d62e3b5abc6 100644
--- a/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py
@@ -120,8 +120,8 @@ def test_json_summary_time_series(session_tools):
assert [output["depth"] for output in outputs] == ["depth.01", "depth.02"]
assert outputs[0]["simulated_time"] < outputs[1]["simulated_time"]
assert summary["simulated_time"] >= outputs[1]["simulated_time"]
- for output in outputs:
- assert output["timestamp"] == timestamp_for(output["simulated_time"])
+ for minute, output in enumerate(outputs, start=1):
+ assert output["timestamp"] == f"{minute} minutes"
assert output["discharge"] is None
assert session_tools.r_info(map=output["depth"], format="json")["rows"] == 5
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_test.py
index 8dcf2f02736..0e6f67a61f3 100644
--- a/raster/r.sim/r.sim.water/tests/r_sim_water_test.py
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_test.py
@@ -529,12 +529,12 @@ def test_duration_affects_time_series_progression(long_slope_session):
"""A longer duration must create more time-series output maps.
With output_step=5, duration=10 produces maps at t=5,10 while
- duration=20 produces maps at t=5,10,15 and, if walkers survive to the
- end, also t=20. The exact count of the longer run is not pinned: whether
- the final t=20 map appears depends on whether the last walkers leave the
- domain before the final step, which varies with seed and platform. We
- assert the robust lower bound (at least three maps) and that the longer
- simulation yields strictly more maps than the shorter one.
+ duration=20 produces maps at t=5,10,15,20. The exact count of the longer
+ run is not pinned: when the last walkers leave the domain early, which
+ varies with seed and platform, the steps after that are not written,
+ only the final t=20 map. We assert the robust lower bound (at least
+ three maps) and that the longer simulation yields strictly more maps
+ than the shorter one.
"""
tools = Tools(session=long_slope_session)
@@ -573,8 +573,8 @@ def test_duration_affects_time_series_progression(long_slope_session):
f"10-min simulation with output_step=5 should produce 2 maps, got {len(maps_10)}"
)
- # 20-min run should produce at least 3 time-series maps (t=5, t=10, t=15;
- # t=20 as well when the last walkers survive to the final step).
+ # 20-min run should produce at least 3 time-series maps (t=5, t=10 and
+ # t=20, and t=15 as well when the last walkers stay past it).
maps_20 = list(tools.g_list(type="raster", pattern="depth_20min*", format="json"))
assert len(maps_20) >= 3, (
f"20-min simulation with output_step=5 should produce at least 3 maps, "
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_time_series_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_time_series_test.py
new file mode 100644
index 00000000000..dbcdd628559
--- /dev/null
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_time_series_test.py
@@ -0,0 +1,169 @@
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+"""Tests of when and under which names r.sim.water writes time series maps"""
+
+import os
+
+import pytest
+
+import grass.script as gs
+from grass.experimental.mapset import TemporaryMapsetSession
+from grass.tools import Tools
+
+# On this slope, the time step is about 10.7 s, so a 1-minute output step
+# does not fall on an iteration and the old schedule, which wrote every
+# (int)(60 / 10.7) = 5 iterations and named maps by the rounded simulated
+# time, wrote dep.04 twice.
+SLOPE = {
+ "elevation": "elevation",
+ "man_value": 0.1,
+ "duration": 10,
+ "random_seed": 1,
+ "nprocs": 1,
+}
+
+
+@pytest.fixture(scope="module")
+def slope_project(tmp_path_factory):
+ """Project with a 50 m by 500 m plane at 10 m rising 0.1 m per row"""
+ project = tmp_path_factory.mktemp("simwe_time_series") / "simwe"
+ gs.create_project(project)
+ with gs.setup.init(project, env=os.environ.copy()) as session:
+ tools = Tools(session=session)
+ tools.g_region(s=0, n=500, w=0, e=50, res=10, flags="s")
+ tools.r_mapcalc(expression="elevation = row() * 0.1")
+ yield session
+
+
+@pytest.fixture
+def slope_tools(slope_project):
+ with (
+ TemporaryMapsetSession(env=slope_project.env) as session,
+ Tools(session=session, capture_stderr=True) as tools,
+ ):
+ yield tools
+
+
+def minutes_of_timestamp(tools, name):
+ """Minutes of a raster timestamp such as '1 minute' or '10 minutes'"""
+ number, unit = tools.r_timestamp(map=name).text.split()
+ assert unit in {"minute", "minutes"}
+ return int(number)
+
+
+def check_series(tools, summary, minutes):
+ """Series maps are named and timestamped by the given minutes.
+
+ Each map is written at the iteration closest to its time, except the
+ last one which is written at the end of the run.
+ """
+ outputs = summary["outputs"]
+ names = [f"dep.{minute:02d}" for minute in minutes]
+ assert [output["depth"] for output in outputs] == names
+ listed = [item["name"] for item in tools.g_list(type="raster", format="json")]
+ assert sorted(listed) == sorted([*names, "elevation"])
+ for minute, output in zip(minutes[:-1], outputs[:-1], strict=True):
+ assert output["timestamp"] == f"{minute} minutes"
+ assert minutes_of_timestamp(tools, output["depth"]) == minute
+ assert abs(output["simulated_time"] - minute * 60) <= summary["time_step"] / 2
+ assert outputs[-1]["timestamp"] == f"{minutes[-1]} minutes"
+ assert minutes_of_timestamp(tools, names[-1]) == minutes[-1]
+ assert outputs[-1]["simulated_time"] == summary["simulated_time"]
+
+
+def test_every_output_step_written_once(slope_tools):
+ """Each output step gets one map named and timestamped by its time"""
+ summary = slope_tools.r_sim_water(
+ **SLOPE,
+ depth="dep",
+ walkers_output="walkers",
+ output_step=1,
+ flags="tp",
+ format="json",
+ ).json
+ check_series(slope_tools, summary, list(range(1, 11)))
+ # Walker maps follow the same schedule with an underscore separator.
+ walkers = [output["walkers"] for output in summary["outputs"]]
+ assert walkers == [f"walkers_{minute:02d}" for minute in range(1, 11)]
+ listed = slope_tools.g_list(type="vector", format="json")
+ assert sorted(item["name"] for item in listed) == walkers
+
+
+def test_final_step_named_by_duration(slope_tools):
+ """Duration which is not a multiple of output_step still ends the series"""
+ summary = slope_tools.r_sim_water(
+ **SLOPE, depth="dep", output_step=3, flags="tp", format="json"
+ ).json
+ check_series(slope_tools, summary, [3, 6, 9, 10])
+
+
+@pytest.mark.parametrize(
+ ("mintimestep", "time_step", "minutes"),
+ [
+ (17.5, 70, [1, 2, 4, 5, 6, 7, 8, 9, 10]),
+ (30, 120, [2, 4, 6, 8, 10]),
+ (45, 180, [3, 6, 9, 10]),
+ ],
+)
+def test_time_step_longer_than_output_step(
+ slope_tools, mintimestep, time_step, minutes
+):
+ """An iteration writes only the output step closest to it
+
+ The simulated time per iteration is the minimum time step times the
+ time coefficient 4. With 70 s, the iteration ending at 210 s is equally
+ close to 3 and 4 minutes and writes the later one. With 120 s and 180 s,
+ the iterations fall exactly on every second and third output step. The
+ old schedule wrote every (int)(60 / time step) = 0 iterations, that is,
+ no maps at all.
+ """
+ result = slope_tools.r_sim_water(
+ **SLOPE,
+ depth="dep",
+ output_step=1,
+ mintimestep=mintimestep,
+ flags="tp",
+ format="json",
+ )
+ summary = result.json
+ assert summary["time_step"] == time_step
+ check_series(slope_tools, summary, minutes)
+ assert "longer than output_step" in result.stderr
+
+
+def test_final_step_after_early_stop(slope_tools):
+ """The series ends with the duration also when all walkers left early
+
+ With a very low Manning's n, all walkers leave the slope before the
+ first output step.
+ """
+ summary = slope_tools.r_sim_water(
+ **(SLOPE | {"man_value": 0.002, "duration": 2}),
+ depth="dep",
+ output_step=1,
+ flags="tp",
+ format="json",
+ ).json
+ assert summary["stopped_early"] is True
+ assert summary["simulated_time"] < 60
+ assert [output["depth"] for output in summary["outputs"]] == ["dep.02"]
+ assert summary["outputs"][0]["timestamp"] == "2 minutes"
+
+
+def test_no_output_step_past_duration(slope_tools):
+ """An output step past the duration is not written
+
+ With a time step of 160 s and a duration of 11 minutes, the last
+ iteration ends at 640 s. The output step closest to it is 12 minutes,
+ which is past the duration.
+ """
+ summary = slope_tools.r_sim_water(
+ **(SLOPE | {"duration": 11}),
+ depth="dep",
+ output_step=3,
+ mintimestep=40,
+ flags="tp",
+ format="json",
+ ).json
+ assert summary["time_step"] == 160
+ check_series(slope_tools, summary, [3, 6, 9, 11])
diff --git a/raster/r.sim/simlib/hydro.c b/raster/r.sim/simlib/hydro.c
index 0c03eb2476d..cfecbbfa1ec 100644
--- a/raster/r.sim/simlib/hydro.c
+++ b/raster/r.sim/simlib/hydro.c
@@ -44,6 +44,8 @@ void main_loop(const Setup *setup, const Geometry *geometry,
int iblock;
double conn = 1.0;
double addac;
+ // Time of the last time series output step written [s]
+ double series_time = 0.;
// nblock is reserved for Monte Carlo replicas. A future
// change will allow nblock > 1, give each replica an
@@ -124,12 +126,6 @@ void main_loop(const Setup *setup, const Geometry *geometry,
for (i = 1; i <= setup->miter;
i++) { /* iteration loop depending on simulation time and deltap */
G_percent(i, setup->miter, 1);
- if (setup->iterout > 0 && i % setup->iterout == 0) {
- /* nfiterw = i / iterout + 10;
- nfiterh = i / iterout + 40; */
- G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d", iblock,
- i, setup->miter, sim->nwalk, sim->nwalka);
- }
if (sim->nwalka == 0 && i > 1) {
summary->stopped_early = true;
@@ -315,13 +311,33 @@ void main_loop(const Setup *setup, const Geometry *geometry,
/* Total remaining walkers for this iteration */
sim->nwalka = nwalka;
+ // Output step j is written at the iteration closest to
+ // j * output_step and named by that time. With a time step
+ // longer than output_step, an iteration is the closest one to
+ // several steps and writes the step closest to it.
+ bool write_series = false;
+ if (settings->ts && settings->iterout > 0) {
+ double iteration_time = simulated_seconds(setup, i);
+ double step_time =
+ settings->iterout *
+ floor(iteration_time / settings->iterout + 0.5);
+ // A step more than half a time step ahead is closer to the
+ // next iteration. No step is past the duration.
+ if (step_time > series_time &&
+ step_time <=
+ iteration_time + 0.5 * time_step_seconds(setup) &&
+ step_time <= settings->timesec) {
+ series_time = step_time;
+ write_series = true;
+ }
+ }
+
/* Changes made by Soeren 8. Mar 2011 to replace the site walker
* output implementation */
/* Save all walkers located within the computational region and with
valid z coordinates */
if (outputs->outwalk != NULL &&
- (i == setup->miter ||
- (setup->iterout > 0 && i % setup->iterout == 0))) {
+ (i == setup->miter || write_series)) {
sim->nstack = 0;
for (lw = 0; lw < sim->nwalk; lw++) {
@@ -352,15 +368,18 @@ void main_loop(const Setup *setup, const Geometry *geometry,
} /* lw loop */
}
- if (settings->ts && setup->iterout > 0 && i % setup->iterout == 0) {
+ if (write_series) {
+ G_debug(2, "iblock=%d i=%d miter=%d nwalk=%d nwalka=%d", iblock,
+ i, setup->miter, sim->nwalk, sim->nwalka);
/* call output for iteration output */
if (outputs->erdep != NULL)
erod(grids->gama, setup, geometry,
grids); /* divergence of gama field */
double itime = simulated_seconds(setup, i);
- int ii = output_data(itime, conn, setup, geometry, settings,
- sim, inputs, outputs, grids, summary);
+ int ii =
+ output_data(itime, series_time, conn, setup, geometry,
+ settings, sim, inputs, outputs, grids, summary);
if (ii != 1)
G_fatal_error(_("Unable to write raster maps"));
}
@@ -453,13 +472,16 @@ void main_loop(const Setup *setup, const Geometry *geometry,
}
/* Write final maps here because we know the last time stamp here */
- if (!settings->ts) {
+ // A time series always ends with a step named by the duration, also when
+ // the duration is not a multiple of output_step or the run stopped early.
+ if (!settings->ts || series_time < settings->timesec) {
// All blocks have completed; gama is the eventual cumulative total,
// so no extrapolation is needed.
conn = 1.0;
double itime = simulated_seconds(setup, summary->iterations_completed);
- int ii = output_data(itime, conn, setup, geometry, settings, sim,
- inputs, outputs, grids, summary);
+ double name_time = settings->ts ? settings->timesec : itime;
+ int ii = output_data(itime, name_time, conn, setup, geometry, settings,
+ sim, inputs, outputs, grids, summary);
if (ii != 1)
G_fatal_error(_("Cannot write raster maps"));
}
diff --git a/raster/r.sim/simlib/input.c b/raster/r.sim/simlib/input.c
index 9987711e71d..ddb6298fbe8 100644
--- a/raster/r.sim/simlib/input.c
+++ b/raster/r.sim/simlib/input.c
@@ -339,10 +339,11 @@ int grad_check(Setup *setup, const Geometry *geometry, const Settings *settings,
(int)(settings->timesec /
(setup->deltap * setup->timec)); /* number of iterations = number
of cells to pass */
- setup->iterout =
- (int)(settings->iterout /
- (setup->deltap * setup->timec)); /* number of cells to pass for
- time series output */
+ if (settings->ts && settings->iterout > 0 &&
+ time_step_seconds(setup) > settings->iterout)
+ G_warning(_("Time step of %.2f s is longer than output_step of %d s, "
+ "some output steps may be skipped"),
+ time_step_seconds(setup), settings->iterout);
fprintf(stderr, "\n");
G_message(_("Min elevation \t= %.2f m\nMax elevation \t= %.2f m\n"), zmin,
diff --git a/raster/r.sim/simlib/output.c b/raster/r.sim/simlib/output.c
index f322ac65221..393e2b12eb7 100644
--- a/raster/r.sim/simlib/output.c
+++ b/raster/r.sim/simlib/output.c
@@ -143,12 +143,14 @@ static void write_history(const char *name, double tt, const Setup *setup,
G_write_raster_timestamp(name, timestamp);
}
-/* conn is the sequential-block extrapolation factor nblock/iblock: scales
+/* tt is the simulated time reached [s] and name_time the time [s] the maps
+ * are named and timestamped by, rounded to minutes.
+ * conn is the sequential-block extrapolation factor nblock/iblock: scales
* the cumulative partial sum in gama into an estimator of the eventual total
* for snapshots taken before all blocks have run. With nblock = 1 (or after
* the final block) conn = 1.0 and the output formulas are unchanged. err is
* written from gammas, which already has conn baked into its accumulator. */
-int output_data(double tt, double conn, const Setup *setup,
+int output_data(double tt, double name_time, double conn, const Setup *setup,
const Geometry *geometry, const Settings *settings,
const Simulation *sim, const Inputs *inputs,
const Outputs *outputs, const Grids *grids, Summary *summary)
@@ -187,7 +189,7 @@ int output_data(double tt, double conn, const Setup *setup,
ndigit = 5;
/* Convert to minutes */
- tt_minutes = (int)(tt / 60. + 0.5);
+ tt_minutes = (int)(name_time / 60. + 0.5);
/* Create timestamp */
snprintf(timestamp_buf, sizeof(timestamp_buf), "%d minutes", tt_minutes);
diff --git a/raster/r.sim/simlib/simlib.h b/raster/r.sim/simlib/simlib.h
index be11d77350b..fe71bea6012 100644
--- a/raster/r.sim/simlib/simlib.h
+++ b/raster/r.sim/simlib/simlib.h
@@ -33,14 +33,13 @@ typedef struct {
double hbeta; // Weighting factor for water flow velocity vector
double hhmax; // Threshold water depth [m]
double frac; // Water diffusion constant
- int iterout; // Time interval for creating output maps [minutes]
- int timesec; // Time how long the simulation runs [minutes]
+ int iterout; // Time interval for creating output maps [seconds]
+ int timesec; // Time how long the simulation runs [seconds]
bool ts; // Time series output
double mintimestep; // Minimum time step for the simulation [seconds]
} Settings;
typedef struct {
- int iterout; // Number of iterations for creating output maps
int miter; // Total number of iterations
double chmean; // Mean Manning's n
double si0; // Mean rainfall excess (or sediment concentration?)
@@ -177,7 +176,7 @@ void main_loop(const Setup *setup, const Geometry *geometry,
const Settings *settings, Simulation *sim,
ObservationPoints *points, const Inputs *inputs,
const Outputs *outputs, Grids *grids, Summary *summary);
-int output_data(double tt, double conn, const Setup *setup,
+int output_data(double tt, double name_time, double conn, const Setup *setup,
const Geometry *geometry, const Settings *settings,
const Simulation *sim, const Inputs *inputs,
const Outputs *outputs, const Grids *grids, Summary *summary);