Climatemachine.jl: implement implicit solvers

[kmdeck]

The soil model (branch em/agl/kd/LandFramework) currently implements Richard's equation for water in soil. We have a test case here: https://github.com/CliMA/ClimateMachine.jl/blob/em/agl/kd/LandFramework/test/Land/Model/haverkamp_test.jl

Issues with the explicit solvers arise for certain parameter choices and hydraulics models.

[kmdeck]

@charleskawczynski @simonbyrne

[charleskawczynski]

Concretely, we tried a first swing at this on this branch (which built ontop of the existing DiffEq PR). The file we were running is test/Numerics/ODESolvers/heat_eq_diff_eq.jl. I see that we were trying Kvaerno3() and Rosenbrock23(). IIRC, these resulted in different errors, one was related to ldiv! not being defined, I think the other was related to ForwardDiff.

[charleskawczynski]

Per @thomasgibson: at the moment we construct our linear solvers here. And, per @ChrisRackauckas the DiffEq linsolve interface is here.

[ChrisRackauckas]

The DiffEq interface for linear solvers is https://docs.sciml.ai/latest/features/linear_nonlinear/#Linear-Solvers:-linsolve-Specification-1 . Basically, an init dispatch lets it establish any caches that it needs, and then it has a function that's essentially f(x,A,b,...) for Ax=b. From the chat I know about:

https://github.com/CliMA/ClimateMachine.jl/tree/master/src/Numerics/SystemSolvers

So I guess two things come to mind:

1. What is the right Jacobian type to use?
2. On that Jacobian type, which linear solver should we use?

[ChrisRackauckas]

Bringing @kanav99 in.

[kmdeck]

bringing in @gagnelandmanna who is working on the soil model as well.

[charleskawczynski]

So I guess two things come to mind:

1. What is the right Jacobian type to use?
2. On that Jacobian type, which linear solver should we use?

Can we use a Jacobian-free method to start with? The Jacobian is really problem-dependent and we'll have more flexibility when we go to the real soil model if we can use a Jacobian-free method--this heat equation is a toy model, but uses the same interface that the real soil model will use.

[ChrisRackauckas]

Yeah, we could do that. That would just be something like TRBDF2(linsolve=LinSolveGMRES()) if you just want to use IterativeSolvers.jl. Maybe that just works, or maybe we have to wrap your GMRES, but either way that wouldn't be difficult. Let's start with that direction then.