From 7f7ce17d73422f9e981308752322c52a7e39d293 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Thu, 5 Aug 2021 10:08:09 -0700 Subject: [PATCH 1/2] Convert to 1-based indexing --- integration_tests/utils/Cases.jl | 23 ++++++++++++----------- src/EDMF_Environment.jl | 19 +++++++++++-------- src/EDMF_Updrafts.jl | 10 ++++++---- src/Grid.jl | 24 ++++++++++++------------ src/ReferenceState.jl | 7 +++---- src/TurbulenceConvection.jl | 2 +- src/Turbulence_PrognosticTKE.jl | 10 +++++----- src/Variables.jl | 20 ++++++++++++-------- src/python_primitives.jl | 18 +++++++----------- src/turbulence_functions.jl | 3 +-- 10 files changed, 70 insertions(+), 66 deletions(-) diff --git a/integration_tests/utils/Cases.jl b/integration_tests/utils/Cases.jl index 58cfb323e..b4771ce8c 100644 --- a/integration_tests/utils/Cases.jl +++ b/integration_tests/utils/Cases.jl @@ -1010,8 +1010,8 @@ function initialize_forcing(self::CasesBase{TRMM_LBA}, Gr::Grid, Ref::ReferenceS A = off_arr(A') self.rad = A # store matrix in self - ind1 = Int(trunc(10.0 / 600.0)) - ind2 = Int(ceil(10.0 / 600.0)) - 1 + ind1 = Int(trunc(10.0 / 600.0)) + 1 + ind2 = Int(ceil(10.0 / 600.0)) @inbounds for k in face_indicies(Gr) if 10 % 600.0 == 0 self.Fo.dTdt[k] = self.rad[ind1, k] @@ -1047,15 +1047,15 @@ end function update_forcing(self::CasesBase{TRMM_LBA}, GMV::GridMeanVariables, TS::TimeStepping) Gr = self.Fo.Gr - ind2 = Int(ceil(TS.t / 600.0)) - ind1 = Int(trunc(TS.t / 600.0)) + ind2 = Int(ceil(TS.t / 600.0)) + 1 + ind1 = Int(trunc(TS.t / 600.0)) + 1 @inbounds for k in face_indicies(Gr) if Gr.z_half[k] >= 22699.48 self.Fo.dTdt[k] = 0.0 else if TS.t < 600.0 # first 10 min use the radiative forcing of t=10min (as in the paper) - self.Fo.dTdt[k] = self.rad[0, k] - elseif TS.t < 21600.0 && ind2 < 36 + self.Fo.dTdt[k] = self.rad[1, k] + elseif TS.t < 21600.0 && ind2 < 37 if TS.t % 600.0 == 0 self.Fo.dTdt[k] = self.rad[ind1, k] else @@ -1064,7 +1064,8 @@ function update_forcing(self::CasesBase{TRMM_LBA}, GMV::GridMeanVariables, TS::T (TS.t - self.rad_time[ind1]) + self.rad[ind1, k] end else - self.Fo.dTdt[k] = self.rad[35, k] + # TODO: remove hard-coded index + self.Fo.dTdt[k] = self.rad[36, k] end end end @@ -1168,8 +1169,8 @@ function update_surface(self::CasesBase{ARM_SGP}, GMV::GridMeanVariables, TS::Ti t_Sur_in = off_arr([0.0, 4.0, 6.5, 7.5, 10.0, 12.5, 14.5]) .* 3600 #LES time is in sec SH = off_arr([-30.0, 90.0, 140.0, 140.0, 100.0, -10, -10]) # W/m^2 LH = off_arr([5.0, 250.0, 450.0, 500.0, 420.0, 180.0, 0.0]) # W/m^2 - self.Sur.shf = pyinterp(off_arr([TS.t]), t_Sur_in, SH)[0] - self.Sur.lhf = pyinterp(off_arr([TS.t]), t_Sur_in, LH)[0] + self.Sur.shf = pyinterp(off_arr([TS.t]), t_Sur_in, SH)[1] + self.Sur.lhf = pyinterp(off_arr([TS.t]), t_Sur_in, LH)[1] # if fluxes vanish bflux vanish and wstar and obukov length are NaNs ## CK +++ I commented out the lines below as I don"t think this is how we want to fix things! # if self.Sur.shf < 1.0 @@ -1188,8 +1189,8 @@ function update_forcing(self::CasesBase{ARM_SGP}, GMV::GridMeanVariables, TS::Ti AT_in = off_arr([0.0, 0.0, 0.0, -0.08, -0.016, -0.016]) ./ 3600.0 # Advective forcing for theta [K/h] converted to [K/sec] RT_in = off_arr([-0.125, 0.0, 0.0, 0.0, 0.0, -0.1]) ./ 3600.0 # Radiative forcing for theta [K/h] converted to [K/sec] Rqt_in = off_arr([0.08, 0.02, 0.04, -0.1, -0.16, -0.3]) ./ 1000.0 ./ 3600.0 # Radiative forcing for qt converted to [kg/kg/sec] - dTdt = pyinterp(off_arr([TS.t]), t_in, AT_in)[0] + pyinterp(off_arr([TS.t]), t_in, RT_in)[0] - dqtdt = pyinterp(off_arr([TS.t]), t_in, Rqt_in)[0] + dTdt = pyinterp(off_arr([TS.t]), t_in, AT_in)[1] + pyinterp(off_arr([TS.t]), t_in, RT_in)[1] + dqtdt = pyinterp(off_arr([TS.t]), t_in, Rqt_in)[1] Gr = self.Fo.Gr @inbounds for k in center_indicies(Gr) if Gr.z_half[k] <= 1000.0 diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index 274c68de6..81f696093 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -1,16 +1,18 @@ function set_bcs(self::EnvironmentVariable, Gr::Grid) - start_low = Gr.gw - 1 - start_high = Gr.nzg - Gr.gw - 1 + start_low = Gr.gw + start_high = Gr.nzg - Gr.gw if self.name == "w" self.values[start_high] = 0.0 self.values[start_low] = 0.0 - @inbounds for k in xrange(1, Gr.gw) + @inbounds for kk in xrange(1, Gr.gw) + k = kk - 1 self.values[start_high + k] = -self.values[start_high - k] self.values[start_low - k] = -self.values[start_low + k] end else - @inbounds for k in xrange(Gr.gw) + @inbounds for kk in xrange(Gr.gw) + k = kk - 1 self.values[start_high + k + 1] = self.values[start_high - k] self.values[start_low - k] = self.values[start_low + 1 + k] end @@ -18,10 +20,11 @@ function set_bcs(self::EnvironmentVariable, Gr::Grid) end function set_bcs(self::EnvironmentVariable_2m, Gr::Grid) - start_low = Gr.gw - 1 - start_high = Gr.nzg - Gr.gw - 1 + start_low = Gr.gw + start_high = Gr.nzg - Gr.gw - @inbounds for k in xrange(Gr.gw) + @inbounds for kk in xrange(Gr.gw) + k = kk - 1 self.values[start_high + k + 1] = self.values[start_high - k] self.values[start_low - k] = self.values[start_low + 1 + k] end @@ -270,7 +273,7 @@ function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVari corr = fmax(fmin(EnvVar.HQTcov.values[k] / fmax(sd_h * sd_q, 1e-13), 1.0), -1.0) # limit sd_q to prevent negative qt_hat - sd_q_lim = (1e-10 - EnvVar.QT.values[k]) / (sqrt2 * abscissas[0]) + sd_q_lim = (1e-10 - EnvVar.QT.values[k]) / (sqrt2 * abscissas[1]) # walking backwards to assure your q_t will not be smaller than 1e-10 # TODO - check # TODO - change 1e-13 and 1e-10 to some epislon diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl index e4ebd9aa1..b932dcbc7 100644 --- a/src/EDMF_Updrafts.jl +++ b/src/EDMF_Updrafts.jl @@ -1,6 +1,6 @@ function set_bcs(self::UpdraftVariable, Gr::Grid) - start_low = Gr.gw - 1 - start_high = Gr.nzg - Gr.gw - 1 + start_low = Gr.gw + start_high = Gr.nzg - Gr.gw n_updrafts = size(self.values)[1] @@ -8,13 +8,15 @@ function set_bcs(self::UpdraftVariable, Gr::Grid) @inbounds for i in xrange(n_updrafts) self.values[i, start_high] = 0.0 self.values[i, start_low] = 0.0 - @inbounds for k in xrange(1, Gr.gw) + @inbounds for kk in xrange(1, Gr.gw) + k = kk - 1 self.values[i, start_high + k] = -self.values[i, start_high - k] self.values[i, start_low - k] = -self.values[i, start_low + k] end end else - @inbounds for k in xrange(Gr.gw) + @inbounds for kk in xrange(Gr.gw) + k = kk - 1 @inbounds for i in xrange(n_updrafts) self.values[i, start_high + k + 1] = self.values[i, start_high - k] self.values[i, start_low - k] = self.values[i, start_low + 1 + k] diff --git a/src/Grid.jl b/src/Grid.jl index 556e6bd9a..5774b13c2 100644 --- a/src/Grid.jl +++ b/src/Grid.jl @@ -24,20 +24,20 @@ struct Grid{A1, FT} nz = namelist["grid"]["nz"] nzg = nz + 2 * gw - cinterior = gw:((nzg - gw) - 1) - finterior = gw:((nzg - gw) - 1) + cinterior = (gw + 1):(nzg - gw) + finterior = (gw + 1):(nzg - gw) # TODO: make cell centers and cell faces different sizes z_half = pyzeros(nz + 2 * gw) z = pyzeros(nz + 2 * gw) - count = 0 - @inbounds for i in xrange(-gw, nz + gw) - z[count] = (i + 1) * dz - z_half[count] = (i + 0.5) * dz + count = 1 + @inbounds for i in range(1 - gw, stop = nz + gw) + z[count] = i * dz + z_half[count] = (i - 0.5) * dz count += 1 end - zmin = z[gw - 1] - zmax = z[nzg - gw - 1] + zmin = z[gw] + zmax = z[nzg - gw] A1 = typeof(z) FT = typeof(zmax) return new{A1, FT}(zmin, zmax, cinterior, finterior, dz, dzi, gw, nz, nzg, z_half, z) @@ -45,10 +45,10 @@ struct Grid{A1, FT} end # Index of the first interior cell above the surface -kc_surface(grid::Grid) = grid.gw -kf_surface(grid::Grid) = grid.gw - 1 -kc_top_of_atmos(grid::Grid) = grid.nzg - grid.gw - 1 -kf_top_of_atmos(grid::Grid) = grid.nzg - grid.gw - 1 +kc_surface(grid::Grid) = grid.gw + 1 +kf_surface(grid::Grid) = grid.gw +kc_top_of_atmos(grid::Grid) = grid.nzg - grid.gw +kf_top_of_atmos(grid::Grid) = grid.nzg - grid.gw is_surface_center(grid::Grid, k::Int) = k == kc_surface(grid) is_toa_center(grid::Grid, k::Int) = k == kc_top_of_atmos(grid) diff --git a/src/ReferenceState.jl b/src/ReferenceState.jl index eda8915dc..dacc85ae8 100644 --- a/src/ReferenceState.jl +++ b/src/ReferenceState.jl @@ -68,7 +68,6 @@ function initialize(self::ReferenceState, grid::Grid, Stats::NetCDFIO_Stats) # Perform the integration # TODO: replace with OrdinaryDiffEq - grid z_span = (grid.zmin, grid.zmax) @show z_span prob = ODEProblem(rhs, logp, z_span) @@ -83,12 +82,12 @@ function initialize(self::ReferenceState, grid::Grid, Stats::NetCDFIO_Stats) # Set boundary conditions (in log-space) by mirroring log-pressure # TODO: re-generalize, is setting the BCs like this correct? - p[0] = p[1] + p[1] = p[2] p[end] = p[end - 1] - p_half[1] = p_half[2] + p_half[2] = p_half[3] p_half[end - 1] = p_half[end - 2] - p_half[0] = p_half[3] + p_half[1] = p_half[4] p_half[end] = p_half[end - 3] p .= exp.(p) diff --git a/src/TurbulenceConvection.jl b/src/TurbulenceConvection.jl index bd43b5cbc..b29d26e1c 100644 --- a/src/TurbulenceConvection.jl +++ b/src/TurbulenceConvection.jl @@ -22,7 +22,7 @@ const CPSGS = CLIMAParameters.Atmos.SubgridScale # For dispatching to inherited class struct BaseCase end -up_sum(vals::OffsetArray) = off_arr(reshape(sum(vals; dims = 1), size(vals, 2))) +up_sum(vals::AbstractArray) = off_arr(reshape(sum(vals; dims = 1), size(vals, 2))) function parse_namelist(namelist, keys...; default = nothing, valid_options = nothing) param = namelist diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index f4e30ef81..77db7427a 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -824,7 +824,7 @@ function get_GMV_CoVar( # -1, indexing phi_e.values[-1] yields the # _last_ value in the array. This is certainly # not intended - if k ≠ 0 + if k ≠ 1 phi_diff = interp2pt(phi_e.values[k - 1] - gmv_phi[k - 1], phi_e.values[k] - gmv_phi[k]) psi_diff = interp2pt(psi_e.values[k - 1] - gmv_psi[k - 1], psi_e.values[k] - gmv_psi[k]) else # just use 0th order approximation @@ -838,7 +838,7 @@ function get_GMV_CoVar( # -1, indexing phi_e.values[-1] yields the # _last_ value in the array. This is certainly # not intended - if k ≠ 0 + if k ≠ 1 phi_diff = interp2pt(phi_u.values[i, k - 1] - gmv_phi[k - 1], phi_u.values[i, k] - gmv_phi[k]) psi_diff = interp2pt(psi_u.values[i, k - 1] - gmv_psi[k - 1], psi_u.values[i, k] - gmv_psi[k]) else # just use 0th order approximation @@ -968,9 +968,9 @@ function compute_updraft_closures(self::EDMF_PrognosticTKE, GMV::GridMeanVariabl @inbounds for k in real_center_indicies(grid) @inbounds for i in xrange(self.n_updrafts) input_p.updraft_top = self.UpdVar.updraft_top[i] - alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])) - 1, 0) + alen = max(length(argwhere(self.UpdVar.Area.values[i, cinterior])), 1) avals = off_arr(self.UpdVar.Area.values[i, cinterior]) - input_p.a_med = Statistics.median(avals[0:alen]) + input_p.a_med = Statistics.median(avals[1:alen]) input.zi = self.UpdVar.cloud_base[i] # entrainment input.buoy_ed_flux = self.EnvVar.TKE.buoy[k] @@ -2092,7 +2092,7 @@ function compute_covariance_interdomain_src( # -1, indexing phi_e.values[-1] yields the # _last_ value in the array. This is certainly # not intended - if k ≠ 0 + if k ≠ 1 phi_diff = interp2pt(phi_u.values[i, k - 1], phi_u.values[i, k]) - interp2pt(phi_e.values[k - 1], phi_e.values[k]) diff --git a/src/Variables.jl b/src/Variables.jl index ed7c0756f..838c3e27a 100644 --- a/src/Variables.jl +++ b/src/Variables.jl @@ -6,11 +6,12 @@ function zero_tendencies(self::VariablePrognostic, Gr::Grid) end function set_bcs(self::VariablePrognostic, Gr::Grid) - start_low = Gr.gw - 1 - start_high = Gr.nzg - Gr.gw - 1 + start_low = Gr.gw + start_high = Gr.nzg - Gr.gw if self.bc == "sym" - @inbounds for k in xrange(Gr.gw) + @inbounds for kk in xrange(Gr.gw) + k = kk - 1 self.values[start_high + k + 1] = self.values[start_high - k] self.values[start_low - k] = self.values[start_low + 1 + k] @@ -30,7 +31,8 @@ function set_bcs(self::VariablePrognostic, Gr::Grid) self.new[start_high] = 0.0 self.new[start_low] = 0.0 - @inbounds for k in xrange(1, Gr.gw) + @inbounds for kk in xrange(1, Gr.gw) + k = kk - 1 self.values[start_high + k] = -self.values[start_high - k] self.values[start_low - k] = -self.values[start_low + k] @@ -46,11 +48,12 @@ function set_bcs(self::VariablePrognostic, Gr::Grid) end function set_bcs(self::VariableDiagnostic, Gr::Grid) - start_low = Gr.gw - 1 - start_high = Gr.nzg - Gr.gw + start_low = Gr.gw + start_high = Gr.nzg - Gr.gw + 1 if self.bc == "sym" - @inbounds for k in xrange(Gr.gw) + @inbounds for kk in xrange(Gr.gw) + k = kk - 1 self.values[start_high + k] = self.values[start_high - 1] self.values[start_low - k] = self.values[start_low + 1] end @@ -58,7 +61,8 @@ function set_bcs(self::VariableDiagnostic, Gr::Grid) else self.values[start_high] = 0.0 self.values[start_low] = 0.0 - @inbounds for k in xrange(1, Gr.gw) + @inbounds for kk in xrange(1, Gr.gw) + k = kk - 1 self.values[start_high + k] = 0.0 #-self.values[start_high - k ] self.values[start_low - k] = 0.0 #-self.values[start_low + k ] end diff --git a/src/python_primitives.jl b/src/python_primitives.jl index d613c5309..7c4229356 100644 --- a/src/python_primitives.jl +++ b/src/python_primitives.jl @@ -7,21 +7,17 @@ fabs(a) = abs(a) argwhere(a) = findall(a .> 0) # Use stop-1 to emulate python ranges -xrange(start, stop, step = 1) = range(start, stop - 1; step = step) +xrange(start, stop, step = 1) = range(start + 1, stop; step = step) xrange(stop) = xrange(0, stop) -pyzeros(n::Int) = OffsetArray(zeros(n), 0:(n - 1)) -pyzeros(m::Int, n::Int) = OffsetArray(zeros(m, n), 0:(m - 1), 0:(n - 1)) -pyones(n::Int) = OffsetArray(ones(n), 0:(n - 1)) -pyones(m::Int, n::Int) = OffsetArray(ones(m, n), 0:(m - 1), 0:(n - 1)) +pyzeros(n::Int) = zeros(n) +pyzeros(m::Int, n::Int) = zeros(m, n) -function off_arr(a::AbstractArray) - dims = ntuple(ndims(a)) do i - 0:(size(a, i) - 1) - end - return OffsetArray(a, dims) -end +pyones(n::Int) = ones(n) +pyones(m::Int, n::Int) = ones(m, n) + +off_arr(a::AbstractArray) = a pow(a::Real, b::Real) = a^b pow(a::AbstractVector, b::AbstractVector) = a .^ b diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl index 794c30196..9186be24e 100644 --- a/src/turbulence_functions.jl +++ b/src/turbulence_functions.jl @@ -262,8 +262,7 @@ function construct_tridiag_diffusion(grid, dt, rho_ae_K_m, rho, ae, a, b, c) end function tridiag_solve(b_rhs, a, b, c) - # Note that `1:end` is zero-based indexing. - A = Tridiagonal(a[1:end], parent(b), c[0:(end - 1)]) + A = Tridiagonal(a[2:end], parent(b), c[1:(end - 1)]) return A \ parent(b_rhs) end From 741f2baab21babe1b661ac5ab700e41d1a57e432 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Fri, 6 Aug 2021 14:38:03 -0700 Subject: [PATCH 2/2] -1 indexing in radiation time --- integration_tests/utils/Cases.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/integration_tests/utils/Cases.jl b/integration_tests/utils/Cases.jl index b4771ce8c..676c2901d 100644 --- a/integration_tests/utils/Cases.jl +++ b/integration_tests/utils/Cases.jl @@ -1017,7 +1017,7 @@ function initialize_forcing(self::CasesBase{TRMM_LBA}, Gr::Grid, Ref::ReferenceS self.Fo.dTdt[k] = self.rad[ind1, k] else self.Fo.dTdt[k] = - (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2 + 1] - self.rad_time[ind1 + 1]) * (10.0) + + (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2] - self.rad_time[ind1]) * (10.0) + self.rad[ind1, k] end end @@ -1060,8 +1060,8 @@ function update_forcing(self::CasesBase{TRMM_LBA}, GMV::GridMeanVariables, TS::T self.Fo.dTdt[k] = self.rad[ind1, k] else self.Fo.dTdt[k] = - (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2] - self.rad_time[ind1]) * - (TS.t - self.rad_time[ind1]) + self.rad[ind1, k] + (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2 - 1] - self.rad_time[ind1 - 1]) * + (TS.t - self.rad_time[ind1 - 1]) + self.rad[ind1, k] end else # TODO: remove hard-coded index