Use bcs and avoid accessing ghost points

integration_tests/ARM_SGP.jl

@@ -11,18 +11,18 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 5.5166749573532592e-01
-best_mse["updraft_area"] = 2.8505841871593282e+02
-best_mse["updraft_w"] = 7.7404844506417481e-01
-best_mse["updraft_qt"] = 5.2600486259435868e+01
-best_mse["updraft_thetal"] = 6.9238569266089868e+01
-best_mse["u_mean"] = 8.7994360629094174e+01
-best_mse["tke_mean"] = 2.9609371778386242e+00
-best_mse["temperature_mean"] = 2.5666937507309872e-04
-best_mse["ql_mean"] = 0.0000000000000000e+00
-best_mse["thetal_mean"] = 2.5573465612699223e-04
-best_mse["Hvar_mean"] = 1.1019193692254991e+02
-best_mse["QTvar_mean"] = 3.7543163816570626e+01
+best_mse["qt_mean"] = 5.8920349209467660e-01
+best_mse["updraft_area"] = 2.9161664375922396e+02
+best_mse["updraft_w"] = 9.2064966152387040e-01
+best_mse["updraft_qt"] = 5.2856934357682029e+01
+best_mse["updraft_thetal"] = 6.9072600122408147e+01
+best_mse["u_mean"] = 8.7994360629094189e+01
+best_mse["tke_mean"] = 2.8644243409080627e+00
+best_mse["temperature_mean"] = 2.7496219242928010e-04
+best_mse["ql_mean"] = 5.0850007235659434e-03
+best_mse["thetal_mean"] = 2.7328104527915207e-04
+best_mse["Hvar_mean"] = 2.2244310672644488e+01
+best_mse["QTvar_mean"] = 2.6090373899600383e+01
 
 @testset "ARM_SGP" begin
     println("Running ARM_SGP...")

integration_tests/Bomex.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 8.2239735839064965e-02
-best_mse["updraft_area"] = 7.0282567418077099e+02
-best_mse["updraft_w"] = 8.7997106745562263e+01
-best_mse["updraft_qt"] = 5.9534668549222083e+00
-best_mse["updraft_thetal"] = 2.3059746455183749e+01
-best_mse["v_mean"] = 1.2344255867930939e+02
-best_mse["u_mean"] = 5.3486981021128493e+01
-best_mse["tke_mean"] = 3.2038984505851317e+01
-best_mse["temperature_mean"] = 3.1019627325766795e-05
-best_mse["ql_mean"] = 1.5630034030829634e+01
-best_mse["thetal_mean"] = 3.1601763199735002e-05
-best_mse["Hvar_mean"] = 3.4265143495072330e+01
-best_mse["QTvar_mean"] = 1.4543323453279942e+01
+best_mse["qt_mean"] = 8.4924728795344143e-02
+best_mse["updraft_area"] = 7.0195368138769777e+02
+best_mse["updraft_w"] = 8.7703050153771670e+01
+best_mse["updraft_qt"] = 5.9194060591321715e+00
+best_mse["updraft_thetal"] = 2.2932547498112353e+01
+best_mse["v_mean"] = 1.2350993295040250e+02
+best_mse["u_mean"] = 5.3487467706554988e+01
+best_mse["tke_mean"] = 3.2337733855931930e+01
+best_mse["temperature_mean"] = 3.1767291118329973e-05
+best_mse["ql_mean"] = 1.4250971600048450e+01
+best_mse["thetal_mean"] = 3.2420072527836003e-05
+best_mse["Hvar_mean"] = 6.2697869671981366e+01
+best_mse["QTvar_mean"] = 2.2154880539077912e+01
 
 @testset "Bomex" begin
     println("Running Bomex...")

integration_tests/DYCOMS_RF01.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 3.8001673417054252e-02
-best_mse["ql_mean"] = 8.4574058349035557e+00
-best_mse["updraft_area"] = 2.2323993812330733e+02
-best_mse["updraft_w"] = 3.4445142582864365e+00
-best_mse["updraft_qt"] = 1.3810735086427255e+00
-best_mse["updraft_thetal"] = 1.2761175177117382e+01
-best_mse["v_mean"] = 4.0030917247593457e+01
-best_mse["u_mean"] = 3.5747549053942961e+01
-best_mse["tke_mean"] = 1.4604034781584005e+01
-best_mse["temperature_mean"] = 3.8670907549444306e-06
-best_mse["thetal_mean"] = 5.7402823217894555e-06
-best_mse["Hvar_mean"] = 8.2258518882332137e+04
-best_mse["QTvar_mean"] = 6.2275442597879110e+03
+best_mse["qt_mean"] = 3.7999769022213942e-02
+best_mse["ql_mean"] = 8.4794189122776036e+00
+best_mse["updraft_area"] = 2.2364090186153348e+02
+best_mse["updraft_w"] = 3.4475661482890558e+00
+best_mse["updraft_qt"] = 1.3839010357798487e+00
+best_mse["updraft_thetal"] = 1.2733570833710051e+01
+best_mse["v_mean"] = 4.0030786550849442e+01
+best_mse["u_mean"] = 3.5747430896153816e+01
+best_mse["tke_mean"] = 1.4611287811858389e+01
+best_mse["temperature_mean"] = 3.8666969165057168e-06
+best_mse["thetal_mean"] = 5.7397807898569714e-06
+best_mse["Hvar_mean"] = 8.4325692999540202e+04
+best_mse["QTvar_mean"] = 6.3010689185932879e+03
 
 
 @testset "DYCOMS_RF01" begin

integration_tests/Nieuwstadt.jl

@@ -11,14 +11,14 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["updraft_area"] = 6.2270686101489321e+02
-best_mse["updraft_w"] = 3.9197923601593011e+01
-best_mse["updraft_thetal"] = 2.9892400636969811e+01
-best_mse["u_mean"] = 7.6422297657642343e+02
-best_mse["tke_mean"] = 6.9178066709527812e+01
-best_mse["temperature_mean"] = 1.6456374707932590e-05
-best_mse["thetal_mean"] = 1.6298071936614467e-05
-best_mse["Hvar_mean"] = 2.0218074864508648e+02
+best_mse["updraft_area"] = 5.0454151695959411e+02
+best_mse["updraft_w"] = 3.6872853179114287e+01
+best_mse["updraft_thetal"] = 2.9768390317748754e+01
+best_mse["u_mean"] = 7.6320346670292111e+02
+best_mse["tke_mean"] = 6.9114000131939633e+01
+best_mse["temperature_mean"] = 1.6216693070306151e-05
+best_mse["thetal_mean"] = 1.6067825070848235e-05
+best_mse["Hvar_mean"] = 2.0364798295577725e+02
 
 @testset "Nieuwstadt" begin
     println("Running Nieuwstadt...")

integration_tests/Rico.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 5.4596292300803628e-01
-best_mse["updraft_area"] = 1.6644272370900098e+03
-best_mse["updraft_w"] = 5.5689731617925463e+02
-best_mse["updraft_qt"] = 2.4446559332424545e+01
-best_mse["updraft_thetal"] = 6.5601697606639760e+01
-best_mse["v_mean"] = 1.0610811546766800e+02
-best_mse["u_mean"] = 1.1387870745418634e+02
-best_mse["tke_mean"] = 9.1707254861508807e+02
-best_mse["temperature_mean"] = 2.4294131149581795e-04
-best_mse["ql_mean"] = 1.3122452320829298e+04
-best_mse["thetal_mean"] = 1.8634680078510443e-04
-best_mse["Hvar_mean"] = 1.9466083203062441e+04
-best_mse["QTvar_mean"] = 1.0274739422376710e+05
+best_mse["qt_mean"] = 5.0707745192217146e-01
+best_mse["updraft_area"] = 1.6283285530615315e+03
+best_mse["updraft_w"] = 5.7196482018497807e+02
+best_mse["updraft_qt"] = 2.4565269303654755e+01
+best_mse["updraft_thetal"] = 6.5505090179439435e+01
+best_mse["v_mean"] = 1.0611825566231954e+02
+best_mse["u_mean"] = 1.1390377464452069e+02
+best_mse["tke_mean"] = 8.9525123248795580e+02
+best_mse["temperature_mean"] = 1.4787364980827869e-04
+best_mse["ql_mean"] = 1.7937532436604911e+02
+best_mse["thetal_mean"] = 1.4098220595492702e-04
+best_mse["Hvar_mean"] = 1.6499437998682095e+03
+best_mse["QTvar_mean"] = 7.3181087637066230e+02
 
 @testset "Rico" begin
     println("Running Rico...")

integration_tests/Soares.jl

@@ -11,16 +11,16 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 2.5843807856554074e-01
-best_mse["updraft_area"] = 8.8039804794002623e+02
-best_mse["updraft_w"] = 2.6632875471066519e+01
-best_mse["updraft_qt"] = 1.0651734450298703e+01
-best_mse["updraft_thetal"] = 2.1622611589959376e+01
-best_mse["u_mean"] = 3.9550058876376834e+03
-best_mse["tke_mean"] = 5.8683619102955831e+01
-best_mse["temperature_mean"] = 1.9035536936412582e-05
-best_mse["thetal_mean"] = 3.1601763199735002e-05
-best_mse["Hvar_mean"] = 2.4252112238315871e+02
+best_mse["qt_mean"] = 2.5795949502235077e-01
+best_mse["updraft_area"] = 9.2319813604790897e+02
+best_mse["updraft_w"] = 2.7143888970271533e+01
+best_mse["updraft_qt"] = 1.0608841949065969e+01
+best_mse["updraft_thetal"] = 2.1532908755155134e+01
+best_mse["u_mean"] = 3.9494974585811301e+03
+best_mse["tke_mean"] = 6.3308887592281366e+01
+best_mse["temperature_mean"] = 1.9251642797449680e-05
+best_mse["thetal_mean"] = 1.9375920993365417e-05
+best_mse["Hvar_mean"] = 2.4476641501276720e+02
 
 
 @testset "Soares" begin

integration_tests/TRMM_LBA.jl

@@ -11,19 +11,19 @@ include(joinpath("utils", "compute_mse.jl"))
 using .NameList
 
 best_mse = OrderedDict()
-best_mse["qt_mean"] = 3.9046069426558216e+00
-best_mse["updraft_area"] = 2.8331099674735910e+04
-best_mse["updraft_w"] = 6.5629758352987142e+02
-best_mse["updraft_qt"] = 2.5012703977470046e+01
-best_mse["updraft_thetal"] = 1.1072811439926345e+02
-best_mse["v_mean"] = 2.9406070727525389e+02
-best_mse["u_mean"] = 1.6903159369983448e+03
-best_mse["tke_mean"] = 2.8961135369556791e+03
-best_mse["temperature_mean"] = 8.3946608562868500e-04
-best_mse["ql_mean"] = 4.3631392987274085e+03
-best_mse["thetal_mean"] = 5.7967137448949413e-04
-best_mse["Hvar_mean"] = 1.3320049142979460e+04
-best_mse["QTvar_mean"] = 1.5870164338554601e+04
+best_mse["qt_mean"] = 3.9066401784413203e+00
+best_mse["updraft_area"] = 2.8333333642947131e+04
+best_mse["updraft_w"] = 6.5662893402745362e+02
+best_mse["updraft_qt"] = 2.5017304239563906e+01
+best_mse["updraft_thetal"] = 1.1072812435859557e+02
+best_mse["v_mean"] = 2.9403959926399466e+02
+best_mse["u_mean"] = 1.6903038672404982e+03
+best_mse["tke_mean"] = 2.8972134698112682e+03
+best_mse["temperature_mean"] = 8.3944025946955543e-04
+best_mse["ql_mean"] = 4.3199525166349194e+03
+best_mse["thetal_mean"] = 5.7974301596273598e-04
+best_mse["Hvar_mean"] = 1.3419697125412104e+04
+best_mse["QTvar_mean"] = 1.6238639092892809e+04
 
 @testset "TRMM_LBA" begin
     println("Running TRMM_LBA...")

src/Operators.jl

@@ -16,6 +16,20 @@ struct Extrapolate end
 function ∇f2c(f, grid::Grid, k::Int)
     return (f[k] - f[k - 1]) * grid.dzi
 end
+function ∇f2c(f_dual::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
+    if is_surface_bc_faces(grid, k)
+        return ∇f2c(f_dual, grid, BottomBCTag(), bottom)
+    elseif is_toa_bc_faces(grid, k)
+        return ∇f2c(f_dual, grid, TopBCTag(), top)
+    else
+        return ∇f2c(f_dual, grid, InteriorTag())
+    end
+end
+∇f2c(f::SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.dzi
+∇f2c(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * grid.dzi
+∇f2c(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value
+∇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 ∇c2f(f, grid::Grid, k::Int)
     return (f[k + 1] - f[k]) * grid.dzi
@@ -45,13 +59,12 @@ end
 ##### ∇(center data)
 #####
 
-c∇(f, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) = c∇(cut(f, k), grid, k; bottom, top)
+c∇(f, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) = c∇(cut(f, grid, k), grid, k; bottom, top)
 
 function c∇(f_cut::SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError())
-    gw = grid.gw
-    if k == gw # NOTE: 0-based indexing (k = 3 for 1-based indexing) bottom boundary
+    if is_surface_bc_centers(grid, k)
         return c∇(f_cut, grid, BottomBCTag(), bottom)
-    elseif k == grid.nzg - gw - 1 # NOTE: 0-based indexing
+    elseif is_toa_bc_centers(grid, k)
         return c∇(f_cut, grid, TopBCTag(), top)
     else
         return c∇(f_cut, grid, InteriorTag())
@@ -65,43 +78,41 @@ function c∇(f::SVector, grid::Grid, ::InteriorTag)
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 function c∇(f::SVector, grid::Grid, ::TopBCTag, bc::SetValue)
-    @assert length(f) == 3
+    @assert length(f) == 2
     # 2fb = cg+ci => cg = 2fb-ci
     f_dual⁺ = SVector(f[2], 2 * bc.value - f[2])
     f_dual⁻ = SVector(f[1], f[2])
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 function c∇(f::SVector, grid::Grid, ::BottomBCTag, bc::SetValue)
-    @assert length(f) == 3
+    @assert length(f) == 2
     # 2fb = cg+ci => cg = 2fb-ci
-    f_dual⁺ = SVector(f[2], f[3])
-    f_dual⁻ = SVector(2 * bc.value - f[2], f[2])
+    f_dual⁺ = SVector(f[1], f[2])
+    f_dual⁻ = SVector(2 * bc.value - f[1], f[1])
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 function c∇(f::SVector, grid::Grid, ::TopBCTag, bc::SetGradient)
-    @assert length(f) == 3
-    f_dual⁺ = SVector(f[2], f[3])
+    @assert length(f) == 2
     f_dual⁻ = SVector(f[1], f[2])
     return (bc.value + ∇_staggered(f_dual⁻, grid)) / 2
 end
 function c∇(f::SVector, grid::Grid, ::BottomBCTag, bc::SetGradient)
-    @assert length(f) == 3
-    f_dual⁺ = SVector(f[2], f[3])
-    f_dual⁻ = SVector(f[1], f[2])
+    @assert length(f) == 2
+    f_dual⁺ = SVector(f[1], f[2])
     return (∇_staggered(f_dual⁺, grid) + bc.value) / 2
 end
 function c∇(f::SVector, grid::Grid, ::TopBCTag, ::Extrapolate)
-    @assert length(f) == 3
+    @assert length(f) == 2
     # 2ci = cg+cii => cg = 2ci-cii. Note: f[3] not used
     f_dual⁺ = SVector(f[2], 2 * f[2] - f[1])
     f_dual⁻ = SVector(f[1], f[2])
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 function c∇(f::SVector, grid::Grid, ::BottomBCTag, ::Extrapolate)
-    @assert length(f) == 3
+    @assert length(f) == 2
     # 2ci = cg+cii => cg = 2ci-cii. Note: f[1] not used
-    f_dual⁺ = SVector(f[2], f[3])
-    f_dual⁻ = SVector(2 * f[2] - f[3], f[2])
+    f_dual⁺ = SVector(f[1], f[2])
+    f_dual⁻ = SVector(2 * f[1] - f[2], f[1])
     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
 end
 
@@ -114,25 +125,48 @@ function ∇_staggered(f::SVector, grid::Grid)
     return (f[2] - f[1]) * grid.dzi
 end
 
-# A 3-point field stencil
-function cut(f::AbstractVector, k::Int)
-    return SVector(f[k - 1], f[k], f[k + 1])
+# A 3-point field stencil for ordinary and updraft variables
+function cut(f::AbstractVector, grid, k::Int)
+    if is_surface_bc_centers(grid, k)
+        return SVector(f[k], f[k + 1])
+    elseif is_toa_bc_centers(grid, k)
+        return SVector(f[k - 1], f[k])
+    else
+        return SVector(f[k - 1], f[k], f[k + 1])
+    end
 end
-
-# A 3-point field stencil for updraft variables
-function cut(f::AbstractMatrix, k::Int, i_up::Int)
-    return SVector(f[i_up, k - 1], f[i_up, k], f[i_up, k + 1])
+function cut(f::AbstractMatrix, grid, k::Int, i_up::Int)
+    if is_surface_bc_centers(grid, k)
+        return SVector(f[i_up, k], f[i_up, k + 1])
+    elseif is_toa_bc_centers(grid, k)
+        return SVector(f[i_up, k - 1], f[i_up, k])
+    else
+        return SVector(f[i_up, k - 1], f[i_up, k], f[i_up, k + 1])
+    end
 end
 
-function ∇_collocated(f::SVector, grid::Grid)
-    @assert length(f) == 3
-    return (f[3] - f[1]) * 0.5 * grid.dzi
+# A 2-point field stencil for ordinary and updraft variables
+function dual_faces(f::AbstractVector, grid, k::Int)
+    if is_surface_bc_faces(grid, k)
+        return SVector(f[k])
+    elseif is_toa_bc_centers(grid, k)
+        return SVector(f[k - 1])
+    else
+        return SVector(f[k], f[k - 1])
+    end
 end
+function dual_faces(f::AbstractMatrix, grid, k::Int, i_up::Int)
+    if is_surface_bc_faces(grid, k)
+        return SVector(f[i_up, k])
+    elseif is_toa_bc_centers(grid, k)
+        return SVector(f[i_up, k - 1])
+    else
+        return SVector(f[i_up, k], f[i_up, k - 1])
+    end
+end
+
+is_surface_bc_centers(grid::Grid, k::Int) = k == grid.gw            # NOTE: 0-based indexing
+is_toa_bc_centers(grid::Grid, k::Int) = k == grid.nzg - grid.gw - 1 # NOTE: 0-based indexing
 
-# TODO: use this implementation
-# function ∇_collocated(f::SVector, grid::Grid)
-#     @assert length(f) == 3
-#     f_dual⁺ = SVector(f[2], f[3])
-#     f_dual⁻ = SVector(f[1], f[2])
-#     return (∇_staggered(f_dual⁺, grid) + ∇_staggered(f_dual⁻, grid)) / 2
-# end
+is_surface_bc_faces(grid::Grid, k::Int) = k == grid.gw            # NOTE: 0-based indexing
+is_toa_bc_faces(grid::Grid, k::Int) = k == grid.nzg - grid.gw - 1 # NOTE: 0-based indexing