Simplify source calculations for NoScatLWRTE solver. (#499)

Ref: earth-system-radiation/rte-rrtmgp#250

perf/benchmark.jl

@@ -49,7 +49,23 @@ println()
 # @suppress_out begin
 include(joinpath(root_dir, "test", "clear_sky_utils.jl"))
 context = ClimaComms.context()
-clear_sky(ClimaComms.context(), TwoStream, TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT; exfiltrate = true)
+
+toler_lw_noscat = Dict(Float64 => Float64(1e-4), Float32 => Float32(0.04))
+toler_lw_2stream = Dict(Float64 => Float64(4.5), Float32 => Float32(4.5))
+toler_sw = Dict(Float64 => Float64(1e-3), Float32 => Float32(0.04))
+
+clear_sky(
+    ClimaComms.context(),
+    TwoStream,
+    TwoStream,
+    TwoStreamLWRTE,
+    TwoStreamSWRTE,
+    VmrGM,
+    FT,
+    toler_lw_2stream,
+    toler_sw;
+    exfiltrate = true,
+)
 # end
 (; slv_lw, slv_sw, as, lookup_sw, lookup_lw) = Infiltrator.exfiltrated
 

src/optics/GrayOpticsKernels.jl

@@ -14,7 +14,7 @@ Computes optical properties for the longwave gray radiation problem.
 function compute_optical_props!(op::OneScalar, as::GrayAtmosphericState, sf::SourceLWNoScat, gcol::Int)
     nlay = AtmosphericStates.get_nlay(as)
     (; p_lay, p_lev, t_lay, t_lev, otp) = as
-    (; lay_source, lev_source_inc, lev_source_dec, sfc_source) = sf
+    (; lay_source, lev_source, sfc_source) = sf
     τ = op.τ
     FT = eltype(τ)
     sbc = FT(RP.Stefan(sf.param_set))
@@ -25,7 +25,8 @@ function compute_optical_props!(op::OneScalar, as::GrayAtmosphericState, sf::Sou
         t_lev_dec = t_lev[1, gcol]
         t_sfc = as.t_sfc[gcol]
         sfc_source[gcol] = sbc * (t_sfc * t_sfc * t_sfc * t_sfc) / FT(π)   # computing sfc_source
-
+        lev_src_inc_prev = FT(0)
+        lev_src_dec_prev = FT(0)
         for glay in 1:nlay
             # compute optical thickness
             p_lev_glayplus1 = p_lev[glay + 1, gcol]
@@ -37,10 +38,18 @@ function compute_optical_props!(op::OneScalar, as::GrayAtmosphericState, sf::Sou
             t_lev_inc = t_lev[glay + 1, gcol]
             t_lay_loc = t_lay[glay, gcol]
             lay_source[glay, gcol] = sbc * (t_lay_loc * t_lay_loc * t_lay_loc * t_lay_loc) / FT(π)   # computing lay_source
-            lev_source_inc[glay, gcol] = sbc * (t_lev_inc * t_lev_inc * t_lev_inc * t_lev_inc) / FT(π)
-            lev_source_dec[glay, gcol] = sbc * (t_lev_dec * t_lev_dec * t_lev_dec * t_lev_dec) / FT(π)
+            lev_src_inc = sbc * (t_lev_inc * t_lev_inc * t_lev_inc * t_lev_inc) / FT(π)
+            lev_src_dec = sbc * (t_lev_dec * t_lev_dec * t_lev_dec * t_lev_dec) / FT(π)
+            if glay == 1
+                lev_source[glay, gcol] = lev_src_dec
+            else
+                lev_source[glay, gcol] = sqrt(lev_src_inc_prev * lev_src_dec)
+            end
+            lev_src_dec_prev = lev_src_dec
+            lev_src_inc_prev = lev_src_inc
             t_lev_dec = t_lev_inc
         end
+        lev_source[nlay + 1, gcol] = lev_src_inc_prev
     end
     return nothing
 end

src/optics/Optics.jl

@@ -143,7 +143,7 @@ Computes optical properties for the longwave problem.
     nlay = AtmosphericStates.get_nlay(as)
     (; vmr) = as
     (; t_planck) = lkp.planck
-    (; lay_source, lev_source_inc, lev_source_dec, sfc_source) = sf
+    (; lay_source, lev_source, sfc_source) = sf
     @inbounds begin
         t_sfc = as.t_sfc[gcol]
         ibnd = lkp.band_data.major_gpt2bnd[igpt]
@@ -152,7 +152,10 @@ Computes optical properties for the longwave problem.
         t_lev_col = view(as.t_lev, :, gcol)
         τ = view(op.τ, :, gcol)
 
+        lev_src_inc_prev = zero(t_sfc)
+        lev_src_dec_prev = zero(t_sfc)
         t_lev_dec = t_lev_col[1]
+
         for glay in 1:nlay
             col_dry, p_lay, t_lay = as_layerdata[1, glay], as_layerdata[2, glay], as_layerdata[3, glay]
             # compute gas optics
@@ -161,13 +164,19 @@ Computes optical properties for the longwave problem.
             t_lev_inc = t_lev_col[glay + 1]
 
             lay_source[glay, gcol] = interp1d(t_lay, t_planck, totplnk) * planckfrac
-            lev_source_inc[glay, gcol] = interp1d(t_lev_inc, t_planck, totplnk) * planckfrac
-            lev_source_dec[glay, gcol] = interp1d(t_lev_dec, t_planck, totplnk) * planckfrac
+            lev_src_inc = interp1d(t_lev_inc, t_planck, totplnk) * planckfrac
+            lev_src_dec = interp1d(t_lev_dec, t_planck, totplnk) * planckfrac
             if glay == 1
                 sfc_source[gcol] = interp1d(t_sfc, t_planck, totplnk) * planckfrac
+                lev_source[glay, gcol] = lev_src_dec
+            else
+                lev_source[glay, gcol] = sqrt(lev_src_inc_prev * lev_src_dec)
             end
+            lev_src_dec_prev = lev_src_dec
+            lev_src_inc_prev = lev_src_inc
             t_lev_dec = t_lev_inc
         end
+        lev_source[nlay + 1, gcol] = lev_src_inc_prev
     end
     return nothing
 end

src/optics/Sources.jl

@@ -24,36 +24,28 @@ and at the surface for no scattering calculations
 
 $(DocStringExtensions.FIELDS)
 """
-struct SourceLWNoScat{S, D, V, PS} <: AbstractSourceLW
+struct SourceLWNoScat{S, D, PS} <: AbstractSourceLW
     "Parameter set"
     param_set::PS
     "Surface source `[W/m2]` `(ncol)`"
     sfc_source::S
-    "storage for `lay_source`, `lev_source_inc` and `lev_source_dec` `(3, nlay, ncol)`"
-    layerdata::D
     "Planck source at layer average temperature `[W/m2]` `(nlay, ncol)`"
-    lay_source::V
-    "Planck source at layer edge in increasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
-    lev_source_inc::V
-    "Planck source at layer edge in decreasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
-    lev_source_dec::V
+    lay_source::D
+    "Planck level source at layer edges `[W/m2]` `(nlay+1, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
+    lev_source::D
 end
 Adapt.@adapt_structure SourceLWNoScat
 
 function SourceLWNoScat(param_set::RP.ARP, ::Type{FT}, ::Type{DA}, nlay::Int, ncol::Int) where {FT <: AbstractFloat, DA}
     sfc_source = DA{FT, 1}(undef, ncol)
-    layerdata = DA{FT, 3}(undef, 3, nlay, ncol)
-    lay_source = view(layerdata, 1, :, :)
-    lev_source_inc = view(layerdata, 2, :, :)
-    lev_source_dec = view(layerdata, 3, :, :)
+    lay_source = DA{FT, 2}(undef, nlay, ncol)
+    lev_source = DA{FT, 2}(undef, nlay + 1, ncol)
 
-    return SourceLWNoScat{typeof(sfc_source), typeof(layerdata), typeof(lay_source), typeof(param_set)}(
+    return SourceLWNoScat{typeof(sfc_source), typeof(lay_source), typeof(param_set)}(
         param_set,
         sfc_source,
-        layerdata,
         lay_source,
-        lev_source_inc,
-        lev_source_dec,
+        lev_source,
     )
 end
 

src/rte/longwave1scalar.jl

@@ -122,7 +122,7 @@ Transport for no-scattering longwave problem.
 )
     # setting references
     (; sfc_source) = src_lw
-    (; lay_source, lev_source_inc, lev_source_dec) = src_lw
+    (; lay_source, lev_source) = src_lw
     (; sfc_emis, inc_flux) = bcs_lw
     (; flux_up, flux_dn) = flux
 
@@ -147,9 +147,10 @@ Transport for no-scattering longwave problem.
     @inbounds while ilev ≥ 1
         τ_loc = τ[ilev, gcol] * Ds[1]
         trans = exp(-τ_loc)
-        lay_src, lev_src_dec = lay_source[ilev, gcol], lev_source_dec[ilev, gcol]
+        lay_src = lay_source[ilev, gcol]
         intensity_dn_ilev =
-            trans * intensity_dn_ilevplus1 + lw_noscat_source_dn(lev_src_dec, lay_src, τ_loc, trans, τ_thresh)
+            trans * intensity_dn_ilevplus1 +
+            lw_noscat_source_dn(lev_source[ilev, gcol], lay_src, τ_loc, trans, τ_thresh)
         intensity_dn_ilevplus1 = intensity_dn_ilev
         flux_dn[ilev, gcol] = intensity_dn_ilev * intensity_to_flux
         ilev -= 1
@@ -165,9 +166,10 @@ Transport for no-scattering longwave problem.
     @inbounds for ilev in 2:(nlay + 1)
         τ_loc = τ[ilev - 1, gcol] * Ds[1]
         trans = exp(-τ_loc)
-        lay_src, lev_src_inc = lay_source[ilev - 1, gcol], lev_source_inc[ilev - 1, gcol]
+        lay_src = lay_source[ilev - 1, gcol]
         intensity_up_ilev =
-            trans * intensity_up_ilevminus1 + lw_noscat_source_up(lev_src_inc, lay_src, τ_loc, trans, τ_thresh)
+            trans * intensity_up_ilevminus1 +
+            lw_noscat_source_up(lev_source[ilev, gcol], lay_src, τ_loc, trans, τ_thresh)
         intensity_up_ilevminus1 = intensity_up_ilev
         flux_up[ilev, gcol] = intensity_up_ilev * intensity_to_flux
     end

test/clear_sky.jl

@@ -4,10 +4,24 @@ include("clear_sky_utils.jl")
 
 context = ClimaComms.context()
 
+toler_lw_noscat = Dict(Float64 => Float64(1e-4), Float32 => Float32(0.05))
+toler_lw_2stream = Dict(Float64 => Float64(4.5), Float32 => Float32(4.5))
+toler_sw = Dict(Float64 => Float64(1e-3), Float32 => Float32(0.04))
+
 @testset "clear-sky: NoScatLWRTE (OneScalar optics) & TwoStreamSWRTE (TwoStream optics)" begin
-    @time clear_sky(context, OneScalar, TwoStream, NoScatLWRTE, TwoStreamSWRTE, VmrGM, FT)
+    @time clear_sky(context, OneScalar, TwoStream, NoScatLWRTE, TwoStreamSWRTE, VmrGM, FT, toler_lw_noscat, toler_sw)
 end
 
 @testset "clear-sky: TwoStreamLWRTE (TwoStream optics) & TwoStreamSWRTE (TwoStream optics)" begin
-    @time clear_sky(context, TwoStream, TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT)
+    @time clear_sky(
+        context,
+        TwoStream,
+        TwoStream,
+        TwoStreamLWRTE,
+        TwoStreamSWRTE,
+        VmrGM,
+        FT,
+        toler_lw_2stream,
+        toler_sw,
+    )
 end

test/clear_sky_utils.jl

@@ -33,7 +33,9 @@ function clear_sky(
     ::Type{SLVLW},
     ::Type{SLVSW},
     ::Type{VMR},
-    ::Type{FT};
+    ::Type{FT},
+    toler_lw,
+    toler_sw;
     exfiltrate = false,
 ) where {FT, OPLW, OPSW, SLVLW, SLVSW, VMR}
     overrides = (; grav = 9.80665, molmass_dryair = 0.028964, molmass_water = 0.018016)
@@ -198,15 +200,12 @@ function clear_sky(
     println("L∞ error in flux_net          = $max_err_flux_net_sw")
     println("L∞ relative error in flux_net = $(max_rel_err_flux_net_sw * 100) %\n")
 
-    toler_lw = Dict(Float64 => Float64(1e-4), Float32 => Float32(0.04))
-    toler_sw = Dict(Float64 => Float64(1e-3), Float32 => Float32(0.04))
-    # TODO: The reference results for the longwave solver are generated using a non-scattering solver with rescaling,
+    # Note: The reference results for the longwave solver are generated using a non-scattering solver with rescaling,
     # which differ from the results generated by the TwoStream currently used. Due to the above difference, 
-    # the longwave tests currently have higher error and are changed to "broken" state.
-    # These will be fixed once the Non-scattering solver with rescaling is added.
-    @test_broken max_err_flux_up_lw ≤ toler_lw[FT]
-    @test_broken max_err_flux_dn_lw ≤ toler_lw[FT]
-    @test_broken max_err_flux_net_lw ≤ toler_lw[FT]
+    # the longwave tests using currently have higher error.
+    @test max_err_flux_up_lw ≤ toler_lw[FT]
+    @test max_err_flux_dn_lw ≤ toler_lw[FT]
+    @test max_err_flux_net_lw ≤ toler_lw[FT]
 
     @test max_err_flux_up_sw ≤ toler_sw[FT]
     @test max_err_flux_dn_sw ≤ toler_sw[FT]

test/runtests.jl

@@ -13,10 +13,16 @@ printstyled("\n\nRRTMGP clear-sky tests \n", color = color1)
 printstyled("==========================\n\n", color = color1)
 @testset "RRTMGP clear-sky tests" begin
     context = ClimaComms.context()
+
+    toler_lw_noscat = Dict(Float64 => Float64(1e-4), Float32 => Float32(0.05))
+    toler_lw_2stream = Dict(Float64 => Float64(4.5), Float32 => Float32(4.5))
+    toler_sw = Dict(Float64 => Float64(1e-3), Float32 => Float32(0.04))
+
     include("clear_sky_utils.jl")
+
     for FT in (Float32, Float64)
-        clear_sky(context, OneScalar, TwoStream, NoScatLWRTE, TwoStreamSWRTE, VmrGM, FT)
-        clear_sky(context, TwoStream, TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT)
+        clear_sky(context, OneScalar, TwoStream, NoScatLWRTE, TwoStreamSWRTE, VmrGM, FT, toler_lw_noscat, toler_sw)
+        clear_sky(context, TwoStream, TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT, toler_lw_2stream, toler_sw)
     end
 end