diff --git a/Project.toml b/Project.toml
index 28478b6dfd55b69be12556064e1491b165630c67..f8d0de1bd7a41aa93d3400de913ca0c99dff3ab9 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
name = "CryoGrid"
uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
authors = ["Brian Groenke <brian.groenke@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.6.5"
+version = "0.6.6"
[deps]
ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl
index 62a31e607aafe6f5b2d1f3f08daa8c4e427cd290..44033241e417d22189e2c36734937f34f4fdcead 100755
--- a/src/CryoGrid.jl
+++ b/src/CryoGrid.jl
@@ -12,8 +12,6 @@ using Reexport
# Common types and methods
export Layer, SubSurface, Top, Bottom
export Process, SubSurfaceProcess, BoundaryProcess, CoupledProcesses, Coupled
-export BoundaryStyle, Dirichlet, Neumann
-export Callback, CallbackStyle
include("types.jl")
export variables, initialcondition!, diagnosticstep!, prognosticstep!, interact!
export boundaryflux, boundaryvalue, criterion, affect!, observe
@@ -28,10 +26,11 @@ include("IO/InputOutput.jl")
include("Diagnostics/Diagnostics.jl")
include("Drivers/Drivers.jl")
+using .Numerics
+export Grid, cells, edges, subgridinds, Δ, volume, area, initializer
using .Utils
+export convert_t, convert_tspan, @sym_str
# Re-exported submodules
-@reexport using .Numerics
-@reexport using .Utils: convert_t, convert_tspan, @sym_str
@reexport using .Physics
@reexport using .Strat
@reexport using .Drivers
diff --git a/src/Numerics/Numerics.jl b/src/Numerics/Numerics.jl
index 8c65d2a399a8d4d785da18f8c72601f1f836a10f..76cc033a307c263374a5b9e89a77613f258bfb63 100644
--- a/src/Numerics/Numerics.jl
+++ b/src/Numerics/Numerics.jl
@@ -38,7 +38,7 @@ struct UnitVolume <: Geometry end
export ∇, Tabulated
include("math.jl")
-export Grid, cells, edges, indexmap, subgridinds, Δ, volume, area
+export Grid, cells, edges, subgridinds, Δ, volume, area
include("grid.jl")
export Profile, profile2array, array2profile
diff --git a/src/Numerics/grid.jl b/src/Numerics/grid.jl
index 025f1b43f0bbf497cf5a0ba89549f6548afba21e..3e3392b154d0cc74189a3c719c8a71dca122b649 100644
--- a/src/Numerics/grid.jl
+++ b/src/Numerics/grid.jl
@@ -53,7 +53,7 @@ ConstructionBase.constructorof(::Type{Grid{S,G,Q,A}}) where {S,G,Q,A} = (geom,va
Base.show(io::IO, grid::Grid{S,G}) where {S,G} = print(io, "Grid{$S}($(grid[1])..$(grid[end])) of length $(length(grid)) with geometry $G")
Base.show(io::IO, ::MIME{Symbol("text/plain")}, grid::Grid) = show(io, grid)
-function subgridinds(grid::Grid{S,G,Q,A}, interval::Interval{L,R,Q}) where {S,G,Q,A,L,R}
+function subgridinds(grid::Grid, interval::Interval{L,R}) where {L,R}
@assert interval.left <= interval.right "Invalid interval: $interval"
vals = values(grid)
# Determine indices which lie in the given interval
diff --git a/src/Numerics/varstates.jl b/src/Numerics/varstates.jl
index 260a22df164065a4379c3ea239a5e7748b7ac016..6ef0c28c8c08a9e910e809ee6bac026ad5559640 100644
--- a/src/Numerics/varstates.jl
+++ b/src/Numerics/varstates.jl
@@ -12,8 +12,8 @@ struct VarStates{names,griddvars,TU,TD,TV,DF,DG}
griddiag::NamedTuple{griddvars,DG} # on-grid non-prognostic variables
gridcache::Dict{ClosedInterval{Int},TD} # grid cache; indices -> subgrid
end
-@generated function getvar(::Val{name}, vs::VarStates{layers,griddvars,TU,TD,TV}, u::TU, du::Union{Nothing,TU}=nothing) where
- {name,layers,griddvars,T,A,pax,TU<:ComponentVector{T,A,Tuple{Axis{pax}}},TD,TV}
+@generated function getvar(::Val{name}, vs::VarStates{layers,griddvars}, u, du=nothing) where {name,layers,griddvars}
+ pax = ComponentArrays.indexmap(first(ComponentArrays.getaxes(u)))
dnames = map(n -> Symbol(:d,n), keys(pax))
if name ∈ griddvars
quote
diff --git a/src/Physics/Boundaries/Boundaries.jl b/src/Physics/Boundaries/Boundaries.jl
index 9e2e23091759d6bd053009eecdfc4e1ec079a2fb..e2e0a49cfdfda9097969005a050675bcb756c17b 100644
--- a/src/Physics/Boundaries/Boundaries.jl
+++ b/src/Physics/Boundaries/Boundaries.jl
@@ -18,74 +18,11 @@ using Unitful
import Flatten: flattenable
-export Constant, Periodic, Bias
export BoundaryEffect, Damping
include("effects.jl")
export Forcing, TimeSeriesForcing, ForcingData
include("forcing.jl")
-
-"""
- struct Constant{S,T} <: BoundaryProcess
-
-Constant boundary condition (of any type/unit) specified by `value`.
-"""
-struct Constant{S,T} <: BoundaryProcess
- value::T
- Constant(::Type{S}, value::T) where {S<:BoundaryStyle,T} = new{S,T}(value)
-end
-ConstructionBase.constructorof(::Type{<:Constant{S}}) where {S} = value -> Constant(S,value)
-boundaryvalue(bc::Constant{S,T},l1,p2,l2,s1,s2) where {S,T} = bc.value
-
-BoundaryStyle(::Type{<:Constant{S}}) where {S} = S()
-
-"""
- struct Periodic{S,T} <: BoundaryProcess
-
-Periodic boundary condition (of any type/unit) specified by `period`, `amplitude`, and `phaseshift`.
-"""
-struct Periodic{S,T} <: BoundaryProcess
- period::Float"s"
- amplitude::T
- phaseshift::T
- Periodic{S}(period::Q, amplitude::T=one(T), phaseshift::T=one(T)) where
- {S<:BoundaryStyle,Q<:TimeQuantity,T} =
- new{S,T}(uconvert(u"s",period) |> dustrip, amplitude, phaseshift)
-end
-
-@inline boundaryvalue(bc::Periodic,l1,p2,l2,s1,s2) = bc.amplitude*sin(Ï€*(1/bc.period)*t + bc.phaseshift)
-
-BoundaryStyle(::Type{<:Periodic{S}}) where {S} = S()
-
-@with_kw struct Bias{P} <: BoundaryProcess
- bias::P = Param(0.0)
-end
-@inline boundaryvalue(bc::Bias,l1,p2,l2,s1,s2) = bc.bias
-
-BoundaryStyle(::Type{<:Bias}) = Dirichlet()
-
-"""
- struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
-
-Represents a composition of two boundary processes, `B1` and `B2`, via an operator `F`.
-A typical use case is combining `Constant` with a forcing-driven boundary process to
-scale or shift the forcing.
-"""
-struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
- op::F
- bc1::B1
- bc2::B2
- function CombinedBoundaryProcess(op::F, bc1::B1, bc2::B2) where {F,B1<:BoundaryProcess,B2<:BoundaryProcess}
- @assert BoundaryStyle(bc1) == BoundaryStyle(bc2) "boundary condition styles (e.g. Dirichlet vs Neumann) must match"
- new{B1,B2,F,typeof(BoundaryStyle(bc1))}(op,bc1,bc2)
- end
-end
-@inline boundaryvalue(cbc::CombinedBoundaryProcess,l1,p2,l2,s1,s2) = cbc.op(cbc.bc1(l1,l2,p2,s1,s2), cbc.bc2(l1,l2,p2,s1,s2))
-variables(top::Top, cbc::CombinedBoundaryProcess) = tuplejoin(variables(top, cbc.bc1), variables(top, cbc.bc2))
-BoundaryStyle(::Type{CombinedBoundaryProcess{B1,B2,F,S}}) where {F,B1,B2,S} = S()
-# Overload arithmetic operators on boundary processes.
-Base.:+(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(+, bc1, bc2)
-Base.:-(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(-, bc1, bc2)
-Base.:*(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(*, bc1, bc2)
-Base.:/(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(/, bc1, bc2)
+export ConstantBC, PeriodicBC, Bias
+include("bc.jl")
end
diff --git a/src/Physics/Boundaries/bc.jl b/src/Physics/Boundaries/bc.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f1d8a3a1db64635ba92c597e23dfc947608a75b4
--- /dev/null
+++ b/src/Physics/Boundaries/bc.jl
@@ -0,0 +1,63 @@
+"""
+ ConstantBC{S,T} <: BoundaryProcess
+
+Constant boundary condition (of any type/unit) specified by `value`.
+"""
+struct ConstantBC{S,T} <: BoundaryProcess
+ value::T
+ ConstantBC(::Type{S}, value::T) where {S<:BoundaryStyle,T} = new{S,T}(value)
+end
+ConstructionBase.constructorof(::Type{<:ConstantBC{S}}) where {S} = value -> ConstantBC(S,value)
+boundaryvalue(bc::ConstantBC{S,T},l1,p2,l2,s1,s2) where {S,T} = bc.value
+
+BoundaryStyle(::Type{<:ConstantBC{S}}) where {S} = S()
+
+"""
+ PeriodicBC{S,T} <: BoundaryProcess
+
+Periodic boundary condition (of any type/unit) specified by `period`, `amplitude`, and `phaseshift`.
+"""
+struct PeriodicBC{S,T} <: BoundaryProcess
+ period::Float"s"
+ amplitude::T
+ phaseshift::T
+ PeriodicBC(::Type{S}, period::Q, amplitude::T=one(T), phaseshift::T=one(T)) where
+ {S<:BoundaryStyle,Q<:TimeQuantity,T} =
+ new{S,T}(uconvert(u"s",period) |> dustrip, amplitude, phaseshift)
+end
+
+@inline boundaryvalue(bc::PeriodicBC,l1,p2,l2,s1,s2) = bc.amplitude*sin(Ï€*(1/bc.period)*t + bc.phaseshift)
+
+BoundaryStyle(::Type{<:PeriodicBC{S}}) where {S} = S()
+
+@with_kw struct Bias{P} <: BoundaryProcess
+ bias::P = Param(0.0)
+end
+@inline boundaryvalue(bc::Bias,l1,p2,l2,s1,s2) = bc.bias
+
+BoundaryStyle(::Type{<:Bias}) = Dirichlet()
+
+"""
+ struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
+
+Represents a composition of two boundary processes, `B1` and `B2`, via an operator `F`.
+A typical use case is combining `ConstantBC` with a forcing-driven boundary process to
+scale or shift the forcing.
+"""
+struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
+ op::F
+ bc1::B1
+ bc2::B2
+ function CombinedBoundaryProcess(op::F, bc1::B1, bc2::B2) where {F,B1<:BoundaryProcess,B2<:BoundaryProcess}
+ @assert BoundaryStyle(bc1) == BoundaryStyle(bc2) "boundary condition styles (e.g. Dirichlet vs Neumann) must match"
+ new{B1,B2,F,typeof(BoundaryStyle(bc1))}(op,bc1,bc2)
+ end
+end
+@inline boundaryvalue(cbc::CombinedBoundaryProcess,l1,p2,l2,s1,s2) = cbc.op(cbc.bc1(l1,l2,p2,s1,s2), cbc.bc2(l1,l2,p2,s1,s2))
+variables(top::Top, cbc::CombinedBoundaryProcess) = tuplejoin(variables(top, cbc.bc1), variables(top, cbc.bc2))
+BoundaryStyle(::Type{CombinedBoundaryProcess{B1,B2,F,S}}) where {F,B1,B2,S} = S()
+# Overload arithmetic operators on boundary processes.
+Base.:+(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(+, bc1, bc2)
+Base.:-(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(-, bc1, bc2)
+Base.:*(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(*, bc1, bc2)
+Base.:/(bc1::BoundaryProcess, bc2::BoundaryProcess) = CombinedBoundaryProcess(/, bc1, bc2)
diff --git a/src/Physics/HeatConduction/heat_bc.jl b/src/Physics/HeatConduction/heat_bc.jl
index 4b793ee31eb64b76d02ba55a3a6c9faa4089bd7a..9483e3f456b8ee4783cc456ce0de2adf3821b9f4 100644
--- a/src/Physics/HeatConduction/heat_bc.jl
+++ b/src/Physics/HeatConduction/heat_bc.jl
@@ -1,9 +1,9 @@
# Boundary condition type aliases
-const ConstantTemp = Constant{Dirichlet,Float"°C"}
-ConstantTemp(value::UFloat"K") = Constant(Dirichlet, dustrip(u"°C", value))
-ConstantTemp(value::UFloat"°C") = Constant(Dirichlet, dustrip(value))
-const GeothermalHeatFlux = Constant{Neumann,Float"J/s/m^2"}
-GeothermalHeatFlux(value::UFloat"J/s/m^2"=0.053xu"J/s/m^2") = Constant(Neumann, dustrip(value))
+const ConstantTemp = ConstantBC{Dirichlet,Float"°C"}
+ConstantTemp(value::UFloat"K") = ConstantBC(Dirichlet, dustrip(u"°C", value))
+ConstantTemp(value::UFloat"°C") = ConstantBC(Dirichlet, dustrip(value))
+const GeothermalHeatFlux = ConstantBC{Neumann,Float"J/s/m^2"}
+GeothermalHeatFlux(value::UFloat"J/s/m^2"=0.053xu"J/s/m^2") = ConstantBC(Neumann, dustrip(value))
struct TemperatureGradient{E,F} <: BoundaryProcess
T::F # temperature forcing
diff --git a/src/Physics/Soils/sfcc_solvers.jl b/src/Physics/Soils/sfcc_solvers.jl
index ac5bc5da88aead824ae2b5f6c1986e6c6f98e326..419da37548f32c133104fb2f28d43a56734f32fa 100644
--- a/src/Physics/Soils/sfcc_solvers.jl
+++ b/src/Physics/Soils/sfcc_solvers.jl
@@ -165,15 +165,17 @@ function initialcondition!(soil::Soil, heat::Heat, sfcc::SFCC{F,∇F,<:SFCCPreSo
Hmax = enthalpy(Tmax, C(Tmax), L, θ(Tmax)),
dH = sfcc.solver.dH,
Hs = Hmin:dH:Hmax;
- θs = [θres]
- Ts = [Tmin]
- for _ in Hs[2:end]
- θᵢ = θs[end]
- Táµ¢ = Ts[end]
+ θs = Vector{eltype(state.θl)}(undef, length(Hs))
+ θs[1] = θres
+ Ts = Vector{eltype(state.T)}(undef, length(Hs))
+ Ts[1] = Tmin
+ for i in 2:length(Hs)
+ θᵢ = θs[i-1]
+ Táµ¢ = Ts[i-1]
dTdH = 1.0 / (C(θᵢ) + dθdT(Tᵢ)*(Tᵢ*cw + L))
dθdH = 1.0 / (1.0/dθdT(Tᵢ)*C(θᵢ)+Tᵢ*cw + L)
- push!(θs, θᵢ + dH*dθdH)
- push!(Ts, Táµ¢ + dH*dTdH)
+ θs[i] = θᵢ + dH*dθdH
+ Ts[i] = Táµ¢ + dH*dTdH
end
sfcc.solver.cache.f = Interpolations.extrapolate(
Interpolations.interpolate((Vector(Hs),), θs, Interpolations.Gridded(Interpolations.Linear())),
diff --git a/src/Presets/Presets.jl b/src/Presets/Presets.jl
index cea6aa8573cd3bb7709f6b336a68b3b26814ec4c..62b27e8289d182fa5618ae47545097d73f75fe38 100644
--- a/src/Presets/Presets.jl
+++ b/src/Presets/Presets.jl
@@ -31,9 +31,9 @@ end
SoilHeatColumn(upperbc::BoundaryProcess, soilprofile::Profile{N,D,<:SoilParameterization}; grid::Grid=DefaultGrid_2cm, freezecurve::F=FreeWater()) where {N,D,F<:FreezeCurve} = SoilHeatColumn(:H, upperbc, soilprofile; grid=grid, freezecurve=freezecurve)
Forcings = (
- Samoylov_ERA5_fitted_daily_1979_2020 = Resource("samoylov_era5_fitted_daily_1979-2020", "json", "https://nextcloud.awi.de/s/Jk274jKGegkkTie/download/samoylov_era5_fitted_daily_1979-2020.json"),
- Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", "json", "https://nextcloud.awi.de/s/F98s5WEo9xMPod7/download"),
- Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", "json", "https://nextcloud.awi.de/s/RSqBtp5sPwkCf45/download"),
+ Samoylov_ERA5_fitted_daily_1979_2020 = Resource("samoylov_era5_fitted_daily_1979-2020", "json", "https://nextcloud.awi.de/s/DSx3BWsfjCdy9MF"),
+ Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044 = Resource("Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044", "json", "https://nextcloud.awi.de/s/dp74KC2ceQKaG43"),
+ Samoylov_ERA_MkL3_CCSM4_long_term = Resource("Samoylov_ERA_MkL3_CCSM4_long_term", "json", "https://nextcloud.awi.de/s/gyoMTy9jpk2pMxL"),
)
Parameters = (
# Faroux et al. doi:10.1109/IGARSS.2007.4422971
diff --git a/src/Strat/state.jl b/src/Strat/state.jl
index d93c2b76ba1e8ee86028dfd0b74ec51f5f5f2f89..3216fe2c8186aa9cf4476937e8eff26e38f77f7b 100644
--- a/src/Strat/state.jl
+++ b/src/Strat/state.jl
@@ -24,7 +24,7 @@ function Base.getproperty(state::LayerState, sym::Symbol)
getproperty(getfield(state, :states), sym)
end
end
-@inline function LayerState(vs::VarStates, zs::NTuple{2,Tz}, u, du, t, ::Val{layername}, ::Val{iip}=Val{inplace}()) where {Tz,layername,iip}
+@inline function LayerState(vs::VarStates, zs::NTuple{2}, u, du, t, ::Val{layername}, ::Val{iip}=Val{inplace}()) where {layername,iip}
z_inds = subgridinds(edges(vs.grid), zs[1]..zs[2])
return LayerState(
_makegrids(Val{layername}(), getproperty(vs.vars, layername), vs, z_inds),
@@ -58,9 +58,15 @@ function Base.getproperty(state::TileState, sym::Symbol)
end
end
@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))
+ layerstates = (
+ quote
+ bounds_i = (ustrip(bounds[$i][1]), ustrip(bounds[$i][2]))
+ LayerState(vs, bounds_i, u, du, t, Val{$(QuoteNode(names[i]))}(), Val{iip}())
+ end
+ for i in 1:length(names)
+ )
quote
- bounds = boundarypairs(zs, vs.grid[end])
+ bounds = boundarypairs(zs, convert(eltype(zs), vs.grid[end]))
return TileState(
vs.grid,
NamedTuple{tuple($(map(QuoteNode,names)...))}(tuple($(layerstates...))),
diff --git a/src/Utils/Utils.jl b/src/Utils/Utils.jl
index 6d0c172f7912a640b47dba63e5b40fc3ea52dc08..514ca0d5835407fae00677a7acb0a30c61bedea0 100644
--- a/src/Utils/Utils.jl
+++ b/src/Utils/Utils.jl
@@ -28,7 +28,7 @@ const DistQuantity{T,U} = Quantity{T,Unitful.ð‹,U} where {T,U<:DistUnit}
const TempUnit{N,A} = Unitful.FreeUnits{N,Unitful.ðš¯,A} where {N,A}
const TempQuantity{T,U} = Quantity{T,Unitful.ðš¯,U} where {T,U<:TempUnit}
const TimeUnit{N,A} = Unitful.FreeUnits{N,Unitful.ð“,A} where {N,A}
-const TimeQuantity{T,U} = Quantity{T,Unitful.ð“,U} where {T,U<:TempUnit}
+const TimeQuantity{T,U} = Quantity{T,Unitful.ð“,U} where {T,U<:TimeUnit}
StructTypes.StructType(::Type{<:Quantity}) = StructTypes.CustomStruct()
StructTypes.lower(value::Quantity) = string(value)
diff --git a/test/Physics/HeatConduction/heat_conduction_tests.jl b/test/Physics/HeatConduction/heat_conduction_tests.jl
index 91b11c129b81f8e483ea6f88bc9567607873da3f..f697daa9b68c9fe1b4900e2618d1c565f3273699 100644
--- a/test/Physics/HeatConduction/heat_conduction_tests.jl
+++ b/test/Physics/HeatConduction/heat_conduction_tests.jl
@@ -32,7 +32,7 @@ include("../../types.jl")
Δk = Δ(x)
sub = TestGroundLayer()
heat = Heat()
- bc = Constant(Dirichlet, 0.0u"°C")
+ bc = ConstantBC(CryoGrid.Dirichlet, 0.0u"°C")
@testset "top: +, bot: -" begin
T₀ = Vector(LinRange(-23,27,length(xc)))u"°C"
∂H = zeros(length(T₀))u"J/s/m^3"
@@ -62,7 +62,7 @@ include("../../types.jl")
@test ∂H[end] < 0.0u"J/s/m^3"
end
@testset "Neumann boundary" begin
- bc = Constant(Neumann, -1.0u"J/m^3")
+ bc = ConstantBC(CryoGrid.Neumann, -1.0u"J/m^3")
T₀ = Vector(LinRange(-23,27,length(xc)))u"°C"
∂H = zeros(length(T₀))u"J/s/m^3"
state = (T=T₀,k=k,dH=∂H,grids=(T=xc,k=x),t=0.0)
@@ -93,7 +93,7 @@ end
f_analytic(x,t) = exp(-t*4Ï€^2)*sin(2.0*Ï€*x)
sub = TestGroundLayer()
heat = Heat()
- bc = Constant(Dirichlet, 0.0u"°C")
+ bc = ConstantBC(CryoGrid.Dirichlet, 0.0u"°C")
function dTdt(T,p,t)
dT = similar(T)u"J/s/m^3"
dT .= zero(eltype(dT))