From 381433aa95d4868648801082e87c603b533754d9 Mon Sep 17 00:00:00 2001
From: Markus Metz <33666869+metzm@users.noreply.github.com>
Date: Fri, 31 Jan 2025 09:22:57 +0100
Subject: [PATCH] r.sim.water: add info to manual, add info about progress
 (#4662)

 * `r.sim.water` is a module with a number of important but not well documented parameters. This PR adds more information about some parameters to the manual.
 * It takes quite some time to run a simulation, another message is added to the progress of the simulation.

This PR tries to add info requested in https://github.com/OSGeo/grass/discussions/4475 to the manual


Co-authored-by: Anna Petrasova <kratochanna@gmail.com>
Co-authored-by: Markus Neteler <neteler@osgeo.org>
---
 raster/r.sim/r.sim.sediment/main.c            | 15 +++++-
 .../r.sim/r.sim.sediment/r.sim.sediment.html  |  4 ++
 raster/r.sim/r.sim.water/main.c               | 15 +++++-
 raster/r.sim/r.sim.water/r.sim.water.html     | 51 +++++++++++++++++++
 raster/r.sim/simlib/hydro.c                   |  3 ++
 raster/r.sim/simlib/input.c                   |  5 ++
 raster/r.sim/simlib/simlib.h                  |  9 ++--
 raster/r.sim/simlib/waterglobs.h              |  1 +
 8 files changed, 97 insertions(+), 6 deletions(-)

diff --git a/raster/r.sim/r.sim.sediment/main.c b/raster/r.sim/r.sim.sediment/main.c
index d0197c2abb1..f49fa8fb880 100644
--- a/raster/r.sim/r.sim.sediment/main.c
+++ b/raster/r.sim/r.sim.sediment/main.c
@@ -237,6 +237,18 @@ int main(int argc, char *argv[])
     parm.niter->description = _("Time used for iterations [minutes]");
     parm.niter->guisection = _("Parameters");
 
+    parm.mintimestep = G_define_option();
+    parm.mintimestep->key = "mintimestep";
+    parm.mintimestep->type = TYPE_DOUBLE;
+    parm.mintimestep->answer = "0.0";
+    parm.mintimestep->required = NO;
+    parm.mintimestep->label =
+        _("Minimum time step for the simulation [seconds]");
+    parm.mintimestep->description =
+        _("A larger minimum time step substantially reduces processing time, "
+          "but at the cost of accuracy");
+    parm.mintimestep->guisection = _("Parameters");
+
     parm.outiter = G_define_option();
     parm.outiter->key = "output_step";
     parm.outiter->type = TYPE_INTEGER;
@@ -286,7 +298,7 @@ int main(int argc, char *argv[])
     parm.threads->answer = NUM_THREADS;
     parm.threads->required = NO;
     parm.threads->description =
-        _("Number of threads which will be used for parallel compute");
+        _("Number of threads which will be used for parallel computation.");
     parm.threads->guisection = _("Parameters");
 
     if (G_parser(argc, argv))
@@ -389,6 +401,7 @@ int main(int argc, char *argv[])
 
     /*      sscanf(parm.nwalk->answer, "%d", &wp.maxwa); */
     sscanf(parm.niter->answer, "%d", &wp.timesec);
+    sscanf(parm.mintimestep->answer, "%lf", &wp.mintimestep);
     sscanf(parm.outiter->answer, "%d", &wp.iterout);
     /*    sscanf(parm.density->answer, "%d", &wp.ldemo); */
     sscanf(parm.diffc->answer, "%lf", &wp.frac);
diff --git a/raster/r.sim/r.sim.sediment/r.sim.sediment.html b/raster/r.sim/r.sim.sediment/r.sim.sediment.html
index caebb9c308c..7f1b80867ca 100644
--- a/raster/r.sim/r.sim.sediment/r.sim.sediment.html
+++ b/raster/r.sim/r.sim.sediment/r.sim.sediment.html
@@ -41,6 +41,10 @@ <h2>DESCRIPTION</h2>
 given by <i>nwalkers</i> should be used.
 <br>
 
+<p>
+Increasing the number of threads with <b>nprocs</b> does not really
+speed up the simulation.
+
 <!--
 <h2>NOTES</h2>
 
diff --git a/raster/r.sim/r.sim.water/main.c b/raster/r.sim/r.sim.water/main.c
index 209be3e33d9..d4a545172a5 100644
--- a/raster/r.sim/r.sim.water/main.c
+++ b/raster/r.sim/r.sim.water/main.c
@@ -239,6 +239,18 @@ int main(int argc, char *argv[])
     parm.niter->description = _("Time used for iterations [minutes]");
     parm.niter->guisection = _("Parameters");
 
+    parm.mintimestep = G_define_option();
+    parm.mintimestep->key = "mintimestep";
+    parm.mintimestep->type = TYPE_DOUBLE;
+    parm.mintimestep->answer = "0.0";
+    parm.mintimestep->required = NO;
+    parm.mintimestep->label =
+        _("Minimum time step for the simulation [seconds]");
+    parm.mintimestep->description =
+        _("A larger minimum time step substantially reduces processing time, "
+          "but at the cost of accuracy");
+    parm.mintimestep->guisection = _("Parameters");
+
     parm.outiter = G_define_option();
     parm.outiter->key = "output_step";
     parm.outiter->type = TYPE_INTEGER;
@@ -320,7 +332,7 @@ int main(int argc, char *argv[])
     parm.threads->answer = NUM_THREADS;
     parm.threads->required = NO;
     parm.threads->description =
-        _("Number of threads which will be used for parallel compute");
+        _("Number of threads which will be used for parallel computation.");
     parm.threads->guisection = _("Parameters");
 
     if (G_parser(argc, argv))
@@ -388,6 +400,7 @@ int main(int argc, char *argv[])
     G_debug(3, "Parsing numeric parameters");
 
     sscanf(parm.niter->answer, "%d", &wp.timesec);
+    sscanf(parm.mintimestep->answer, "%lf", &wp.mintimestep);
     sscanf(parm.outiter->answer, "%d", &wp.iterout);
     sscanf(parm.diffc->answer, "%lf", &wp.frac);
     sscanf(parm.hmax->answer, "%lf", &wp.hhmax);
diff --git a/raster/r.sim/r.sim.water/r.sim.water.html b/raster/r.sim/r.sim.water/r.sim.water.html
index af54dea4490..e3965c27008 100644
--- a/raster/r.sim/r.sim.water/r.sim.water.html
+++ b/raster/r.sim/r.sim.water/r.sim.water.html
@@ -25,6 +25,19 @@ <h2>DESCRIPTION</h2>
 can be used to compute partial derivates of the predefined flow using its direction given
 by aspect and slope.
 
+<p>
+The equations are
+
+<div class="code"><pre>
+dx = tan(slope) * cos(aspect)
+</pre></div>
+
+and
+
+<div class="code"><pre>
+dy = tan(slope) * sin(aspect)
+</pre></div>
+
 <p>
 <div align="center" style="margin: 10px;">
 <img style="margin: 0.5em;" src="r_sim_water.png" alt="r.sim.water generated depth map"><br>
@@ -99,6 +112,13 @@ <h2>DESCRIPTION</h2>
 diffusion term increases as given by <b>halpha</b> and advection term
 (direction of flow) is given as "prevailing" direction of flow computed
 as average of flow directions from the previous <b>hbeta</b> number of grid cells.
+The model tries to keep water "shallow" with maximum shallow water
+depth defined by <b>hmax</b> default 0.3 meters. However, water depths
+much higher than <b>hmax</b> can be observed if water accumulates in
+natural sinks or river beds. Depending on the area of interest and the
+used digital elevation model, <b>hmax</b>, <b>halpha</b> and
+<b>hbeta</b> might need to be adjusted in order to deal realistically
+with elevation depressions or obstacles.
 
 <h2>NOTES</h2>
 
@@ -146,6 +166,37 @@ <h2>NOTES</h2>
 are useful both for everyday exploratory work using a desktop computer and
 for large, cutting-edge applications using high performance computing.
 
+<p>
+<b>Suggested Manning's n for surface roughness</b><br>
+from <a href="https://baharmon.github.io/hydrology-in-grass">https://baharmon.github.io/hydrology-in-grass</a>
+<p>
+<table border="1">
+ <tr><th>NLCD Landcover Category</th><th>Manning’s n value</th></tr>
+ <tr><td>Open Water</td><td>0.001</td></tr>
+ <tr><td>Developed, Open Space</td><td>0.0404</td></tr>
+ <tr><td>Developed, Low Intensity</td><td>0.0678</td></tr>
+ <tr><td>Developed, Medium Intensity</td><td>0.0678</td></tr>
+ <tr><td>Developed, High Intensity</td><td>0.0404</td></tr>
+ <tr><td>Barren Land</td><td>0.0113</td></tr>
+ <tr><td>Deciduous Forest</td><td>0.36</td></tr>
+ <tr><td>Evergreen Forest</td><td>0.32</td></tr>
+ <tr><td>Mixed Forest</td><td>0.4</td></tr>
+ <tr><td>Shrub/Scrub</td><td>0.4</td></tr>
+ <tr><td>Grassland/Herbaceuous</td><td>0.368</td></tr>
+ <tr><td>Pasture/Hay</td><td>0.325</td></tr>
+ <tr><td>Cultivated Crops</td><td>0.325</td></tr>
+ <tr><td>Woody Wetlands</td><td>0.086</td></tr>
+ <tr><td>Emergent Herbaceuous Wetlands</td><td>0.1825</td></tr>
+</table>
+
+<p>
+The <a href="https://www.usgs.gov/centers/eros/science/annual-nlcd-science-product-user-guide">NLCD user guide</a>
+provides more information about the different NLCD classes.
+
+<p>
+Increasing the number of threads with <b>nprocs</b> does not really
+speed up the simulation.
+
 <h2>EXAMPLE</h2>
 
 Using the North Carolina full sample dataset:
diff --git a/raster/r.sim/simlib/hydro.c b/raster/r.sim/simlib/hydro.c
index eb7649763e9..d8eea788641 100644
--- a/raster/r.sim/simlib/hydro.c
+++ b/raster/r.sim/simlib/hydro.c
@@ -104,6 +104,7 @@ int iterout, mx2o, my2o;
 int miter, nwalka;
 double timec;
 int ts, timesec;
+double mintimestep;
 
 double rain_val;
 double manin_val;
@@ -156,6 +157,8 @@ void main_loop(void)
     for (iblock = 1; iblock <= nblock; iblock++) {
         /* ++icoub; */
 
+        G_message(_("Processing block %d of %d"), iblock, nblock);
+
         lw = 0;
         walkwe = 0.;
         barea = stepx * stepy;
diff --git a/raster/r.sim/simlib/input.c b/raster/r.sim/simlib/input.c
index e9f41b88709..44a948d0244 100644
--- a/raster/r.sim/simlib/input.c
+++ b/raster/r.sim/simlib/input.c
@@ -95,6 +95,7 @@ void WaterParams_init(struct WaterParams *wp)
     wp->timec = 0;
     wp->ts = 0;
     wp->timesec = 0;
+    wp->mintimestep = 0;
 
     wp->rain_val = 0;
     wp->manin_val = 0;
@@ -201,6 +202,7 @@ void init_library_globals(struct WaterParams *wp)
     timec = wp->timec;
     ts = wp->ts;
     timesec = wp->timesec;
+    mintimestep = wp->mintimestep;
 
     rain_val = wp->rain_val;
     manin_val = wp->manin_val;
@@ -508,6 +510,9 @@ int grad_check(void)
         deltaw = 0.8 / (sigmax * vmax);          /*time step for sediment */
     deltap = 0.25 * sqrt(stepx * stepy) / vmean; /*time step for water */
 
+    if (deltap < mintimestep)
+        deltap = mintimestep;
+
     if (deltaw > deltap)
         timec = 4.;
     else
diff --git a/raster/r.sim/simlib/simlib.h b/raster/r.sim/simlib/simlib.h
index 4a0fe592bb7..2b651bb5503 100644
--- a/raster/r.sim/simlib/simlib.h
+++ b/raster/r.sim/simlib/simlib.h
@@ -35,6 +35,7 @@ struct WaterParams {
     int miter, nwalka;
     double timec;
     int ts, timesec;
+    double mintimestep;
 
     double rain_val;
     double manin_val;
@@ -86,10 +87,10 @@ void free_walkers(void);
 
 struct options {
     struct Option *elevin, *dxin, *dyin, *rain, *infil, *traps, *manin,
-        *observation, *depth, *disch, *err, *outwalk, *nwalk, *niter, *outiter,
-        *density, *diffc, *hmax, *halpha, *hbeta, *wdepth, *detin, *tranin,
-        *tauin, *tc, *et, *conc, *flux, *erdep, *rainval, *maninval, *infilval,
-        *logfile, *seed, *threads;
+        *observation, *depth, *disch, *err, *outwalk, *nwalk, *niter,
+        *mintimestep, *outiter, *density, *diffc, *hmax, *halpha, *hbeta,
+        *wdepth, *detin, *tranin, *tauin, *tc, *et, *conc, *flux, *erdep,
+        *rainval, *maninval, *infilval, *logfile, *seed, *threads;
 };
 
 struct flags {
diff --git a/raster/r.sim/simlib/waterglobs.h b/raster/r.sim/simlib/waterglobs.h
index 7bf9fab74f2..10acce39c5b 100644
--- a/raster/r.sim/simlib/waterglobs.h
+++ b/raster/r.sim/simlib/waterglobs.h
@@ -111,6 +111,7 @@ extern int iterout, mx2o, my2o;
 extern int miter, nwalka;
 extern double timec;
 extern int ts, timesec;
+extern double mintimestep;
 
 extern double rain_val;
 extern double manin_val;