speed up of computing grids

[mbakker7]

I have been testing the speed of computing a grid with the current implementation (also in comments at the latest workshop where users reported that a string of wells is much slower than a string of line-sinks, which it shouldn't be). I compared two models to compute a grid of 100 by 100:

1. Single layer, semi-confined, steady, 100 wells on a row
2. Single layer, semi-confined, steady, 10 line-sinks along same row (one line-sink is computationally at least the same as the 8 wells used when doing numerical integration)
3. Single layer, semi-confined, steady, 100a line-sinks

Results:

wells, parallel=False: 5.44 s
wells, parallel=True: 3.07 s
10 line-sinks, parallel=False: 0.51 s
10 line-sinks, parallel=True: 2.47 s !!
100 line-sinks: parallel=False: 4.56 s
100 line-sinks, parallel=True: 3.32 s

Two things are surprising: 100 wells is 10 times slower than 10 line-sinks.
parallel=True doesn't work at all for 10 line-sinks, but does work a bit for 100 line-sinks. And 100 line-sinks ar about 10 times slower than 10 line-sinks, as you expect.

Then I did the same thing, but for a 10 layer model, again with semi-confined top. Results:

wells, parallel=False: 5.43 s
wells, parallel=True: 3.16 s
10 line-sinks, parallel=False: 0.52 s
10 line-sinks, parallel=True: 2.49 s !!
100 line-sinks: parallel=False: 4.54 s
100 line-sinks, parallel=True: 3.49 s

So this is almost the same. That is surprising, as we are computing things now for 10 layers rather than 1 layer. This may be, because all new lambda values are really small. So the potential goes to zero quickly. Next, I made the vertical conductivity really small, so that all new lambda values are substantial. New results:

wells, parallel=False: 5.26 s
wells, parallel=True: 3.2 s
10 line-sinks, parallel=False: 0.72 s
10 line-sinks, parallel=True: 2.6 s !!
100 line-sinks: parallel=False: 6.25 s
100 line-sinks, parallel=True: 3.79 s

So, a bit of an increase for the line-sinks, but none for the wells. Then again, the wells are not really set to zero beyond r > 8*lab now. But they probably should.

[dbrakenhoff]

I did some testing as well, and I found that for models with many wells (600 in my model), the current implementation of parallel does not speed things up at all. It is barely faster than the single core run and often slower (depending on the number of workers). My feeling is the numba implementation will bypass this issue, but I haven't tested it yet. But definitely something we should look into.

I did a profile run of the head computation for my model and I guess computing the potential in Python loops causes it to be slow? It takes 0.05 s for 600 potinf computations, then 0.093 s for summing the potinfs. The element.potential function takes 0.178s to perform the multiplications and the slicing/reshaping.

         21771 function calls (21768 primitive calls) in 0.404 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      600    0.178    0.000    0.327    0.001 element.py:100(potential)
      601    0.093    0.000    0.093    0.000 {method 'reduce' of 'numpy.ufunc' objects}
        1    0.073    0.073    0.397    0.397 model.py:212(potential)
      600    0.051    0.000    0.054    0.000 well.py:92(potinf)
      601    0.002    0.000    0.095    0.000 fromnumeric.py:69(_wrapreduction)
     1201    0.001    0.000    0.001    0.000 {built-in method numpy.zeros}
      601    0.001    0.000    0.096    0.000 fromnumeric.py:2338(sum)
    15000    0.001    0.000    0.001    0.000 {built-in method builtins.abs}
        2    0.001    0.000    0.006    0.003 {method '__exit__' of 'sqlite3.Connection' objects}
      601    0.001    0.000    0.001    0.000 {method 'reshape' of 'numpy.ndarray' objects}
      609    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}
      601    0.000    0.000    0.000    0.000 {method 'items' of 'dict' objects}
        1    0.000    0.000    0.000    0.000 {method 'execute' of 'sqlite3.Connection' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
      2/1    0.000    0.000    0.403    0.403 {built-in method builtins.exec}
      601    0.000    0.000    0.000    0.000 fromnumeric.py:2333(_sum_dispatcher)
        1    0.000    0.000    0.006    0.006 __init__.py:244(fun)
      2/1    0.000    0.000    0.397    0.397 model.py:345(head)
        1    0.000    0.000    0.000    0.000 aquifer.py:202(potential_to_head)

Another issue that popped up is that the number of workers that I choose for the parallel run isn't always optimal. Some users in the workshop got better results with parallel=4 for example. So it is worth trying a few values to see what works best. And yes for small workflows, parallel=True will lead to performance degradation because the overhead of starting N workers for <=N tasks is really not worth it. It might be worth checking the length of the task list vs the number of workers so that we don't spawn too many workers for the job.

Finally, the current parallel implementation is very very memory hungry, since i send the entire model along with each task (🙈). That can be greatly reduced by adding an initializer that creates a single global model variable for each worker that all tasks (e.g. head computations) can then use.