From fc7407bb2f5b6d6940b20921d4c2bd2cb115efb5 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Fri, 30 Jul 2021 08:28:21 -0700 Subject: [PATCH] Add basic operators --- Project.toml | 1 + src/Forcing.jl | 8 ++++---- src/Operators.jl | 35 +++++++++++++++++++++++++++++++++ src/Turbulence.jl | 8 ++++---- src/TurbulenceConvection.jl | 2 ++ src/Turbulence_PrognosticTKE.jl | 8 +++++--- 6 files changed, 51 insertions(+), 11 deletions(-) create mode 100644 src/Operators.jl diff --git a/Project.toml b/Project.toml index 39ea9cd692..0f0bc54fd2 100644 --- a/Project.toml +++ b/Project.toml @@ -19,6 +19,7 @@ OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" diff --git a/src/Forcing.jl b/src/Forcing.jl index b429187510..9638e0f06a 100644 --- a/src/Forcing.jl +++ b/src/Forcing.jl @@ -58,8 +58,8 @@ function update(self::ForcingBase{ForcingStandard}, GMV::GridMeanVariables) if self.apply_subsidence @inbounds for k in real_center_indicies(self.Gr) # Apply large-scale subsidence tendencies - GMV.H.tendencies[k] -= (GMV.H.values[k + 1] - GMV.H.values[k]) * self.Gr.dzi * self.subsidence[k] - GMV.QT.tendencies[k] -= (GMV.QT.values[k + 1] - GMV.QT.values[k]) * self.Gr.dzi * self.subsidence[k] + GMV.H.tendencies[k] -= ∇_upwind(GMV.H.values, self.Gr, k) * self.subsidence[k] + GMV.QT.tendencies[k] -= ∇_upwind(GMV.QT.values, self.Gr, k) * self.subsidence[k] end end if self.apply_coriolis @@ -85,8 +85,8 @@ function update(self::ForcingBase{ForcingDYCOMS_RF01}, GMV::GridMeanVariables) @inbounds for k in real_center_indicies(self.Gr) GMV.QT.tendencies[k] += self.dqtdt[k] # Apply large-scale subsidence tendencies - GMV.H.tendencies[k] -= (GMV.H.values[k + 1] - GMV.H.values[k]) * self.Gr.dzi * self.subsidence[k] - GMV.QT.tendencies[k] -= (GMV.QT.values[k + 1] - GMV.QT.values[k]) * self.Gr.dzi * self.subsidence[k] + GMV.H.tendencies[k] -= ∇_upwind(GMV.H.values, self.Gr, k) * self.subsidence[k] + GMV.QT.tendencies[k] -= ∇_upwind(GMV.QT.values, self.Gr, k) * self.subsidence[k] end if self.apply_coriolis diff --git a/src/Operators.jl b/src/Operators.jl new file mode 100644 index 0000000000..d3a67e3826 --- /dev/null +++ b/src/Operators.jl @@ -0,0 +1,35 @@ + + +function ∇_collocated(f, grid::Grid, k::Int) + return (f[k+1] - f[k-1]) / 2 * grid.dzi +end + +function ∇_collocated(f::SVector, grid::Grid) + @assert length(f) == 3 + return (∇_staggered(f[2:3], grid) - ∇_staggered(f[1:2], grid)) * grid.dzi +end + +function ∇_staggered(f, grid::Grid) + @assert length(f) == 2 + return (f[2] - f[1]) * grid.dzi +end + +function ∇f2c(f, grid::Grid, k::Int) + return (f[k] - f[k-1]) * grid.dzi +end + +function ∇c2f(f, grid::Grid, k::Int) + return (f[k+1] - f[k]) * grid.dzi +end + +function ∇_upwind(f, grid::Grid, k::Int) + return (f[k+1] - f[k]) * grid.dzi +end + +function interpc2f(f, grid::Grid, k::Int) + return 0.5*(f[k+1] + f[k]) +end + +function cut(f, k) + return SVector(f[k-1], f[k], f[k+1]) +end diff --git a/src/Turbulence.jl b/src/Turbulence.jl index ceb918620b..7843ccb6c1 100644 --- a/src/Turbulence.jl +++ b/src/Turbulence.jl @@ -27,7 +27,7 @@ end # Update the diagnosis of the inversion height, using the maximum temperature gradient method function update_inversion(self::ParameterizationBase, GMV::GridMeanVariables, option) theta_rho = center_field(self.Gr) - maxgrad = 0.0 + ∇θ_liq_max = 0.0 k_fi = first_center(self.Gr) @inbounds for k in real_center_indicies(self.Gr) @@ -46,9 +46,9 @@ function update_inversion(self::ParameterizationBase, GMV::GridMeanVariables, op elseif option == "thetal_maxgrad" @inbounds for k in real_center_indicies(self.Gr) - grad = (GMV.THL.values[k + 1] - GMV.THL.values[k]) * self.Gr.dzi - if grad > maxgrad - maxgrad = grad + ∇θ_liq = ∇_upwind(GMV.THL.values, self.Gr, k) + if ∇θ_liq > ∇θ_liq_max + ∇θ_liq_max = ∇θ_liq self.zi = self.Gr.z[k] end end diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl index 778f99d7e1..ad021dd299 100644 --- a/src/TurbulenceConvection.jl +++ b/src/TurbulenceConvection.jl @@ -2,6 +2,7 @@ module TurbulenceConvection import Dierckx using StatsBase +using StaticArrays import Statistics using PoissonRandom: pois_rand import LambertW @@ -33,6 +34,7 @@ include("python_primitives.jl") include("parameters.jl") include("Grid.jl") include("Fields.jl") +include("Operators.jl") include("NetCDFIO.jl") include("ReferenceState.jl") include("types.jl") diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 4090398697..a7d3415a24 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -411,9 +411,9 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl # Buoyancy-shear-subdomain exchange-dissipation TKE equilibrium scale shear2 = - pow((GMV.U.values[k + 1] - GMV.U.values[k - 1]) * 0.5 * dzi, 2) + - pow((GMV.V.values[k + 1] - GMV.V.values[k - 1]) * 0.5 * dzi, 2) + - pow((self.EnvVar.W.values[k] - self.EnvVar.W.values[k - 1]) * dzi, 2) + pow(∇_collocated(GMV.U.values, grid, k), 2) + + pow(∇_collocated(GMV.V.values, grid, k), 2) + + pow(∇f2c(self.EnvVar.W.values, grid, k), 2) qt_dry = self.EnvThermo.qt_dry[k] th_dry = self.EnvThermo.th_dry[k] t_cloudy = self.EnvThermo.t_cloudy[k] @@ -428,6 +428,8 @@ function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariabl grad_qt_plus = (self.EnvVar.QT.values[k + 1] - self.EnvVar.QT.values[k]) * dzi grad_thl = interp2pt(grad_thl_low, grad_thl_plus) grad_qt = interp2pt(grad_qt_low, grad_qt_plus) + # grad_thl = ∇_collocated(self.EnvVar.THL.values, grid, k) + # grad_qt = ∇_collocated(self.EnvVar.QT.values, grid, k) # g/theta_ref prefactor = g * (Rd / ref_state.alpha0_half[k] / ref_state.p0_half[k]) * exner_c(ref_state.p0_half[k]) d_buoy_thetal_dry = prefactor * (1.0 + (eps_vi - 1.0) * qt_dry)