Skip to content
Snippets Groups Projects
Commit 97473876 authored by Brian Groenke's avatar Brian Groenke
Browse files

Move CFL to Drivers module and remove Callbacks module

parent 17d251e9
No related branches found
No related tags found
1 merge request!56Fix bug in temperature gradient
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
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
...@@ -28,7 +28,6 @@ include("Land/Land.jl") ...@@ -28,7 +28,6 @@ include("Land/Land.jl")
include("IO/InputOutput.jl") include("IO/InputOutput.jl")
include("Diagnostics/Diagnostics.jl") include("Diagnostics/Diagnostics.jl")
include("Drivers/Drivers.jl") include("Drivers/Drivers.jl")
include("Callbacks/Callbacks.jl")
using .Utils using .Utils
# Re-exported submodules # Re-exported submodules
...@@ -40,7 +39,6 @@ using .Utils ...@@ -40,7 +39,6 @@ using .Utils
@reexport using .Drivers @reexport using .Drivers
@reexport using .InputOutput @reexport using .InputOutput
@reexport using .Diagnostics @reexport using .Diagnostics
@reexport using .Callbacks
# Re-exported packages # Re-exported packages
@reexport using Dates: Dates, DateTime @reexport using Dates: Dates, DateTime
......
...@@ -20,7 +20,7 @@ using Unitful ...@@ -20,7 +20,7 @@ using Unitful
import CryoGrid.Land: Tile, Stratigraphy, StratComponent import CryoGrid.Land: Tile, Stratigraphy, StratComponent
export DefaultJac, TridiagJac, JacobianStyle, HeatOnlyLandModel export DefaultJac, TridiagJac, JacobianStyle, HeatOnlyTile
""" """
JacobianStyle JacobianStyle
...@@ -39,11 +39,13 @@ JacobianStyle(::Type{<:Tile}) = DefaultJac() ...@@ -39,11 +39,13 @@ JacobianStyle(::Type{<:Tile}) = DefaultJac()
# Auto-detect Jacobian sparsity for problems with one or more heat-only layers. # 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! # 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 :( # 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} const HeatOnlyTile = Tile{<:Stratigraphy{N,<:Tuple{TTop,Vararg{<:Union{<:StratComponent{<:SubSurface, <:CompoundProcess{<:Tuple{<:Heat}}},TBot}}}}} where {N,TTop,TBot}
JacobianStyle(::Type{<:HeatOnlyLandModel}) = TridiagJac() JacobianStyle(::Type{<:HeatOnlyTile}) = TridiagJac()
# DiffEq/SciML driver # DiffEq/SciML driver
export CryoGridProblem export CryoGridProblem
include("diffeq.jl") include("diffeq.jl")
export CFLStepLimiter
include("courant_step.jl")
end end
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment