Merge branch 'bugfix/tempgrad' into 'master'

Fix bug in temperature gradient

See merge request sparcs/cryogrid/cryogridjulia!56

Project.toml

 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <brian.groenke@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.5.2"
+version = "0.5.3"
 
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"

src/Callbacks/Callbacks.jl

-module Callbacks
-
-import CryoGrid
-
-using CryoGrid.Drivers: HeatOnlyLandModel
-using CryoGrid.Numerics
-using CryoGrid.Land: Tile, getvar
-using CryoGrid.Utils
-
-using DiffEqCallbacks
-using IfElse
-
-"""
-CryoGridCallbackFunction{TState,TSetup}(setup, state)
-
-Helper type for defining callbacks on CryoGrid models. Given a Tile and some additional user-defined state type
-`TState`, the user can provide dispatches for `CryoGridCallbackFunction{TState}` that satisfy the relevant `DifferentialEquations.jl`
-callback function signature. For example:
-
-```julia
-struct MyState
-    # some state variables
-    ...
-end
-function (fn::CryoGridCallbackFunction{MyState})(u,p,t)
-    ...
-end
-function MyCallback(setup::Tile)
-    state = MyState(...)
-    fn = CryoGridCallbackFunction(setup, state)
-    # create and return SciML callback here
-end
-```
-"""
-struct CryoGridCallbackFunction{TState,TSetup}
-    setup::TSetup
-    state::TState
-    CryoGridCallbackFunction(setup::Tile, state::TState) where TState = new{TState, typeof(setup)}(setup, state)
-end
-(fn::CryoGridCallbackFunction)(u,p,t) = error("no method dispatch provided for callback function $(typeof(fn))")
-
-export CryoGridCallbackFunction
-
-include("courant_step.jl")
-
-end

src/Callbacks/courant_step.jl

-struct CFLHeatState{A}
-    Δt::A
-    default_dt::Float64
-end
-function (fn::CryoGridCallbackFunction{<:CFLHeatState{<:AbstractArray}})(u,p,t)
-    let Δt = fn.state.Δt,
-        Δx = dustrip(Δ(fn.setup.grid)),
-        Ceff = getvar(:Ceff, fn.setup, u), # apparent heat capacity
-        kc = getvar(:kc, fn.setup, u); # thermal cond. at grid centers
-        @. Δt = Utils.adstrip(Δx^2 * Ceff / kc)
-        Δt_min = minimum(Δt)
-        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, fn.state.default_dt)
-    end
-end
-function (fn::CryoGridCallbackFunction{<:CFLHeatState{Nothing}})(u,p,t)
-    let Δx = dustrip(Δ(fn.setup.grid)),
-        Ceff = getvar(:Ceff, fn.setup, u), # apparent heat capacity
-        kc = getvar(:kc, fn.setup, u); # thermal cond. at grid centers
-        Δt = Utils.adstrip(Δx^2 * Ceff / kc)
-        Δt_min = minimum(Δt)
-        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, fn.state.default_dt)
-    end
-end
-function CFLStepLimiter(setup::HeatOnlyLandModel; courant_number=1//2, default_dt=60.0, iip=true)
-    state = iip ? CFLHeatState(zero(dustrip(Δ(setup.grid))), default_dt) : CFLHeatState(nothing, default_dt)
-    fn = CryoGridCallbackFunction(setup, state)
-    StepsizeLimiter(fn; safety_factor=courant_number, max_step=true)
-end
-
-export CFLStepLimiter
\ No newline at end of file

src/CryoGrid.jl

@@ -28,7 +28,6 @@ include("Land/Land.jl")
 include("IO/InputOutput.jl")
 include("Diagnostics/Diagnostics.jl")
 include("Drivers/Drivers.jl")
-include("Callbacks/Callbacks.jl")
 
 using .Utils
 # Re-exported submodules
@@ -40,7 +39,6 @@ using .Utils
 @reexport using .Drivers
 @reexport using .InputOutput
 @reexport using .Diagnostics
-@reexport using .Callbacks
 
 # Re-exported packages
 @reexport using Dates: Dates, DateTime

src/Drivers/Drivers.jl

@@ -20,7 +20,7 @@ using Unitful
 
 import CryoGrid.Land: Tile, Stratigraphy, StratComponent
 
-export DefaultJac, TridiagJac, JacobianStyle, HeatOnlyLandModel
+export DefaultJac, TridiagJac, JacobianStyle, HeatOnlyTile
 
 """
     JacobianStyle
@@ -39,11 +39,13 @@ JacobianStyle(::Type{<:Tile}) = DefaultJac()
 # Auto-detect Jacobian sparsity for problems with one or more heat-only layers.
 # Note: This assumes that the processes/forcings on the boundary layers do not violate the tridiagonal structure!
 # Unfortunately, the Stratigraphy type signature is a bit nasty to work with :(
-const HeatOnlyLandModel = Tile{<:Stratigraphy{N,<:Tuple{TTop,Vararg{<:Union{<:StratComponent{<:SubSurface, <:CompoundProcess{<:Tuple{<:Heat}}},TBot}}}}} where {N,TTop,TBot}
-JacobianStyle(::Type{<:HeatOnlyLandModel}) = TridiagJac()
+const HeatOnlyTile = Tile{<:Stratigraphy{N,<:Tuple{TTop,Vararg{<:Union{<:StratComponent{<:SubSurface, <:CompoundProcess{<:Tuple{<:Heat}}},TBot}}}}} where {N,TTop,TBot}
+JacobianStyle(::Type{<:HeatOnlyTile}) = TridiagJac()
 
 # DiffEq/SciML driver
 export CryoGridProblem
 include("diffeq.jl")
+export CFLStepLimiter
+include("courant_step.jl")
 
 end

src/Drivers/courant_step.jl

+struct CFLStepLimiter{TTile,A}
+    tile::TTile
+    Δt::A
+    default_dt::Float64
+end
+function (cfl::CFLStepLimiter{<:Tile,<:AbstractArray})(u,p,t)
+    let Δt = cfl.Δt,
+        Δx = dustrip(Δ(cfl.tile.grid)),
+        Ceff = getvar(:Ceff, cfl.tile, u), # apparent heat capacity
+        kc = getvar(:kc, cfl.tile, u); # thermal cond. at grid centers
+        @. Δt = Utils.adstrip(Δx^2 * Ceff / kc)
+        Δt_min = minimum(Δt)
+        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, cfl.default_dt)
+    end
+end
+function (cfl::CFLStepLimiter{<:Tile,Nothing})(u,p,t)
+    let Δt = cfl.Δt,
+        Δx = dustrip(Δ(cfl.tile.grid)),
+        Ceff = getvar(:Ceff, cfl.tile, u), # apparent heat capacity
+        kc = getvar(:kc, cfl.tile, u); # thermal cond. at grid centers
+        Δt = Utils.adstrip(Δx^2 * Ceff / kc)
+        Δt_min = minimum(Δt)
+        IfElse.ifelse(isfinite(Δt_min) && Δt_min > 0, Δt_min, cfl.default_dt)
+    end
+end
+function CFLStepLimiter(tile::HeatOnlyTile; courant_number=1//2, default_dt=60.0, iip=true)
+    cfl = iip ? CFLStepLimiter(tile, zero(dustrip(Δ(setup.grid))), default_dt) : CFLHeatState(tile, nothing, default_dt)
+    StepsizeLimiter(cfl; safety_factor=courant_number, max_step=true)
+end

src/Land/Land.jl

@@ -32,11 +32,11 @@ export StratComponent, componentname, copmonenttypes, components, boundaries
 export top, bottom, subsurface
 include("stratigraphy.jl")
 
-export LandModelState, LayerState
+export TileState, LayerState
 include("state.jl")
 
 export Tile, withaxes, getstate
-include("model.jl")
+include("tile.jl")
 
 export ParameterVector, LinearTrend, parameters
 include("params.jl")

src/Land/state.jl

@@ -36,32 +36,32 @@ end
 end
 
 """
-    LandModelState{iip,TGrid,TStates,Tt,names}
+    TileState{iip,TGrid,TStates,Tt,names}
 
 Represents the instantaneous state of a CryoGrid `Tile`.
 """
-struct LandModelState{iip,TGrid,TStates,Tt,names}
+struct TileState{iip,TGrid,TStates,Tt,names}
     grid::TGrid
     states::NamedTuple{names,TStates}
     t::Tt
-    LandModelState(grid::TGrid, states::NamedTuple{names,TS}, t::Tt, ::Val{iip}=Val{inplace}()) where
+    TileState(grid::TGrid, states::NamedTuple{names,TS}, t::Tt, ::Val{iip}=Val{inplace}()) where
         {TGrid<:Numerics.AbstractDiscretization,TS<:Tuple{Vararg{<:LayerState}},Tt,names,iip} =
             new{iip,TGrid,TS,Tt,names}(grid, states, t)
 end
-Base.getindex(state::LandModelState, sym::Symbol) = Base.getproperty(state, sym)
-Base.getindex(state::LandModelState, i::Int) = state.states[i]
-function Base.getproperty(state::LandModelState, sym::Symbol)
+Base.getindex(state::TileState, sym::Symbol) = Base.getproperty(state, sym)
+Base.getindex(state::TileState, i::Int) = state.states[i]
+function Base.getproperty(state::TileState, sym::Symbol)
     return if sym ∈ (:grid,:states,:t)
         getfield(state, sym)
     else
         getproperty(getfield(state, :states), sym)
     end
 end
-@inline @generated function LandModelState(vs::VarStates{names}, zs::NTuple, u=copy(vs.uproto), du=similar(vs.uproto), t=0.0, ::Val{iip}=Val{inplace}()) where {names,iip}
+@inline @generated function TileState(vs::VarStates{names}, zs::NTuple, u=copy(vs.uproto), du=similar(vs.uproto), t=0.0, ::Val{iip}=Val{inplace}()) where {names,iip}
     layerstates = (:(LayerState(vs, (ustrip(bounds[$i][1]), ustrip(bounds[$i][2])), u, du, t, Val{$(QuoteNode(names[i]))}(), Val{iip}())) for i in 1:length(names))
     quote
         bounds = boundaryintervals(zs, vs.grid[end])
-        return LandModelState(
+        return TileState(
             vs.grid,
             NamedTuple{tuple($(map(QuoteNode,names)...))}(tuple($(layerstates...))),
             t,

src/Land/model.jl → src/Land/tile.jl

@@ -130,7 +130,7 @@ is only executed during compilation and will not appear in the compiled version.
     _du .= zero(eltype(_du))
     du = ComponentArray(_du, getaxes(model.state.uproto))
     u = ComponentArray(_u, getaxes(model.state.uproto))
-    state = LandModelState(model.state, boundaries(strat), u, du, t, Val{inplace}())
+    state = TileState(model.state, boundaries(strat), u, du, t, Val{inplace}())
     end push!(expr.args)
     # Initialize variables for all layers
     for i in 1:N
@@ -208,7 +208,7 @@ initialcondition!(model::Tile, tspan::NTuple{2,DateTime}, p::AbstractVector, arg
     du = zero(similar(model.state.uproto, eltype(p)))
     u = zero(similar(model.state.uproto, eltype(p)))
     strat = Flatten.reconstruct(model.strat, p, ModelParameters.SELECT, ModelParameters.IGNORE)
-    state = LandModelState(model.state, boundaries(strat), u, du, tspan[1], Val{iip}())
+    state = TileState(model.state, boundaries(strat), u, du, tspan[1], Val{iip}())
     end push!(expr.args)
     # Call initializers
     for i in 1:N
@@ -315,7 +315,7 @@ Gets the
 function getstate(model::Tile{TStrat,TGrid,TStates,iip}, _u, _du, t) where {TStrat,TGrid,TStates,iip}
     du = ComponentArray(_du, getaxes(model.state.uproto))
     u = ComponentArray(_u, getaxes(model.state.uproto))
-    return LandModelState(model.strat, model.state, u, du, t, Val{iip}())
+    return TileState(model.strat, model.state, u, du, t, Val{iip}())
 end
 """
 Collects and validates all declared variables (`Var`s) for the given strat component.

src/Physics/HeatConduction/heat_bc.jl

@@ -12,7 +12,7 @@ struct TemperatureGradient{E,F} <: BoundaryProcess
 end
 BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
 BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
-@inline boundaryvalue(bc::TemperatureGradient,l2,p2,s1,s2) where {F} = bc.T(s1.t)
+@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) where {F} = bc.T(s1.t)
 
 @with_kw struct NFactor{W,S} <: BoundaryEffect
     winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0