r.sim: fix infiltration implementation

raster/r.sim/r.sim.water/testsuite/test_r_sim_water.py

@@ -19,6 +19,7 @@ class TestRSimWater(TestCase):
     reference_discharge_default = "discharge_default"
     reference_depth_complex = "depth_complex"
     reference_discharge_complex = "discharge_complex"
+    reference_depth_infil = "depth_infil"
 
     @classmethod
     def setUpClass(cls):
@@ -46,7 +47,7 @@ def setUpClass(cls):
         )
         cls.runModule(
             "r.mapcalc",
-            expression=f"{cls.infil} = if (x() < 636500 && y() > 223500, 0.001, 0)",
+            expression=f"{cls.infil} = if (x() > 636500 && y() > 223500, 10, 0)",
         )
         cls.runModule(
             "r.unpack",
@@ -58,6 +59,11 @@ def setUpClass(cls):
             input="data/discharge_complex.pack",
             output=cls.reference_discharge_complex,
         )
+        cls.runModule(
+            "r.unpack",
+            input="data/depth_infil.pack",
+            output=cls.reference_depth_infil,
+        )
 
     @classmethod
     def tearDownClass(cls):
@@ -154,6 +160,28 @@ def test_complex(self):
             precision="0.000001",
         )
 
+    def test_infiltration(self):
+        """Test r.sim.water execution with infiltration"""
+        # Run the r.sim.water simulation
+        self.assertModule(
+            "r.sim.water",
+            elevation=self.elevation,
+            dx=self.dx,
+            dy=self.dy,
+            rain=self.infil,
+            infil=self.infil,
+            depth=self.depth,
+            niterations=30,
+            random_seed=1,
+        )
+
+        # Assert that the output rasters exist
+        self.assertRasterExists(self.depth)
+        # Assert that the output rasters are the same
+        self.assertRastersEqual(
+            self.depth, reference=self.reference_depth_infil, precision="0.000001"
+        )
+
 
 @unittest.skip("runs too long")
 class TestRSimWaterLarge(TestCase):

raster/r.sim/simlib/hydro.c

@@ -180,23 +180,26 @@ void main_loop(const Setup *setup, const Geometry *geometry,
                         if (grids->zz[k][l] != UNDEF) {
                             if (grids->inf[k][l] !=
                                 UNDEF) { /* infiltration part */
-                                if (grids->inf[k][l] - grids->si[k][l] > 0.) {
 
-                                    double decr = pow(
-                                        addac * sim->w[lw].m,
+                                double decr =
+                                    pow(addac * sim->w[lw].m,
                                         3. / 5.); /* decreasing factor in m */
-                                    if (grids->inf[k][l] > decr) {
-                                        grids->inf[k][l] -=
-                                            decr; /* decrease infilt. in cell
-                                                     and eliminate the walker */
+                                if (grids->inf[k][l] > decr) {
+                                    grids->inf[k][l] -=
+                                        decr; /* decrease infilt. in cell
+                                                    and eliminate the walker */
+                                    sim->w[lw].m = 0.;
+                                    continue;
+                                }
+                                else {
+                                    sim->w[lw].m -= sim->rwalk *
+                                                    grids->inf[k][l] /
+                                                    setup->sisum;
+                                    grids->inf[k][l] = 0.;
+                                    // eliminate walker
+                                    if (sim->w[lw].m < 0.) {
                                         sim->w[lw].m = 0.;
-                                    }
-                                    else {
-                                        sim->w[lw].m -=
-                                            pow(grids->inf[k][l], 5. / 3.) /
-                                            addac; /* use just proportional part
-                                                      of the walker weight */
-                                        grids->inf[k][l] = 0.;
+                                        continue;
                                     }
                                 }
                             }

raster/r.sim/simlib/input.c

@@ -291,9 +291,6 @@ int grad_check(Setup *setup, const Geometry *geometry, const Settings *settings,
             } /* DEFined area */
         }
     }
-    if (grids->inf != NULL && smax < infmax)
-        G_warning(_("Infiltration exceeds the rainfall rate everywhere! No "
-                    "overland flow."));
     if (n == 0) {
         G_fatal_error(_("No values in the elevation raster."));
     }
@@ -372,8 +369,6 @@ int grad_check(Setup *setup, const Geometry *geometry, const Settings *settings,
                 /*if(v1[k][l]*v1[k][l]+v2[k][l]*v2[k][l] > cellsize, warning,
                  *napocitaj ak viac ako 10%a*/
                 /* THIS IS CORRECT SOLUTION currently commented out */
-                if (grids->inf)
-                    grids->inf[k][l] *= settings->timesec;
                 if (inputs->wdepth)
                     grids->gama[k][l] = 0.;
                 if (outputs->et) {