initial full soil imp solver

[juliasloan25]

Purpose

Adds a simple implicit tendency for the EnergyHydrology model, which only treats the liquid water content implicitly (using the functions already defined for RichardsModel)

closes #538

Requirements for a model with implicit and explicit prognostic variables

Changes in src

• all prognostic variables must have an entry (@name(var), @name(var)) => block in the Jacobian. The block should be the negative identity matrix for variables that are stepped explicitly, and a tridiagonal matrix for variables that are stepped implicitly
• in the explicit tendency, we need to reset dY.model .= 0 for implicitly-stepped variables
• in the implicit tendency, we need to reset dY.model .= 0 for explicitly-stepped variables

  Changes in driver

• in addition to exp_tendency!, create imp_tendency!, tendency_jacobian!, and jac_kwargs
• change the timestepper from an explicit to an imex method (e.g. RK4 -> ARS111 or ARS343)
• change the type of the ODE algorithm from CTS.ExplicitAlgorithm to CTS.IMEXAlgorithm, and add the NewtonsMethod argument
• when constructing the problem, pass an ODE function for both the explicit and implicit tendencies (e.g. add SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...) to CTS.ClimaODEFunction call)

Next steps (separate PRs)

• [ ] add implicit tendency for energy dT/dI
• [ ] include boundary condition contributions
• [ ] include contribution from drhoe_int / dtheta_l

Content

• [x] move jacobian struct and related implicit timestepping functions into src/Soil
• [x] add imp_tendency function for EnergyHydrology
• [x] update all instances of EnergyHydrology to use implicit tendency (docs, tests)
• [x] use implicit solver with SoilCanopyModel
  • add jacobian terms for canopy, soilco2 prognostic variables in ImplicitEquationJacobian
  • set implicit tendency for canopy, soilco2 prognostic variables to 0 in imp_tendency! (replace do_nothing_imp_tendency ? --> update compute_imp_tendency (and exp_tendency defaults to set dY .= 0)
  • add make_jac_tendency function for AbstractLandModel - calling functions of component models, analogous to ALM update_aux
• [x] add plots from docs

Check run output

docs

• [x] soil/canopy tutorial

  

• [x] coarse sand evaporation

  

• [x] Gilat Loess evaporation

  

• [x] freezing front

  

• [x] energy hydrology

  

• [x] sublimation

  

experiments (visually inspected)

• [x] vaira
  • plots (ET, SHF, temperature) look smoother now! (probably from timestepper, dt changes)
  • slower because dt decreased by factor of 2
• [x] ozark
  • plots are nearly identical
• [x] niwot
  • plots are nearly identical
• [x] harvard
  • plots look identical
• [x] ozark PFT
  • plots are nearly identical
• [ ] ozark conservation
  • energy conservation fractional error shows same pattern, but larger by ~10^13 orders of magnitude
  • water conservation fractional error larger by ~10^13 orders of magnitude
• [x] soil/canopy lsm performance
  • solve takes longer now

---

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

[kmdeck]

for the freezing front test, maybe we can try a higher order IMEX stepper? E.g. ARS343

It may be that using explicit euler as the explicit stepper is not sufficient or requires a very small step, while when we were using RK4 we could use a bigger step.

[juliasloan25]

for the freezing front test, maybe we can try a higher order IMEX stepper? E.g. ARS343

It may be that using explicit euler as the explicit stepper is not sufficient or requires a very small step, while when we were using RK4 we could use a bigger step.

Even using ARS343 with a timestep as small as 0.1s (was 60s when explicit), all prognostic variables are NaN at the second timestep. Maybe we do need to add more terms to our jacobian approximation?

[kmdeck]

for the freezing front test, maybe we can try a higher order IMEX stepper? E.g. ARS343 It may be that using explicit euler as the explicit stepper is not sufficient or requires a very small step, while when we were using RK4 we could use a bigger step.

Even using ARS343 with a timestep as small as 0.1s (was 60s when explicit), all prognostic variables are NaN at the second timestep. Maybe we do need to add more terms to our jacobian approximation?

Let's think about this a bit more. It seems odd that treating RHS terms all explicitly is OK, but moving one to be treated implicitly (which should only making things more stable) and leaving the others explicit causes errors. Do you want to meet this afternoon to try some debugging?

[juliasloan25]

for the freezing front test, maybe we can try a higher order IMEX stepper? E.g. ARS343 It may be that using explicit euler as the explicit stepper is not sufficient or requires a very small step, while when we were using RK4 we could use a bigger step.

Even using ARS343 with a timestep as small as 0.1s (was 60s when explicit), all prognostic variables are NaN at the second timestep. Maybe we do need to add more terms to our jacobian approximation?

Let's think about this a bit more. It seems odd that treating RHS terms all explicitly is OK, but moving one to be treated implicitly (which should only making things more stable) and leaving the others explicit causes errors. Do you want to meet this afternoon to try some debugging?

Yeah, good point. Want to meet at 1pm?

[kmdeck]

Yeah, good point. Want to meet at 1pm?

sure!