add option to use thomas algorithm within SemiImplicitModel

[pat-schmitt]

Here I add an option to the SemiImplicitModel to use a strapped down version of scipy.linalg.solve_banded, which replicates the thomas algorithm. This means we are not checking that the inputs are valid (no NaNs and finite) and we automatically select the correct fortran function for our input type. All this is normally happening inside of scipy.linalg.solve_banded at every call

I quickly tested the speed up with this:

from oggm import cfg, utils, workflow, DEFAULT_BASE_URL, tasks
from oggm.core.flowline import SemiImplicitModel
from functools import partial
import time
import numpy as np

cfg.initialize()
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='semi_implicit_speedup', reset=True)
rgi_ids = ['RGI60-11.00897']
gdirs = workflow.init_glacier_directories(rgi_ids,prepro_base_url=DEFAULT_BASE_URL,
                                          from_prepro_level=4,prepro_border=80)

gdir = gdirs[0]

SemiImplicitModelThomas = partial(SemiImplicitModel, use_thomas=True)

for nyears in [10, 100, 1000]:
    start_time = time.time()
    tasks.run_random_climate(gdir, y0=2000, seed=0, evolution_model=SemiImplicitModel,
                             nyears=nyears,
                             output_filesuffix='_original')
    time_orig = time.time() - start_time
    
    start_time = time.time()
    tasks.run_random_climate(gdir, y0=2000, seed=0, evolution_model=SemiImplicitModelThomas,
                             nyears=nyears,
                             output_filesuffix='_thomas')
    time_thomas = time.time() - start_time

    print(f"nyears: {nyears}")
    print(f"original: {time_orig:.2f} s")
    print(f"thomas: {time_thomas:.2f} s")
    print(f"thomas time is {100 / time_orig * time_thomas:.0f} % of original time")

    # small test everything works as expected, for production should add a test with fl_diagnositcs
    ds_orig = utils.compile_run_output(gdir, input_filesuffix='_original')
    ds_thomas = utils.compile_run_output(gdir, input_filesuffix='_thomas')
    for key in ds_orig:
        np.testing.assert_allclose(ds_orig[key], ds_thomas[key])

And found that the new implementation needs between 75% to 85% of the original time. Not much but for regional to global simulations this will add up.

For discussion if it is worth to implement properly (adding more tests that the model output is not effected).

• Tests added/passed
• Fully documented
• Entry in whats-new.rst

[fmaussion]

I think that it does not need to be an option unless for testing? For me a few GCM simulations for a few regions (say, alps and antarctica) to check that results are unchanged would be enough as a test.

If results and error rates are the same, then this should simply be the new default?