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

• move jacobian struct and related implicit timestepping functions into src/Soil
• add imp_tendency function for EnergyHydrology
• update all instances of EnergyHydrology to use implicit tendency (docs, tests)
• 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
• add plots from docs

Check run output

docs

• soil/canopy tutorial

• coarse sand evaporation

• Gilat Loess evaporation

• freezing front

• energy hydrology

• sublimation

experiments (visually inspected)

• vaira
  • plots (ET, SHF, temperature) look smoother now! (probably from timestepper, dt changes)
  • slower because dt decreased by factor of 2
• ozark
  • plots are nearly identical
• niwot
  • plots are nearly identical
• harvard
  • plots look identical
• 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
• 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!

[juliasloan25]

Interestingly, this file was already using ARS111 and an IMEXAlgorithm (lines 18-32 of this file)

[kmdeck]

oh! that's totally it. We were using explicit euler before, but the only option for that is IMEX ars111 (with imp_tendency! doing nothing). there is no pure explicit euler option in ClimaTimestppers. this matters just because the expected energy and water changes we compute are computed assuming explicit euler timesteps.

like, if Var(t+dt)= Var(t) + function(fluxes(t)) * dt, we can check conservation by comparing Var(t+dt) - Var(t) = function(fluxes(t)) * dt (but over many steps - and just expect roundoff error)

but if the timestepper changes, this is not the case. if it was implicit, we could do
Var(t+dt) - Var(t) = function(fluxes(t+dt)) * dt

basically im not really sure how to compute the "truth" for higher order imex schemes.

[kmdeck]

maybe we can ask Dennis. We mostly want to verify that d(energy in land) = energy in - energy out over the whole sim.

I think the important thing is that this is not growing. so I am inclined to think this is OK, we just are computing what we want

[juliasloan25]

Ohh I see! Would it be possible to use the equation you gave for the explicit case for explicitly-stepped variables, and the implicit one for implicitly-stepped variables? Or does that not make sense because we have to consider both variables and fluxes?

[kmdeck]

im not sure! maybe we can discuss more tomorrow? Im not worried about this really for this PR though... we maybe need to rethink that conversation test in general

[juliasloan25]

happy to discuss tomorrow! I opened an issue for this too so we can keep track of it #662

[kmdeck]

Should this do nothing? or should it return minus identity matrix?

or do we never rely on this default (and never will)? if so, maybe we can note that this is only a stub and must be extended or ? not sure.

[juliasloan25]

I don't think it needs to return -I. I know that we had previously set the jacobian blocks for explicitly-stepped variables to -I in update_jacobian!, but I realized that we just can initialize these as -I in our ImplicitEquationJacobian constructor. Then we just do nothing for these blocks in update_jacobian!.

And now that I'm writing this out, I notice that update_jacobian! is now identical for RichardsModel and EnergyHydrology with the jacobian approximation we're using right now. I can move it into src/Soil and change the type dispatch to use one method for both, or we could just keep it split up since we'll be changing the EnergyHydrology version as we iterate on the jacobian anyway. I'm fine with either.

I'll add a note to the generic update_jacobian! to explain that it needs to be extended

[kmdeck]

I think this looks great as a starting point! Mostly comments about keeping track of our to-do list afterwards