From 0c26658696447e51516f7befb0b20631fc74ca0c Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Wed, 15 Dec 2021 22:21:40 +0100
Subject: [PATCH] Optimize freeW freeze curve and related functions

---
 src/Physics/HeatConduction/HeatConduction.jl |  1 +
 src/Physics/HeatConduction/heat.jl           | 12 +++++++++---
 src/Physics/HeatConduction/heat_bc.jl        |  2 +-
 src/Physics/HeatConduction/soil/soilheat.jl  | 12 ++++++++----
 4 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index 1db0bf5e..067fd43a 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -11,6 +11,7 @@ using CryoGrid.Numerics: nonlineardiffusion!, harmonicmean!, harmonicmean, heavi
 using CryoGrid.Layers: Soil, porosity
 using CryoGrid.Utils
 
+using Base: @propagate_inbounds
 using DimensionalData
 using IfElse
 using Interpolations: Linear, Flat
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index 79c07c58..ab1f8b46 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -201,8 +201,14 @@ total water content (θw), and liquid water content (θl).
             I_c*(H/Lθ) + I_t
         end
     end
-    @inbounds @. state.θl = freezethaw(state.H, heat.L, state.θw)*state.θw
-    heatcapacity!(layer, heat, state) # update heat capacity, C
-    @inbounds @. state.T = enthalpyinv(state.H, state.C, heat.L, state.θw)
+    @inbounds for i in 1:length(state.H)
+        let θw = state.θw[i],
+            H = state.H[i],
+            L = heat.L;
+            state.θl[i] = θw*freezethaw(H, L, θw)
+            state.C[i] = heatcapacity(layer, heat, state, i) # update heat capacity, C
+            state.T[i] = enthalpyinv(H, state.C[i], L, θw)
+        end
+    end
     return nothing
 end
diff --git a/src/Physics/HeatConduction/heat_bc.jl b/src/Physics/HeatConduction/heat_bc.jl
index 9983809f..e7ea2fdf 100644
--- a/src/Physics/HeatConduction/heat_bc.jl
+++ b/src/Physics/HeatConduction/heat_bc.jl
@@ -12,7 +12,7 @@ struct TemperatureGradient{E,F} <: BoundaryProcess
 end
 BoundaryStyle(::Type{<:TemperatureGradient}) = Dirichlet()
 BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
-@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) where {F} = bc.T(s1.t)
+@inline boundaryvalue(bc::TemperatureGradient, l1, p2, l2, s1, s2) = bc.T(s1.t)
 
 @with_kw struct NFactor{W,S} <: BoundaryEffect
     winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
diff --git a/src/Physics/HeatConduction/soil/soilheat.jl b/src/Physics/HeatConduction/soil/soilheat.jl
index e1e08992..8d200fab 100644
--- a/src/Physics/HeatConduction/soil/soilheat.jl
+++ b/src/Physics/HeatConduction/soil/soilheat.jl
@@ -1,3 +1,4 @@
+# Thermal properties
 @inline function heatcapacity(soil::Soil, totalwater, liquidwater, mineral, organic)
     @unpack cw,co,cm,ca,ci = soil.hc
     let air = 1.0 - totalwater - mineral - organic,
@@ -14,13 +15,16 @@ end
         (liq*kw^0.5 + ice*ki^0.5 + mineral*km^0.5 + organic*ko^0.5 + air*ka^0.5)^2
     end
 end
-@inline function heatcapacity!(soil::Soil, ::Heat, state)
-    @inbounds @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
+@inline @propagate_inbounds heatcapacity(soil::Soil, ::Heat, state, i) = heatcapacity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
+@inline @propagate_inbounds function heatcapacity!(soil::Soil, ::Heat, state)
+    @. state.C = heatcapacity(soil, state.θw, state.θl, state.θm, state.θo)
 end
-@inline function thermalconductivity!(soil::Soil, ::Heat, state)
-    @inbounds @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
+@inline @propagate_inbounds thermalconductivity(soil::Soil, ::Heat, state, i) = thermalconductivity(soil, state.θw[i], state.θl[i], state.θm[i], state.θo[i])
+@inline @propagate_inbounds function thermalconductivity!(soil::Soil, ::Heat, state)
+    @. state.kc = thermalconductivity(soil, state.θw, state.θl, state.θm, state.θo)
 end
 
+# SFCC
 include("sfcc.jl")
 
 """
-- 
GitLab