From 767b96447e12576ea55384628c038e037b8923c7 Mon Sep 17 00:00:00 2001 From: Moritz Langer Date: Thu, 23 Feb 2023 15:49:32 +0100 Subject: [PATCH 01/17] lake scheme --- src/Physics/Lakes/Lakes.jl | 80 ++++++++++++++++++++++++++++++++++++++ src/Physics/Physics.jl | 2 + 2 files changed, 82 insertions(+) create mode 100644 src/Physics/Lakes/Lakes.jl diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl new file mode 100644 index 00000000..a5a981b7 --- /dev/null +++ b/src/Physics/Lakes/Lakes.jl @@ -0,0 +1,80 @@ +module Lakes +import CryoGrid +using CryoGrid.Physics +using CryoGrid.Utils +using CryoGrid.Numerics +using Unitful + +export Lake + +abstract type LakeParameterization <: CryoGrid.Parameterization end + +struct SimpleLakeScheme <: LakeParameterization end + +Base.@kwdef struct Lake{Tpara<:LakeParameterization,Tsp,Tproc,Tprop} <: CryoGrid.SubSurface{Tproc} + para::Tpara = SimpleLakeScheme() + prop::Tprop = CryoGrid.ThermalProperties() + sp::Tsp = nothing + proc::Tproc +end + +Lake(proc::Tproc; kwargs...) where {Tproc} = Lake(;proc, kwargs...) + +CryoGrid.variables(lake::Lake, heat::CryoGrid.HeatBalance) = ( + CryoGrid.variables(heat)..., + Diagnostic(:ρ_w, Scalar, u"kg*m^-3", domain=0..Inf, desc = "density of water with temperature"), + Diagnostic(:T_ub, Scalar, u"°C"), +) + +function CryoGrid.initialcondition!(lake::Lake, heat::HeatBalance, state) + L = heat.prop.L + # initialize liquid water content based on temperature + @inbounds for i in 1:length(state.T) + θwi = Hydrology.watercontent(lake, state, i) + state.θw[i] = ifelse(state.T[i] > 0.0, θwi, 0.0) + state.C[i] = heatcapacity(lake, heat, volumetricfractions(lake, state, i)...) + state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i]) + end +end + +function CryoGrid.diagnosticstep!(sub::Lake, heat::HeatBalance, 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.prognosticstep!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.Enthalpy}, 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, + stop, + slake +) + slake.T_ub = CryoGrid.boundaryvalue(bc, top, heat, lake, stop, slake) + # boundary flux + slake.jH[1] += CryoGrid.boundaryflux(bc, top, heat, lake, stop, slake) + return nothing +end + + +end \ No newline at end of file diff --git a/src/Physics/Physics.jl b/src/Physics/Physics.jl index 95824833..d78695b7 100644 --- a/src/Physics/Physics.jl +++ b/src/Physics/Physics.jl @@ -59,6 +59,7 @@ include("Snow/Snow.jl") include("Soils/Soils.jl") include("SEB/SEB.jl") include("Sources/Sources.jl") +include("Lakes/Lakes.jl") @reexport using .Boundaries @reexport using .Heat @@ -67,5 +68,6 @@ include("Sources/Sources.jl") @reexport using .Soils @reexport using .SEB @reexport using .Sources +@reexport using .Lakes end \ No newline at end of file From 08432b47abf88f8040f65ccdf52a1db7f7a3e5ce Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Thu, 23 Feb 2023 16:10:06 +0100 Subject: [PATCH 02/17] Initial fix for import issues --- src/Physics/Lakes/Lakes.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index a5a981b7..57042d56 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -1,8 +1,12 @@ module Lakes + import CryoGrid -using CryoGrid.Physics + +using CryoGrid +using CryoGrid.Physics.Heat using CryoGrid.Utils using CryoGrid.Numerics + using Unitful export Lake @@ -20,7 +24,7 @@ end Lake(proc::Tproc; kwargs...) where {Tproc} = Lake(;proc, kwargs...) -CryoGrid.variables(lake::Lake, heat::CryoGrid.HeatBalance) = ( +CryoGrid.variables(lake::Lake, heat::HeatBalance) = ( CryoGrid.variables(heat)..., Diagnostic(:ρ_w, Scalar, u"kg*m^-3", domain=0..Inf, desc = "density of water with temperature"), Diagnostic(:T_ub, Scalar, u"°C"), From d437d935be0e902d01fd153c52aac4285868de36 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Fri, 24 Feb 2023 15:05:01 +0100 Subject: [PATCH 03/17] Re-add Lakes module inclusion post-merge --- src/Physics/Physics.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Physics/Physics.jl b/src/Physics/Physics.jl index 9a7dbbeb..0704ecb7 100644 --- a/src/Physics/Physics.jl +++ b/src/Physics/Physics.jl @@ -42,3 +42,5 @@ include("SEB/SEB.jl") @reexport using .SEB include("Sources/Sources.jl") @reexport using .Sources +include("Lakes/Lakes.jl") +@reexport using .Lakes From 3fae713ee978a53e5bfa276e75a9a45cef446441 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Fri, 24 Feb 2023 15:17:48 +0100 Subject: [PATCH 04/17] Add example script for lakes --- examples/heat_freeW_samoylov_lakes.jl | 51 +++++++++++++++++++++++++++ src/Physics/Lakes/Lakes.jl | 7 +--- 2 files changed, 52 insertions(+), 6 deletions(-) create mode 100644 examples/heat_freeW_samoylov_lakes.jl diff --git a/examples/heat_freeW_samoylov_lakes.jl b/examples/heat_freeW_samoylov_lakes.jl new file mode 100644 index 00000000..649d790d --- /dev/null +++ b/examples/heat_freeW_samoylov_lakes.jl @@ -0,0 +1,51 @@ +using CryoGrid +using Plots + +forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C"); +# use air temperature as upper boundary forcing; +tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); +T0 = values(tair[tspan[1]])[1] +tempprofile = TemperatureProfile( + 0.0u"m" => T0, + 1.0u"m" => -8.0u"°C", + 20.0u"m" => -10u"°C", + 1000.0u"m" => 1.0u"°C" +) +# soil profile: depth => (excess ice, natural porosity, saturation, organic fraction) +soilprofile = SoilProfile( + 0.0u"m" => HomogeneousMixture(por=0.80,sat=1.0,org=0.75), + 0.1u"m" => HomogeneousMixture(por=0.80,sat=1.0,org=0.25), + 0.4u"m" => HomogeneousMixture(por=0.55,sat=1.0,org=0.25), + 3.0u"m" => HomogeneousMixture(por=0.50,sat=1.0,org=0.0), + 10.0u"m" => HomogeneousMixture(por=0.30,sat=1.0,org=0.0), + 100.0u"m" => HomogeneousMixture(por=0.10,sat=0.1,org=0.0), +); +initT = initializer(:T, tempprofile) +# Heat diffusion with temperature as prognostic state variable +op = Heat.Diffusion(:H) +# Construct soil layers +soil_layers = map(enumerate(soilprofile)) do (i, (depth, para)) + name = Symbol(:soil, i) + depth => name => Soil(HeatBalance(op, dtlim=CryoGrid.MaxDelta(50u"kJ")), para=para) +end +# Build stratigraphy: +# @Stratigraphy macro lets us list multiple subsurface layers +strat = @Stratigraphy( + -2.0*u"m" => Top(TemperatureGradient(tair)), + -2.0u"m" => Lake(), + soil_layers..., + 1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")), +); +# define time span +tspan = (DateTime(2010,10,30),DateTime(2011,10,30)) +u0, du0 = initialcondition!(tile, tspan) +# construct CryoGridProblem with tile, initial condition, and timespan; +# we disable the default timestep limiter since we will use an adaptive solver. +prob = CryoGridProblem(tile, u0, tspan, savevars=(:T,)) +@info "Running model" +out = @time solve(prob, Euler(), dt=300.0, saveat=24*3600.0, progress=true) |> CryoGridOutput; +# Plot it! +zs = [1,10,20,30,50,100,200,500,1000]u"cm" +cg = Plots.cgrad(:copper,rev=true); +plot(out.H[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150) +plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150) diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 57042d56..411ceffb 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -1,14 +1,9 @@ module Lakes -import CryoGrid - using CryoGrid -using CryoGrid.Physics.Heat using CryoGrid.Utils using CryoGrid.Numerics -using Unitful - export Lake abstract type LakeParameterization <: CryoGrid.Parameterization end @@ -74,7 +69,7 @@ function CryoGrid.interact!( stop, slake ) - slake.T_ub = CryoGrid.boundaryvalue(bc, top, heat, lake, stop, slake) + @setscalar slake.T_ub = CryoGrid.boundaryvalue(bc, top, heat, lake, stop, slake) # boundary flux slake.jH[1] += CryoGrid.boundaryflux(bc, top, heat, lake, stop, slake) return nothing From ed881361986cd0e40be4e234521e709a63c8a5af Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Fri, 24 Feb 2023 15:34:36 +0100 Subject: [PATCH 05/17] Add necessary material properties dispatches for Lake --- examples/heat_freeW_samoylov_lakes.jl | 16 ++++++++++------ src/Physics/Lakes/Lakes.jl | 10 +++++++--- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/examples/heat_freeW_samoylov_lakes.jl b/examples/heat_freeW_samoylov_lakes.jl index 649d790d..04534f65 100644 --- a/examples/heat_freeW_samoylov_lakes.jl +++ b/examples/heat_freeW_samoylov_lakes.jl @@ -2,6 +2,8 @@ using CryoGrid using Plots forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C"); +# define time span +tspan = (DateTime(2010,10,30),DateTime(2011,10,30)) # use air temperature as upper boundary forcing; tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); T0 = values(tair[tspan[1]])[1] @@ -22,26 +24,28 @@ soilprofile = SoilProfile( ); initT = initializer(:T, tempprofile) # Heat diffusion with temperature as prognostic state variable -op = Heat.Diffusion(:H) # Construct soil layers soil_layers = map(enumerate(soilprofile)) do (i, (depth, para)) name = Symbol(:soil, i) - depth => name => Soil(HeatBalance(op, dtlim=CryoGrid.MaxDelta(50u"kJ")), para=para) + depth => name => Soil(HeatBalance(), para=para) end # Build stratigraphy: # @Stratigraphy macro lets us list multiple subsurface layers strat = @Stratigraphy( -2.0*u"m" => Top(TemperatureGradient(tair)), - -2.0u"m" => Lake(), + -2.0u"m" => :lake => Lake(HeatBalance()), soil_layers..., - 1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")), + 998.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")), ); -# define time span -tspan = (DateTime(2010,10,30),DateTime(2011,10,30)) +# shift grid up by 2 m +modelgrid = Grid(CryoGrid.Presets.DefaultGrid_5cm .- 2.0u"m") +tile = Tile(strat, modelgrid, initT) u0, du0 = initialcondition!(tile, tspan) # construct CryoGridProblem with tile, initial condition, and timespan; # we disable the default timestep limiter since we will use an adaptive solver. prob = CryoGridProblem(tile, u0, tspan, savevars=(:T,)) +# test step function (good to do debugging here) +tile(du0, u0, prob.p, prob.tspan[1]) @info "Running model" out = @time solve(prob, Euler(), dt=300.0, saveat=24*3600.0, progress=true) |> CryoGridOutput; # Plot it! diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 411ceffb..e3dce3ca 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -19,6 +19,11 @@ end Lake(proc::Tproc; kwargs...) where {Tproc} = Lake(;proc, kwargs...) +# Material properties +Heat.thermalproperties(lake::Lake) = lake.prop +Hydrology.watercontent(::Lake, state) = 1.0 +CryoGrid.volumetricfractions(::Lake, state, i) = (state.θw[i], 1 - state.θw[i], 0.0) + CryoGrid.variables(lake::Lake, heat::HeatBalance) = ( CryoGrid.variables(heat)..., Diagnostic(:ρ_w, Scalar, u"kg*m^-3", domain=0..Inf, desc = "density of water with temperature"), @@ -42,9 +47,9 @@ function CryoGrid.diagnosticstep!(sub::Lake, heat::HeatBalance, state) Heat.freezethaw!(sub, heat, state) # force unfrozen water to Tair - if all(state.Θw .>= 1.) + if all(state.θw .>= 1.) state.T .= state.T_ub - state.H .= Heat.enthalpy(state.T, state.C, heat.prop.L, state.Θw) + state.H .= Heat.enthalpy.(state.T, state.C, heat.prop.L, state.θw) end # Update thermal conductivity @@ -52,7 +57,6 @@ function CryoGrid.diagnosticstep!(sub::Lake, heat::HeatBalance, state) return nothing end - function CryoGrid.prognosticstep!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.Enthalpy}, state) Δk = Δ(state.grids.k) # cell sizes ΔT = Δ(state.grids.T) # midpoint distances From 03843d504e35ea2db63ab22459c9f1a12ba10dde Mon Sep 17 00:00:00 2001 From: Moritz Langer Date: Fri, 3 Mar 2023 18:36:42 +0100 Subject: [PATCH 06/17] lake thing Co-authored-by: Brian Groenke --- examples/heat_freeW_lake_lite_implicit.jl | 56 +++++++++++++++++++++++ src/Physics/Lakes/Lakes.jl | 45 ++++++++++++++++++ 2 files changed, 101 insertions(+) create mode 100644 examples/heat_freeW_lake_lite_implicit.jl diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl new file mode 100644 index 00000000..fc4a35a2 --- /dev/null +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -0,0 +1,56 @@ +using CryoGrid +using CryoGrid.LiteImplicit +using Plots + +forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term, :Tair => u"°C", :Ptot => u"mm"); +# use air temperature as upper boundary forcing +tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); +# snowdepth = TimeSeriesForcing(ustrip.(forcings.data.Dsn), forcings.timestamps, :Dsn); +soilprofile = SoilProfile( + 0.0u"m" => HomogeneousMixture(por=0.80,sat=0.9,org=0.75), #(θwi=0.80,θm=0.05,θo=0.15,ϕ=0.80), + 0.1u"m" => HomogeneousMixture(por=0.80,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.80), + 0.4u"m" => HomogeneousMixture(por=0.55,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55), + 3.0u"m" => HomogeneousMixture(por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50), + 10.0u"m" => HomogeneousMixture(por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30), +) +tempprofile_linear = TemperatureProfile( + 0.0u"m" => -30.0u"°C", + 10.0u"m" => -10.0u"°C", + 1000.0u"m" => 10.2u"°C" +) +modelgrid = Grid(CryoGrid.Presets.DefaultGrid_2cm .- 2.0u"m") +z_top = -2.0u"m" +z_sub = map(knot -> knot.depth, soilprofile) +z_bot = modelgrid[end] +upperbc = TemperatureGradient(tair, NFactor()) +initT = initializer(:T, tempprofile_linear) +@info "Building stratigraphy" +heatop = Heat.EnthalpyImplicit() +strat = @Stratigraphy( + z_top => Top(upperbc), + z_top => :lake => Lake(HeatBalance(heatop)), + z_sub[1] => :topsoil1 => Soil(HeatBalance(heatop), para=soilprofile[1].value), + z_sub[2] => :topsoil2 => Soil(HeatBalance(heatop), para=soilprofile[2].value), + z_sub[3] => :sediment1 => Soil(HeatBalance(heatop), para=soilprofile[3].value), + z_sub[4] => :sediment2 => Soil(HeatBalance(heatop), para=soilprofile[4].value), + z_sub[5] => :sediment3 => Soil(HeatBalance(heatop), para=soilprofile[5].value), + 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(2015,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,)) +lake_state = getstate(:lake, tile, u0, du0, prob.tspan[1]) +# debug one step +@run tile(du0, u0, prob.p, prob.tspan[1]) +@info "Running model" +sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600) +out = CryoGridOutput(sol) + +# Plot the results +zs = [5,10,15,20,25,30,40,50,100,500,1000,5000]u"cm" +cg = Plots.cgrad(:copper,rev=true); +plot(out.T[Z(Near(zs))] |> ustrip, color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", leg=false, dpi=150) diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index e3dce3ca..9477305d 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -26,6 +26,7 @@ CryoGrid.volumetricfractions(::Lake, state, i) = (state.θw[i], 1 - state.θw[i] 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"), ) @@ -57,6 +58,41 @@ function CryoGrid.diagnosticstep!(sub::Lake, heat::HeatBalance, state) return nothing end +function CryoGrid.diagnosticstep!( + sub::Lake, + heat::HeatBalanceImplicit, + state +) + Heat.resetfluxes!(sub, heat, state) + # Evaluate freeze/thaw processes + Heat.freezethaw!(sub, heat, state) + # Update thermal conductivity + Heat.thermalconductivity!(sub, heat, state) + isthawed = true + @inbounds for i in eachindex(state.θw) + isthawed = isthawed && state.θw[i] ≈ 1.0 + end + I_t = Float64(isthawed) + # Compute diffusion coefficients + an = state.DT_an + as = state.DT_as + ap = state.DT_ap + k = state.k + dx = Δ(cells(state.grid)) + dxp = Δ(state.grid) + k_inner = @view k[2:end-1] + dxpn = @view dxp[1:end-1] + dxps = @view dxp[2:end] + @. an[2:end] = k_inner / dx / dxpn + @. as[1:end-1] = (k_inner / dx / dxps)*I_t + @. ap[1:end-1] += as[1:end-1] + @. ap[2:end] += an[2:end] + @. an[2:end] *= I_t + return nothing +end + +CryoGrid.prognosticstep!(::Lake, ::HeatBalanceImplicit, state) = nothing + function CryoGrid.prognosticstep!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.Enthalpy}, state) Δk = Δ(state.grids.k) # cell sizes ΔT = Δ(state.grids.T) # midpoint distances @@ -65,6 +101,15 @@ function CryoGrid.prognosticstep!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.Ent return nothing end +function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake) + Δk = CryoGrid.thickness(lake, slake, first) + jH_top = boundaryflux(bc, top, heat, lake, stop, slake) + k = slake.k[1] + slake.DT_bp[1] += jH_top / Δk + slake.DT_ap[1] += Heat._ap(CryoGrid.BoundaryStyle(bc), k, Δk) + return nothing +end + function CryoGrid.interact!( top::Top, bc::HeatBC, From 3c09ab14aba0b4a8550553f70a4efab0ee0c927c Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Sat, 4 Mar 2023 16:17:54 +0100 Subject: [PATCH 07/17] WIP --- examples/heat_freeW_lake_lite_implicit.jl | 7 +++++-- src/Physics/Lakes/Lakes.jl | 6 +++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index fc4a35a2..e79590f3 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -18,7 +18,7 @@ tempprofile_linear = TemperatureProfile( 10.0u"m" => -10.0u"°C", 1000.0u"m" => 10.2u"°C" ) -modelgrid = Grid(CryoGrid.Presets.DefaultGrid_2cm .- 2.0u"m") +modelgrid = Grid(vcat(-2.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm)) z_top = -2.0u"m" z_sub = map(knot -> knot.depth, soilprofile) z_bot = modelgrid[end] @@ -45,7 +45,10 @@ u0, du0 = @time initialcondition!(tile, tspan); prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:θw,:T,)) lake_state = getstate(:lake, tile, u0, du0, prob.tspan[1]) # debug one step -@run tile(du0, u0, prob.p, prob.tspan[1]) +CryoGrid.debug(true) +tile(du0, u0, prob.p, prob.tspan[1]) +integrator = init(prob, LiteImplicitEuler(), dt=24*3600) +step!(integrator) @info "Running model" sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600) out = CryoGridOutput(sol) diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 9477305d..5591b060 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -72,7 +72,7 @@ function CryoGrid.diagnosticstep!( @inbounds for i in eachindex(state.θw) isthawed = isthawed && state.θw[i] ≈ 1.0 end - I_t = Float64(isthawed) + I_f = 1 - Float64(isthawed) # Compute diffusion coefficients an = state.DT_an as = state.DT_as @@ -84,10 +84,10 @@ function CryoGrid.diagnosticstep!( dxpn = @view dxp[1:end-1] dxps = @view dxp[2:end] @. an[2:end] = k_inner / dx / dxpn - @. as[1:end-1] = (k_inner / dx / dxps)*I_t + @. as[1:end-1] = (k_inner / dx / dxps)*I_f @. ap[1:end-1] += as[1:end-1] @. ap[2:end] += an[2:end] - @. an[2:end] *= I_t + @. an[2:end] *= I_f return nothing end From b8e71b859d5764d3af2dd57fd69cc21c2e8a27d0 Mon Sep 17 00:00:00 2001 From: Moritz Langer Date: Mon, 2 Oct 2023 13:01:41 +0200 Subject: [PATCH 08/17] have no idea --- examples/heat_freeW_lake_lite_implicit.jl | 24 +++++++++++++++++------ src/Physics/Lakes/Lakes.jl | 18 +++++++++++++++++ 2 files changed, 36 insertions(+), 6 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index e79590f3..caad1078 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -2,6 +2,8 @@ using CryoGrid using CryoGrid.LiteImplicit using Plots +CryoGrid.debug(true) + forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term, :Tair => u"°C", :Ptot => u"mm"); # use air temperature as upper boundary forcing tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); @@ -18,7 +20,11 @@ tempprofile_linear = TemperatureProfile( 10.0u"m" => -10.0u"°C", 1000.0u"m" => 10.2u"°C" ) +<<<<<<< HEAD modelgrid = Grid(vcat(-2.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm)) +======= +modelgrid = CryoGrid.Presets.DefaultGrid_2cm +>>>>>>> have no idea z_top = -2.0u"m" z_sub = map(knot -> knot.depth, soilprofile) z_bot = modelgrid[end] @@ -28,12 +34,12 @@ initT = initializer(:T, tempprofile_linear) heatop = Heat.EnthalpyImplicit() strat = @Stratigraphy( z_top => Top(upperbc), - z_top => :lake => Lake(HeatBalance(heatop)), - z_sub[1] => :topsoil1 => Soil(HeatBalance(heatop), para=soilprofile[1].value), - z_sub[2] => :topsoil2 => Soil(HeatBalance(heatop), para=soilprofile[2].value), - z_sub[3] => :sediment1 => Soil(HeatBalance(heatop), para=soilprofile[3].value), - z_sub[4] => :sediment2 => Soil(HeatBalance(heatop), para=soilprofile[4].value), - z_sub[5] => :sediment3 => Soil(HeatBalance(heatop), para=soilprofile[5].value), + z_sub[1] => :lake => Lake(HeatBalance(heatop)), + #z_sub[1] => :topsoil1 => Soil(HeatBalance(heatop), para=soilprofile[1].value), + #z_sub[2] => :topsoil2 => Soil(HeatBalance(heatop), para=soilprofile[2].value), + #z_sub[3] => :sediment1 => Soil(HeatBalance(heatop), para=soilprofile[3].value), + #z_sub[4] => :sediment2 => Soil(HeatBalance(heatop), para=soilprofile[4].value), + #z_sub[5] => :sediment3 => Soil(HeatBalance(heatop), para=soilprofile[5].value), z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2")) ); @info "Building tile" @@ -45,10 +51,16 @@ u0, du0 = @time initialcondition!(tile, tspan); prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:θw,:T,)) lake_state = getstate(:lake, tile, u0, du0, prob.tspan[1]) # debug one step +<<<<<<< HEAD CryoGrid.debug(true) tile(du0, u0, prob.p, prob.tspan[1]) integrator = init(prob, LiteImplicitEuler(), dt=24*3600) step!(integrator) +======= +@run tile(du0, u0, prob.p, prob.tspan[1]) +integrator = @time init(prob, LiteImplicitEuler(), dt=24*3600) +@run step!(integrator) +>>>>>>> have no idea @info "Running model" sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600) out = CryoGridOutput(sol) diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 5591b060..f47caabf 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -68,11 +68,21 @@ function CryoGrid.diagnosticstep!( Heat.freezethaw!(sub, heat, state) # Update thermal conductivity Heat.thermalconductivity!(sub, heat, state) +<<<<<<< HEAD isthawed = true @inbounds for i in eachindex(state.θw) isthawed = isthawed && state.θw[i] ≈ 1.0 end I_f = 1 - Float64(isthawed) +======= + #isthawed = true + #@inbounds for i in eachindex(state.θw) + # isthawed = isthawed && state.θw[i] ≈ 1.0 + #end + + isthawed = state.θw .≈ 1.0 + +>>>>>>> have no idea # Compute diffusion coefficients an = state.DT_an as = state.DT_as @@ -84,10 +94,18 @@ function CryoGrid.diagnosticstep!( dxpn = @view dxp[1:end-1] dxps = @view dxp[2:end] @. an[2:end] = k_inner / dx / dxpn +<<<<<<< HEAD @. as[1:end-1] = (k_inner / dx / dxps)*I_f @. ap[1:end-1] += as[1:end-1] @. ap[2:end] += an[2:end] @. an[2:end] *= I_f +======= + @. as[1:end-1] = (k_inner / dx / dxps) + @. as[isthawed] = 0. + @. ap[1:end-1] += as[1:end-1] + @. ap[2:end] += an[2:end] + @. an[isthawed] *= 0. +>>>>>>> have no idea return nothing end From b33c3245fcd4cdf29019eb9a96dc0933ca118bcc Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 2 Oct 2023 13:31:33 +0200 Subject: [PATCH 09/17] Update Lake impl to be compatible with master --- examples/heat_freeW_lake_lite_implicit.jl | 22 +++++++++++----------- src/Physics/Lakes/Lakes.jl | 17 +++++++++-------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index fa43ec13..0ad263d5 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -9,11 +9,11 @@ forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_t tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); # snowdepth = TimeSeriesForcing(ustrip.(forcings.data.Dsn), forcings.timestamps, :Dsn); soilprofile = SoilProfile( - 0.0u"m" => HomogeneousMixture(por=0.80,sat=0.9,org=0.75), #(θwi=0.80,θm=0.05,θo=0.15,ϕ=0.80), - 0.1u"m" => HomogeneousMixture(por=0.80,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.80), - 0.4u"m" => HomogeneousMixture(por=0.55,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55), - 3.0u"m" => HomogeneousMixture(por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50), - 10.0u"m" => HomogeneousMixture(por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30), + 0.0u"m" => MineralOrganic(por=0.80,sat=0.9,org=0.75), #(θwi=0.80,θm=0.05,θo=0.15,ϕ=0.80), + 0.1u"m" => MineralOrganic(por=0.80,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.80), + 0.4u"m" => MineralOrganic(por=0.55,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55), + 3.0u"m" => MineralOrganic(por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50), + 10.0u"m" => MineralOrganic(por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30), ) tempprofile_linear = TemperatureProfile( 0.0u"m" => -30.0u"°C", @@ -30,12 +30,12 @@ initT = initializer(:T, tempprofile_linear) heatop = Heat.EnthalpyImplicit() strat = @Stratigraphy( z_top => Top(upperbc), - z_sub[1] => :lake => Lake(HeatBalance(heatop)), - #z_sub[1] => :topsoil1 => Soil(HeatBalance(heatop), para=soilprofile[1].value), - #z_sub[2] => :topsoil2 => Soil(HeatBalance(heatop), para=soilprofile[2].value), - #z_sub[3] => :sediment1 => Soil(HeatBalance(heatop), para=soilprofile[3].value), - #z_sub[4] => :sediment2 => Soil(HeatBalance(heatop), para=soilprofile[4].value), - #z_sub[5] => :sediment3 => Soil(HeatBalance(heatop), para=soilprofile[5].value), + z_sub[1] => Lake(HeatBalance(heatop)), + #z_sub[1] => Ground(HeatBalance(heatop), para=soilprofile[1].value), + #z_sub[2] => Ground(HeatBalance(heatop), para=soilprofile[2].value), + #z_sub[3] => Ground(HeatBalance(heatop), para=soilprofile[3].value), + #z_sub[4] => Ground(HeatBalance(heatop), para=soilprofile[4].value), + #z_sub[5] => Ground(HeatBalance(heatop), para=soilprofile[5].value), z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2")) ); @info "Building tile" diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index dae9156a..6de935cc 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -8,19 +8,20 @@ export Lake abstract type LakeParameterization <: CryoGrid.Parameterization end -struct SimpleLakeScheme <: LakeParameterization end +Base.@kwdef struct SimpleLakeScheme{Thp} <: LakeParameterization + heat::Thp = ThermalProperties() +end -Base.@kwdef struct Lake{Tpara<:LakeParameterization,Tsp,Tproc,Tprop} <: CryoGrid.SubSurface{Tproc} +Base.@kwdef struct Lake{Tpara<:LakeParameterization,Theat<:HeatBalance,Taux} <: CryoGrid.SubSurface para::Tpara = SimpleLakeScheme() - prop::Tprop = CryoGrid.ThermalProperties() - sp::Tsp = nothing - proc::Tproc + heat::Theat = HeatBalance() + aux::Taux = nothing end Lake(proc::Tproc; kwargs...) where {Tproc} = Lake(;proc, kwargs...) # Material properties -Heat.thermalproperties(lake::Lake) = lake.prop +Heat.thermalproperties(lake::Lake) = lake.para.heat Hydrology.watercontent(::Lake, state) = 1.0 CryoGrid.volumetricfractions(::Lake, state, i) = (state.θw[i], 1 - state.θw[i], 0.0) @@ -93,7 +94,7 @@ end CryoGrid.computefluxes!(::Lake, ::HeatBalanceImplicit, state) = nothing -function CryoGrid.computefluxes!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.Enthalpy}, state) +function CryoGrid.computefluxes!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.EnthalpyBased}, 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 @@ -118,7 +119,7 @@ function CryoGrid.interact!( stop, slake ) - @setscalar slake.T_ub = CryoGrid.boundaryvalue(bc, top, heat, lake, stop, slake) + @setscalar slake.T_ub = stop.T_ub # boundary flux slake.jH[1] += CryoGrid.boundaryflux(bc, top, heat, lake, stop, slake) return nothing From a24b00bb09decb576226eca0036c68aacb8a0e49 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 2 Oct 2023 13:35:21 +0200 Subject: [PATCH 10/17] Fix out-of-date forcing setup in lake example --- examples/heat_freeW_lake_lite_implicit.jl | 9 +++------ src/Physics/Lakes/Lakes.jl | 2 -- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index 0ad263d5..3293e2d5 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -4,10 +4,7 @@ using Plots CryoGrid.debug(true) -forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term, :Tair => u"°C", :Ptot => u"mm"); -# use air temperature as upper boundary forcing -tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); -# snowdepth = TimeSeriesForcing(ustrip.(forcings.data.Dsn), forcings.timestamps, :Dsn); +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), #(θwi=0.80,θm=0.05,θo=0.15,ϕ=0.80), 0.1u"m" => MineralOrganic(por=0.80,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.80), @@ -24,13 +21,13 @@ modelgrid = Grid(vcat(-2.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_ z_top = -2.0u"m" z_sub = map(knot -> knot.depth, soilprofile) z_bot = modelgrid[end] -upperbc = TemperatureGradient(tair, NFactor()) +upperbc = TemperatureGradient(forcings.Tair, NFactor()) initT = initializer(:T, tempprofile_linear) @info "Building stratigraphy" heatop = Heat.EnthalpyImplicit() strat = @Stratigraphy( z_top => Top(upperbc), - z_sub[1] => Lake(HeatBalance(heatop)), + z_sub[1] => Lake(heat=HeatBalance(heatop)), #z_sub[1] => Ground(HeatBalance(heatop), para=soilprofile[1].value), #z_sub[2] => Ground(HeatBalance(heatop), para=soilprofile[2].value), #z_sub[3] => Ground(HeatBalance(heatop), para=soilprofile[3].value), diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 6de935cc..b2d1b9cf 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -18,8 +18,6 @@ Base.@kwdef struct Lake{Tpara<:LakeParameterization,Theat<:HeatBalance,Taux} <: aux::Taux = nothing end -Lake(proc::Tproc; kwargs...) where {Tproc} = Lake(;proc, kwargs...) - # Material properties Heat.thermalproperties(lake::Lake) = lake.para.heat Hydrology.watercontent(::Lake, state) = 1.0 From 3888266bfd3451a4d615bbc77f05a90dc9c191ca Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 2 Oct 2023 13:37:44 +0200 Subject: [PATCH 11/17] Fix another out of date reference --- src/Physics/Lakes/Lakes.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index b2d1b9cf..509c924b 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -20,7 +20,11 @@ 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) = ( @@ -105,7 +109,7 @@ function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceI jH_top = boundaryflux(bc, top, heat, lake, stop, slake) k = slake.k[1] slake.DT_bp[1] += jH_top / Δk - slake.DT_ap[1] += Heat._ap(CryoGrid.BoundaryStyle(bc), k, Δk) + slake.DT_ap[1] += Heat._ap(CryoGrid.BCKind(bc), k, Δk) return nothing end From d404d9fd5afe94f4cab443b52233852ccfa8a258 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 2 Oct 2023 14:11:33 +0200 Subject: [PATCH 12/17] Fix another bug in lake example --- examples/heat_freeW_lake_lite_implicit.jl | 28 +++++++++-------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index 3293e2d5..9f19a048 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -6,11 +6,11 @@ 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), #(θwi=0.80,θm=0.05,θo=0.15,ϕ=0.80), - 0.1u"m" => MineralOrganic(por=0.80,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.80), - 0.4u"m" => MineralOrganic(por=0.55,sat=0.9,org=0.25), #(θwi=0.80,θm=0.15,θo=0.05,ϕ=0.55), - 3.0u"m" => MineralOrganic(por=0.50,sat=1.0,org=0.0), #(θwi=0.50,θm=0.50,θo=0.0,ϕ=0.50), - 10.0u"m" => MineralOrganic(por=0.30,sat=1.0,org=0.0), #(θwi=0.30,θm=0.70,θo=0.0,ϕ=0.30), + 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( 0.0u"m" => -30.0u"°C", @@ -27,12 +27,8 @@ initT = initializer(:T, tempprofile_linear) heatop = Heat.EnthalpyImplicit() strat = @Stratigraphy( z_top => Top(upperbc), - z_sub[1] => Lake(heat=HeatBalance(heatop)), - #z_sub[1] => Ground(HeatBalance(heatop), para=soilprofile[1].value), - #z_sub[2] => Ground(HeatBalance(heatop), para=soilprofile[2].value), - #z_sub[3] => Ground(HeatBalance(heatop), para=soilprofile[3].value), - #z_sub[4] => Ground(HeatBalance(heatop), para=soilprofile[4].value), - #z_sub[5] => Ground(HeatBalance(heatop), para=soilprofile[5].value), + -2.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" @@ -42,17 +38,15 @@ tspan = (DateTime(2010,12,30), DateTime(2015,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,)) -lake_state = getstate(:lake, tile, u0, du0, prob.tspan[1]) -# debug one step -CryoGrid.debug(true) -tile(du0, u0, prob.p, prob.tspan[1]) +# set up integrator integrator = init(prob, LiteImplicitEuler(), dt=24*3600) -step!(integrator) +# debug one step +@run step!(integrator) @info "Running model" sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600) out = CryoGridOutput(sol) # Plot the results -zs = [5,10,15,20,25,30,40,50,100,500,1000,5000]u"cm" +zs = [-200,-180,-160,-120,-100,-50,-10,-5,-1,1,5,11]u"cm" cg = Plots.cgrad(:copper,rev=true); plot(out.T[Z(Near(zs))] |> ustrip, color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", leg=false, dpi=150) From 726e672d8e4adc2dd2b5d453f96e82d5fd54dae5 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Mon, 2 Oct 2023 14:11:56 +0200 Subject: [PATCH 13/17] Remove old example script --- examples/heat_freeW_samoylov_lakes.jl | 55 --------------------------- 1 file changed, 55 deletions(-) delete mode 100644 examples/heat_freeW_samoylov_lakes.jl diff --git a/examples/heat_freeW_samoylov_lakes.jl b/examples/heat_freeW_samoylov_lakes.jl deleted file mode 100644 index 04534f65..00000000 --- a/examples/heat_freeW_samoylov_lakes.jl +++ /dev/null @@ -1,55 +0,0 @@ -using CryoGrid -using Plots - -forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C"); -# define time span -tspan = (DateTime(2010,10,30),DateTime(2011,10,30)) -# use air temperature as upper boundary forcing; -tair = TimeSeriesForcing(forcings.data.Tair, forcings.timestamps, :Tair); -T0 = values(tair[tspan[1]])[1] -tempprofile = TemperatureProfile( - 0.0u"m" => T0, - 1.0u"m" => -8.0u"°C", - 20.0u"m" => -10u"°C", - 1000.0u"m" => 1.0u"°C" -) -# soil profile: depth => (excess ice, natural porosity, saturation, organic fraction) -soilprofile = SoilProfile( - 0.0u"m" => HomogeneousMixture(por=0.80,sat=1.0,org=0.75), - 0.1u"m" => HomogeneousMixture(por=0.80,sat=1.0,org=0.25), - 0.4u"m" => HomogeneousMixture(por=0.55,sat=1.0,org=0.25), - 3.0u"m" => HomogeneousMixture(por=0.50,sat=1.0,org=0.0), - 10.0u"m" => HomogeneousMixture(por=0.30,sat=1.0,org=0.0), - 100.0u"m" => HomogeneousMixture(por=0.10,sat=0.1,org=0.0), -); -initT = initializer(:T, tempprofile) -# Heat diffusion with temperature as prognostic state variable -# Construct soil layers -soil_layers = map(enumerate(soilprofile)) do (i, (depth, para)) - name = Symbol(:soil, i) - depth => name => Soil(HeatBalance(), para=para) -end -# Build stratigraphy: -# @Stratigraphy macro lets us list multiple subsurface layers -strat = @Stratigraphy( - -2.0*u"m" => Top(TemperatureGradient(tair)), - -2.0u"m" => :lake => Lake(HeatBalance()), - soil_layers..., - 998.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2")), -); -# shift grid up by 2 m -modelgrid = Grid(CryoGrid.Presets.DefaultGrid_5cm .- 2.0u"m") -tile = Tile(strat, modelgrid, initT) -u0, du0 = initialcondition!(tile, tspan) -# construct CryoGridProblem with tile, initial condition, and timespan; -# we disable the default timestep limiter since we will use an adaptive solver. -prob = CryoGridProblem(tile, u0, tspan, savevars=(:T,)) -# test step function (good to do debugging here) -tile(du0, u0, prob.p, prob.tspan[1]) -@info "Running model" -out = @time solve(prob, Euler(), dt=300.0, saveat=24*3600.0, progress=true) |> CryoGridOutput; -# Plot it! -zs = [1,10,20,30,50,100,200,500,1000]u"cm" -cg = Plots.cgrad(:copper,rev=true); -plot(out.H[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Enthalpy", leg=false, dpi=150) -plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150) From 05a0f21bc80e435b0d409f8d3432b62d6f9582f8 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 3 Oct 2023 10:15:16 +0200 Subject: [PATCH 14/17] Fix issues with lake impl --- examples/heat_freeW_lake_lite_implicit.jl | 16 ++-- src/Physics/Heat/heat_implicit.jl | 54 +++++------ src/Physics/Lakes/Lakes.jl | 107 ++++++++++++---------- src/Solvers/LiteImplicit/cglite_step.jl | 2 +- 4 files changed, 99 insertions(+), 80 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index 9f19a048..7e4ca0c7 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -13,7 +13,8 @@ soilprofile = SoilProfile( 10.0u"m" => MineralOrganic(por=0.30,sat=1.0,org=0.0), ) tempprofile_linear = TemperatureProfile( - 0.0u"m" => -30.0u"°C", + -2.0u"m" => 0.0u"°C", + 0.0u"m" => 0.0u"°C", 10.0u"m" => -10.0u"°C", 1000.0u"m" => 10.2u"°C" ) @@ -21,7 +22,7 @@ modelgrid = Grid(vcat(-2.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_ z_top = -2.0u"m" z_sub = map(knot -> knot.depth, soilprofile) z_bot = modelgrid[end] -upperbc = TemperatureGradient(forcings.Tair, NFactor()) +upperbc = TemperatureGradient(forcings.Tair, NFactor(nf=0.5)) initT = initializer(:T, tempprofile_linear) @info "Building stratigraphy" heatop = Heat.EnthalpyImplicit() @@ -34,19 +35,22 @@ strat = @Stratigraphy( @info "Building tile" tile = @time Tile(strat, modelgrid, initT) # define time span, 5 years -tspan = (DateTime(2010,12,30), DateTime(2015,12,30)) +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 -@run step!(integrator) +step!(integrator) @info "Running model" sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600) out = CryoGridOutput(sol) # Plot the results -zs = [-200,-180,-160,-120,-100,-50,-10,-5,-1,1,5,11]u"cm" +zs = [-200,-150.0,-100.0,-50.0,-1.0,3.0]u"cm" cg = Plots.cgrad(:copper,rev=true); -plot(out.T[Z(Near(zs))] |> ustrip, color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", title="", leg=false, dpi=150) +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)) + +Plots.plotly() diff --git a/src/Physics/Heat/heat_implicit.jl b/src/Physics/Heat/heat_implicit.jl index ec60894e..d8a74795 100644 --- a/src/Physics/Heat/heat_implicit.jl +++ b/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,14 +71,13 @@ 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) @@ -84,7 +85,7 @@ function CryoGrid.interact!(sub::SubSurface, heat::HeatBalanceImplicit, bot::Bot 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 diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 509c924b..816ca979 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -18,6 +18,19 @@ Base.@kwdef struct Lake{Tpara<:LakeParameterization,Theat<:HeatBalance,Taux} <: 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)) + break + end + end + end + return ubc_idx +end + # Material properties Heat.thermalproperties(lake::Lake) = lake.para.heat @@ -32,6 +45,7 @@ CryoGrid.variables(lake::Lake, heat::HeatBalance) = ( # 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.initialcondition!(lake::Lake, heat::HeatBalance, state) @@ -45,7 +59,50 @@ function CryoGrid.initialcondition!(lake::Lake, heat::HeatBalance, state) end end -function CryoGrid.updatestate!(sub::Lake, heat::HeatBalance, state) +function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake) + jH_top = boundaryflux(bc, top, heat, lake, stop, slake) + T_ub = slake.T_ub[1] = getscalar(stop.T_ub) + ubc_idx = get_upper_boundary_index(T_ub, slake.θw) + # 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) + # first lake cell + bp[1] += jH_top / dxp[1] + ap[1] += Heat.apbc(CryoGrid.BCKind(bc), k[1], dxp[1]) + # deeper lake cells + @inbounds for i in 2:ubc_idx + bp[i] += jH_top / dxp[i] + ap[i] -= as[i] + as[i] = zero(eltype(as)) + an[i] = zero(eltype(an)) + ap[i] += Heat.apbc(CryoGrid.BCKind(bc), k[i], dxp[i], dx[i-1]) + 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) + # 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 + ubc_idx = get_upper_boundary_index(slake.T_ub[1], slake.θw) + 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) @@ -61,42 +118,7 @@ function CryoGrid.updatestate!(sub::Lake, heat::HeatBalance, state) return nothing end -function CryoGrid.updatestate!( - sub::Lake, - heat::HeatBalanceImplicit, - state -) - Heat.resetfluxes!(sub, heat, state) - # Evaluate freeze/thaw processes - Heat.freezethaw!(sub, heat, state) - # Update thermal conductivity - Heat.thermalconductivity!(sub, heat, state) - isthawed = true - @inbounds for i in eachindex(state.θw) - isthawed = isthawed && state.θw[i] ≈ 1.0 - end - I_f = 1 - Float64(isthawed) - # Compute diffusion coefficients - an = state.DT_an - as = state.DT_as - ap = state.DT_ap - k = state.k - dx = Δ(cells(state.grid)) - dxp = Δ(state.grid) - k_inner = @view k[2:end-1] - dxpn = @view dxp[1:end-1] - dxps = @view dxp[2:end] - @. an[2:end] = k_inner / dx / dxpn - @. as[1:end-1] = (k_inner / dx / dxps)*I_f - @. ap[1:end-1] += as[1:end-1] - @. ap[2:end] += an[2:end] - @. an[2:end] *= I_f - return nothing -end - -CryoGrid.computefluxes!(::Lake, ::HeatBalanceImplicit, state) = nothing - -function CryoGrid.computefluxes!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.EnthalpyBased}, state) +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 @@ -104,20 +126,12 @@ function CryoGrid.computefluxes!(::Lake, ::HeatBalance{<:FreezeCurve,<:Heat.Enth return nothing end -function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake) - Δk = CryoGrid.thickness(lake, slake, first) - jH_top = boundaryflux(bc, top, heat, lake, stop, slake) - k = slake.k[1] - slake.DT_bp[1] += jH_top / Δk - slake.DT_ap[1] += Heat._ap(CryoGrid.BCKind(bc), k, Δk) - return nothing -end function CryoGrid.interact!( top::Top, bc::HeatBC, lake::Lake, - heat::HeatBalance, + heat::HeatBalance{FreeWater,<:Heat.MOLEnthalpy}, stop, slake ) @@ -127,5 +141,4 @@ function CryoGrid.interact!( return nothing end - end \ No newline at end of file diff --git a/src/Solvers/LiteImplicit/cglite_step.jl b/src/Solvers/LiteImplicit/cglite_step.jl index 4119861c..22e8e857 100644 --- a/src/Solvers/LiteImplicit/cglite_step.jl +++ b/src/Solvers/LiteImplicit/cglite_step.jl @@ -9,7 +9,7 @@ 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) integrator.t = t From e9a5d045f4baaf9aa4a4312e1793660828bc214b Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 3 Oct 2023 10:17:28 +0200 Subject: [PATCH 15/17] Split Lakes module files --- src/Physics/Lakes/Lakes.jl | 136 +------------------------------ src/Physics/Lakes/lake_simple.jl | 134 ++++++++++++++++++++++++++++++ 2 files changed, 135 insertions(+), 135 deletions(-) create mode 100644 src/Physics/Lakes/lake_simple.jl diff --git a/src/Physics/Lakes/Lakes.jl b/src/Physics/Lakes/Lakes.jl index 816ca979..ebb59181 100644 --- a/src/Physics/Lakes/Lakes.jl +++ b/src/Physics/Lakes/Lakes.jl @@ -5,140 +5,6 @@ using CryoGrid.Utils using CryoGrid.Numerics export Lake - -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)) - break - end - end - end - return ubc_idx -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.initialcondition!(lake::Lake, heat::HeatBalance, state) - L = heat.prop.L - # initialize liquid water content based on temperature - @inbounds for i in 1:length(state.T) - θwi = Hydrology.watercontent(lake, state, i) - state.θw[i] = ifelse(state.T[i] > 0.0, θwi, 0.0) - state.C[i] = heatcapacity(lake, heat, volumetricfractions(lake, state, i)...) - state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i]) - end -end - -function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake) - jH_top = boundaryflux(bc, top, heat, lake, stop, slake) - T_ub = slake.T_ub[1] = getscalar(stop.T_ub) - ubc_idx = get_upper_boundary_index(T_ub, slake.θw) - # 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) - # first lake cell - bp[1] += jH_top / dxp[1] - ap[1] += Heat.apbc(CryoGrid.BCKind(bc), k[1], dxp[1]) - # deeper lake cells - @inbounds for i in 2:ubc_idx - bp[i] += jH_top / dxp[i] - ap[i] -= as[i] - as[i] = zero(eltype(as)) - an[i] = zero(eltype(an)) - ap[i] += Heat.apbc(CryoGrid.BCKind(bc), k[i], dxp[i], dx[i-1]) - 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) - # 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 - ubc_idx = get_upper_boundary_index(slake.T_ub[1], slake.θw) - 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 +include("lake_simple.jl") end \ No newline at end of file diff --git a/src/Physics/Lakes/lake_simple.jl b/src/Physics/Lakes/lake_simple.jl new file mode 100644 index 00000000..4626f8c1 --- /dev/null +++ b/src/Physics/Lakes/lake_simple.jl @@ -0,0 +1,134 @@ +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)) + break + end + end + end + return ubc_idx +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.initialcondition!(lake::Lake, heat::HeatBalance, state) + L = heat.prop.L + # initialize liquid water content based on temperature + @inbounds for i in 1:length(state.T) + θwi = Hydrology.watercontent(lake, state, i) + state.θw[i] = ifelse(state.T[i] > 0.0, θwi, 0.0) + state.C[i] = heatcapacity(lake, heat, volumetricfractions(lake, state, i)...) + state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i]) + end +end + +function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake) + jH_top = boundaryflux(bc, top, heat, lake, stop, slake) + T_ub = slake.T_ub[1] = getscalar(stop.T_ub) + ubc_idx = get_upper_boundary_index(T_ub, slake.θw) + # 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) + # first lake cell + bp[1] += jH_top / dxp[1] + ap[1] += Heat.apbc(CryoGrid.BCKind(bc), k[1], dxp[1]) + # deeper lake cells + @inbounds for i in 2:ubc_idx + bp[i] += jH_top / dxp[i] + ap[i] -= as[i] + as[i] = zero(eltype(as)) + an[i] = zero(eltype(an)) + ap[i] += Heat.apbc(CryoGrid.BCKind(bc), k[i], dxp[i], dx[i-1]) + 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) + # 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 + ubc_idx = get_upper_boundary_index(slake.T_ub[1], slake.θw) + 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 \ No newline at end of file From 6d4d96dfc1b34ed6146bfeab29c449b1138050b6 Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 3 Oct 2023 14:47:56 +0200 Subject: [PATCH 16/17] More fixes for simple lake scheme --- examples/heat_freeW_lake_lite_implicit.jl | 27 +++++++++++++++++------ src/Physics/Lakes/lake_simple.jl | 25 +++++++++++---------- 2 files changed, 33 insertions(+), 19 deletions(-) diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index 7e4ca0c7..f9b6dcca 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -2,6 +2,8 @@ using CryoGrid using CryoGrid.LiteImplicit using Plots +Plots.plotly() + CryoGrid.debug(true) forcings = loadforcings(CryoGrid.Presets.Forcings.Samoylov_ERA_MkL3_CCSM4_long_term); @@ -18,8 +20,8 @@ tempprofile_linear = TemperatureProfile( 10.0u"m" => -10.0u"°C", 1000.0u"m" => 10.2u"°C" ) -modelgrid = Grid(vcat(-2.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm)) -z_top = -2.0u"m" +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)) @@ -28,7 +30,7 @@ initT = initializer(:T, tempprofile_linear) heatop = Heat.EnthalpyImplicit() strat = @Stratigraphy( z_top => Top(upperbc), - -2.0u"m" => Lake(heat=HeatBalance(heatop)), + -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")) ); @@ -43,14 +45,25 @@ prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600.0, savevars=(:θw,:T,)) 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) out = CryoGridOutput(sol) # Plot the results -zs = [-200,-150.0,-100.0,-50.0,-1.0,3.0]u"cm" +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)) - -Plots.plotly() +plot!(convert_t.(dims(out.T, Ti)), t -> forcings.Tair.(t), c=:blue, linestyle=:dash) diff --git a/src/Physics/Lakes/lake_simple.jl b/src/Physics/Lakes/lake_simple.jl index 4626f8c1..5002d5ff 100644 --- a/src/Physics/Lakes/lake_simple.jl +++ b/src/Physics/Lakes/lake_simple.jl @@ -16,11 +16,11 @@ function get_upper_boundary_index(T_ub, θw) @inbounds for i in eachindex(θw) ubc_idx = i if θw[i] < one(eltype(θw)) - break + return ubc_idx end end end - return ubc_idx + return ubc_idx+1 end # Material properties @@ -52,7 +52,6 @@ function CryoGrid.initialcondition!(lake::Lake, heat::HeatBalance, state) end function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceImplicit, stop, slake) - jH_top = boundaryflux(bc, top, heat, lake, stop, slake) T_ub = slake.T_ub[1] = getscalar(stop.T_ub) ubc_idx = get_upper_boundary_index(T_ub, slake.θw) # get variables @@ -63,16 +62,18 @@ function CryoGrid.interact!(top::Top, bc::HeatBC, lake::Lake, heat::HeatBalanceI k = slake.k dx = Δ(cells(slake.grid)) dxp = Δ(slake.grid) - # first lake cell - bp[1] += jH_top / dxp[1] - ap[1] += Heat.apbc(CryoGrid.BCKind(bc), k[1], dxp[1]) - # deeper lake cells - @inbounds for i in 2:ubc_idx - bp[i] += jH_top / dxp[i] - ap[i] -= as[i] + # 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)) - ap[i] += Heat.apbc(CryoGrid.BCKind(bc), k[i], dxp[i], dx[i-1]) end return nothing end @@ -89,7 +90,7 @@ function CryoGrid.interact!(lake::Lake, ::HeatBalanceImplicit, sub::SubSurface, harmonicmean(k₁, k₂, Δ₁, Δ₂) end ubc_idx = get_upper_boundary_index(slake.T_ub[1], slake.θw) - slake.DT_ap[end] += slake.DT_as[end] = (ubc_idx < length(slake.θw))*k / Δz / Δk₁ + 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 From ee99fc70fafc488f5f369b554d04126d816508bf Mon Sep 17 00:00:00 2001 From: Brian Groenke Date: Tue, 3 Oct 2023 16:08:56 +0200 Subject: [PATCH 17/17] Add diagnosticstep! for discrete time state updates --- docs/src/manual/architecture.md | 4 ++-- examples/heat_freeW_lake_lite_implicit.jl | 2 +- src/CryoGrid.jl | 2 +- src/Physics/Lakes/lake_simple.jl | 22 +++++++++++----------- src/Solvers/LiteImplicit/cglite_step.jl | 3 +++ src/Tiles/stratigraphy.jl | 6 ++++++ src/Tiles/tile.jl | 16 +++++++++------- src/methods.jl | 8 ++++++++ 8 files changed, 41 insertions(+), 22 deletions(-) diff --git a/docs/src/manual/architecture.md b/docs/src/manual/architecture.md index 47299de8..dfaba343 100644 --- a/docs/src/manual/architecture.md +++ b/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: diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl index f9b6dcca..8f65e3b3 100644 --- a/examples/heat_freeW_lake_lite_implicit.jl +++ b/examples/heat_freeW_lake_lite_implicit.jl @@ -59,7 +59,7 @@ state.top.T_ub state.lake.T @info "Running model" -sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600) +sol = @time solve(prob, LiteImplicitEuler(), dt=24*3600.0) out = CryoGridOutput(sol) # Plot the results diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl index a3d80b85..16a747b2 100755 --- a/src/CryoGrid.jl +++ b/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") diff --git a/src/Physics/Lakes/lake_simple.jl b/src/Physics/Lakes/lake_simple.jl index 5002d5ff..59bc825e 100644 --- a/src/Physics/Lakes/lake_simple.jl +++ b/src/Physics/Lakes/lake_simple.jl @@ -40,20 +40,20 @@ CryoGrid.variables(lake::Lake, heat::HeatBalance) = ( Diagnostic(:ubc_idx, Scalar, NoUnits, Int), ) -function CryoGrid.initialcondition!(lake::Lake, heat::HeatBalance, state) - L = heat.prop.L - # initialize liquid water content based on temperature - @inbounds for i in 1:length(state.T) - θwi = Hydrology.watercontent(lake, state, i) - state.θw[i] = ifelse(state.T[i] > 0.0, θwi, 0.0) - state.C[i] = heatcapacity(lake, heat, volumetricfractions(lake, state, i)...) - state.H[i] = enthalpy(state.T[i], state.C[i], L, state.θw[i]) - end +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 = get_upper_boundary_index(T_ub, slake.θw) + ubc_idx = Int(getscalar(slake.ubc_idx)) # get variables an = slake.DT_an as = slake.DT_as @@ -81,6 +81,7 @@ function CryoGrid.interact!(lake::Lake, ::HeatBalanceImplicit, sub::SubSurface, Δ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], @@ -89,7 +90,6 @@ function CryoGrid.interact!(lake::Lake, ::HeatBalanceImplicit, sub::SubSurface, Δ₂ = Δk₂[1]; harmonicmean(k₁, k₂, Δ₁, Δ₂) end - ubc_idx = get_upper_boundary_index(slake.T_ub[1], slake.θw) 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 diff --git a/src/Solvers/LiteImplicit/cglite_step.jl b/src/Solvers/LiteImplicit/cglite_step.jl index 22e8e857..2e74a256 100644 --- a/src/Solvers/LiteImplicit/cglite_step.jl +++ b/src/Solvers/LiteImplicit/cglite_step.jl @@ -12,6 +12,9 @@ function CryoGrid.perform_step!(integrator::CGLiteIntegrator) # 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 diff --git a/src/Tiles/stratigraphy.jl b/src/Tiles/stratigraphy.jl index 02fa921a..d62f0221 100644 --- a/src/Tiles/stratigraphy.jl +++ b/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) diff --git a/src/Tiles/tile.jl b/src/Tiles/tile.jl index 1b28d6c1..44dec465 100644 --- a/src/Tiles/tile.jl +++ b/src/Tiles/tile.jl @@ -15,9 +15,9 @@ abstract type AbstractTile{iip} end (tile::AbstractTile{true})(du,u,p,t,dt=1.0) (tile::AbstractTile{false})(u,p,t,dt=1.0) -Invokes the corresponding `evaluate!` function to compute the time derivative du/dt. +Invokes the corresponding `prognostic!` function to compute the time derivative du/dt. """ -(tile::AbstractTile{true})(du, u, p, t, dt=1.0) = evaluate!(tile,du,u,p,t,dt) +(tile::AbstractTile{true})(du, u, p, t, dt=1.0) = prognostic!(tile,du,u,p,t,dt) (tile::AbstractTile{false})(u, p, t, dt=1.0) = evaluate(tile,u,p,t,dt) """ @@ -28,11 +28,11 @@ Returns true if `tile` uses in-place mode, false if out-of-place. isinplace(::AbstractTile{iip}) where {iip} = iip """ - evaluate!(::T,du,u,p,t) where {T<:AbstractTile} + prognostic!(::T,du,u,p,t) where {T<:AbstractTile} In-place update function for tile `T`. Computes du/dt and stores the result in `du`. """ -evaluate!(::T, du, u, p, t, dt=1.0) where {T<:AbstractTile} = error("no implementation of in-place evaluate! for $T") +prognostic!(::T, du, u, p, t, dt=1.0) where {T<:AbstractTile} = error("no implementation of in-place prognostic! for $T") """ evaluate(::T,u,p,t) where {T<:AbstractTile} @@ -143,7 +143,7 @@ function Tiles.Tile(integrator::SciMLBase.DEIntegrator) end """ - evaluate!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents} + prognostic!(_tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, p, t) where {TStrat,TGrid,TStates,TInits,TEvents} Time derivative step function (i.e. du/dt) for any arbitrary Tile. Specialized code is generated and compiled on the fly via the @generated macro to ensure type stability. The generated code updates each layer in the stratigraphy @@ -155,7 +155,7 @@ interact!(layer[i], ..., layer[i+1], ...) computefluxes!(layer[i], ...) ``` """ -function evaluate!( +function prognostic!( _tile::Tile{TStrat,TGrid,TStates,TInits,TEvents,true}, _du, _u, @@ -181,6 +181,8 @@ function evaluate!( return nothing end +CryoGrid.diagnosticstep!(tile::Tile, state) = diagnosticstep!(tile.strat, state) + """ timestep(_tile::Tile, _du, _u, p, t) @@ -222,7 +224,7 @@ function CryoGrid.initialcondition!(tile::Tile, tspan::NTuple{2,Float64}, p=noth state = TileState(tile.state, zs, u, du, t0, 1.0, Val{iip}()) CryoGrid.initialcondition!(strat, state, tile.inits) # evaluate initial time derivative - evaluate!(tile, du, u, p, t0) + prognostic!(tile, du, u, p, t0) return u, du end diff --git a/src/methods.jl b/src/methods.jl index 708e3245..fff4659f 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -127,6 +127,14 @@ actual chosen timestep will depend on the integrator being used and other user c timestep(layer::Layer, state) = timestep(layer, processes(layer), state) timestep(::Layer, ::Process, state) = Inf +""" + diagnosticstep!(layer::Layer, state) + +Optionally performs discontinuous/discrete-time updates to the layer state. Should return `true` +if the prognostic state was modified and `false` otherwise. Defaults to returning `false`. +""" +diagnosticstep!(layer::Layer, state) = false + """ parameterize(x::T) where {T} parameterize(x::Unitful.AbstractQuantity; props...)