Skip to content

Commit

Permalink
Use more ClimaCore operators, rm upwind gradient
Browse files Browse the repository at this point in the history
  • Loading branch information
@charleskawczynski
charleskawczynski committed Dec 3, 2021
1 parent d1d8ac6 commit 7e65abe
Show file tree
Hide file tree
Showing 4 changed files with 31 additions and 37 deletions.
4 changes: 4 additions & 0 deletions integration_tests/utils/main.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,11 @@ cent_aux_vars_edmf(FT, n_up) = (;
KM = FT(0),
KH = FT(0),
mixing_length = FT(0),
m_entr_detr = FT(0),
∇m_entr_detr = FT(0),
θ_virt = FT(0),
∇θ_virt = FT(0),
∇Ri_bulk = FT(0),
massflux_tendency_h = FT(0),
massflux_tendency_qt = FT(0),
diffusive_tendency_h = FT(0),
Expand Down
18 changes: 0 additions & 18 deletions src/Operators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,24 +36,6 @@ end
∇c2f(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[1] - bc.value) * grid.Δzi * 2.0
∇c2f(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value

c∇_upwind(f, grid::Grid, k; bottom = NoBCGivenError(), top = NoBCGivenError()) =
c∇_upwind(ccut_upwind(f, grid, k), grid, k; bottom, top)

function c∇_upwind(f_dual::SA.SVector, grid::Grid, k; bottom = NoBCGivenError(), top = NoBCGivenError())
if is_surface_center(grid, k)
return c∇_upwind(f_dual, grid, BottomBCTag(), bottom)
elseif is_toa_center(grid, k)
return c∇_upwind(f_dual, grid, TopBCTag(), top)
else
return c∇_upwind(f_dual, grid, InteriorTag())
end
end
c∇_upwind(f::SA.SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.Δzi
c∇_upwind(f::SA.SVector, grid::Grid, ::TopBCTag, bc::FreeBoundary) = (f[2] - f[1]) * grid.Δzi
c∇_upwind(f::SA.SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value
c∇_upwind(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[1] - bc.value) * (grid.Δzi * 2)
c∇_upwind(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value

# Used when traversing cell faces

interpc2f(f, grid::Grid, k::CCO.PlusHalf; bottom = NoBCGivenError(), top = NoBCGivenError()) =
Expand Down
14 changes: 10 additions & 4 deletions src/turbulence_functions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,15 +16,20 @@ function get_inversion(grid::Grid, state, param_set, Ri_bulk_crit)
u = center_prog_grid_mean(state).u
v = center_prog_grid_mean(state).v
Ri_bulk = center_aux_bulk(state).Ri
∇θ_virt = center_aux_turbconv(state).∇θ_virt
∇Ri_bulk = center_aux_turbconv(state).∇Ri_bulk
θ_virt_b = θ_virt[kc_surf]
z_c = grid.zc
∇c = CCO.DivergenceF2C()
wvec = CC.Geometry.WVector

# test if we need to look at the free convective limit
if (u[kc_surf]^2 + v[kc_surf]^2) <= 0.01
kmask = map(k -> (k, θ_virt[k] > θ_virt_b), real_center_indices(grid))
k_star = first(kmask[findlast(km -> km[2], kmask)])
∇θ_virt = c∇_upwind(θ_virt, grid, k_star; bottom = SetGradient(0), top = FreeBoundary())
h = (θ_virt_b - θ_virt[k_star - 1]) / ∇θ_virt + z_c[k_star - 1]
LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(θ_virt[kc_surf]))
@. ∇θ_virt = ∇c(wvec(LB(θ_virt)))
h = (θ_virt_b - θ_virt[k_star - 1]) / ∇θ_virt[k_star] + z_c[k_star - 1]
else
Ri_bulk_fn(k) = g * (θ_virt[k] - θ_virt_b) * z_c[k] / θ_virt_b / (u[k] * u[k] + v[k] * v[k])

Expand All @@ -33,8 +38,9 @@ function get_inversion(grid::Grid, state, param_set, Ri_bulk_crit)
end
kmask = map(k -> (k, Ri_bulk_fn(k) > Ri_bulk_crit), real_center_indices(grid))
k_star = first(kmask[findlast(km -> km[2], kmask)])
∇Ri_bulk = c∇_upwind(Ri_bulk, grid, k_star; bottom = SetGradient(0), top = FreeBoundary())
h = (Ri_bulk_crit - Ri_bulk[k_star - 1]) / ∇Ri_bulk + z_c[k_star - 1]
LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(Ri_bulk[kc_surf]))
@. ∇Ri_bulk = ∇c(wvec(LB(Ri_bulk)))
h = (Ri_bulk_crit - Ri_bulk[k_star - 1]) / ∇Ri_bulk[k_star] + z_c[k_star - 1]
end

return h
Expand Down
32 changes: 17 additions & 15 deletions src/update_aux.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,9 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para
prog_up_f = face_prog_updrafts(state)
aux_en_unsat = aux_en.unsat
aux_en_sat = aux_en.sat
m_entr_detr = aux_tc.m_entr_detr
∇m_entr_detr = aux_tc.∇m_entr_detr
wvec = CC.Geometry.WVector

#####
##### center variables
Expand Down Expand Up @@ -238,20 +241,20 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para
##### compute_updraft_closures
#####

@inbounds for k in real_center_indices(grid)
@inbounds for i in 1:N_up
Ic = CCO.InterpolateF2C()
∇c = CCO.DivergenceF2C()
LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(FT(0)))

@inbounds for i in 1:N_up
# compute ∇m at cell centers
a_up = aux_up[i].area
w_up = aux_up_f[i].w
w_gm = prog_gm_f.w
@. m_entr_detr = a_up * (Ic(w_up) - Ic(w_gm))
@. ∇m_entr_detr = ∇c(wvec(LB(m_entr_detr)))
@inbounds for k in real_center_indices(grid)
# entrainment
if aux_up[i].area[k] > 0.0
# compute ∇m at cell centers
a_up_c = aux_up[i].area[k]
w_up_c = interpf2c(aux_up_f[i].w, grid, k)
w_gm_c = interpf2c(prog_gm_f.w, grid, k)
m = a_up_c * (w_up_c - w_gm_c)
a_up_cut = ccut_upwind(aux_up[i].area, grid, k)
w_up_cut = daul_f2c_upwind(aux_up_f[i].w, grid, k)
w_gm_cut = daul_f2c_upwind(prog_gm_f.w, grid, k)
m_cut = a_up_cut .* (w_up_cut .- w_gm_cut)
∇m = FT(c∇_upwind(m_cut, grid, k; bottom = SetValue(0), top = FreeBoundary()))

w_min = 0.001

Expand All @@ -267,8 +270,8 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para
b_up = aux_up[i].buoy[k],
b_en = aux_en.buoy[k],
tke = aux_en.tke[k],
dMdz = ∇m,
M = m,
dMdz = ∇m_entr_detr[k],
M = m_entr_detr[k],
a_up = aux_up[i].area[k],
a_en = aux_en.area[k],
R_up = edmf.pressure_plume_spacing[i],
Expand Down Expand Up @@ -298,7 +301,6 @@ function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, para
end
end

wvec = CC.Geometry.WVector
@inbounds for i in 1:N_up

# pressure
Expand Down

0 comments on commit 7e65abe

Please sign in to comment.