Implemented template for LES eddy viscosity models.

Author: b-fg

src/Metrics.jl

@@ -9,6 +9,13 @@ Base.@propagate_inbounds @fastmath function permute(f,i)
     f(j,k)-f(k,j)
 end
 ×(a,b) = fSV(i->permute((j,k)->a[j]*b[k],i),3)
+@fastmath @inline function dot(a,b)
+    init=zero(eltype(a))
+    @inbounds for ij in eachindex(a)
+     init += a[ij] * b[ij]
+    end
+    return init
+end
 
 """
     ke(I::CartesianIndex,u,U=0)

src/util.jl

@@ -258,4 +258,37 @@ function interp(x::SVector{D}, varr::AbstractArray) where {D}
     # Shift to align with each staggered grid component and interpolate
     @inline shift(i) = SVector{D}(ifelse(i==j,0.5,0.0) for j in 1:D)
     return SVector{D}(interp(x+shift(i),@view(varr[..,i])) for i in 1:D)
+end
+
+# Turbulence modelling utils
+"""
+    sgs(flow, t; νₜ, S, C=0.16)
+
+Implements a user-defined function `udf` to model subgrid-scale LES stresses based on the Boussinesq approximation
+    τᵃᵢⱼ = τʳᵢⱼ - (1/3)τʳₖₖδᵢⱼ = -2νₜS̅ᵢⱼ
+where
+            ▁▁▁▁
+    τʳᵢⱼ =  uᵢuⱼ - u̅ᵢu̅ⱼ
+
+and we add -∂ⱼ(τᵃᵢⱼ) to the RHS as a body force (the isotropic part of the tensor is automatically modelled by the pressure gradient term).
+Users need to define the turbulent viscosity function `νₜ` and pass it as a keyword argument to this function together with rate-of-strain
+    tensor array buffer `S` and model constant `C`. For example, the standard Smagorinsky–Lilly model for the sub-grid scale stresses
+
+    νₜ = (CΔ)²|√(2S̅ᵢⱼS̅ᵢⱼ)|
+
+could be defined as
+    `smagorinsky(I::CartesianIndex{m} where m, S; C=0.16) = @views C^2*sqrt(2dot(S[I,:,:],S[I,:,:]))`
+and passed into `sim_step!` as a keyword argument:
+    `sim_step!(sim, ...; udf=sgs, νₜ=smagorinsky, S, C)`
+"""
+function sgs!(flow, t; νₜ, S, C=0.16)
+    N,n = size_u(flow.u)
+    @loop S[I,:,:] .= WaterLily.S(I,flow.u) over I ∈ inside(flow.σ)
+    for i ∈ 1:n, j ∈ 1:n
+        WaterLily.@loop (
+            flow.σ[I] = -νₜ(I,S;C)*∂(j,CI(I,i),flow.u);
+            flow.f[I,i] += flow.σ[I];
+        ) over I ∈ inside_u(N,j)
+        WaterLily.@loop flow.f[I-δ(j,I),i] -= flow.σ[I] over I ∈ WaterLily.inside_u(N,j)
+    end
 end