Merge branch 'dev_moritz' into 'master'

Add preliminary implementation of Lake scheme from CryoGridLite

See merge request cryogrid/cryogridjulia!116

docs/src/manual/architecture.md

@@ -46,7 +46,7 @@ The `Stratigraphy` is simply a `Tuple` of layers in ascending order of depth (i.
 
 Each `SubSurface` layer in the stratigraphy will typically consist of one or more `Process`es as fields on the layer `struct` which should then be explicitly declared via a dispatch of the [`processes`](@ref) method. The [`variables`](@ref) and [`events`](@ref) methods similarly declare state variables and events respectively that should be defined for any given configuration of the layer.
 
-The `Tile` constructor collects all of the relevant state variables declared by `variables` and discretizes them according to the given [`DiscretizationStrategy`](@ref). The resulting state vectors are initialized in the forward-diff compatible [`StateVars`](@ref) cache. On each invocation of the `Tile` (i.e. the [`evaluate!`](@ref) method), the current [`TileState`](@ref) is constructed from the current prognostic state variable `u`, parameter vector `p`, and time step `t`. The `TileState` consists of named [`LayerState`](@ref)s which provide layer-local `view`s of each state variable array, i.e. only grid cells within the layer boundaries are included.
+The `Tile` constructor collects all of the relevant state variables declared by `variables` and discretizes them according to the given [`DiscretizationStrategy`](@ref). The resulting state vectors are initialized in the forward-diff compatible [`StateVars`](@ref) cache. On each invocation of the `Tile` (i.e. the [`prognostic!`](@ref) method), the current [`TileState`](@ref) is constructed from the current prognostic state variable `u`, parameter vector `p`, and time step `t`. The `TileState` consists of named [`LayerState`](@ref)s which provide layer-local `view`s of each state variable array, i.e. only grid cells within the layer boundaries are included.
 
 ## Control flow
 
@@ -56,7 +56,7 @@ The `CryoGrid` module defines three primary methods that can be used to implemen
 2. [`interact!`](@ref) defines interactions between adjacent layers in the stratigraphy, including fluxes over the layer boundary.
 3. [`computefluxes!`](@ref) computes all internal fluxes (and the divergence thereof) within each layer, after boundary fluxes are taken into account by `interact!`.
 
-Layer and/or process specific implementations of each of these methods can generally assume that the previous methods have already been invoked by the caller (it is the responsibility of the calling code to ensure that this is the case). This is, for example, the order in which these methods will be invoked by [`evaluate!(du, u, p t)`](@ref).
+Layer and/or process specific implementations of each of these methods can generally assume that the previous methods have already been invoked by the caller (it is the responsibility of the calling code to ensure that this is the case). This is, for example, the order in which these methods will be invoked by [`prognostic!(du, u, p t)`](@ref).
 
 Note that, due to the nature of multiple dispatch, the execution path (i.e. with respect to the actual source code) of any given model configuration will typically be quite nonlinear and may span multiple source files depending on where the matching method dispatches are defined. Users may find the `which` provided by Julia (and correspondingly the `@which` macro from `InteractiveUtils`) useful in figuring out where executing code is located. For example:
 

examples/heat_freeW_lake_lite_implicit.jl

@@ -0,0 +1,69 @@
+using CryoGrid
+using CryoGrid.LiteImplicit
+using Plots
+
+Plots.plotly()
+
+CryoGrid.debug(true)
+
+forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term);
+soilprofile = SoilProfile(
+    0.0u"m" => MineralOrganic(por=0.80,sat=0.9,org=0.75),
+    0.1u"m" => MineralOrganic(por=0.80,sat=1.0,org=0.25),
+    0.4u"m" => MineralOrganic(por=0.55,sat=1.0,org=0.25),
+    3.0u"m" => MineralOrganic(por=0.50,sat=1.0,org=0.0),
+    10.0u"m" => MineralOrganic(por=0.30,sat=1.0,org=0.0),
+)
+tempprofile_linear = TemperatureProfile(
+    -2.0u"m" => 0.0u"°C",
+    0.0u"m" => 0.0u"°C",
+    10.0u"m" => -10.0u"°C", 
+    1000.0u"m" => 10.2u"°C"
+)
+modelgrid = Grid(vcat(-1.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm))
+z_top = -1.0u"m"
+z_sub = map(knot -> knot.depth, soilprofile)
+z_bot = modelgrid[end]
+upperbc = TemperatureGradient(forcings.Tair, NFactor(nf=0.5))
+initT = initializer(:T, tempprofile_linear)
+@info "Building stratigraphy"
+heatop = Heat.EnthalpyImplicit()
+strat = @Stratigraphy(
+    z_top => Top(upperbc),
+    -1.0u"m" => Lake(heat=HeatBalance(heatop)),
+    0.0u"m" => Ground(soilprofile[1].value, heat=HeatBalance(heatop)),
+    z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
+);
+@info "Building tile"
+tile = @time Tile(strat, modelgrid, initT)
+# define time span, 5 years
+tspan = (DateTime(2010,12,30), DateTime(2012,12,30))
+tspan_sol = convert_tspan(tspan)
+u0, du0 = @time initialcondition!(tile, tspan);
+prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:θw,:T,))
+# set up integrator
+integrator = init(prob, LiteImplicitEuler(), dt=24*3600)
+# debug one step
+step!(integrator)
+
+for i in integrator
+    state = getstate(integrator)
+    if state.top.T_ub[1] > 0.0
+        break
+    end
+end
+
+step!(integrator)
+state = getstate(integrator)
+state.top.T_ub
+state.lake.T
+
+@info "Running model"
+sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600.0)
+out = CryoGridOutput(sol)
+
+# Plot the results
+zs = [-100,-50.0,-10.0,11.0]u"cm"
+cg = Plots.cgrad(:copper,rev=true);
+plot(convert_t.(dims(out.T, Ti)), Array(out.T[Z(Near(zs))])' |> ustrip, color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", dpi=150)
+plot!(convert_t.(dims(out.T, Ti)), t -> forcings.Tair.(t), c=:blue, linestyle=:dash)

src/CryoGrid.jl

@@ -53,7 +53,7 @@ include("variables.jl")
 export BCKind, Volume, FixedVolume, DiagnosticVolume, PrognosticVolume
 include("traits.jl")
 
-export initialcondition!, updatestate!, interact!, interactmaybe!, computefluxes!, resetfluxes!
+export initialcondition!, updatestate!, interact!, interactmaybe!, computefluxes!, resetfluxes!, diagnosticstep!
 export variables, processes, initializers, timestep, isactive, caninteract
 export boundaryflux, boundaryvalue, criterion, criterion!, trigger!
 include("methods.jl")

src/Physics/Heat/heat_implicit.jl

@@ -3,6 +3,26 @@ Type alias for the implicit enthalpy formulation of HeatBalance.
 """
 const HeatBalanceImplicit{Tfc} = HeatBalance{Tfc,<:EnthalpyImplicit} where {Tfc<:FreezeCurve}
 
+apbc(::Dirichlet, k, Δk) = k / (Δk^2/2)
+apbc(::Dirichlet, k, Δk, Δx) = k / Δx / Δk
+apbc(::Neumann, k, Δk) = 0
+
+function prefactors!(ap, an, as, k, dx, dxp)
+    # loop over grid cells
+    @inbounds for i in eachindex(dxp)
+        if i == 1
+            as[1] = k[2] / dx[1] / dxp[1]
+        elseif i == length(dxp)
+            an[end] = k[end-1] / dx[end] / dxp[end]
+        else
+            an[i] = k[i] / dx[i-1] / dxp[i]
+            as[i] = k[i+1] / dx[i] / dxp[i]
+        end
+        ap[i] = an[i] + as[i]
+    end
+    return nothing
+end
+
 CryoGrid.variables(::HeatBalanceImplicit) = (
     Prognostic(:H, OnGrid(Cells), u"J/m^3"),
     Diagnostic(:T, OnGrid(Cells), u"°C"),
@@ -12,12 +32,13 @@ CryoGrid.variables(::HeatBalanceImplicit) = (
     Diagnostic(:k, OnGrid(Edges), u"W/m/K"),
     Diagnostic(:kc, OnGrid(Cells), u"W/m/K"),
     Diagnostic(:θw, OnGrid(Cells), domain=0..1),
-    # coefficients and cache variables for diffusion operator
+    # coefficients and cache variables for implicit diffusion operator
     Diagnostic(:DT_an, OnGrid(Cells)),
     Diagnostic(:DT_as, OnGrid(Cells)),
     Diagnostic(:DT_ap, OnGrid(Cells)),
     Diagnostic(:DT_bp, OnGrid(Cells)),
 )
+
 function CryoGrid.updatestate!(
     sub::SubSurface,
     heat::HeatBalanceImplicit,
@@ -34,29 +55,10 @@ function CryoGrid.updatestate!(
     k = state.k
     dx = Δ(cells(state.grid))
     dxp = Δ(state.grid)
-    # loop over grid cells
-    @inbounds for i in eachindex(dxp)
-        if i == 1
-            as[1] = k[2] / dx[1] / dxp[1]
-        elseif i == length(dxp)
-            an[end] = k[end-1] / dx[end] / dxp[end]
-        else
-            an[i] = k[i] / dx[i-1] / dxp[i]
-            as[i] = k[i+1] / dx[i] / dxp[i]
-        end
-    end
-    @. ap = an + as
-    # k_inner = @view k[2:end-1]
-    # dxn = @view dx[1:end-1]
-    # dxs = @view dx[2:end]
-    # dxpn = @view dxp[1:end-1]
-    # dxps = @view dxp[2:end]
-    # @. an[2:end] = k_inner / dxn / dxpn
-    # @. as[1:end-1] = k_inner / dxs / dxps
-    # @. ap[1:end-1] += as[1:end-1]
-    # @. ap[2:end] += an[2:end]
+    prefactors!(ap, an, as, k, dx, dxp)
     return nothing
 end
+
 function CryoGrid.boundaryflux(::Dirichlet, bc::HeatBC, top::Top, heat::HeatBalanceImplicit, sub::SubSurface, stop, ssub)
     T_ub = boundaryvalue(bc, stop)
     k = ssub.k[1]
@@ -69,22 +71,21 @@ function CryoGrid.boundaryflux(::Dirichlet, bc::HeatBC, bot::Bottom, heat::HeatB
     Δk = CryoGrid.thickness(sub, ssub, last)
     return 2*T_lb*k / Δk
 end
-_ap(::Dirichlet, k, Δk) = k / (Δk^2/2)
-_ap(::Neumann, k, Δk) = 0
+
 function CryoGrid.interact!(top::Top, bc::HeatBC, sub::SubSurface, heat::HeatBalanceImplicit, stop, ssub)
     Δk = CryoGrid.thickness(sub, ssub, first)
     jH_top = boundaryflux(bc, top, heat, sub, stop, ssub)
     k = ssub.k[1]
     ssub.DT_bp[1] += jH_top / Δk
-    ssub.DT_ap[1] += _ap(CryoGrid.BCKind(bc), k, Δk)
+    ssub.DT_ap[1] += apbc(CryoGrid.BCKind(bc), k, Δk)
     return nothing
 end
 function CryoGrid.interact!(sub::SubSurface, heat::HeatBalanceImplicit, bot::Bottom, bc::HeatBC, ssub, sbot)
     Δk = CryoGrid.thickness(sub, ssub, last)
     jH_bot = boundaryflux(bc, bot, heat, sub, sbot, ssub)
     k = ssub.k[1]
     ssub.DT_bp[end] += jH_bot / Δk
-    ssub.DT_ap[end] += _ap(CryoGrid.BCKind(bc), k, Δk)
+    ssub.DT_ap[end] += apbc(CryoGrid.BCKind(bc), k, Δk)
     return nothing
 end
 function CryoGrid.interact!(sub1::SubSurface, ::HeatBalanceImplicit, sub2::SubSurface, ::HeatBalanceImplicit, s1, s2)
@@ -103,6 +104,7 @@ function CryoGrid.interact!(sub1::SubSurface, ::HeatBalanceImplicit, sub2::SubSu
     s2.DT_ap[1] += s2.DT_an[1] = k / Δz / Δk₂
     return nothing
 end
+
 # do nothing in computefluxes!
 CryoGrid.computefluxes!(::SubSurface, ::HeatBalanceImplicit, state) = nothing
 

src/Physics/Lakes/Lakes.jl

@@ -0,0 +1,10 @@
+module Lakes
+
+using CryoGrid
+using CryoGrid.Utils
+using CryoGrid.Numerics
+
+export Lake
+include("lake_simple.jl")
+
+end

src/Physics/Lakes/lake_simple.jl

@@ -0,0 +1,135 @@
+abstract type LakeParameterization <: CryoGrid.Parameterization end
+
+Base.@kwdef struct SimpleLakeScheme{Thp} <: LakeParameterization
+    heat::Thp = ThermalProperties()
+end
+
+Base.@kwdef struct Lake{Tpara<:LakeParameterization,Theat<:HeatBalance,Taux} <: CryoGrid.SubSurface
+    para::Tpara = SimpleLakeScheme()
+    heat::Theat = HeatBalance()
+    aux::Taux = nothing
+end
+
+function get_upper_boundary_index(T_ub, θw)
+    ubc_idx = 1
+    if T_ub >= zero(T_ub)
+        @inbounds for i in eachindex(θw)
+            ubc_idx = i
+            if θw[i] < one(eltype(θw))
+                return ubc_idx
+            end
+        end
+    end
+    return ubc_idx+1
+end
+
+# Material properties
+Heat.thermalproperties(lake::Lake) = lake.para.heat
+
+Hydrology.watercontent(::Lake, state) = 1.0
+
+CryoGrid.processes(lake::Lake) = lake.heat
+
+CryoGrid.volumetricfractions(::Lake, state, i) = (state.θw[i], 1 - state.θw[i], 0.0)
+
+CryoGrid.variables(lake::Lake, heat::HeatBalance) = (
+    CryoGrid.variables(heat)...,
+    # Diagnostic(:I_t, OnGrid(Cells), desc="Indicator variable for is thawed (1 or 0)."),
+    Diagnostic(:ρ_w, Scalar, u"kg*m^-3", domain=0..Inf, desc = "density of water with temperature"),
+    Diagnostic(:T_ub, Scalar, u"°C"),
+    Diagnostic(:ubc_idx, Scalar, NoUnits, Int),
+)
+
+function CryoGrid.diagnosticstep!(::Lake, state)
+    T_ub = getscalar(state.T_ub)
+    @setscalar state.ubc_idx = get_upper_boundary_index(T_ub, state.θw)
+    return false
+end
+
+function CryoGrid.initialcondition!(lake::Lake, state)
+    initialcondition!(lake, lake.heat, state)
+    diagnosticstep!(lake, state)
+end
+
+function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake)
+    T_ub = slake.T_ub[1] = getscalar(stop.T_ub)
+    ubc_idx = Int(getscalar(slake.ubc_idx))
+    # get variables
+    an = slake.DT_an
+    as = slake.DT_as
+    ap = slake.DT_ap
+    bp = slake.DT_bp
+    k = slake.k
+    dx = Δ(cells(slake.grid))
+    dxp = Δ(slake.grid)
+    # outer lake cell
+    bp[1] = T_ub*k[1] / (dxp[1]/2) / dxp[1]
+    ap[1] = k[1] / (dxp[1]/2) / dxp[1]
+    if ubc_idx > 1
+        as[1] = an[1] = zero(eltype(as))
+    end
+    # inner lake cells
+    @inbounds for i in 2:ubc_idx-1
+        bp[i] = T_ub*k[i] / dx[i-1] / dxp[i]
+        ap[i] = an[i]
+        as[i] = zero(eltype(as))
+        an[i] = zero(eltype(an))
+    end
+    return nothing
+end
+function CryoGrid.interact!(lake::Lake, ::HeatBalanceImplicit, sub::SubSurface, ::HeatBalanceImplicit, slake, ssub)
+    Δk₁ = CryoGrid.thickness(lake, slake, last)
+    Δk₂ = CryoGrid.thickness(sub, ssub, first)
+    Δz = CryoGrid.midpoint(sub, ssub, first) - CryoGrid.midpoint(lake, slake, last)
+    ubc_idx = Int(getscalar(slake.ubc_idx))
+    # thermal conductivity between cells
+    k = slake.k[end] = ssub.k[1] =
+        @inbounds let k₁ = slake.kc[end],
+            k₂ = ssub.kc[1],
+            Δ₁ = Δk₁[end],
+            Δ₂ = Δk₂[1];
+            harmonicmean(k₁, k₂, Δ₁, Δ₂)
+        end
+    slake.DT_ap[end] += slake.DT_as[end] = (ubc_idx <= length(slake.θw))*k / Δz / Δk₁
+    ssub.DT_ap[1] += ssub.DT_an[1] = k / Δz / Δk₂
+    return nothing
+end
+
+function CryoGrid.updatestate!(sub::Lake, heat::HeatBalance{FreeWater,<:Heat.MOLEnthalpy}, state)
+    Heat.resetfluxes!(sub, heat, state)
+    # Evaluate freeze/thaw processes
+    Heat.freezethaw!(sub, heat, state)
+
+    # force unfrozen water to Tair
+    if all(state.θw .>= 1.)
+        state.T .= state.T_ub
+        state.H .= Heat.enthalpy.(state.T, state.C, heat.prop.L, state.θw)
+    end
+
+    # Update thermal conductivity
+    Heat.thermalconductivity!(sub, heat, state)
+    return nothing
+end
+
+function CryoGrid.computefluxes!(::Lake, ::HeatBalance{FreeWater,<:Heat.MOLEnthalpy}, state)
+    Δk = Δ(state.grids.k) # cell sizes
+    ΔT = Δ(state.grids.T) # midpoint distances
+    # compute internal fluxes and non-linear diffusion assuming boundary fluxes have been set
+    Numerics.nonlineardiffusion!(state.∂H∂t, state.jH, state.T, ΔT, state.k, Δk)
+    return nothing
+end
+
+
+function CryoGrid.interact!(
+    top::Top,
+    bc::HeatBC,
+    lake::Lake,
+    heat::HeatBalance{FreeWater,<:Heat.MOLEnthalpy},
+    stop,
+    slake
+)
+    @setscalar slake.T_ub = stop.T_ub
+    # boundary flux
+    slake.jH[1] += CryoGrid.boundaryflux(bc, top, heat, lake, stop, slake)
+    return nothing
+end

src/Physics/Physics.jl

@@ -40,3 +40,5 @@ include("Surface/Surface.jl")
 @reexport using .Surface
 include("Sources/Sources.jl")
 @reexport using .Sources
+include("Lakes/Lakes.jl")
+@reexport using .Lakes

src/Solvers/LiteImplicit/cglite_step.jl

@@ -9,9 +9,12 @@ function CryoGrid.perform_step!(integrator::CGLiteIntegrator)
     t = t₀ + dt
     tile = Tiles.resolve(Tile(integrator.sol.prob.f), u, p, t)
     # explicit update, if necessary
-    explicit_step!(integrator, tile, du, u, p, t)
+    # explicit_step!(integrator, tile, du, u, p, t)
     # implicit update for energy state
     implicit_step!(integrator, tile, du, u, p, t)
+    # diagnostic udpate
+    CryoGrid.diagnosticstep!(tile, getstate(integrator))
+    # update t and step number
     integrator.t = t
     integrator.step += 1
     return nothing

src/Tiles/stratigraphy.jl

@@ -146,6 +146,12 @@ function CryoGrid.updatestate!(strat::Stratigraphy, state)
     end
 end
 
+function CryoGrid.diagnosticstep!(strat::Stratigraphy, state)
+    fastiterate(namedlayers(strat)) do named_layer
+        CryoGrid.diagnosticstep!(named_layer.val, getproperty(state, nameof(named_layer)))
+    end
+end
+
 """
     interact!(strat::Stratigraphy, state)