Add simple two-stage nfactor impl

src/Processes/HeatConduction/heat_bc.jl

-struct TemperatureGradient{F} <: BoundaryProcess
+# Boundary condition type aliases
+const ConstantTemp = Constant{Dirichlet,Float"°C"}
+ConstantTemp(value::UFloat"K") = Constant{Dirichlet}(dustrip(u"°C", value))
+ConstantTemp(value::UFloat"°C") = Constant{Dirichlet}(dustrip(value))
+const GeothermalHeatFlux = Constant{Neumann,Float"J/s/m^2"}
+GeothermalHeatFlux(value::UFloat"J/s/m^2"=0.053xu"J/s/m^2") = Constant{Neumann}(dustrip(value))
+
+abstract type TempGradTransform end
+struct NoTransform <: TempGradTransform end
+struct StationaryNFactor{nfw,nfs} <: TempGradTransform end
+struct TwoStageNFactor{nfw1,nfw2,nfs1,nfs2,nftc} <: TempGradTransform end
+
+struct TemperatureGradient{T,F} <: BoundaryProcess
     forcing::F
-    TemperatureGradient(forcing::Forcing{Float"°C"}) = new{typeof(forcing)}(forcing)
+    transform::T
+    TemperatureGradient(forcing::Forcing{Float"°C"}, transform::T=NoTransform()) where {T<:TempGradTransform} = new{T,typeof(forcing)}(forcing, transform)
 end
 
-@inline (bc::TemperatureGradient)(t) = bc.forcing(t)
-@inline (bc::TemperatureGradient)(l1,l2,p2,s1,s2) = bc(s1.t)
-
-BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
+NFactor(type::Symbol=:stationary, name::Symbol=:nf) = NFactor(Val{type}(), name)
+NFactor(::Val{:stationary}, name::Symbol=:nf) = StationaryNFactor{Symbol(name,:w),Symbol(name,:s)}()
+NFactor(::Val{:twostage}, name::Symbol=:nf) = TwoStageNFactor{Symbol(name,:w1),Symbol(name,:w2),Symbol(name,:s1),Symbol(name,:s2),Symbol(name,:tc)}()
 
-struct NFactor{F,nf1,nf2,nt} <: BoundaryProcess
-    tgrad::TemperatureGradient{F}
-    factor1::Float64
-    factor2::Float64
-    threshold::Float64
-    NFactor(tgrad::TemperatureGradient{F}, factor1::Float64=0.5, factor2::Float64=1.0, threshold::Float64=0.0, name::Symbol=:n) where {F} =
-        new{F,Symbol(name,:_factor1),Symbol(name,:_factor2),Symbol(name,:_thresh)}(tgrad, factor1, factor2, threshold)
+@inline (bc::TemperatureGradient)(l1,l2,p2,s1,s2) where {F} = bc(s1.t)
+@inline function (bc::TemperatureGradient{StationaryNFactor{nfw,nfs}})(l1,l2,p2,s1,s2) where {nfw,nfs}
+    let nf_winter = s1.params[nfw] |> getscalar,
+        nf_summer = s1.params[nfs] |> getscalar,
+        Tair = bc.forcing(s1.t),
+        factor_eff = (Tair <= zero(Tair))*nf_winter + (Tair > zero(Tair))*nf_summer;
+        Tbc = factor_eff*Tair
+    end
 end
-
-@inline function (bc::NFactor{F,nf1,nf2,nt})(l1,l2,p2,s1,s2) where {F,nf1,nf2,nt}
-    let factor1 = s1.params[nf1] |> getscalar,
-        factor2 = s1.params[nf2] |> getscalar,
-        thresh = s1.params[nt] |> getscalar,
-        Tair = bc.tgrad(l1,l2,p2,s1,s2),
-        check = Tair - thresh,
-        factor_eff = (check <= zero(check))*factor1 + (check > zero(check))*factor2;
-        factor_eff*Tair
+@inline function (bc::TemperatureGradient{TwoStageNFactor{nfw1,nfw2,nfs1,nfs2,nftc}})(l1,l2,p2,s1,s2) where {nfw1,nfw2,nfs1,nfs2,nftc}
+    let nf_winter1 = s1.params[nfw1] |> getscalar,
+        nf_winter2 = s1.params[nfw2] |> getscalar,
+        nf_summer1 = s1.params[nfs1] |> getscalar,
+        nf_summer2 = s1.params[nfs2] |> getscalar,
+        nf_tchange = s1.params[nftc] |> getscalar,
+        tcheck = s1.t >= nf_tchange,
+        Tair = bc.forcing(s1.t),
+        factor_eff = (Tair <= zero(Tair))*(1-tcheck)*nf_winter1 +
+            (Tair > zero(Tair))*(1-tcheck)*nf_summer1 +
+            (Tair <= zero(Tair))*tcheck*nf_winter2 +
+            (Tair > zero(Tair))*tcheck*nf_summer2
+        Tbc = factor_eff*Tair
     end
 end
 
-variables(::Top, bc::NFactor{F,nf1,nf2,nt}) where {F,nf1,nf2,nt} = (Parameter(nf1, bc.factor1, 0..1), Parameter(nf2, bc.factor2, 0..1), Parameter(nt, bc.threshold))
-BoundaryStyle(::Type{<:NFactor}) = Dirichlet()
+variables(::Top, bc::TemperatureGradient{StationaryNFactor{nfw,nfs}}) where {nfw,nfs} = (Parameter(nfw, 1.0, 0..1), Parameter(nfs, 1.0, 0..1))
+variables(::Top, bc::TemperatureGradient{TwoStageNFactor{nfw1,nfw2,nfs1,nfs2,nftc}}) where {nfw1,nfw2,nfs1,nfs2,nftc} = (
+    (Parameter(name, 1.0, 0..1) for name in (nfw1,nfw2,nfs1,nfs2))...,
+    Parameter(nftc, 0.0, 0..Inf)
+)
 
-# Boundary condition type aliases
-const ConstantTemp = Constant{Dirichlet,Float"°C"}
-ConstantTemp(value::UFloat"K") = Constant{Dirichlet}(dustrip(u"°C", value))
-ConstantTemp(value::UFloat"°C") = Constant{Dirichlet}(dustrip(value))
-const GeothermalHeatFlux = Constant{Neumann,Float"J/s/m^2"}
-GeothermalHeatFlux(value::UFloat"J/s/m^2"=0.053xu"J/s/m^2") = Constant{Neumann}(dustrip(value))
+BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()

test/Processes/HeatConduction/heat_conduction_tests.jl

@@ -65,22 +65,28 @@ end
 	@testset "n-factors" begin
 		ts = DateTime(2010,1,1):Hour(1):DateTime(2010,1,1,4)
 		forcing = TimeSeriesForcing([1.0,0.5,-0.5,-1.0,0.1], ts, :Tair)
-		nfactor = NFactor(TemperatureGradient(forcing), 0.5, 1.0, 0.0, :n)
-		vars = variables(Top(), nfactor)
-		@test length(vars) == 3
-		@test getscalar(vars[1].default_value) == 0.5
+		tgrad = TemperatureGradient(forcing, NFactor(:stationary, :nf))
+		vars = variables(Top(), tgrad)
+		@test length(vars) == 2
+		@test getscalar(vars[1].default_value) == 1.0
 		@test getscalar(vars[2].default_value) == 1.0
-		@test getscalar(vars[3].default_value) == 0.0
 		sub = TestGroundLayer()
 		heat = Heat{:H}()
-		# standard zero threshold
-		f1(t) = nfactor(Top(),sub,heat,(t=t, params=(n_factor1=0.5, n_factor2=1.0, n_thresh=0.0)), (t=t,))
+		# stationary case
+		f1(t) = tgrad(Top(),sub,heat,(t=t, params=(nfw=0.5, nfs=1.0)), (t=t,))
 		Tres = f1.(Dates.datetime2epochms.(ts)./1000.0)
 		@test all(Tres .≈ [1.0,0.5,-0.25,-0.5,0.1])
-		# test with non-zero threshold
-		f2(t) = nfactor(Top(),sub,heat,(t=t, params=(n_factor1=0.5, n_factor2=1.0, n_thresh=0.5)), (t=t,))
-		Tres = f2.(Dates.datetime2epochms.(ts)./1000.0)
-		@test all(Tres .≈ [1.0,0.25,-0.25,-0.5,0.05])
+		# two-stage non-stationary case
+		tgrad = TemperatureGradient(forcing, NFactor(:twostage, :nf))
+		vars = variables(Top(), tgrad)
+		@test length(vars) == 5
+		@test all([getscalar(vars[i].default_value) == 1.0 for i in 1:4])
+		@test getscalar(vars[5].default_value) == 0.0
+		tchange = Dates.datetime2epochms(ts[end-1])./1000
+		f2(t) = tgrad(Top(),sub,heat,(t=t, params=(nfw1=0.5, nfw2=0.2, nfs1=1.0, nfs2=0.9, nftc=tchange)), (t=t,))
+		Tres = f2.(Dates.datetime2epochms.(ts)./1000)
+		println(Tres)
+		@test all(Tres .≈ [1.0,0.5,-0.25,-0.2,0.09])
 	end
 end
 @testset "Fourier solution" begin