From 771a660b601131d3f3a1e6b68acfd7f53c319b8f Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Sat, 28 Aug 2021 21:53:33 -0700 Subject: [PATCH 1/2] =?UTF-8?q?Rename=20=E2=88=87=5Fonesided=20->=20c?= =?UTF-8?q?=E2=88=87=5Fonesided?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Forcing.jl | 12 ++++++------ src/Operators.jl | 20 ++++++++++---------- src/Turbulence_PrognosticTKE.jl | 2 +- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/Forcing.jl b/src/Forcing.jl index 3c39e0ff36..f82f142a2a 100644 --- a/src/Forcing.jl +++ b/src/Forcing.jl @@ -50,8 +50,8 @@ function update(self::ForcingBase{ForcingStandard}, GMV::GridMeanVariables) # Apply large-scale subsidence tendencies H_cut = ccut_onesided(GMV.H.values, grid, k) q_tot_cut = ccut_onesided(GMV.QT.values, grid, k) - ∇H = ∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) - ∇q_tot = ∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇H = c∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇q_tot = c∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) GMV.H.tendencies[k] -= ∇H * self.subsidence[k] GMV.QT.tendencies[k] -= ∇q_tot * self.subsidence[k] end @@ -79,8 +79,8 @@ function update(self::ForcingBase{ForcingDYCOMS_RF01}, GMV::GridMeanVariables) @inbounds for k in real_center_indicies(grid) H_cut = ccut_onesided(GMV.H.values, grid, k) q_tot_cut = ccut_onesided(GMV.QT.values, grid, k) - ∇H = ∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) - ∇q_tot = ∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇H = c∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇q_tot = c∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) GMV.QT.tendencies[k] += self.dqtdt[k] # Apply large-scale subsidence tendencies @@ -140,8 +140,8 @@ function update(self::ForcingBase{ForcingLES}, GMV::GridMeanVariables) H_cut = ccut_onesided(GMV.H.values, grid, k) q_tot_cut = ccut_onesided(GMV.QT.values, grid, k) - ∇H = ∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) - ∇q_tot = ∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇H = c∇_onesided(H_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇q_tot = c∇_onesided(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) GMV_H_subsidence_k = -∇H * self.subsidence[k] GMV_QT_subsidence_k = -∇q_tot * self.subsidence[k] else diff --git a/src/Operators.jl b/src/Operators.jl index 61c8635074..271bf348fa 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -35,26 +35,26 @@ end ∇f2c(f::SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[2] - bc.value) * grid.dzi ∇f2c(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value -function ∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) +function c∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) if is_surface_center(grid, k) - return ∇_onesided(f_dual, grid, BottomBCTag(), bottom) + return c∇_onesided(f_dual, grid, BottomBCTag(), bottom) elseif is_toa_center(grid, k) - return ∇_onesided(f_dual, grid, TopBCTag(), top) + return c∇_onesided(f_dual, grid, TopBCTag(), top) else - return ∇_onesided(f_dual, grid, InteriorTag()) + return c∇_onesided(f_dual, grid, InteriorTag()) end end -∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi -∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2) +c∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi +c∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2) # TODO: this is a crud approximation, as we're specifying what should be the derivative # at the boundary, and we're taking this as the derivative at the first interior at the # top of the domain. -∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value -∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info +c∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value +c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info # TODO: this is a crud approximation, as we're specifying what should be the derivative # at the boundary, and we're taking this as the derivative at the first interior at the # top of the domain. -∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value +c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value # Used when traversing cell faces @@ -111,7 +111,7 @@ function upwind_advection_area(ρ0_half::Vector{Float64}, a_up::Vector{Float64}, ρ_0_cut = ccut_onesided(ρ0_half, grid, k) a_up_cut = ccut_onesided(a_up, grid, k) m_cut = ρ_0_cut .* a_up_cut .* w_up_cut - ∇m = ∇_onesided(m_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇m = c∇_onesided(m_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) # TODO: Why are we dividing by ρ0_half[k + 1]? return -∇m / ρ0_half[k + 1] end diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 2a967f625e..cfcf6f5c93 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -1837,7 +1837,7 @@ function update_inversion(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, opti @inbounds for k in real_center_indicies(grid) ∇θ_liq_cut = ccut_onesided(GMV.H.values, grid, k) - ∇θ_liq = ∇_onesided(∇θ_liq_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + ∇θ_liq = c∇_onesided(∇θ_liq_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) if ∇θ_liq > ∇θ_liq_max ∇θ_liq_max = ∇θ_liq self.zi = grid.z[k] From 81125fd0a8f8c3dff358a781d769965b02e9c22c Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Sat, 28 Aug 2021 22:44:29 -0700 Subject: [PATCH 2/2] Use bcs with operators in upwind_advection_velocity --- src/Operators.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/Operators.jl b/src/Operators.jl index 271bf348fa..091915c7da 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -56,6 +56,21 @@ c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - # top of the domain. c∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value +function f∇_onesided(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) + if is_surface_face(grid, k) + return f∇_onesided(f_dual, grid, BottomBCTag(), bottom) + elseif is_toa_face(grid, k) + return f∇_onesided(f_dual, grid, TopBCTag(), top) + else + return f∇_onesided(f_dual, grid, InteriorTag()) + end +end +f∇_onesided(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi +f∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * (grid.dzi / 2) +f∇_onesided(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value +f∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.dzi # don't use BC info +f∇_onesided(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value + # Used when traversing cell faces interpc2f(f, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) = @@ -124,6 +139,16 @@ function upwind_advection_velocity(ρ0::Vector{Float64}, a_up::Vector{Float64}, return adv end +function upwind_advection_velocity(ρ0::Vector{Float64}, a_up::Vector{Float64}, w_up::Vector{Float64}, grid, k) + a_dual = fdaul_onesided(a_up, grid, k) + ρ_0_dual = dual_faces(ρ0, grid, k) + w_up_dual = dual_faces(w_up, grid, k) + adv_dual = a_dual .* ρ_0_dual .* w_up_dual .* w_up_dual + ∇ρaw = f∇_onesided(adv_dual, grid, k; bottom = FreeBoundary(), top = SetGradient(0)) + return ∇ρaw +end + + function upwind_advection_scalar( ρ0_half::Vector{Float64}, a_up::Vector{Float64},