r.sim: Write time series maps at their output step times

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

@@ -90,7 +90,15 @@ <h2>DESCRIPTION</h2>
 the time series flag <b>-t</b> and <b>output_step</b> 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).
+(e.g. 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 <b>output_step</b>, 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 <b>duration</b> which hold the state
+at the end of the run, also when the duration is not a multiple of <b>output_step</b>
+or when the simulation stopped early.
 <!-- TODO add example of registering outputs into temporal framework when  #2146 is fixed -->
 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
@@ -209,12 +217,12 @@ <h3>Run summary</h3>
  <tr><td><code>mean_source_rate</code></td><td>Mean rainfall excess</td><td>m/s</td></tr>
  <tr><td><code>mean_infiltration</code></td><td>Mean infiltration rate, 0 without infiltration input</td><td>m/s</td></tr>
  <tr><td><code>threads</code></td><td>Threads used for the computation</td><td>count</td></tr>
- <tr><td><code>outputs</code></td><td>One entry per set of written maps: one per <b>output_step</b> with <b>-t</b>, otherwise a single entry</td><td></td></tr>
+ <tr><td><code>outputs</code></td><td>One entry per set of written maps: with <b>-t</b>, one per written output step, the last one named by <b>duration</b>, otherwise a single entry</td><td></td></tr>
 </table>
 
 <p>
-Each entry of <code>outputs</code> contains the <code>simulated_time</code> (s) and <code>timestamp</code>
-of the written maps, the number of <code>walkers_remaining</code> at that time, and
+Each entry of <code>outputs</code> contains the <code>simulated_time</code> (s) when the maps
+were written, their <code>timestamp</code>, the number of <code>walkers_remaining</code> at that time, and
 the names of the <code>depth</code>, <code>discharge</code>, <code>error</code> and <code>walkers</code> maps, or
 <code>null</code> for maps which were not requested.
 

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:
 

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
 

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, "

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])

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"));
     }

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,