Implement OGGM's flowline model in Julia

[JordiBolibar]

A first step towards OGGM-ODINN compatibility will be to implement OGGM's flowline model in Julia.

A first version is more or less working, but it needs to be revised and optimized (i.e. enforce in-place operations). The next step would be to make the interface more versatile, in order to take gdirs and multiple mass balance models as an entity task.

[facusapienza21]

Although not efficient, the code seems to be implemented correctly. We can observe this from a comparison between the output of the Python code provided in the OGGM tutorial and the Julia code available in 1D_SIA.jl:

This was obtain with the following specifications for the ODE solver from DifferentialEquations.jl.

iceflow_sol = solve(iceflow_prob,
                            alg=Tsit5(),
                            adaptive=true,
                            reltol=1e-6, 
                            save_everystep=false, 
                            progress=true, 
                            progress_steps = 10)

The total running time is 183s. This is far from being efficient, since the Python code does the job in 5s. We also observe wiggles in the solution, which are caused by instabilities. The source of these is still unknown. I think this can be caused by problems in the boundary conditions. In order to exactly reproduce the Python solver we can use the Euler() solver and specify the timestep:

    iceflow_sol = solve(iceflow_prob,
                        alg=Euler(),
                        dt= 10 * sec_in_day,
                        reltol=1e-6, save_everystep=false, 
                        progress=true, 
                        progress_steps = 10)

and this is the solution we observed by t=52.5, just before the solution unstabilizes. I write a manual explicit Euler scheme (same as in Python) and the solver explodes at the same time, so I don't think this problem is coming from DifferentialEquations.jl.

It's still not clear what is the reason for this. It is hard to believe there is an obvious typo in the code since the first scenario shows similar solution (ignoring the wiggles), but there must be something there causing this unstabilities. My guess is that there is something with the boundary conditions that is not working, but I don't understand why is this is the case since the Python code is working.

I also tried implicit schemes, but this doesn't solve the issue.

[fmaussion]

I'm still unsure why we can't simply start with a "pure" translation, i.e. the same as the python code but written in Julia, and then see from here. @pat-schmitt also has an implicit implementation somewhere.

Having a "one to one" translation to the code would allow to compare performances of python, python+Numba (which is blasingly fast), Julia.

Even if we don't use the fancy julia solvers, if julia is really that fast, then converting the entire core of OGGM (dynamics + mass-balance models) would be a very achievable target (much better than trying to numba-ize the OGGM code which is impossible without huge refactoring)

[JordiBolibar]

Today I worked on a new raw version mimicking the exact type of explicit solver from the original Python script. The code runs really fast, but I'm still having some weird instabilities related to the diffusivity. I still haven't been able to isolate the issue. Maybe @facusapienza21 can take a look and go on from there?

I think the issue might come from some boundary conditions or some index problem (which are different in Python and Julia). Someone else taking a look at my code might be a good idea, since my brain is pretty fried this week.

[facusapienza21]

Ok, I finally found the problem with the code. You can see the details in #3 (@JordiBolibar , I made a PR since I don't have permission for this repo).

I also this opportunity to benchmark a little bit the different methods when running the same Explicit Euler scheme with same parameters:

1. Python code: 4.27s
2. Julia code with manual Euler method (pure translation from Python, as @fmaussion suggested): 0.99s
3. Julia code with Euler method using DifferentialEquations.jl: 2.06s

A pure translation from the Python to the Julia code shows a x4 factor of speed up in the code. I am pretty sure that this last option using DifferentialEquations.jl can be improved using the same tricks that @JordiBolibar used for the manual Euler method and also by using a more cleaver solver. It is also worth mentioning that one of the advantages of Julia here is that we can play more with the solvers to get better tradeoffs between accuracy and computational speed. By specifying the solver and the desired level of accuracy, we can rely more on the solver without using a fixed timestep.

[JordiBolibar]

Thanks Facu for having taken a look. I have integrated your PR, and I have some additional improvements to make it faster. The new benchmarking shows the following results:

• Python: 4.27s
• Julia: 0.39s

This is roughly a 11x improvement in speed. Not bad. I'm sure it can be further improved doing some more tricks, since I can still see around 700MB of memory being allocated.

As Facu mentioned, next step will be to implement this using DifferentialEquations.jl in an efficient manner.

[JordiBolibar]

After a short discussion in the Julia Discourse, I managed to implement some more tricks, further boosting the performance: dc69d61

New benchmark:

• Python: 4.27s
• Julia: 0.14s

We are now at a 30x speed improvement. Not bad, I think we can stop here for now 😄.

[fmaussion]

Just because this has now made it to twitter 😉 this is not the flowline model used by OGGM, this is a toy model from Karthaus and Oerlemans 1997.

I still expect Julia to be order of magnitudes faster nevertheless, also on the "real" model.