diff --git a/src/Callbacks/Callbacks.jl b/src/Callbacks/Callbacks.jl deleted file mode 100644 index 0fe659f8b4e9cd62136f24b69747caa1bad3cd48..0000000000000000000000000000000000000000 --- a/src/Callbacks/Callbacks.jl +++ /dev/null @@ -1,46 +0,0 @@ -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 diff --git a/src/Callbacks/courant_step.jl b/src/Callbacks/courant_step.jl deleted file mode 100644 index 4cd09b778e61b53d6b63c12d8a7ae1ff0412cbde..0000000000000000000000000000000000000000 --- a/src/Callbacks/courant_step.jl +++ /dev/null @@ -1,30 +0,0 @@ -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 diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl index ce2b3234d02d7ce0b685922e97bafff7807741db..2925ebd6c55452d15a3a603528251cbfe0c9254e 100755 --- a/src/CryoGrid.jl +++ b/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 diff --git a/src/Drivers/Drivers.jl b/src/Drivers/Drivers.jl index 3894bd45496922a31a7faa0a45c567cce105af2e..8aa39dbb34f799cd9a7cd9821f10d077a2989c45 100644 --- a/src/Drivers/Drivers.jl +++ b/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 diff --git a/src/Drivers/courant_step.jl b/src/Drivers/courant_step.jl new file mode 100644 index 0000000000000000000000000000000000000000..f9e609a7edc46bd1a1634c6bcc30a6ee2733ead6 --- /dev/null +++ b/src/Drivers/courant_step.jl @@ -0,0 +1,29 @@ +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