With @facusapienza we’re trying to find a good solver to solve a 2D nonlinear ice flow PDE with mass addition and subtraction (i.e. heat equation with source). We are encountering a big difference in solving it with and without source. Here are the main differences:

Without source:

- PDE: A
`Ralston()`

solver provided the best performance, with a benchmarked reference time (simulation of 5 years) of around 15s. - UDE: A
`ROCK4()`

solver provided the best performance. Since there is AD, the same solver as for the PDE. didn’t give the best results The reference time is around 30s.

With source:

- PDE: This one was hard to find, after trying over 20 solvers, I encountered the best results with a
`RDPK3Sp35()`

solver (5-stage, third order low-storage scheme with embedded error estimator, optimized for compressible fluid mechanics). The reference time is around 80s. - UDE:
**Here comes the problem**. So far I have not found any good solver for this case,`RDPK3Sp35()`

and similar solvers seem to be the ones which perform the best, but simulations are awfully slow, with reference times over 1h (totally not useable). The simulations seem correct, since I’m not seeing any`NaN`

values. I think the problem comes from the fact that sources change the surface gradient quite a lot, and since the diffusivity is a nonlinear function of the surface gradient (4th power), this changes drastically the timestep required for convergence.

Does anybody have any suggestions based on their experience on this? The benchmark we’re using can be found here. I don’t have any self-contained MWE, since what we’re trying to optimize is the full model itself for scientific simulations. @ChrisRackauckas might surely have some insights on this. Thanks in advance!