O2.3.6 Implicit solver for full soil model

[kmdeck]

Purpose

We need an implicit solver for the soil model (both for Richards equation alone, and the full soil system), and possibly in the future for other models as well. Currently, the implicit solver for Richards equation is implemented and in use. The full soil implicit solver is still in progress.

Components

• A way to specify separate explicit and implicit tendencies
  • Richards equation: https://github.com/CliMA/ClimaLand.jl/pull/191
• A way to set tolerances and other convergence requirements (max iters, etc)
  • Included by using ClimaTimeSteppers: https://github.com/CliMA/ClimaLand.jl/pull/218/files
• A way to specify Jacobians, and how often to update them; Function to compute Richards model Jacobian (derivation here)
  • Richards equation: https://github.com/CliMA/ClimaLand.jl/pull/190

Inputs

• For Richards equation: https://www.overleaf.com/9252289794ykktrwjxqvcj. We will be implementing modified Picard iterations.

Results and deliverables/Validation

• [x] For Richards equation: plots of mass conservation and solution error vs timestep for modified Picard iterations, for constant flux and state BC test runs. Comparison to an explicit solver of the same order - obtain max dt the explicit solver runs at. - in our buildkite
• [x] Richards equation: unit test of Jacobian computed by the code with the analytic expression
• [x] Make sure the above validation results are preserved with new CTS interface
• [x] For full soil model : plots of mass conservation and solution error vs timestep for modified Picard iterations, for constant flux, moisture state BC, and for AtmosDrivenBC. Comparison to an explicit solver of the same order - obtain max dt the explicit solver runs at.
• [x] Unit tests of jacobian
• [x] flexible implementation that allows for LSM models - with many component models - to also have Jacobians specified. Start with soil/canopy model

Initial task breakdown

A preliminary list of PRs and a preliminary timeline of PRs, milestones, and key results.

• [X] Have reviewers review equations (Tapio)
• [X] Add method computing Jacobian for approximate Newton (modified Picard). Make sure A and b update every X timesteps, where X is chosen at the top-most level.
• [x] Simulations with the implicit solver using ClimaTimesteppers. comparison of approx Newton to explicit runs
• [X] State varying boundary conditions treated implicitly
• [x] Implement in general fashion #190 #191 #218
• [x] update to new ClimaTimesteppers interface #528
• [x] Add Jacobian for full soil model; repeat tests #538
• [x] Add Jacobian for an LSM model; repeat tests #538
• [ ] Time dependent BC TBD

Producers

@juliasloan25 (design and implementation) @kmdeck @dennisYatunin (design, implementation support)

Reviewers

Katherine, Dennis - PRs Tapio - SDI review (@tapios)

SDI Changelog

• 4 March 2024: SDI updated for clarity and to more accurately reflect current status (by @juliasloan25)
• 2 Oct 2024: issues related to implicit timestepping but not required for this milestone have been moved to #809

### step RichardsModel implicitly
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/224
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/268
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/528
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/541
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/647

### step EnergyHydrology implicitly
- [ ] #538
- [ ] #660
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/663

### QA
- [x] write tutorial demonstrating how to use IMEX model (done - EnergyHydrologyModel example)

### canopy implicit timestepping
- [ ] https://github.com/CliMA/ClimaLand.jl/issues/696

[tapios]

LGTM.

We should keep the tridiagonal solver in one central place that we do not have to duplicate it (so that we update it in one place, e.g., for GPU performance, this percolates through).

[kmdeck]

@dennisYatunin could we move the thomas algorithm function to ClimaTimesteppers?

also, if the matrix passed is the identity, will the algorithm spend time inverting it, or will it just return b?

[dennisYatunin]

Actually, ClimaCore might be a slightly better location than ClimaTimeSteppers. We're currently computing Jacobian blocks as ClimaCore Fields, then copying them into arrays to run the Thomas Algorithm. It would probably be better to have a function in ClimaCore that can directly invert a Jacobian block Field (I left a note in ClimaCore last year to add a function like this).

If the Jacobian matrix is I (or some constant times I, which is encoded by the same struct), the matrix inversion will be ldiv!(x, I, b), which is defined in LinearAlgebra to just set x .= b (scaled by a constant, if applicable).

[kmdeck]

If the Jacobian matrix is I (or some constant times I, which is encoded by the same struct), the matrix inversion will be ldiv!(x, I, b), which is defined in LinearAlgebra to just set x .= b (scaled by a constant, if applicable).

excellent!

[juliasloan25]

completed - see #809 for followup