From 5c3277338525a00a023f6e0f99ae23062b11eada Mon Sep 17 00:00:00 2001 From: kmdeck Date: Tue, 9 May 2023 08:16:52 -0700 Subject: [PATCH] jacobian functions added richards tridiagonal w and bc added unit tests and multi column example update test toml doc strings --- Project.toml | 1 + docs/src/APIs/SharedUtilities.md | 4 + docs/src/APIs/Soil.md | 7 + experiments/Manifest.toml | 4 +- src/ClimaLSM.jl | 1 + src/SharedUtilities/boundary_conditions.jl | 4 +- src/SharedUtilities/implicit_tendencies.jl | 72 +++++++ src/SharedUtilities/utils.jl | 7 + src/Soil/Soil.jl | 5 +- src/Soil/boundary_conditions.jl | 86 ++++++++- src/Soil/energy_hydrology.jl | 28 +-- src/Soil/rre.jl | 158 +++++++++++++++- src/Soil/soil_hydrology_parameterizations.jl | 20 +- test/Project.toml | 1 + test/implicit_timestepping/richards_model.jl | 189 +++++++++++++++++++ test/runtests.jl | 1 + 16 files changed, 565 insertions(+), 23 deletions(-) create mode 100644 src/SharedUtilities/implicit_tendencies.jl create mode 100644 test/implicit_timestepping/richards_model.jl diff --git a/Project.toml b/Project.toml index 103c053047..0aa7115dfa 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SurfaceFluxes = "49b00bb7-8bd4-4f2b-b78c-51cd0450215f" diff --git a/docs/src/APIs/SharedUtilities.md b/docs/src/APIs/SharedUtilities.md index d081521246..f201c8f1f6 100644 --- a/docs/src/APIs/SharedUtilities.md +++ b/docs/src/APIs/SharedUtilities.md @@ -48,6 +48,10 @@ ClimaLSM.BottomBoundary ClimaLSM.boundary_flux ClimaLSM.diffusive_flux ClimaLSM.get_Δz +ClimaLSM.make_tendency_jacobian +ClimaLSM.make_update_jacobian +ClimaLSM.∂tendencyBC∂Y +ClimaLSM.AbstractTridiagonalW ``` ## Drivers diff --git a/docs/src/APIs/Soil.md b/docs/src/APIs/Soil.md index 0595a17177..840d981734 100644 --- a/docs/src/APIs/Soil.md +++ b/docs/src/APIs/Soil.md @@ -36,6 +36,7 @@ ClimaLSM.Soil.temperature_from_ρe_int ClimaLSM.Soil.thermal_conductivity ClimaLSM.Soil.phase_change_source ClimaLSM.Soil.thermal_time +ClimaLSM.Soil.dψdϑ ``` ## Soil Parameters @@ -56,4 +57,10 @@ ClimaLSM.Soil.FreeDrainage ClimaLSM.Soil.AtmosDrivenFluxBC ClimaLSM.Soil.AbstractSoilSource ClimaLSM.Soil.PhaseChange +``` + +## Soil Jacobian Structures + +```@docs +ClimaLSM.Soil.RichardsTridiagonalW ``` \ No newline at end of file diff --git a/experiments/Manifest.toml b/experiments/Manifest.toml index ab2bbb9b87..00342177aa 100644 --- a/experiments/Manifest.toml +++ b/experiments/Manifest.toml @@ -193,10 +193,10 @@ uuid = "d934ef94-cdd4-4710-83d6-720549644b70" version = "0.3.7" [[deps.ClimaLSM]] -deps = ["ArtifactWrappers", "CLIMAParameters", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "Dates", "DocStringExtensions", "Insolation", "IntervalSets", "JLD2", "NCDatasets", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnPack"] +deps = ["ArtifactWrappers", "CLIMAParameters", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "Dates", "DocStringExtensions", "Insolation", "IntervalSets", "JLD2", "LinearAlgebra", "NCDatasets", "StaticArrays", "SurfaceFluxes", "Thermodynamics", "UnPack"] path = ".." uuid = "7884a58f-fab6-4fd0-82bb-ecfedb2d8430" -version = "0.2.2" +version = "0.2.3" [[deps.CloseOpenIntervals]] deps = ["ArrayInterface", "Static"] diff --git a/src/ClimaLSM.jl b/src/ClimaLSM.jl index 3c46e425d1..618c685ad1 100644 --- a/src/ClimaLSM.jl +++ b/src/ClimaLSM.jl @@ -16,6 +16,7 @@ include("SharedUtilities/drivers.jl") include("SharedUtilities/utils.jl") include("SharedUtilities/boundary_conditions.jl") include("SharedUtilities/sources.jl") +include("SharedUtilities/implicit_tendencies.jl") include("Bucket/Bucket.jl") export make_interactions_update_aux export domain diff --git a/src/SharedUtilities/boundary_conditions.jl b/src/SharedUtilities/boundary_conditions.jl index 6029eebdea..4e0ad8a48f 100644 --- a/src/SharedUtilities/boundary_conditions.jl +++ b/src/SharedUtilities/boundary_conditions.jl @@ -17,13 +17,13 @@ convenient conditions. abstract type AbstractBC end """ - AbstractBoundary{} + AbstractBoundary An abstract type to indicate which boundary we are doing calculations for. Currently, we support the top boundary (TopBoundary) and bottom boundary (BottomBoundary). """ -abstract type AbstractBoundary{} end +abstract type AbstractBoundary end """ diff --git a/src/SharedUtilities/implicit_tendencies.jl b/src/SharedUtilities/implicit_tendencies.jl new file mode 100644 index 0000000000..9096caee6a --- /dev/null +++ b/src/SharedUtilities/implicit_tendencies.jl @@ -0,0 +1,72 @@ +export make_tendency_jacobian, + make_update_jacobian, AbstractTridiagonalW, ∂tendencyBC∂Y + + +""" + make_tendency_jacobian(model::AbstractModel) + +Creates and returns a function which updates the auxiliary +variables `p` in place and then updates the entries of the +Jacobian matrix `W` for the `model` in place. + +The default is that no updates are required, no implicit tendency is +present, and hence the timestepping is entirely explicit. + +Note that the returned function `tendency_jacobian!` should be +used as `Wfact!` in `ClimaTimeSteppers.jl` and `OrdinaryDiffEq.jl`. +""" +function make_tendency_jacobian(model::AbstractModel) + update_aux! = make_update_aux(model) + update_jacobian! = make_update_jacobian(model) + function tendency_jacobian!(W, Y, p, dtγ, t) + update_aux!(p, Y, t) + update_jacobian!(W, Y, p, dtγ, t) + end + return tendency_jacobian! +end + +""" + make_update_jacobian(model::AbstractModel) + +Creates and returns a function which updates the entries +of the Jacobian matrix `W` in place. + +If the implicit tendency function is given by +`T!(dY, Y, p, t) = make_implicit_tendency(model)`, the Jacobian +should be given by `W_{i,j}! = ∂T!_i/∂Y_j`, where `Y_j` is the +`j-th` state variable +and `T!_i` is the implicit tendency of the `i-th` state variable. + +The default is that no updates are required, no implicit tendency is +present, and hence the timestepping is entirely explicit. +""" +function make_update_jacobian(model::AbstractModel) + function update_jacobian!(W, Y, p, dtγ, t) end + return update_jacobian +end + +""" + ∂tendencyBC∂Y(::AbstractModel, + ::AbstractBC, + ::AbstractBoundary, + _...)::Union{ClimaCore.Fields.FieldVector, Nothing} + +A function stub which returns the derivative of the implicit tendency +term of the `model` arising from the boundary condition, +with respect to the state Y. +""" +function ∂tendencyBC∂Y( + ::AbstractModel, + ::AbstractBC, + ::AbstractBoundary, + _..., +)::Union{ClimaCore.Fields.FieldVector, Nothing} end + +""" + AbstractTridiagonalW + +An abstract type for tridiagonal Jacobian matrices. +""" +abstract type AbstractTridiagonalW end + +Base.similar(w::AbstractTridiagonalW) = w diff --git a/src/SharedUtilities/utils.jl b/src/SharedUtilities/utils.jl index a579f33514..4cc424c4c3 100644 --- a/src/SharedUtilities/utils.jl +++ b/src/SharedUtilities/utils.jl @@ -10,3 +10,10 @@ function heaviside(x::FT)::FT where {FT} return FT(0.0) end end + +""" + to_scalar_coefs(vector_coefs) + +Helper function used in computing tendencies of vertical diffusion terms. +""" +to_scalar_coefs(vector_coefs) = map(vector_coef -> vector_coef.u₃, vector_coefs) diff --git a/src/Soil/Soil.jl b/src/Soil/Soil.jl index 043300700e..e0e83a0cbf 100644 --- a/src/Soil/Soil.jl +++ b/src/Soil/Soil.jl @@ -56,18 +56,20 @@ for interactions between components. =# using ClimaLSM -using UnPack using DocStringExtensions +using LinearAlgebra using ClimaCore import ..Parameters as LSMP import ClimaCore: Fields, Operators, Geometry, Spaces using Thermodynamics import ClimaLSM.Domains: Column, HybridBox, SphericalShell +using ClimaLSM: AbstractTridiagonalW import ClimaLSM: AbstractModel, make_update_aux, make_rhs, + make_update_jacobian, prognostic_vars, auxiliary_vars, name, @@ -84,6 +86,7 @@ import ClimaLSM: surface_height export RichardsModel, RichardsParameters, + RichardsTridiagonalW, EnergyHydrology, EnergyHydrologyParameters, AbstractSoilModel, diff --git a/src/Soil/boundary_conditions.jl b/src/Soil/boundary_conditions.jl index f20602d8af..a0d0edd64c 100644 --- a/src/Soil/boundary_conditions.jl +++ b/src/Soil/boundary_conditions.jl @@ -159,7 +159,7 @@ function ClimaLSM.boundary_flux( ψ_c = Fields.level(p.soil.ψ, p_len) # Calculate pressure head using boundary condition - @unpack vg_α, vg_n, vg_m, θ_r, ν, S_s = params + (; vg_α, vg_n, vg_m, θ_r, ν, S_s) = params θ_bc = rre_bc.bc(p, t) ψ_bc = @. pressure_head(vg_α, vg_n, vg_m, θ_r, θ_bc, ν, S_s) @@ -186,7 +186,7 @@ function ClimaLSM.boundary_flux( ψ_c = Fields.level(p.soil.ψ, 1) # Calculate pressure head using boundary condition - @unpack vg_α, vg_n, vg_m, θ_r, ν, S_s = params + (; vg_α, vg_n, vg_m, θ_r, ν, S_s) = params θ_bc = rre_bc.bc(p, t) ψ_bc = @. pressure_head(vg_α, vg_n, vg_m, θ_r, θ_bc, ν, S_s) @@ -270,3 +270,85 @@ function soil_boundary_fluxes(bc::NamedTuple, boundary, model, Y, Δz, p, t) return ClimaLSM.boundary_flux(bc.water, boundary, Δz, p, t, params), ClimaLSM.boundary_flux(bc.heat, boundary, Δz, p, t, params) end + +""" + ClimaLSM.∂tendencyBC∂Y( + model::RichardsModel, + ::MoistureStateBC, + boundary::ClimaLSM.TopBoundary, + Δz, + Y, + p, + t, +) + +Computes and returns the derivative of the part of the +implicit tendency in the top layer, due to the boundary +condition, with respect to the state variable in the top layer. + +For a diffusion equation like Richards equation with a single state +variable, this is given by +`∂T_N∂Y_N = [-∂/∂z(∂F_bc/∂Y_N)]_N`, where `N` indicates the top +layer cell index. +""" +function ClimaLSM.∂tendencyBC∂Y( + model::RichardsModel, + ::MoistureStateBC, + boundary::ClimaLSM.TopBoundary, + Δz, + Y, + p, + t, +) + (; ν, vg_α, vg_n, vg_m, S_s, θ_r) = model.parameters + fs = ClimaLSM.Domains.obtain_face_space(model.domain.space) + face_len = ClimaCore.Utilities.PlusHalf(ClimaCore.Spaces.nlevels(fs) - 1) + interpc2f_op = Operators.InterpolateC2F( + bottom = Operators.Extrapolate(), + top = Operators.Extrapolate(), + ) + K = Fields.level(interpc2f_op.(p.soil.K), face_len) + ϑ_l = Fields.level(interpc2f_op.(Y.soil.ϑ_l), face_len) + return ClimaCore.Fields.FieldVector(; + :soil => (; + :ϑ_l => @. -K / Δz * dψdϑ(ϑ_l, ν, θ_r, vg_α, vg_n, vg_m, S_s) / + (2 * Δz) + ), + ) +end + +""" + ClimaLSM.∂tendencyBC∂Y( + ::AbstractSoilModel, + ::AbstractSoilBC, + boundary::ClimaLSM.TopBoundary, + Δz, + Y, + p, + t, +) + +A default method which computes and returns the zero for the +derivative of the part of the +implicit tendency in the top layer, due to the boundary +condition, with respect to the state variable in the top layer. + +For a diffusion equation like Richards equation with a single state +variable, this is given by +`∂T_N∂Y_N = [-∂/∂z(∂F_bc/∂Y_N)]_N`, where `N` indicates the top +layer cell index. + +If `F_bc` can be approximated as independent of `Y_N`, the derivative +is zero. +""" +function ClimaLSM.∂tendencyBC∂Y( + ::AbstractSoilModel, + ::AbstractSoilBC, + boundary::ClimaLSM.TopBoundary, + Δz, + Y, + p, + t, +) + return ClimaCore.Fields.FieldVector(; :soil => (; :ϑ_l => zeros(axes(Δz)))) +end diff --git a/src/Soil/energy_hydrology.jl b/src/Soil/energy_hydrology.jl index 95bf2b23d5..94e4776c1a 100644 --- a/src/Soil/energy_hydrology.jl +++ b/src/Soil/energy_hydrology.jl @@ -318,18 +318,20 @@ This has been written so as to work with Differential Equations.jl. """ function ClimaLSM.make_update_aux(model::EnergyHydrology) function update_aux!(p, Y, t) - @unpack ν, - vg_α, - vg_n, - vg_m, - K_sat, - S_s, - θ_r, - Ω, - γ, - γT_ref, - κ_sat_frozen, - κ_sat_unfrozen = model.parameters + (; + ν, + vg_α, + vg_n, + vg_m, + K_sat, + S_s, + θ_r, + Ω, + γ, + γT_ref, + κ_sat_frozen, + κ_sat_unfrozen, + ) = model.parameters p.soil.θ_l = volumetric_liquid_fraction.(Y.soil.ϑ_l, ν .- Y.soil.θ_i) p.soil.κ = @@ -406,7 +408,7 @@ function ClimaLSM.source!( model, ) where {FT} params = model.parameters - @unpack ν, earth_param_set = params + (; ν, earth_param_set) = params _ρ_l = FT(LSMP.ρ_cloud_liq(earth_param_set)) _ρ_i = FT(LSMP.ρ_cloud_ice(earth_param_set)) ρc = volumetric_heat_capacity.(p.soil.θ_l, Y.soil.θ_i, Ref(params)) diff --git a/src/Soil/rre.jl b/src/Soil/rre.jl index 39be3cf035..a8734b58ce 100644 --- a/src/Soil/rre.jl +++ b/src/Soil/rre.jl @@ -87,7 +87,7 @@ This has been written so as to work with Differential Equations.jl. """ function ClimaLSM.make_rhs(model::RichardsModel) function rhs!(dY, Y, p, t) - @unpack ν, vg_α, vg_n, vg_m, K_sat, S_s, θ_r = model.parameters + (; ν, vg_α, vg_n, vg_m, K_sat, S_s, θ_r) = model.parameters z = ClimaCore.Fields.coordinate_field(model.domain.space).z Δz_top, Δz_bottom = get_Δz(z) @@ -204,7 +204,7 @@ This has been written so as to work with Differential Equations.jl. """ function ClimaLSM.make_update_aux(model::RichardsModel) function update_aux!(p, Y, t) - @unpack ν, vg_α, vg_n, vg_m, K_sat, S_s, θ_r = model.parameters + (; ν, vg_α, vg_n, vg_m, K_sat, S_s, θ_r) = model.parameters @. p.soil.K = hydraulic_conductivity( K_sat, vg_m, @@ -214,3 +214,157 @@ function ClimaLSM.make_update_aux(model::RichardsModel) end return update_aux! end + + +""" + RichardsTridiagonalW{R, J, W, T} <: ClimaLSM.AbstractTridiagonalW + +A struct containing the necessary information for constructing a tridiagonal +Jacobian matrix (`W`) for solving Richards equation, treating only the vertical +diffusion term implicitly. + +Note that the diagonal, upper diagonal, and lower diagonal entry values +are stored in this struct and updated in place. +$(DocStringExtensions.FIELDS) +""" +struct RichardsTridiagonalW{R, J, JA, T, A} <: ClimaLSM.AbstractTridiagonalW + "Reference to dtγ, which is specified by the ODE solver" + dtγ_ref::R + "Diagonal entries of the Jacobian stored as a ClimaCore.Fields.Field" + ∂ϑₜ∂ϑ::J + "Array of tridiagonal matrices containing W for each column" + W_column_arrays::JA + "An allocated cache used to evaluate ldiv!" + temp1::T + "An allocated cache used to evaluate ldiv!" + temp2::T + "A flag indicating whether this struct is used to compute Wfact_t or Wfact" + transform::Bool + "A pre-allocated cache storing ones on the face space" + ones_face_space::A +end + +""" + RichardsTridiagonalW( + Y::ClimaCore.Fields.FieldVector; + transform::Bool = false +) + +Outer constructor for the RichardsTridiagonalW Jacobian +matrix struct. + +Initializes all variables to zeros. +""" +function RichardsTridiagonalW( + Y::ClimaCore.Fields.FieldVector; + transform::Bool = false, +) + FT = eltype(Y.soil.ϑ_l) + center_space = axes(Y.soil.ϑ_l) + N = Spaces.nlevels(center_space) + + tridiag_type = Operators.StencilCoefs{-1, 1, NTuple{3, FT}} + ∂ϑₜ∂ϑ = Fields.Field(tridiag_type, center_space) + + ∂ϑₜ∂ϑ.coefs[1] .= FT(0) + ∂ϑₜ∂ϑ.coefs[2] .= FT(0) + ∂ϑₜ∂ϑ.coefs[3] .= FT(0) + + W_column_arrays = [ + LinearAlgebra.Tridiagonal( + zeros(FT, N - 1) .+ FT(0), + zeros(FT, N) .+ FT(0), + zeros(FT, N - 1) .+ FT(0), + ) for _ in 1:Threads.nthreads() + ] + dtγ_ref = Ref(FT(0)) + temp1 = similar(Y.soil.ϑ_l) + temp1 .= FT(0) + temp2 = similar(Y.soil.ϑ_l) + temp2 .= FT(0) + + face_space = ClimaLSM.Domains.obtain_face_space(center_space) + ones_face_space = ones(face_space) + + return RichardsTridiagonalW( + dtγ_ref, + ∂ϑₜ∂ϑ, + W_column_arrays, + temp1, + temp2, + transform, + ones_face_space, + ) +end + + +""" + ClimaLSM.make_update_jacobian(model::RichardsModel) + +Creates and returns the update_jacobian! function for RichardsModel. + +Using this Jacobian with a backwards Euler timestepper is equivalent +to using the modified Picard scheme of Celia et al. (1990). +""" +function ClimaLSM.make_update_jacobian(model::RichardsModel) + function update_jacobian!(W::RichardsTridiagonalW, Y, p, dtγ, t) + FT = eltype(Y.soil.ϑ_l) + (; dtγ_ref, ∂ϑₜ∂ϑ) = W + dtγ_ref[] = dtγ + (; ν, vg_α, vg_n, vg_m, S_s, θ_r) = model.parameters + + divf2c_op = Operators.DivergenceF2C( + top = Operators.SetValue(Geometry.WVector.(FT(0))), + bottom = Operators.SetValue(Geometry.WVector.(FT(0))), + ) + divf2c_stencil = Operators.Operator2Stencil(divf2c_op) + gradc2f_op = Operators.GradientC2F( + top = Operators.SetGradient(Geometry.WVector.(FT(0))), + bottom = Operators.SetGradient(Geometry.WVector.(FT(0))), + ) + gradc2f_stencil = Operators.Operator2Stencil(gradc2f_op) + interpc2f_op = Operators.InterpolateC2F( + bottom = Operators.Extrapolate(), + top = Operators.Extrapolate(), + ) + compose = Operators.ComposeStencils() + + @. ∂ϑₜ∂ϑ = compose( + divf2c_stencil(Geometry.Covariant3Vector(W.ones_face_space)), + ( + interpc2f_op(p.soil.K) * ClimaLSM.to_scalar_coefs( + gradc2f_stencil( + ClimaLSM.Soil.dψdϑ( + Y.soil.ϑ_l, + ν, + θ_r, + vg_α, + vg_n, + vg_m, + S_s, + ), + ), + ) + ), + ) + # Hardcoded for single column: FIX! + z = Fields.coordinate_field(axes(Y.soil.ϑ_l)).z + Δz_top, Δz_bottom = get_Δz(z) + ∂T_bc∂YN = ClimaLSM.∂tendencyBC∂Y( + model, + model.boundary_conditions.top.water, + ClimaLSM.TopBoundary(), + Δz_top, + Y, + p, + t, + ) + #TODO: allocate space in W? See how final implementation of stencils with boundaries works out + N = ClimaCore.Spaces.nlevels(axes(Y.soil.ϑ_l)) + parent(ClimaCore.Fields.level(∂ϑₜ∂ϑ.coefs.:2, N)) .= + parent(ClimaCore.Fields.level(∂ϑₜ∂ϑ.coefs.:2, N)) .+ + parent(∂T_bc∂YN.soil.ϑ_l) + + end + return update_jacobian! +end diff --git a/src/Soil/soil_hydrology_parameterizations.jl b/src/Soil/soil_hydrology_parameterizations.jl index 5ce91fc5bb..f77c8ae51e 100644 --- a/src/Soil/soil_hydrology_parameterizations.jl +++ b/src/Soil/soil_hydrology_parameterizations.jl @@ -5,7 +5,25 @@ export volumetric_liquid_fraction, pressure_head, hydraulic_conductivity, impedance_factor, - viscosity_factor + viscosity_factor, + dψdϑ + +""" + dψdϑ(ϑ, ν, θ_r, vg_α, vg_n, vg_m, S_s) + +Computes and returns the derivative of the pressure head +with respect to ϑ for the van Genuchten formulation. +""" +function dψdϑ(ϑ, ν, θ_r, vg_α, vg_n, vg_m, S_s) + S = effective_saturation(ν, ϑ, θ_r) + if S < 1.0 + return 1.0 / (vg_α * vg_m * vg_n) / (ν - θ_r) * + (S^(-1 / vg_m) - 1)^(1 / vg_n - 1) * + S^(-1 / vg_m - 1) + else + return 1.0 / S_s + end +end """ volumetric_liquid_fraction(ϑ_l::FT, ν_eff::FT) where {FT} diff --git a/test/Project.toml b/test/Project.toml index c6e56cac17..a846531b3b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,7 @@ ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884" ClimaLSM = "7884a58f-fab6-4fd0-82bb-ecfedb2d8430" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Insolation = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/test/implicit_timestepping/richards_model.jl b/test/implicit_timestepping/richards_model.jl new file mode 100644 index 0000000000..6883afcd05 --- /dev/null +++ b/test/implicit_timestepping/richards_model.jl @@ -0,0 +1,189 @@ +using Test +using LinearAlgebra +using ClimaCore +import CLIMAParameters as CP +using ClimaLSM +using ClimaLSM.Domains: Column, HybridBox +using ClimaLSM.Soil + +import ClimaLSM +import ClimaLSM.Parameters as LSMP + +include(joinpath(pkgdir(ClimaLSM), "parameters", "create_parameters.jl")) + +@testset "Richards Jacobian entries, Moisture BC" begin + FT = Float64 + ν = FT(0.495) + K_sat = FT(0.0443 / 3600 / 100) # m/s + S_s = FT(1e-3) #inverse meters + vg_n = FT(1.43) + vg_α = FT(2.6) # inverse meters + vg_m = FT(1) - FT(1) / vg_n + θ_r = FT(0.124) + zmax = FT(0) + zmin = FT(-1.5) + nelems = 150 + soil_domains = [ + Column(; zlim = (zmin, zmax), nelements = nelems), + HybridBox(; + xlim = (0.0, 1.0), + ylim = (0.0, 1.0), + zlim = (zmin, zmax), + nelements = (1, 1, nelems), + npolynomial = 3, + ), + ] + top_state_bc = MoistureStateBC((p, t) -> eltype(t)(ν - 1e-3)) + bot_flux_bc = FreeDrainage() + sources = () + boundary_states = + (; top = (water = top_state_bc,), bottom = (water = bot_flux_bc,)) + params = Soil.RichardsParameters{FT}(ν, vg_α, vg_n, vg_m, K_sat, S_s, θ_r) + + for domain in soil_domains + soil = Soil.RichardsModel{FT}(; + parameters = params, + domain = domain, + boundary_conditions = boundary_states, + sources = sources, + ) + + Y, p, coords = initialize(soil) + Y.soil.ϑ_l .= FT(0.24) + W = RichardsTridiagonalW(Y) + Wfact! = make_tendency_jacobian(soil) + dtγ = FT(1.0) + Wfact!(W, Y, p, dtγ, FT(0.0)) + + K_ic = hydraulic_conductivity( + K_sat, + vg_m, + effective_saturation(ν, FT(0.24), θ_r), + ) + dz = FT(0.01) + dψdϑ_ic = dψdϑ(FT(0.24), ν, θ_r, vg_α, vg_n, vg_m, S_s) + + @test all(parent(W.temp2) .== FT(0.0)) + @test all(parent(W.temp2) .== FT(0.0)) + @test W.transform == false + @test typeof(W.W_column_arrays) <: + Vector{LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}} + @test length(W.W_column_arrays) == 1 + if typeof(domain) <: Column + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:1)[2:end] .≈ + parent(W.∂ϑₜ∂ϑ.coefs.:3)[1:(end - 1)], + ) + @test parent(W.∂ϑₜ∂ϑ.coefs.:1)[1] == FT(0) + @test parent(W.∂ϑₜ∂ϑ.coefs.:3)[end] == FT(0) + @test all(parent(W.∂ϑₜ∂ϑ.coefs.:1)[2:end] .≈ K_ic / dz^2 * dψdϑ_ic) + @test parent(W.∂ϑₜ∂ϑ.coefs.:2)[1] .≈ -K_ic / dz^2 * dψdϑ_ic + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[2:(end - 1)] .≈ + -2 * K_ic / dz^2 * dψdϑ_ic, + ) + @test parent(W.∂ϑₜ∂ϑ.coefs.:2)[end] .≈ + -K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic + elseif typeof(domain) <: HybridBox + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:1)[2:end, :, 1, 1, 1] .≈ + parent(W.∂ϑₜ∂ϑ.coefs.:3)[1:(end - 1), :, 1, 1, 1], + ) + @test all(parent(W.∂ϑₜ∂ϑ.coefs.:1)[1, :, 1, 1, 1] .== FT(0)) + @test all(parent(W.∂ϑₜ∂ϑ.coefs.:3)[end, :, 1, 1, 1] .== FT(0)) + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:1)[2:end, :, 1, 1, 1] .≈ + K_ic / dz^2 * dψdϑ_ic, + ) + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[1, :, 1, 1, 1] .≈ + -K_ic / dz^2 * dψdϑ_ic, + ) + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[2:(end - 1), :, 1, 1, 1] .≈ + -2 * K_ic / dz^2 * dψdϑ_ic, + ) + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[end, :, 1, 1, 1] .≈ + -K_ic / dz^2 * dψdϑ_ic - K_ic / (dz * dz / 2) * dψdϑ_ic, + ) + end + + end +end + + +@testset "Richards Jacobian entries, Flux BC" begin + FT = Float64 + ν = FT(0.495) + K_sat = FT(0.0443 / 3600 / 100) # m/s + S_s = FT(1e-3) #inverse meters + vg_n = FT(1.43) + vg_α = FT(2.6) # inverse meters + vg_m = FT(1) - FT(1) / vg_n + θ_r = FT(0.124) + zmax = FT(0) + zmin = FT(-1.5) + nelems = 150 + soil_domains = [ + Column(; zlim = (zmin, zmax), nelements = nelems), + HybridBox(; + xlim = (0.0, 1.0), + ylim = (0.0, 1.0), + zlim = (zmin, zmax), + nelements = (1, 1, nelems), + npolynomial = 3, + ), + ] + top_flux_bc = FluxBC((p, t) -> eltype(t)(-K_sat)) + bot_flux_bc = FreeDrainage() + sources = () + boundary_states = + (; top = (water = top_flux_bc,), bottom = (water = bot_flux_bc,)) + params = Soil.RichardsParameters{FT}(ν, vg_α, vg_n, vg_m, K_sat, S_s, θ_r) + for domain in soil_domains + soil = Soil.RichardsModel{FT}(; + parameters = params, + domain = domain, + boundary_conditions = boundary_states, + sources = sources, + ) + + Y, p, coords = initialize(soil) + Y.soil.ϑ_l .= FT(0.24) + W = RichardsTridiagonalW(Y) + Wfact! = make_tendency_jacobian(soil) + dtγ = FT(1.0) + Wfact!(W, Y, p, dtγ, FT(0.0)) + + K_ic = hydraulic_conductivity( + K_sat, + vg_m, + effective_saturation(ν, FT(0.24), θ_r), + ) + dz = FT(0.01) + dψdϑ_ic = dψdϑ(FT(0.24), ν, θ_r, vg_α, vg_n, vg_m, S_s) + if typeof(domain) <: Column + @test parent(W.∂ϑₜ∂ϑ.coefs.:2)[1] .≈ -K_ic / dz^2 * dψdϑ_ic + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[2:(end - 1)] .≈ + -2 * K_ic / dz^2 * dψdϑ_ic, + ) + @test parent(W.∂ϑₜ∂ϑ.coefs.:2)[end] .≈ -K_ic / dz^2 * dψdϑ_ic + elseif typeof(domain) <: HybridBox + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[1, :, 1, 1, 1] .≈ + -K_ic / dz^2 * dψdϑ_ic, + ) + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[2:(end - 1), :, 1, 1, 1] .≈ + -2 * K_ic / dz^2 * dψdϑ_ic, + ) + @test all( + parent(W.∂ϑₜ∂ϑ.coefs.:2)[end, :, 1, 1, 1] .≈ + -K_ic / dz^2 * dψdϑ_ic, + ) + end + + end +end diff --git a/test/runtests.jl b/test/runtests.jl index fa4527d735..c4692005ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ include("./Snow/parameterizations.jl") +include("./implicit_timestepping/richards_model.jl") include("./Vegetation/test_bigleaf_parameterizations.jl") include("./Vegetation/canopy_model.jl") include("./Soil/Biogeochemistry/co2_parameterizations.jl")