Good solver for nonlinear diffusion PDE with source?

[JordiBolibar]

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!

[ChrisRackauckas]

@JordiBolibar

• 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.

Do you have a good eigenvalue estimate? ROCK4 with a eigen est seems to perform much better on reverse passes than without an eigen estimate.

I can’t run things right now, but a few points to see:

1. Is computing the sparsity pattern possible here? Did you try the preconditioned GMRES from https://diffeq.sciml.ai/dev/tutorials/advanced_ode_example/ ? That’s pretty much always going to be best with scaling unless you get away with ROCK
2. Did you try quadrature adjoint with stiff solvers?

@JordiBolibar

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.

No it’s that this method doesn’t scale to large equations, and your equation with interpolating adjoint includes all of your parameters.