Merge branch 'refactor/layer-naming' into 'master'

Rework layer name handling in stratigraphy

See merge request cryogrid/cryogridjulia!115

Project.toml

@@ -1,7 +1,7 @@
 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <[email protected]>", "Jan Nitzbon <[email protected]>", "Moritz Langer <[email protected]>"]
-version = "0.19.2"
+version = "0.19.3"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

docs/make.jl

@@ -36,7 +36,8 @@ rm(examples_output_dir, recursive=true, force=true)
 # recreate directory
 mkpath(examples_output_dir)
 # generate example docs from scripts
-example_docfiles = map(filter(∉(["Manifest.toml", "Project.toml"]), readdir(examples_dir))) do f
+ignored_files = ["Manifest.toml", "Project.toml", "08_heat_sfcc_richardseq_samoylov.jl"]
+example_docfiles = map(filter(∉(ignored_files), readdir(examples_dir))) do f
        infile = joinpath(examples_dir, f)
        @info "Generating docpage for example script $infile and writing to directory $examples_output_dir"
        Literate.markdown(infile, examples_output_dir, execute=false, documenter=true)

docs/src/manual/overview.md

@@ -55,7 +55,7 @@ julia> out.T
 Notice that, in the example above, it is types such as `Ground`, `HeatBalance`, `DallAmico`, etc. that specify which components the model should use. These components are defined by adding method dispatches to the [CryoGrid interface](@ref toplevel) methods. State variables are declared via the [`variables`](@ref) method, e.g:
 
 ```julia
-variables(soil::Soil, heat::HeatBalance{<:Enthalpy}) = (
+variables(soil::Soil, heat::HeatBalance{<:EnthalpyBased}) = (
     Prognostic(:H, OnGrid(Cells), u"J/m^3"),
     Diagnostic(:T, OnGrid(Cells), u"°C"),
     Diagnostic(:C, OnGrid(Cells), u"J//K*/m^3"),
@@ -74,7 +74,7 @@ The real work finally happens in [`updatestate!`](@ref) and [`computefluxes!`](@
 We can take as an example the implementation of `computefluxes!` for enthalpy-based heat conduction (note that `jH` is a diagnostic variable representing the energy flux over each cell edge):
 
 ```julia
-function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:Enthalpy}, state)
+function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:EnthalpyBased}, state)
     Δk = Δ(state.grid) # cell sizes
     ΔT = Δ(cells(state.grid)) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set

examples/02_heat_sfcc_constantbc.jl

@@ -30,7 +30,7 @@ Plots.plot(-2.0u"°C":0.01u"K":0.0u"°C", sfcc)
 # Enthalpy form of the heat transfer operator (i.e. prognostic :H). In this case, this is equivalent to
 # the shorthand `SoilHeatTile(:H, ...)`. However, it's worth demonstrating how the operator can be explicitly
 # specified.
-heatop = Heat.EnthalpyForm(SFCCPreSolver())
+heatop = Heat.MOLEnthalpy(SFCCPreSolver())
 initT = initializer(:T, tempprofile)
 tile = CryoGrid.Presets.SoilHeatTile(
     heatop,

examples/03_heat_sfcc_samoylov.jl

@@ -4,7 +4,7 @@
 # forced using air temperatures from Samoylov Island. We use the
 # SFCC formulation of Painter and Karra (2014). For the purpose
 # of demonstration, we use the apparent heat capacity form of the
-# heat equation in this example (i.e. [`Heat.TemperatureForm`](@ref)).
+# heat equation in this example (i.e. [`Heat.MOLTemperature`](@ref)).
 
 using CryoGrid
 

examples/04_heat_freeW_snow_samoylov.jl

@@ -29,11 +29,11 @@ snowpack = Snowpack(
 strat = @Stratigraphy(
     z_top => Top(upperbc),
     z_top => :snowpack => snowpack,
-    z_sub[1] => :topsoil1 => Ground(soilprofile[1].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
-    z_sub[2] => :topsoil2 => Ground(soilprofile[2].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
-    z_sub[3] => :sediment1 => Ground(soilprofile[3].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
-    z_sub[4] => :sediment2 => Ground(soilprofile[4].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
-    z_sub[5] => :sediment3 => Ground(soilprofile[5].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    z_sub[1] => Ground(soilprofile[1].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    z_sub[2] => Ground(soilprofile[2].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    z_sub[3] => Ground(soilprofile[3].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    z_sub[4] => Ground(soilprofile[4].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    z_sub[5] => Ground(soilprofile[5].value, heat=HeatBalance(), water=WaterBalance(BucketScheme())),
     z_bot => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
 );
 modelgrid = CryoGrid.Presets.DefaultGrid_5cm

examples/05_heat_freeW_bucketW_samoylov.jl

@@ -9,15 +9,15 @@ forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2
 _, tempprofile = CryoGrid.Presets.SamoylovDefault;
 initT = initializer(:T, tempprofile)
 initsat = initializer(:sat, (l,state) -> state.sat .= l.para.sat);
-upperbc = WaterHeatBC(SurfaceWaterBalance(rainfall=forcings.rainfall), TemperatureGradient(forcings.Tair, NFactor(0.5,0.9)));
+upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureGradient(forcings.Tair, NFactor(0.5,0.9)));
 
 # The `@Stratigraphy` macro is just a small convenience that automatically wraps the three subsurface layers in a tuple.
 # It would be equivalent to use the `Stratigraphy` constructor directly and wrap these layers in a tuple or list.
 strat = @Stratigraphy(
     0.0u"m" => Top(upperbc),
-    0.0u"m" => :topsoil => Ground(MineralOrganic(por=0.80,sat=0.5,org=0.75), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
-    0.2u"m" => :subsoil => Ground(MineralOrganic(por=0.40,sat=0.75,org=0.10), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
-    2.0u"m" => :substrat => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.5,org=0.75), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    0.2u"m" => :middle => Ground(MineralOrganic(por=0.40,sat=0.75,org=0.10), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
+    2.0u"m" => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")),
 );
 modelgrid = CryoGrid.Presets.DefaultGrid_2cm

examples/06_heat_freeW_seb_snow_bucketW_samoylov.jl

@@ -30,12 +30,12 @@ water = WaterBalance(BucketScheme(), DampedET())
 ## build stratigraphy
 strat = @Stratigraphy(
     -z => Top(upperbc), 
-    -z => :snowpack => Snowpack(heat=HeatBalance(), water=water),
-    soilprofile[1].depth => :soil1 => Ground(soilprofile[1].value; heat, water),
-    soilprofile[2].depth => :soil2 => Ground(soilprofile[2].value; heat, water),
-    soilprofile[3].depth => :soil3 => Ground(soilprofile[3].value; heat, water),
-    soilprofile[4].depth => :soil4 => Ground(soilprofile[4].value; heat, water),
-    soilprofile[5].depth => :soil5 => Ground(soilprofile[5].value; heat, water),
+    -z => Snowpack(heat=HeatBalance(), water=water),
+    soilprofile[1].depth => Ground(soilprofile[1].value; heat, water),
+    soilprofile[2].depth => Ground(soilprofile[2].value; heat, water),
+    soilprofile[3].depth => Ground(soilprofile[3].value; heat, water),
+    soilprofile[4].depth => Ground(soilprofile[4].value; heat, water),
+    soilprofile[5].depth => Ground(soilprofile[5].value; heat, water),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")),
 );
 ## create Tile
@@ -72,8 +72,17 @@ Plots.plot(ustrip.(out.T[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]',
 # Saturation:
 Plots.plot(ustrip.(out.sat[Z(Near(zs))]), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Soil saturation", leg=false, size=(800,500), dpi=150)
 
+# Runoff
+Plots.plot(ustrip.(out.top.runoff), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Total runoff", leg=false, size=(800,500), dpi=150)
+
+# Evapotranspiration
+Plots.plot(ustrip.(out.top.ET), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Total ET", leg=false, size=(800,500), dpi=150)
+
 # Snow depth:
 Plots.plot(ustrip.(out.snowpack.dsn), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Snow depth", leg=false, size=(800,500), dpi=150)
 
 # Integrated ground heat flux:
 Plots.plot(ustrip.(cumsum(out.top.Qg, dims=2)), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Integrated ground heat flux", leg=false, size=(800,500), dpi=150)
+
+# Integratoed ground latent heat flux:
+Plots.plot(ustrip.(cumsum(out.top.Qe, dims=2)), color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Integrated ground heat flux", leg=false, size=(800,500), dpi=150)

examples/07_heat_freeW_lite_implicit.jl

@@ -23,23 +23,23 @@ freezecurve = FreeWater()
 # freezecurve = PainterKarra(swrc=VanGenuchten(α=0.1, n=1.8))
 strat = @Stratigraphy(
     z_top => Top(upperbc),
-    0.0u"m" => :topsoil1 => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.75), heat=HeatBalance(heatop; freezecurve)),
-    0.1u"m" => :topsoil2 => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)),
-    0.4u"m" => :sediment1 => Ground(MineralOrganic(por=0.55,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)),
-    3.0u"m" => :sediment2 => Ground(MineralOrganic(por=0.50,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)),
-    10.0u"m" => :sediment3 => Ground(MineralOrganic(por=0.30,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)),
+    0.0u"m" => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.75), heat=HeatBalance(heatop; freezecurve)),
+    0.1u"m" => Ground(MineralOrganic(por=0.80,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)),
+    0.4u"m" => Ground(MineralOrganic(por=0.55,sat=1.0,org=0.25), heat=HeatBalance(heatop; freezecurve)),
+    3.0u"m" => Ground(MineralOrganic(por=0.50,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)),
+    10.0u"m" => Ground(MineralOrganic(por=0.30,sat=1.0,org=0.0), heat=HeatBalance(heatop; freezecurve)),
     z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
 modelgrid = CryoGrid.Presets.DefaultGrid_2cm
 tile = Tile(strat, modelgrid, initT);
 
 # Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
 # Here we specify a time span of 5 years.
-tspan = (DateTime(2010,12,30), DateTime(2011,12,30))
+tspan = (DateTime(2010,1,1), DateTime(2015,1,1))
 tspan_sol = convert_tspan(tspan)
 u0, du0 = initialcondition!(tile, tspan);
 prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:T,), step_limiter=nothing)
-sol = @time solve(prob, LiteImplicitEuler(), dt=3*3600.0)
+sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600.0)
 out = CryoGridOutput(sol)
 
 # Plot the results!
@@ -51,6 +51,7 @@ Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="
 # CryoGridLite can also be embedded into integrators from OrdinaryDiffEq.jl via the `NLCGLite` nonlinear solver interface.
 # Note that these sovers generally will not be faster (in execution time) but may be more stable in some cases. Adaptive timestepping can be employed by
 # removing the `adaptive=false` argument.
+using OrdinaryDiffEq
 sol2 = @time solve(prob, ImplicitEuler(nlsolve=NLCGLite(max_iter=1000)), adaptive=false, dt=24*3600.0)
 ## 3rd order additive scheme from Kennedy and Alan 2001
 sol3 = @time solve(prob, KenCarp3(nlsolve=NLCGLite(max_iter=1000)), adaptive=false, dt=24*3600.0)

examples/08_heat_sfcc_richardseq_samoylov.jl

@@ -19,31 +19,32 @@ initT = initializer(:T, tempprofile);
 # Here we conigure the water retention curve and freeze curve. The van Genuchten parameters coorespond to that
 # which would be reasonable for a silty soil.
 swrc = VanGenuchten(α=0.1, n=1.8)
-sfcc = PainterKarra(ω=0.0, swrc=swrc)
+sfcc = PainterKarra(ω=0.0; swrc)
 waterflow = RichardsEq(swrc=swrc);
 
 # We use the enthalpy-based heat diffusion with high accuracy Newton-based solver for inverse enthalpy mapping
-heatop = Heat.EnthalpyForm(SFCCNewtonSolver())
-upperbc = WaterHeatBC(SurfaceWaterBalance(rainfall=forcings.rainfall), TemperatureGradient(forcings.Tair, NFactor(nf=0.6, nt=0.9)));
+heatop = Heat.MOLEnthalpy(SFCCNewtonSolver())
+upperbc = WaterHeatBC(SurfaceWaterBalance(forcings), TemperatureGradient(forcings.Tair, NFactor(nf=0.6, nt=0.9)));
 
 # We will use a simple stratigraphy with three subsurface soil layers.
 # Note that the @Stratigraphy macro lets us list multiple subsurface layers without wrapping them in a tuple.
+heat = HeatBalance(heatop, freezecurve=sfcc, advection=false)
+water = WaterBalance(RichardsEq(; swrc))
 strat = @Stratigraphy(
     -2.0u"m" => Top(upperbc),
-    0.0u"m" => :topsoil => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))),
-    0.2u"m" => :subsoil => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))),
-    2.0u"m" => :substrat => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(heatop, freezecurve=sfcc), water=WaterBalance(RichardsEq(;swrc))),
+    0.0u"m" => Ground(MineralOrganic(por=0.80,sat=0.7,org=0.75); heat, water),
+    0.2u"m" => Ground(MineralOrganic(por=0.40,sat=0.8,org=0.10); heat, water),
+    2.0u"m" => Ground(MineralOrganic(por=0.10,sat=1.0,org=0.0); heat, water),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
 grid = CryoGrid.Presets.DefaultGrid_2cm
 tile = Tile(strat, grid, initT);
 u0, du0 = initialcondition!(tile, tspan)
 prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw));
+integrator = init(prob, Euler(), dt=60.0)
 
-# This is currently somewhat slow since the integrator must take very small time steps during the thawed season;
-# expect it to take about 3-5 minutes per year on a typical workstation/laptop.
-# Minor speed-ups might be possible by tweaking the dt limiters or by using the `SFCCPreSolver` for the freeze curve.
-integrator = init(prob, CGEuler())
+using BenchmarkTools
+@btime $tile($du0, $u0, $prob.p, $prob.tspan[1])
 
 # Here we take just one step to check if it's working.
 step!(integrator)
@@ -52,9 +53,12 @@ step!(integrator)
 # We can use the `getstate` function to construct the current Tile state from the integrator.
 # We then check that all water fluxes are near zero since we're starting in frozen conditions.
 state = getstate(integrator);
-@assert all(isapprox.(0.0, state.topsoil.jw, atol=1e-14))
+@assert all(isapprox.(0.0, state.ground1.jw, atol=1e-14))
 
-# Run the integrator forward in time until the end of the tspan:
+# Run the integrator forward in time until the end of the tspan;
+# Note that this is currently somewhat slow since the integrator must take very small time steps during the thawed season;
+# expect it to take about 3-5 minutes per year on a typical workstation/laptop.
+# Minor speed-ups might be possible by tweaking the dt limiters or by using the `SFCCPreSolver` for the freeze curve.
 @time while integrator.t < prob.tspan[end]
     @assert all(isfinite.(integrator.u))
     @assert all(0 .<= integrator.u.sat .<= 1)
@@ -65,6 +69,8 @@ state = getstate(integrator);
 end;
 out = CryoGridOutput(integrator.sol)
 
+@run step!(integrator)
+
 # Check mass conservation...
 water_added = values(sum(upreferred.(forcings.rainfall.(tspan[1]:Hour(3):tspan[2]).*u"m/s".*3u"hr")))[1]
 water_mass = Diagnostics.integrate(out.θwi, tile.grid)

examples/09_heat_sfcc_salt_constantbc.jl

@@ -11,7 +11,7 @@ sfcc = DallAmicoSalt(swrc=VanGenuchten(α=4.06, n=2.03))
 upperbc = SaltHeatBC(SaltGradient(benthicSalt=900.0, surfaceState=0), ConstantTemperature(0.0))
 strat = @Stratigraphy(
     0.0u"m" => Top(upperbc),
-    0.0u"m" => :sediment => SalineGround(heat=HeatBalance(:T, freezecurve=sfcc), salt=SaltMassBalance()),
+    0.0u"m" => SalineGround(heat=HeatBalance(:T, freezecurve=sfcc), salt=SaltMassBalance()),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
 modelgrid = CryoGrid.Presets.DefaultGrid_10cm

src/Physics/Heat/heat_conduction.jl

@@ -9,7 +9,7 @@ function heatflux(T₁, T₂, k₁, k₂, Δ₁, Δ₂, z₁, z₂)
 end
 
 # Free water freeze curve
-@propagate_inbounds function enthalpyinv(sub::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy}, state, i)
+@propagate_inbounds function enthalpyinv(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased}, state, i)
     θwi = Hydrology.watercontent(sub, state, i)
     H = state.H[i]
     L = heat.prop.L
@@ -35,14 +35,14 @@ end
 end
 
 """
-    freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy}, state)
+    freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased}, state)
 
 Implementation of "free water" freezing characteristic for any subsurface layer.
 Assumes that `state` contains at least temperature (T), enthalpy (H), heat capacity (C),
 total water content (θwi), and liquid water content (θw).
 """
 freezethaw!(sub, state) = freezethaw!(sub, processes(sub), state)
-function freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy}, state)
+function freezethaw!(sub::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased}, state)
     @inbounds for i in 1:length(state.H)
         # update T, θw, C
         state.T[i], state.θw[i], state.C[i] = enthalpyinv(sub, heat, state, i)
@@ -68,15 +68,15 @@ heatvariables(::HeatBalance) = (
 """
 Variable definitions for heat conduction (enthalpy) on any SubSurface layer.
 """
-CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:Enthalpy}) = (
+CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:EnthalpyBased}) = (
     Prognostic(:H, OnGrid(Cells), u"J/m^3"),
     Diagnostic(:T, OnGrid(Cells), u"°C"),
     heatvariables(heat)...,
 )
 """
 Variable definitions for heat conduction (temperature).
 """
-CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:Temperature}) = (
+CryoGrid.variables(heat::HeatBalance{<:FreezeCurve,<:TemperatureBased}) = (
     Prognostic(:T, OnGrid(Cells), u"°C"),
     Diagnostic(:H, OnGrid(Cells), u"J/m^3"),
     Diagnostic(:∂H∂t, OnGrid(Cells), u"W/m^3"),
@@ -107,7 +107,7 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalance, sub2::SubSurface, :
     return nothing
 end
 
-function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:Enthalpy}, state)
+function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:EnthalpyBased}, state)
     Δk = Δ(state.grid) # cell sizes
     ΔT = Δ(cells(state.grid)) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set;
@@ -118,7 +118,7 @@ function CryoGrid.computefluxes!(::SubSurface, ::HeatBalance{<:FreezeCurve,<:Ent
     divergence!(state.∂H∂t, state.jH, Δk)
     return nothing
 end
-function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:FreezeCurve,<:Temperature}, state)
+function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:FreezeCurve,<:TemperatureBased}, state)
     Δk = Δ(state.grid) # cell sizes
     ΔT = Δ(cells(state.grid)) # midpoint distances
     # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
@@ -134,7 +134,7 @@ function CryoGrid.computefluxes!(sub::SubSurface, ::HeatBalance{<:FreezeCurve,<:
 end
 
 # CFL not defined for free-water freeze curve
-CryoGrid.timestep(::SubSurface, heat::HeatBalance{FreeWater,<:Enthalpy,<:CryoGrid.CFL}, state) = Inf
+CryoGrid.timestep(::SubSurface, heat::HeatBalance{FreeWater,<:EnthalpyBased,<:CryoGrid.CFL}, state) = Inf
 """
     timestep(::SubSurface, ::HeatBalance{Tfc,THeatOp,CFL}, state) where {Tfc,THeatOp}
 
@@ -143,10 +143,10 @@ defined as: Δt_max = u*Δx^2, where`u` is the "characteristic velocity" which h
 is taken to be the diffusivity: `∂H∂T / kc`.
 """
 function CryoGrid.timestep(::SubSurface, heat::HeatBalance{Tfc,THeatOp,<:CryoGrid.CFL}, state) where {Tfc,THeatOp}
-    derivative(::Enthalpy, state) = state.∂H∂t
-    derivative(::Temperature, state) = state.∂T∂t
-    prognostic(::Enthalpy, state) = state.H
-    prognostic(::Temperature, state) = state.T
+    derivative(::EnthalpyBased, state) = state.∂H∂t
+    derivative(::TemperatureBased, state) = state.∂T∂t
+    prognostic(::EnthalpyBased, state) = state.H
+    prognostic(::TemperatureBased, state) = state.T
     Δx = Δ(state.grid)
     dtmax = Inf
     @inbounds for i in eachindex(Δx)

src/Physics/Heat/heat_methods.jl

@@ -76,6 +76,6 @@ heatcapacity(op::HeatOperator) = op.hc
 Retreives the nonlinear solver for the enthalpy/temperature constitutive relation, if defined.
 """
 fcsolver(::HeatOperator) = nothing
-fcsolver(op::Enthalpy) = op.fcsolver
+fcsolver(op::EnthalpyBased) = op.fcsolver
 fcsolver(heat::HeatBalance) = fcsolver(heat.op)
 fcsolver(::Nothing) = nothing