update implicit solver interface

[juliasloan25]

Purpose

Update implicit solver interface for RichardsModel to use ClimaCore.MatrixFields.

Comparing the buildkite timing to the most recent main build, the RichardsModels runs see a ~25% speedup from this PR.

closes #528
#602 has to be merged first for update to ClimaCore v0.14

To Do

• GPU compatibility - requires fixes for land implicit solver ClimaCore.jl#1684 and ClimaCore patch release
• compare allocations/performance to previous implementation

Content

• delete thomas_algorithm
  handling multiple vars
• the full soil case will contain a variable that isn't stepped implicitly (theta_i), so we can't just loop over everything in the state
  • when extending, will add if statement, so that we use tridiagonal for theta_l and heat, and identity for theta_i
• boundary conditions
• allocate space in cache for dtopBCdY (extend boundary_vars, boundary_var_domain_names, boundary_var_types, etc for MoistureStateBC)
• update cache var with expression currently in dtendencyBCdY (except multiply by dz factor)
• update dtendencyBCdY: update in place, rename to set_dfluxBCdY!, multiply by dz factor
• update dtheta_resdtheta setting: add if/else branch where if default use expression currently in code; if p.soil.dtopBCdtheta exists, add dtopBCdtheta term before div applied (use Operators.SetBoundaryOperators, then construct Bidiagonal(?)Matrix of zeros with this BC)

Status 4/17/24

The implicit solver has been updated and works on CPU. It doesn't work on GPU because of a couple of problems. First, one of the operators produces a non-isbits object. This is because the space is not being stripped out of the input object, which should happen via a call to strip_space. We need to extend adapt_structure for the SetBoundaryOperator type (and perhaps other operator types too), but even with that the solver fails on GPU.

It seems like the second bug may be from FT not being inferred correctly when it crosses the function boundary from make_update_jacobian to update_jacobian!

---

• I have read and checked the items on the review checklist.

[kmdeck]

I wonder if this should have a clearer name, or if we need this function at all. We can get rid of this function if we move the update_aux! and update_boundary_fluxes! function calls that make_tendency_jacobian does before it calls update jacobian into update_jacobian.

[juliasloan25]

Good point! I think it makes more sense to get rid of make_tendency_jacobian, and just have make_update_jacobian

[juliasloan25]

I realized that now we're updating the aux and boundary fluxes in exp_tendency!, imp_tendency!, and update_jacobian!. I think this is fine for now, but we might want to avoid recomputing these multiple times per timestep. Dennis mentioned that ClimaTimeSteppers has the capability to take in an update cache function, which will be called once per timestep (Atmos leverages this to call set_precomputed_quantities! as few times as possible)., so we could later on switch to using that as a performance improvement

[juliasloan25]

After removing tendency_jacobian!, we can't directly compare the jacobian entries in test/shared_utilities/implicit_timestepping/richards_model.jl, since update_jacobian! also calls update_aux and boundary_fluxes. I wonder if it's still worth keeping separate functions for cleaner structure and more granular testing

[kmdeck]

This looks good! I left a couple of minor comments. nice job @juliasloan25 @dennisYatunin !!

[juliasloan25]

failing because of a change introduced in CliMA/ClimaCore.jl#1653 - need to debug this

[juliasloan25]

failing because of a change introduced in CliMA/ClimaCore.jl#1653 - need to debug this

fixed in CliMA/ClimaCore.jl#1750

[kmdeck]

looks great! thank you!!