Skip to content

Commit

Permalink
Merge #177
Browse files Browse the repository at this point in the history 
177: Use bcs with operators in upwind_advection_velocity r=charleskawczynski a=charleskawczynski

This PR:
 - Makes upwind_advection_velocity use bcs
 - Renames `∇_onesided` -> `c∇_onesided`
 - defines `f∇_onesided`

Co-authored-by: Charles Kawczynski <kawczynski.charles@gmail.com>
  • Loading branch information
@bors @charleskawczynski
bors[bot] and charleskawczynski committed Aug 29, 2021
2 parents e51a08c + 81125fd commit 56e6387
Show file tree
Hide file tree
Showing 3 changed files with 42 additions and 17 deletions.
12 changes: 6 additions & 6 deletions src/Forcing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
45 changes: 35 additions & 10 deletions src/Operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,26 +35,41 @@ 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

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

Expand Down Expand Up @@ -111,7 +126,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
Expand All @@ -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},
Expand Down
2 changes: 1 addition & 1 deletion src/Turbulence_PrognosticTKE.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down

0 comments on commit 56e6387

Please sign in to comment.