Free threading?

[tyralla]

We are currently creating a model that covers a huge river basin (Rhein) in a detailed spatial resolution (2.5 km). We want to calibrate it automatically "in one piece" based on the Multiscale Parameter Regionalization, which will hardly be possible without some kind of parallelisation. Now that HydPy runs on Python 3.13, the first Python version that provides a (still experimental) build that supports free threading, we asked ourselves if this is (already) the right way to introduce parallelism into HydPy.

(Parallelising HydPy-MPR based on multi-threading or multi-processing could eventually be simpler but would be a less general solution.)

For a start, I installed Python 3.13.1-t and tried to pip-install our requirements. Unfortunately, netCDF4 (version 4-1.7.2) does not seem to be ready:

      reading from setup.cfg...
          HDF5_DIR environment variable not set, checking some standard locations ..
...
        File "<string>", line 277, in <module>
        File "<string>", line 226, in _populate_hdf5_info
      ValueError: did not find HDF5 headers
      [end of output]

I could not find any information directly related to netCDF4. However, I found this issue related to H5Py, so there generally seems to be some movement. I also tried to install h5netcdf (an alternative to netCDF4 that we already considered when starting to add NetCDF features to HydPy), but this also did not work (version 1.4.1):

      running build_ext
      Loading library to get build settings and version: hdf5.dll
      error: Unable to load dependency HDF5, ensure HDF5 is installed properly
      on sys.platform='win32' with platform.machine()='AMD64'
      Library dirs checked: []
      error: Could not find module 'hdf5.dll' (or one of its dependencies). Try using the full path with constructor syntax.

So, at least out of the box, Python 3.13-t and NetCDF support do not currently go hand in hand.

In the next step, I removed netCDF4 from our requirements. The installation failed due to problems with building 'msgpack._cmsgpack' extension. Strangely, Visual Studio reports some severe errors (e.g. syntax errors):

      msgpack/_cmsgpack.c(2343): error C2146: Syntaxfehler: Fehlendes ")" vor Bezeichner "vc"
      msgpack/_cmsgpack.c(2343): error C2081: "__pyx_vectorcallfunc": Name in der formalen Parameterliste ist ung\x81ltig
      msgpack/_cmsgpack.c(2343): error C2061: Syntaxfehler: Bezeichner "vc"
      msgpack/_cmsgpack.c(2343): error C2059: Syntaxfehler: ";"
      msgpack/_cmsgpack.c(2343): error C2059: Syntaxfehler: ","
      msgpack/_cmsgpack.c(2343): error C2059: Syntaxfehler: ")"
      msgpack/_cmsgpack.c(6941): warning C4996: 'Py_OptimizeFlag': deprecated in 3.12
      msgpack/_cmsgpack.c(8188): warning C4244: "Funktion": Konvertierung von "Py_ssize_t" in "unsigned int", m”glicher Datenverlust
      msgpack/_cmsgpack.c(8343): warning C4244: "Funktion": Konvertierung von "long" in "char", m”glicher Datenverlust
      msgpack/_cmsgpack.c(8493): warning C4244: "Funktion": Konvertierung von "Py_ssize_t" in "unsigned int", m”glicher Datenverlust
      msgpack/_cmsgpack.c(9882): warning C4244: "Funktion": Konvertierung von "__int64" in "unsigned int", m”glicher Datenverlust
      msgpack/_cmsgpack.c(10124): warning C4244: "Funktion": Konvertierung von "__int64" in "unsigned int", m”glicher Datenverlust
      msgpack/_cmsgpack.c(10386): warning C4244: "Funktion": Konvertierung von "Py_ssize_t" in "unsigned int", m”glicher Datenverlust
      msgpack/_cmsgpack.c(21391): error C2146: Syntaxfehler: Fehlendes ")" vor Bezeichner "vc"
      msgpack/_cmsgpack.c(21391): error C2081: "__pyx_vectorcallfunc": Name in der formalen Parameterliste ist ung\x81ltig
      msgpack/_cmsgpack.c(21391): error C2061: Syntaxfehler: Bezeichner "vc"
      msgpack/_cmsgpack.c(21391): error C2059: Syntaxfehler: ";"
      msgpack/_cmsgpack.c(21391): error C2059: Syntaxfehler: ","
      msgpack/_cmsgpack.c(21391): error C2059: Syntaxfehler: ")"
      msgpack/_cmsgpack.c(21436): error C2146: Syntaxfehler: Fehlendes ")" vor Bezeichner "vc"
      msgpack/_cmsgpack.c(21436): error C2081: "__pyx_vectorcallfunc": Name in der formalen Parameterliste ist ung\x81ltig
      msgpack/_cmsgpack.c(21436): error C2061: Syntaxfehler: Bezeichner "vc"
      msgpack/_cmsgpack.c(21436): error C2059: Syntaxfehler: ";"
      msgpack/_cmsgpack.c(21436): error C2059: Syntaxfehler: ","
      msgpack/_cmsgpack.c(21436): error C2059: Syntaxfehler: ")"
      msgpack/_cmsgpack.c(22125): error C2065: "__pyx_vectorcallfunc": nichtdeklarierter Bezeichner
      msgpack/_cmsgpack.c(22125): error C2146: Syntaxfehler: Fehlendes ";" vor Bezeichner "vc"
      msgpack/_cmsgpack.c(22125): error C2065: "vc": nichtdeklarierter Bezeichner
      msgpack/_cmsgpack.c(22125): warning C4047: "=": Anzahl der Dereferenzierungen bei "int" und "vectorcallfunc" unterschiedlich
      msgpack/_cmsgpack.c(22126): error C2065: "vc": nichtdeklarierter Bezeichner
      msgpack/_cmsgpack.c(22128): warning C4013: "__Pyx_PyVectorcall_FastCallDict" undefiniert; Annahme: extern mit R\x81ckgabetyp int
      msgpack/_cmsgpack.c(22128): error C2065: "vc": nichtdeklarierter Bezeichner
      msgpack/_cmsgpack.c(22128): warning C4047: "return": Anzahl der Dereferenzierungen bei "PyObject *" und "int" unterschiedlich

I did not even know what msgpack is. According to johnnydep, it is a dependency of cachecontrol, which is a dependency of lastversion. However, lastversion is "only" a helper tool for making the Installer for Windows, so it is irrelevant for a possible parallelisation of a separate version of HydPy. If I also remove latestversion from requirements.txt, pip install finishes successfully.

Here is a tracking issue that documents the free threading support of some critical Python site packages. Numpy should be ready since version 2.1 (so, no compatibility with Arcpy environments, but this is also not overly important). For matplotlib and pandas, there are also proper releases on PyPI. However, for Cython and Scipy, there are not. According to this issue, one would have to deal with Cython's master branch. Maybe it is similar to Scipy; here is the preliminary documentation on this topic.

Scipy should not concern us too much now because we do not use its functionalities during simulation runs (only, for example, when preparing a simulation run). Cython is definitely important, so I started a build with Cython 3.0.11*. python prepare_build.py finished successfully. python -m build results in similar msgpack-like errors:

hydpy\cythons\autogen\annutils.c(2765): error C2146: Syntaxfehler: Fehlendes ")" vor Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(2765): error C2081: "__pyx_vectorcallfunc": Name in der formalen Parameterliste ist ungültig
hydpy\cythons\autogen\annutils.c(2765): error C2061: Syntaxfehler: Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(2765): error C2059: Syntaxfehler: ";"
hydpy\cythons\autogen\annutils.c(2765): error C2059: Syntaxfehler: ","
hydpy\cythons\autogen\annutils.c(2765): error C2059: Syntaxfehler: ")"
hydpy\cythons\autogen\annutils.c(7612): warning C4996: 'Py_OptimizeFlag': deprecated in 3.12
hydpy\cythons\autogen\annutils.c(12342): warning C4996: 'Py_OptimizeFlag': deprecated in 3.12
hydpy\cythons\autogen\annutils.c(29791): error C2146: Syntaxfehler: Fehlendes ")" vor Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(29791): error C2081: "__pyx_vectorcallfunc": Name in der formalen Parameterliste ist ungültig
hydpy\cythons\autogen\annutils.c(29791): error C2061: Syntaxfehler: Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(29791): error C2059: Syntaxfehler: ";"
hydpy\cythons\autogen\annutils.c(29791): error C2059: Syntaxfehler: ","
hydpy\cythons\autogen\annutils.c(29791): error C2059: Syntaxfehler: ")"
hydpy\cythons\autogen\annutils.c(29836): error C2146: Syntaxfehler: Fehlendes ")" vor Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(29836): error C2081: "__pyx_vectorcallfunc": Name in der formalen Parameterliste ist ungültig
hydpy\cythons\autogen\annutils.c(29836): error C2061: Syntaxfehler: Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(29836): error C2059: Syntaxfehler: ";"
hydpy\cythons\autogen\annutils.c(29836): error C2059: Syntaxfehler: ","
hydpy\cythons\autogen\annutils.c(29836): error C2059: Syntaxfehler: ")"
hydpy\cythons\autogen\annutils.c(30525): error C2065: "__pyx_vectorcallfunc": nichtdeklarierter Bezeichner
hydpy\cythons\autogen\annutils.c(30525): error C2146: Syntaxfehler: Fehlendes ";" vor Bezeichner "vc"
hydpy\cythons\autogen\annutils.c(30525): error C2065: "vc": nichtdeklarierter Bezeichner
hydpy\cythons\autogen\annutils.c(30525): warning C4047: "=": Anzahl der Dereferenzierungen bei "int" und "vectorcallfunc" unterschiedlich
hydpy\cythons\autogen\annutils.c(30526): error C2065: "vc": nichtdeklarierter Bezeichner
hydpy\cythons\autogen\annutils.c(30528): warning C4013: "__Pyx_PyVectorcall_FastCallDict" undefiniert; Annahme: extern mit Rückgabetyp int
hydpy\cythons\autogen\annutils.c(30528): error C2065: "vc": nichtdeklarierter Bezeichner
hydpy\cythons\autogen\annutils.c(30528): warning C4047: "return": Anzahl der Dereferenzierungen bei "PyObject *" und "int" unterschiedlich

So, the next step necessary seems to install the master branch Cython version, which, based on this instruction, is easy. Running python prepare_build.py with the current version of Cython 3.1.0a1 succeeds but warns as follows:

<frozen importlib._bootstrap>:488: RuntimeWarning: The global interpreter lock (GIL) has been enabled to load module 'Cython.Utils', which has not declared that it can run safely without the GIL. To override this behaviour and keep the GIL disabled (at your own risk), run with PYTHON_GIL=0 or -Xgil=0.

python -m build still fails. According to this comment, it is because we have to build numpy with the previously built cython version. Numpy's documentation informs us that this requires more effort.

I stop here. It appears that trying to get everything running and then maintaining everything in a fast-changing environment is too much work and we should wait at least for an official release of Cython 3.1 and a fitting release of Numpy. However, implementing HydPy-internal parallelisation differently now, given that free threading seems reachable soon, also does not seem right. So, maybe we should postpone it and "just" try to parallelise HydPy-MPR for now.

@BGWKlein @GernotBelger

[tyralla]

I also quickly checked for the possibility of working with multiple interpreters to achieve parallelization. This feature is expected to land in Python 3.14 fully, but this site package should help to at least test it now. However, numpy's support is still incomplete.

If I remember correctly, we do not really use Numpy's functionalities during simulation runs. During simulations, we are only accessing and setting the values of numpy arrays through the typed memoryviews interface (if this is the right way to say it). So, numpy limitations must not necessarily mean a problem for us. However, it again seems a little bit early to follow this route...

[tyralla]

Oh, I forgot to consider one long-planned opportunity. We always took great care to add the nogil option to many (or even all?) simulation-related methods of our Cython extension classes. This was always a means of not using Python features with much overhead by accident. Also, we thought about parallelising parts of our models via Cython's prange (which did not seem very efficient in some preliminary experiments, if I remember correctly). Packing multiple model instances into one big prange loop was also discussed but appeared very complicated. However, maybe combining the threading feature of "normal" Python builds with actually releasing the GIL within the main calls to our Cython extension classes helps? (I don't know why I never thought about it...)

Here is a very simple preliminary experiment based on our HydPy-H-Lahn example project (with a simulation period of 10 years).

This is how the core of our central simulate method looks:

            with cm:
                for idx in printtools.progressbar(range(idx_start, idx_end)):
                    for func in methodorder:
                        func(idx)

I changed it to:

            from concurrent.futures import as_completed, ThreadPoolExecutor
            with cm, ThreadPoolExecutor(max_workers=4) as executor:
                for idx in printtools.progressbar(range(idx_start, idx_end)):
                    futures = [executor.submit(func, idx) for func in methodorder]
                    for future in as_completed(futures):
                        pass

I set max_workers to four because the Lahn catchment consists of four big subbasins. I have little experience with ThreadPoolExecutor, so there might be much to optimise.

Additionally, I put the body of the simulate method of all Cython extension classes (that represent hydrological models) into a with nogil block. There are some other methods which play some role, but this should cover at least 90% (maybe even 99%) of the total simulation time.

On the first try, I ran the example as is. Introducing threading increased the simulation time from 1 to 2 seconds.

In a second attempt, I set the value of the parameter sclass (number of snow classes) to 1,000. Then, threading decreases the simulation time from 5 to 4 seconds.

With 5,000 snow classes, Threading decreases simulation time from 30 to 12 seconds.

It should be mentioned that the simulation results are wrong because I did not care for the correct execution order in this preliminary example.

So, somehow, HydPy already has a "free threading" functionality thanks to using Cython's nogil option, which could be activated without too much effort. However, it is still quite unclear what performance increases we can expect in realistic applications (or if we can expect them at all).

[tyralla]

Gernot advised to repeat the last experiment based on Queue instead of ThreadPoolExecutor. The code is only a little more complicated and likely easier extendible for our purposes:

            import threading
            import queue

            class Worker(threading.Thread):
                self.queue_ = queue.Queue

                def __init__(self, queue_: queue.Queue, *args, **kwargs) -> None:
                    super().__init__(*args, **kwargs)
                    self.queue_ = queue_

                def run(self) -> None:
                    while True:
                        try:
                            method, idx = self.queue_.get()
                        except queue.ShutDown:
                            return
                        method(idx)
                        self.queue_.task_done()

            queue_ = queue.Queue()
            for _ in range(4):
                Worker(queue_=queue_).start()

            for idx in printtools.progressbar(range(idx_start, idx_end)):
                for method_or_methods in methodorder:
                    if isinstance(method_or_methods, list):
                        for method in method_or_methods:
                            queue_.put((method, idx))
                        queue_.join()
                    else:
                        method_or_methods(idx)

            queue_.shutdown()

The computation time now decreases (compared to the non-threading approach) in all cases (0.7 --> 0.5, 4.8 --> 2.6, 30.0 --> 13.5). Compared the the ThreadPoolExecutor example, the Queue approach is a little bit slower for the "biggest" model. However, we should not take the last comparison too seriously. The Queue approach already includes a small improvement as it only tries to parallelise the simulate methods (hence the preliminary isinstance check).

[tyralla]

I noted that the above threading did not access as much CPU power as the number of workers (4) suggests. I think this is due to the different computation times of the four parallelised model instances. So here are the results for setting RecStep to 1e4, which affects all model instances in the same manner:

• no threading: 3.4 s (100 %)
• Pool, 4 workers: 2.4 s (71 %)
• Queue, 4 workers: 1.3 s (39 %)

And here for 1e5:

• no threading: 25.9 s (100 %)
• Pool, 4 workers: 9.5 s (37 %)
• Queue, 4 workers: 8.3 s (32 %)

Quite close to the optimum of 25 %. But we must find ways to group our method calls in huge, parallelisable chunks.

[tyralla]

All of the above threading examples lead to wrong results. Here is the first approach which ensures correct results (at least I hope so). Therefore, I derived a Queue subclass that is partly tailor-suited to the Lahn example project:

       import threading
            import queue

            class Queue(queue.Queue):

                def __init__(self) -> None:
                    super().__init__()
                    self.upstream2downstream = {
                        "land_dill_assl": "stream_dill_assl_lahn_leun",
                        "land_lahn_marb": "stream_lahn_marb_lahn_leun",
                        "land_lahn_leun": "stream_lahn_leun_lahn_kalk",
                        "stream_dill_assl_lahn_leun": "stream_lahn_leun_lahn_kalk",
                        "stream_lahn_marb_lahn_leun": "stream_lahn_leun_lahn_kalk",
                    }
                    self.starters = {
                        "land_dill_assl",
                        "land_lahn_marb",
                        "land_lahn_leun",
                        "land_lahn_kalk",
                    }
                    self.dependencies = {
                        "stream_dill_assl_lahn_leun": 1,
                        "stream_lahn_marb_lahn_leun": 1,
                        "stream_lahn_leun_lahn_kalk": 3,
                    }

                def register(self, method_name: tuple[Callable, str], idx: int) -> None:
                    self.idx = idx
                    self.waiting = {}
                    for method, name in method_name:
                        if name in self.starters:
                            self.put((method, name, idx))
                        else:
                            self.waiting[name] = [method, self.dependencies[name]]

                def task_done(self, name: str) -> None:
                    name = self.upstream2downstream.get(name)
                    if name is not None:
                        method, nmb = self.waiting[name]
                        if nmb == 1:
                            self.put((method, name, self.idx))
                        else:
                            self.waiting[name][1] -= 1
                    super().task_done()

            class Worker(threading.Thread):
                self.queue_ = Queue

                def __init__(self, queue_: Queue, *args, **kwargs) -> None:
                    super().__init__(*args, **kwargs)
                    self.queue_ = queue_

                def run(self) -> None:
                    while True:
                        try:
                            method, name, idx = self.queue_.get()
                        except queue.ShutDown:
                            return
                        method(idx)
                        self.queue_.task_done(name)

            queue_ = Queue()
            for _ in range(4):
                Worker(queue_=queue_).start()

            for idx in printtools.progressbar(range(idx_start, idx_end)):
                for method_or_methods in methodorder:
                    if isinstance(method_or_methods, list):
                        queue_.register(method_or_methods, idx)
                        queue_.join()
                    else:
                        method_or_methods(idx)

            queue_.shutdown()

There might be scope for further optimisation. Still, the additional computation times compared to the non-tread safe example above are not too big:

RecStep unmodified:

• no threading: 0.74 s (100 %)
• Queue, 4 workers: 0.69 s (93 %)

RecStep = 1e4:

• no threading: 3.4 s (100 %)
• Queue, 4 workers: 1.5 s (44 %)

RecStep = 1e5:

• no threading: 25.9 s (100 %)
• Queue, 4 workers: 8.5 s (33 %)

[tyralla]

Unfortunately, the above Queue approach is not absolutely thread-safe because multiple model instances (handled by different elements) can modify the discharge value handled by a single node at the same time. For HydPy-H-HBV96:

class Pass_Q_V1(modeltools.Method):
    r"""Update the outlet link sequence."""

    REQUIREDSEQUENCES = (hland_fluxes.QT,)
    RESULTSEQUENCES = (hland_outlets.Q,)

    @staticmethod
    def __call__(model: modeltools.Model) -> None:
        flu = model.sequences.fluxes.fastaccess
        out = model.sequences.outlets.fastaccess
        out.q[0] += flu.qt

Here, out.q[0] += flu.q means that the node's original double is queried and modified via a C pointer. Compared to the other calculations, this takes little time. Hence, the corresponding errors should occur rarely and did not affect the above simulation results. However, they can happen, so we need to deal with them. I see three solutions:

1. Disallow multithreading if multiple models modify the same node's value. (Simple but inconvenient for users.)
2. Ensure that models that modify the same node's value are run sequentially. (It could be complicated to decide automatically which is the best running order regarding computational efficiency. In worst cases, this would disallow any parallelization.)
3. Reacquire the GIL when calling certain methods like update_outlets, update_senders, update_outputs, and possibly others. (Reacquiring the GIL and re-releasing takes additional time. Also, we would have to double-check that no method of any model that is not "officially" allowed to do so modifies the values of nodes.)
4. Add an individual lock to all model instances that modify the same node's value. (This is similar to the third suggestion. It could be more efficient, but it is definitely more complicated. First, we would have to implement it based on C. Second, it would only work if all used models were cythonized. Otherwise, one would have to fall back to solution 3 anyhow.)

At first glance, solution 3 seems like the first way to try. Here are the results:

RecStep unmodified:

• no threading: 0.74 s (100 %)
• Queue, 4 workers:
  • reacquire GIL no: 0.69 s (93 %)
  • reacquire GIL yes: 0.71 s (96 %)

RecStep = 1e4:

• no threading: 3.4 s (100 %)
• Queue, 4 workers:
  • reacquire GIL no: 1.5 s (44 %)
  • reacquire GIL yes: 1.5 s (44 %)

RecStep = 1e5:

• no threading: 25.9 s (100 %)
• Queue, 4 workers:
  • reacquire GIL no: 8.5 s (33 %)
  • reacquire GIL yes: 8.5 s (33 %)

I double-checked the slight increase in computation time in the unmodified example (RecStep not increased), and it seems real. However, it is so tiny that I think it is the way to go for now.

[tyralla]

As a first step to tidy up a little, I added the option threads, explained by "The number of additional threads to be opened during a simulation run.". Zero means "no multi-threading" (everything works as before). I also thought about nthreads or (nmbthreads or numberthreads ). When finally deciding on a name, we should take into account that we might also add another option later that allows configuring further details of how the threading should work (threadingmode or so).

[tyralla]

I did some refactoring and applied the preliminary state of our efforts to a finer model of the Lahn basin (1008 raster cells) with 12 additional threads. At the very first glance, the results seem correct, but the simulation time increases from around 8.5 minutes to 11.5 minutes. This was to be expected because each raster cell consists only of one hydrological response unit, so each individual job for a worker is shorter than in the previously discussed model, where each of the four subbasins has at least ten hydrological response units. (Note that the raster-cell model uses HydPy-L-DD instead of HydPy-H-HBV96, which might make a difference, too.)

[tyralla]

As feared, it really seems we need some kind of chunking to achieve relevant speed-ups in real-world applications. I will use the terms "spatial chunking" and "temporal chunking" in the following. I start with spatial chunking because this idea allows us to build a prototype in pure Python that uses the cythonized models without further modifications (no need to access the time series handled by independent node instances).

My first idea is to build a layer of spatial chunks. Assume we want to work with, at most, ten chunks at a time. So, the first layer consists of all headwater catchment models (with no dependencies), which are grouped into ten evenly sized chunks. The second layer consists of the uppermost river segments (that do have one headwater as a dependency). One can first process all models of the second layer, afterwards the models of the second layer, and so on. (Note that uppermost river segments with two headwater catchments would be members of the third layer but could be members of the second one - that is a detail for later.)

The simple core of this idea is that we need a way to group methods. Therefore, I implemented the class Chunk (quick and dirty):

class Chunk:

    methods: list[BoundMethod]
    upstreams: list[str]

    def __init__(self, methods: list[BoundMethod], upstreams: list[str]) -> None:
        self.methods = methods
        self.upstreams = upstreams

    def __call__(self, idx: int) -> None:
        for method in self.methods:
            method(idx)

I also introduced BoundMethod as a type alias for Callable[[int], None], of which Chunk is also a structural subtype. This means we can now pass Chunk instances to our workers that process multiple methods instead of passing individual methods to the workers. By allowing both ways simultaneously, we must still preserve the old "upstream to downstream" mechanism so that each first-layer chunk counts as a dependency for each second-layer chunk.

Running the current implementation (under the optimistic assumption that everything is working correctly) should give a first impression of whether pure-Python chunking (with the GIL very frequently enabled and disabled) brings any benefits for huge models. I checked it for our four-subbasin model and our 1008-grid cell model (the latter with a simulation period of only one year to save time):

• 4 subbasins
  • original RecStep
    • no threading: 0.74 s
    • threading (4 workers)
      • no chunking: 0.76 s
      • chunking: 0.79 s
  • RecStep 1e4
  • no threading: 3.5 s
  • threading (4 workers)
    • no chunking: 1.5 s
    • chunking: 1.6 s
  • RecStep 1e5
    • no threading: 26.1 s
    • threading (4 workers)
      • no chunking: 8.5 s
      • chunking: 8.6 s
• 1008 grid cells
  • no threading: 4.5 s
  • threading (12 workers)
    • no chunking: 10.5 s
    • chunking: 9.0 s

The first example (the four-subbasin model) shows that the parallelization works in principle. The second example shows a slight improvement over the not chunked threading approach, suggesting that we can, in fact, skip a little overhead by building chunks. However, we are far from any improvements compared to a non-threaded simulation run in the latter example.

[tyralla]

I added a small Cython prototype to our Python prototype to see if this brings relevant speed-ups. The general idea is that our Python-based Chunk class now triggers a Cython-based Chunk class:

# python
class Chunk:
    
    _chunk: threadingutils.Chunk
    upstreams: list[str]
    

    def __init__(self, methods: list[BoundMethod], upstreams: list[str]) -> None:
        self._chunk = threadingutils.Chunk(methods)
        self.upstreams = upstreams

    def __call__(self, idx: int) -> None:
        self._chunk.simulate(idx)

# cython
cdef class Chunk:

    def __init__(self, models: Sequence[Model]) -> None:
        free(self.models)
        self.number = len(models)
        size = self.number * cython.sizeof(cython.pointer(PyObject))
        self.models = <PyObject**>malloc(size)
        for i, model in enumerate(models):
            self.models[i] = <PyObject*>model.cymodel

    cpdef void simulate(self, int idx):
        cdef int j
        with nogil:
            for j in range(self.number):
                (<BaseInterface>self.models[j]).simulate(idx)

    def __dealloc__(self) -> None:
        free(self.models)

Therefore, I added the simulate method to our BaseInterface (which is a little dirty because BaseInterface deals with controlling submodels, but this is just for testing for now).

This approach should reduce threading overhead because, potentially, the GIL can now be permanently released during the processing of all methods grouped in one chunk. I even enforced this by not reacquiring the GIL when passing data to node instances. Here are the results:

• 1008 grid cells
  • no threading: 4.5 s
  • threading (12 workers)
    • no chunking: 10.5 s
    • chunking
      • Python-based: 9.0 s
      • Cython-based
        • reacquire GIL when necessary: 7.0 s
        • disable GIL permanently (risky): 5.7 s

So, the decreases in computation time have the expected tendencies but are much smaller than hoped for.

[tyralla]

I performed a very quick-and-dirty experiment to get a first impression of what is achievable with temporal chunking. Therefore, I added a simulate_period method to HydPy and our pure-Python and cythonized models:

# HydPy
    @printtools.print_progress
    def simulate_period(self) -> None:
        idx_start, idx_end = hydpy.pub.timegrids.simindices
        if hydpy.pub.options.threads == 0:
            for method in self._determine_methodorder_part2(multithreading=False):
                method(idx_start, idx_end)
        else:
            methods_names2 = self._determine_methodorder_part2(multithreading=True)
            queue_ = Queue.from_methods(
                elements=self.elements,
                ordered_names=tuple(name for _, name, _ in methods_names2),
            )
            for _ in range(hydpy.pub.options.threads):
                Worker(queue_=queue_).start()
            queue_.register(methods_names2, idx_start, idx_end)
            queue_.join()
            queue_.shutdown()

# cythonized model
    cpdef inline void simulate_period(self, numpy.int64_t idx0, numpy.int64_t idx1)  noexcept nogil:
        cdef numpy.int64_t idx
        with nogil:
            for idx in range(idx0, idx1):
                self.simulate(idx)

I excluded all other model methods from the loop in HydPy.simulate in the treading mode for simplicity and the non-threading mode for comparability. Hence, the simulated output is currently broken, and the following performance gains due to threading might overestimate the actually achievable gains a little:

• 4 subbasins (10 years)
  • RecStep unmodified
    • no threading: 0.71 s (100 %)
    • Queue, temporal chunking, 4 workers: 0.20 s (28 %)
  • RecStep = 1e4
    • no threading: 3.4 s (100 %)
    • Queue, 4 workers: 0.9 s (26 %)
  • RecStep = 1e5
    • no threading: 26.0 s (100 %)
    • Queue, 4 workers: 7.0 s (27 %)
• 1008 raster cells
  • 1 year
    • no threading
      • reacquire GIL yes: 2.2 s (100 %)
      • reacquire GIL no
        • release GIL twice: 2.1 s (100 %)
        • release GIL once: 2.1 s (100 %)
    • Queue, temporal chunking, 12 workers
      • reacquire GIL yes: 2.2 s (100 %)
      • reacquire GIL no
        • release GIL twice: 0.8 s (38 %)
        • release GIL once: 0.9 s (42 %)
  • 70 years
    - no threading
    • reacquire GIL yes: 192 s (100 %)
    • reacquire GIL no
      • release GIL twice: 185 s (100 %)
      • release GIL once: 187 s (100 %)
        - Queue, temporal chunking, 12 workers
    • reacquire GIL yes: 153 s (80 %)
    • reacquire GIL no
      • release GIL twice: 60 s (32 %)
      • release GIL once: 60 s (32 %)

Note that the simulation times of the 1008 raster cell model are increased because the routing models receive nan values as input. I think this is due to the iterative search for the reference water level of HydPy-Musk-MCT.

"release GIL twice" means with nogil is used in Model.simulate_period and Model.simulate, "release GIL once" means it is used only in Model.simulate_period.

"reacquire GIL" is again related to the tread-safety of passing results to node instances.

So, for long simulation periods and models where downstream locations do not affect upstream locations, temporal chunking seems very promising. However, it comes with a few more complications and possible limitations compared to spatial chunking, so we should keep on investigating both.

[tyralla]

All previous threading examples rely on Python features. Another approach is to use Cython's prange function, which uses OpenMP for executing the "iterations" of a loop in parallel threads.

Again, The following experiment is very quick and dirty and only supposed to give a first impression of what could be eventually achievable. Like before, only the models' simulate / simulate_period methods are called. All other methods are simply ignored. Additionally, I did not read any meteorological data, so it's all just a big processing of nan values.

I did not implement any chunking as prange seems to do it on its own in a way that is configurable via schedule and some related parameters. At first, I stuck to the default, which I think is "dynamic" (no chunking). In a second attempt, I passed the argument "static" (one chunk per thread).

I still use the Cython extension class Chunk, even if it does no chunking by itself. Its modified simulate_period works in principle as discussed above and includes temporal chunking. Its new simulate_period_stepwise method refrains from doing temporal chunking for comparison:

    cpdef void simulate_period(self, int idx_start, int idx_end):
        cdef int j
        for j in prange(self.number, nogil=True, schedule="static"):
            (<BaseInterface>self.models[j]).simulate_period(idx_start, idx_end)

    cpdef void simulate_period_stepwise(self, int idx_start, int idx_end):
        cdef int i, j
        for i in range(idx_start, idx_end):
            for j in prange(self.number, nogil=True, schedule="static"):
                (<BaseInterface>self.models[j]).simulate(i)

The results:

• 4 subbasins (10 years)
  • RecStep unmodified
    • no threading: 0.70 s (100 %)
    • OpenMP, 4 workers
      • simulate_period: 0.20 s (29 %)
      • simulate_period_stepwise: 0.26 s (37 %)
  • RecStep = 1e4
    • no threading: 3.5 s (100 %)
    • OpenMP, 4 workers
      • simulate_period: 0.9 s (26 %)
      • simulate_period_stepwise: 1.2 s (34 %)
  • RecStep = 1e5
    • no threading: 25.9 s (100 %)
    • OpenMP, 4 workers
      • simulate_period: 7,0 s (27 %)
      • simulate_period_stepwise: 9.5 s (37 %)
• 1008 raster cells
  • 1 year
    • no threading: 5.1 s
    • OpenMP, number of workers automatically selected
      • simulate_period
        • dynamic: 1.0 s (20 %)
        • static: 1.1 s (22 %)
      • simulate_period_stepwise
        • dynamic: 1.4 s (27 %)
        • static 1.4 s (27 %)
  • 70 years
    • no threading: 358.5 s
    • OpenMP, number of workers automatically selected
      • simulate_period
        • dynamic: 70.7 s (20 %)
        • static: 70.2 s (20 %)
      • simulate_period_stepwise
        • dynamic: 96.7 s (27 %)
        • static: 96.7 (27 %)

A little better than before, but still no final breakthrough...

[tyralla]

I successively increased the number of hydrological response units (HRUs) to get an impression of how the size of individual model instances affects the threading efficiency. All results are for the 1008 raster cell model using Cythons prange (OpenMP). Water is still not routed downstream and meteorological time series are not read.

Between 1 and 100 HRUs, threading efficiency increases (presumably) due to a higher ratio of model calculations and threading overhead. For more than 100 HRUs, threading efficiency sinks again (presumably) due to the increased model sizes. This is different from the 4 subbasin model, where we increased the model calculation workload by increasing the number of iterations only, not by also increasing the size of the related arrays.

This result holds for simulate_period and simulate_period_stepwise. Hence, it seems that temporal chunking cannot improve the situation.

If it is all about limited cache sizes, this experiment might reveal the highest efficiency we can hope for (at least for my computer). I am relatively sure we cannot decrease the size of our models without losing flexibility (e.g., by making parameters scalar) or transparency (e.g., by not storing internal fluxes or states in arrays).

[tyralla]

I repeated some of the last experiments with the same general setting but with meteorological input data available in RAM. The up-and-down pattern in threading efficiency with an increasing number of hydrological response units is very similar: