Feature/performance

.gitignore

@@ -34,3 +34,9 @@ lisflood_optimization_report.md
 lisflood_optimization_report_v2.md
 # version file auto-generated by setuptools_scm (do not track)
 src/lisflood/_version.py
+# output folders generated by the test suite / example runs
+tests/data/**/out/
+# generated profiling artifacts 
+profiling/
+# debug variants of source modules (not part of the package)
+*_debug.py

.zenodo.json

@@ -81,9 +81,9 @@
       "orcid": "0000-0002-2430-3634"
     },
     {
-        "name": "Jensen, Laura",
-        "affiliation": "GFZ Helmholtz Centre for Geosciences",
-        "orcid": "0000-0002-8318-791X"
+      "name": "Jensen, Laura",
+      "affiliation": "GFZ Helmholtz Centre for Geosciences",
+      "orcid": "0000-0002-8318-791X"
     },
     {
       "name": "Morrison, Oisin",

docs/3_step2_preparing-setting-file/index.md

@@ -757,6 +757,25 @@ Within the 'lfoptions' element of the settings file, each option is defined usin
 + **OutputMapsChunks**: this option is used to write output maps every X steps (default 1).  
 + **OutputMapsDataType**: this option sets the output data type and may take the following values: "float64" or "float32" (default float64). This option applies to all output maps.
 
+### Options to control parallelization
+
+LISFLOOD uses several independent parallelization mechanisms. They all affect run time only, never results: model outputs are identical regardless of the values below. Each option accepts a positive integer (number of threads/processes), or `0` / `"all"` / `"auto"` meaning "use all available CPU cores".
+
+These are `lfuser`/`lfbinding` variables (add them like any other `textvar`), for example:
+
+```xml
+<textvar name="numCPUs_parallelNumba" value="0"/>
+<textvar name="numCPUs_parallelNumexpr" value="auto"/>
+<textvar name="numCPUs_BLAS" value="auto"/>
+<textvar name="numCPUs_soilInit" value="0"/>
+```
+
++ **numCPUs_parallelNumba**: number of threads used by the Numba just-in-time compiled kernels (the kinematic-wave / Muskingum-Cunge routing and the soil/vegetation loop). Value `0` uses all cores (or the `NUMBA_NUM_THREADS` environment variable if set). On small catchments, thread fork/join overhead can make multi-threaded execution slower than serial, so LISFLOOD auto-tunes this to 1 (serial) for small domains when left on `0`; large domains benefit from more threads.
++ **numCPUs_parallelNumexpr**: number of threads used by the `numexpr` library (vectorized array maths in the kinematic-wave routing). Default `auto`: a domain-size-aware choice that runs serial on small catchments (where `numexpr`'s pool gives no speedup and oversubscribes against Numba/BLAS) and uses all cores on large/global domains (where it is measurably faster). You can force a value: `1` (serial), `0`/`all` (all cores), or a specific number. The large/small cutoff is the same valid-pixel threshold used by the Numba auto-tune. **Set this to `1` under an outer parallel driver** (e.g. [lisflood-calibration](https://github.com/ec-jrc/lisflood-calibration), one LISFLOOD per DEAP individual) so each instance stays single-core.
++ **numCPUs_BLAS**: number of threads used by the BLAS/OpenMP backend of NumPy/SciPy (e.g. the linear algebra in the soil-hydraulics fitting during initialization). Default `auto`, with the same domain-size-aware behaviour and overrides as `numCPUs_parallelNumexpr`; set it to `1` under an outer parallel driver.
++ **numCPUs_soilInit**: number of worker **processes** used to solve the soil layer-2 saturation during a cold start. This is a per-pixel SciPy solve that can dominate initialization time on large/global domains (turning it from serial minutes/hours into a few minutes with many cores). This is process-based parallelism, independent of the thread options above. Default `0` (use all available cores); the solve is bit-identical regardless of the worker count.
+    - **Set `numCPUs_soilInit` to `1` when an outer tool already runs many LISFLOOD instances in parallel** — most importantly [lisflood-calibration](https://github.com/ec-jrc/lisflood-calibration), which runs one LISFLOOD per DEAP individual and expects each instance to use a single core. As a safeguard, this parallelization automatically disables itself (falls back to serial) when LISFLOOD detects it is running inside a forked/daemon worker, so an accidental high value cannot cause oversubscription under calibration.
+
 ### Reference settings file
 In order to facilitate the preparation of the settings file, a complete example is provided [here](https://github.com/ec-jrc/lisflood-code/tree/master/src/lisfloodSettings_reference.xml). The user is encouraged to update the paths, the names of the maps and of the tables in the provided template.
 

pyproject.toml

@@ -49,6 +49,7 @@ test = [
     "pytest>=6.2.5",
     "pytest-cov>=3.0.0",
     "pytest-mock>=3.6.1",
+    "pytest-xdist>=3.0.0",
     "coverage>=6.0",
 ]
 

pytest.ini

@@ -14,4 +14,15 @@ env =
     D:NUMEXPR_MAX_THREADS=12
 ;    D:PCRASTER_NR_WORKER_THREADS=4 gives errors with pcraster.multicore - need changes in many hydro modules
 ; run pytest tests/ -m "slow" before release !!!
-addopts = -ra -x -l --cov=lisflood --cov-config=.coveragerc -m "slow or not slow" -vv
+; Parallel test options (default): run the suite across CPU cores with pytest-xdist
+; while still collecting coverage (pytest-cov combines per-worker coverage into one
+; report). "-n auto" uses all available cores; "--dist loadscope" keeps every test
+; in a class on the same worker, which is required for safe parallelism here because
+; tests reuse per-class output dirs and different files share 'out/' leaf names
+; (cross-worker collisions are avoided by the per-worker output isolation in
+; tests/test_utils.py). "-x" is omitted because stop-on-first-failure is racy under
+; xdist. This runs both locally and in CI at full speed with coverage.
+addopts = -ra -l --cov=lisflood --cov-config=.coveragerc -m "slow or not slow" -vv -n auto --dist loadscope
+; Serial fallback (uncomment this line and comment the parallel one above) - useful
+; for debugging a single test or when pytest-xdist is unavailable; keeps -x.
+; addopts = -ra -x -l --cov=lisflood --cov-config=.coveragerc -m "slow or not slow" -vv

src/lisflood/Lisflood_initial.py

@@ -21,20 +21,14 @@
 from collections import OrderedDict
 
 import numpy as np
-try:
-    from numba import njit, set_num_threads, get_num_threads, config as numba_config
-    numba_full = True
-except ImportError:
-    import numba
-    print('WARNING! Numba version is probably too old. Only limited functionalities will be available for Numba')
-    print('Numba version: {}'.format(numba.__version__))
-    from numba import njit, config as numba_config
-    numba_full = False
+from numba import njit
+# Thread-pool governance (numba/numexpr/BLAS) lives in parallelization.py.
 
 
 from .global_modules.settings import CutMap, LisSettings, NetCDFMetadata, EPICSettings, MaskInfo
 from .global_modules.zusatz import DynamicModel
 from .global_modules.add1 import loadsetclone, mapattrNetCDF
+from .global_modules.parallelization import configure_parallelism
 from .hydrological_modules.miscInitial import miscInitial
 
 from .hydrological_modules.readmeteo import readmeteo
@@ -98,11 +92,16 @@ def __init__(self):
         flags = self.settings.flags
         report_steps = self.settings.report_steps
 
-        # set the maximum number of threads that numba should use (now used in soilloop only)
-        num_threads = int(binding["numCPUs_parallelNumba"])
-        if (num_threads>0):
-            if num_threads<=numba_config.NUMBA_NUM_THREADS:
-                set_num_threads(num_threads)
+        # Centrally govern the numba, numexpr and BLAS/OpenMP thread pools to
+        # avoid oversubscription. Nested thread pools each sized to the host
+        # core count cause heavy context-switching and can make runs slower
+        # than serial, especially on small/narrow domains. Thread counts affect
+        # performance only, not results.
+        num_pixels = int(self.maskinfo.info.mapC[0]) if self.maskinfo.info.mapC else None
+        # verbose=False: the effective thread counts are reported once, in the
+        # LisfloodRunInfo startup banner (see LisfloodRunInfo in settings.py),
+        # to keep the startup output in a single coherent block.
+        configure_parallelism(binding, num_pixels=num_pixels)
 
         # Mapping of vegetation types to land use fractions (and the other way around)
         ##global VEGETATION_LANDUSE, LANDUSE_VEGETATION, PRESCRIBED_VEGETATION, PRESCRIBED_LAI

src/lisflood/global_modules/add1.py

@@ -291,25 +291,39 @@ def compressArray(map, pcr=True, name=None, force_load_with_nans=False):
 
 def decompress(map):
     maskinfo = MaskInfo.instance()
-    dmap = maskinfo.info.maskall.copy()
-    dmap[~maskinfo.info.maskflat] = map[:]
-    dmap = dmap.reshape(maskinfo.info.shape)
     # check if integer map (like outlets, lakes etc)
     try:
         checkint = str(map.dtype)
     except:
         checkint = None
 
-    if checkint in ("int16", "int32", "int64"):
-        dmap[dmap.mask] = -9999
-        map = numpy2pcr(Nominal, dmap, -9999)
-    elif checkint == "int8":
+    if checkint == "int8":
+        # Rare, special-cased branch: the original relied on np.ma semantics
+        # (dmap[dmap < 0] excludes masked cells, leaving their uninitialized
+        # masked_all .data untouched). Keep the exact original implementation
+        # here to preserve behaviour bit-for-bit; it is not on any hot path
+        # (compressArray returns float, so decompress almost always hits the
+        # Scalar branch below).
+        dmap = maskinfo.info.maskall.copy()
+        dmap[~maskinfo.info.maskflat] = map[:]
+        dmap = dmap.reshape(maskinfo.info.shape)
         dmap[dmap < 0] = -9999
-        map = numpy2pcr(Nominal, dmap, -9999)
-    else:
-        dmap[dmap.mask] = -9999
-        map = numpy2pcr(Scalar, dmap, -9999)
-    return map
+        return numpy2pcr(Nominal, dmap, -9999)
+
+    # Fast path (float / int16 / int32 / int64): build the full-domain array
+    # with a plain numpy buffer instead of copying the masked `maskall` and
+    # operating on it through the np.ma machinery. `maskall` is float64 and the
+    # result handed to numpy2pcr is simply: valid cells = map values,
+    # outside-domain cells = -9999 (the missing-value marker). This avoids the
+    # full-domain masked-array copy and the np.ma wrapped-method overhead that
+    # dominate decompress on large/global domains. Verified bit-identical to the
+    # original masked-array path for all these dtype branches.
+    dmap = np.full(maskinfo.info.shapeflat, -9999.0, dtype=np.float64)
+    dmap[~maskinfo.info.maskflat] = map[:]
+    dmap = dmap.reshape(maskinfo.info.shape)
+    if checkint in ("int16", "int32", "int64"):
+        return numpy2pcr(Nominal, dmap, -9999)
+    return numpy2pcr(Scalar, dmap, -9999)
 
 
 def makenumpy(map):
@@ -322,6 +336,61 @@ def makenumpy(map):
         return map
 
 
+# Cache of parsed 2-column lookup tables: {abspath: {id(int): value(float)}}.
+_LOOKUP_TABLE_CACHE = {}
+
+
+def _read_lookup_table(table_path):
+    """Parse a PCRaster-style 2-column 'key value' lookup table into an
+    {int(key): float(value)} dict. Only the simple exact-match form used by the
+    reservoir/lake parameter tables is supported (no range/interval rows).
+    Cached by absolute path."""
+    key = os.path.abspath(str(table_path))
+    cached = _LOOKUP_TABLE_CACHE.get(key)
+    if cached is not None:
+        return cached
+    mapping = {}
+    with open(str(table_path)) as f:
+        for line in f:
+            parts = line.split()
+            if len(parts) < 2:
+                continue
+            # keys are integer point ids (reservoir/lake ids); skip non-numeric
+            # (e.g. accidental header lines) rather than crash.
+            try:
+                k = int(float(parts[0]))
+                v = float(parts[1])
+            except ValueError:
+                continue
+            mapping[k] = v
+    _LOOKUP_TABLE_CACHE[key] = mapping
+    return mapping
+
+
+def lookup_table_values(table_path, point_ids):
+    """Vectorized equivalent of
+        np.compress(sites>0, compressArray(pcraster.lookupscalar(table, siteMapPcr)))
+    for a simple exact-match 2-column lookup table.
+
+    For each id in `point_ids` (the compressed point ids, e.g.
+    ReservoirSitesCC / LakeSitesCC) return the table value for that id, or NaN if
+    the id is absent from the table -- exactly as pcraster.lookupscalar returns a
+    missing value (which compressArray turns into NaN). This replaces a full-grid
+    PCRaster table lookup (~30 s/call on the global domain) with an O(n_points)
+    dict lookup.
+    """
+    mapping = _read_lookup_table(table_path)
+    ids = np.asarray(point_ids)
+    out = np.empty(ids.shape[0], dtype=float)
+    for i in range(ids.shape[0]):
+        out[i] = mapping.get(int(ids[i]), np.nan)
+    # Cast to float32 then back, matching pcraster.lookupscalar's Scalar (float32)
+    # precision. Without this, the full float64 table values introduce ~1e-7
+    # per-value differences (visible as ~1e-4 in discharge after amplification
+    # through reservoir/lake routing), breaking bit-identical output.
+    return out.astype(np.float32).astype(np.float64)
+
+
 def loadmap(*args, **kwargs):
     settings = LisSettings.instance()
     binding = settings.binding

src/lisflood/global_modules/netcdf.py

@@ -271,10 +271,11 @@ def load_next_chunk(self):
         # handle (NetCDF4DataStore.__del__ -> close). Concurrent read/close on a
         # non-threadsafe HDF5 corrupts its global state and aborts the process
         # (Fatal Python error: Aborted / SIGABRT). Reads were already serialized
-        # by xarray's HDF5 lock, so running synchronously removes the crash with
-        # negligible performance impact.
-        with dask.config.set(scheduler='synchronous'):
-            self.dataset_chunk = chunk.load()  # triggers xarray computation
+        # by xarray's HDF5 lock, so running synchronously removes the crash.
+        # TODO: check the synchronous issue, as actually this slows down the loading of chunks
+        #       so it is commented right now
+        #with dask.config.set(scheduler='synchronous'):
+        self.dataset_chunk = chunk.load()  # triggers xarray computation
 
     def __getitem__(self, step):
 

src/lisflood/global_modules/output.py

@@ -649,23 +649,32 @@ def initial(self):
         # output for single column eg mapmaximum
         self.var.Tss = {}
 
+        # Resolve each distinct output-point location ("where": Gauges, Sites,
+        # LakeSites, Catchments, ...) only ONCE. The resolved map is used read-only by
+        # TimeoutputTimeseries.__init__, so sharing one object per `where` is safe.
+        outpoints_by_where = {}
+
         for tss in report_time_serie_act:
             where = report_time_serie_act[tss].where
-            outpoints = binding[where]
-            if where == "Catchments":
-                outpoints = decompress(outpoints)
+            if where in outpoints_by_where:
+                outpoints = outpoints_by_where[where]
             else:
-                coord = binding[where].split()  # could be gauges, sites, lakeSites etc.
-                if len(coord) % 2 == 0:
-                    outpoints = valuecell(self.var.MaskMap, coord, outpoints)
+                outpoints = binding[where]
+                if where == "Catchments":
+                    outpoints = decompress(outpoints)
                 else:
-                    try:
-                        outpoints = loadmap(where, pcr=True)
-                        outpoints = ifthen(outpoints != 0, outpoints)
-                        # this is necessary if netcdf maps are loaded !! otherwise strange dis.tss
-                    except Exception as e:
-                        msg = "Setting output points\n {}".format(str(e))
-                        raise LisfloodFileError(outpoints, msg)
+                    coord = binding[where].split()  # could be gauges, sites, lakeSites etc.
+                    if len(coord) % 2 == 0:
+                        outpoints = valuecell(self.var.MaskMap, coord, outpoints)
+                    else:
+                        try:
+                            outpoints = loadmap(where, pcr=True)
+                            outpoints = ifthen(outpoints != 0, outpoints)
+                            # this is necessary if netcdf maps are loaded !! otherwise strange dis.tss
+                        except Exception as e:
+                            msg = "Setting output points\n {}".format(str(e))
+                            raise LisfloodFileError(outpoints, msg)
+                outpoints_by_where[where] = outpoints
 
             if option['MonteCarlo']:
                 if os.path.exists(os.path.split(binding[tss])[0]):
@@ -716,7 +725,20 @@ def dynamic(self):
                 if how == 'total':
                     changed = compressArray(catchmenttotal(decompress(eval(what)) * self.var.PixelAreaPcr, self.var.Ldd) * self.var.InvUpArea)
                     what = 'changed'
-                self.var.Tss[tss].sample(decompress(eval(what)))
+                tss_obj = self.var.Tss[tss]
+                value_1d = eval(what)
+                # Fast path: when every gauge is a single cell, the inherited
+                # .sample()'s per-zone areaaverage equals the gauge-cell value,
+                # so read it straight from the 1D array and skip the
+                # decompress -> numpy2pcr full-map conversion (~50-136 ms/call on
+                # a global grid). Only for Scalar series (the areaaverage branch);
+                # a 1D numpy array implies Scalar/Directional here. Falls back to
+                # the exact PCRaster path otherwise (e.g. multi-cell gauge zones,
+                # non-spatial gauges, or non-ndarray values).
+                if getattr(tss_obj, '_compressed_sampling_ok', False) and isinstance(value_1d, np.ndarray):
+                    tss_obj.sample_compressed(value_1d)
+                else:
+                    tss_obj.sample(decompress(value_1d))
 
         # ************************************************************
         # ***** WRITING RESULTS: MAPS   ******************************