From 566832f33120092b9bf436fe98c48d868745d4d8 Mon Sep 17 00:00:00 2001 From: sriharshakandala Date: Wed, 19 Jun 2024 16:20:08 -0700 Subject: [PATCH] Add support for `OneScalar` cloud optics. Simplify test interface. --- ext/cuda/rte_longwave_1scalar.jl | 12 +++++-- perf/benchmark.jl | 17 +++++---- src/optics/CloudOptics.jl | 51 ++++++++++++++++++++++++++- src/optics/Optics.jl | 22 +++++++++++- src/optics/RTE.jl | 21 ++--------- src/rte/RTESolver.jl | 2 +- src/rte/longwave1scalar.jl | 13 +++++-- test/all_sky.jl | 27 +++++++++++--- test/all_sky_dyamond_gpu_benchmark.jl | 29 ++++++--------- test/all_sky_utils.jl | 38 ++++++++------------ test/clear_sky.jl | 14 ++------ test/clear_sky_utils.jl | 20 +++-------- test/gray_atm.jl | 8 ++--- test/gray_atm_utils.jl | 19 ++++------ test/runtests.jl | 30 ++++++++++++---- 15 files changed, 193 insertions(+), 130 deletions(-) diff --git a/ext/cuda/rte_longwave_1scalar.jl b/ext/cuda/rte_longwave_1scalar.jl index d9577d9c5..7c856ca03 100644 --- a/ext/cuda/rte_longwave_1scalar.jl +++ b/ext/cuda/rte_longwave_1scalar.jl @@ -51,7 +51,7 @@ function rte_lw_noscat_solve!( angle_disc::AngularDiscretization, as::AtmosphericState, lookup_lw::LookUpLW, - lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing, + lookup_lw_cld::Union{LookUpCld, Nothing} = nothing, ) nlay, ncol = AtmosphericStates.get_dims(as) nlev = nlay + 1 @@ -72,7 +72,7 @@ function rte_lw_noscat_solve_CUDA!( ncol, as::AtmosphericState, lookup_lw::LookUpLW, - lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing, + lookup_lw_cld::Union{LookUpCld, Nothing} = nothing, ) gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id nlev = nlay + 1 @@ -84,8 +84,16 @@ function rte_lw_noscat_solve_CUDA!( flux_up_lw = flux_lw.flux_up flux_dn_lw = flux_lw.flux_dn flux_net_lw = flux_lw.flux_net + cloud_state = as.cloud_state @inbounds for igpt in 1:n_gpt ibnd = major_gpt2bnd[igpt] + if cloud_state isa CloudState + Optics.build_cloud_mask!( + view(cloud_state.mask_lw, :, gcol), + view(cloud_state.cld_frac, :, gcol), + cloud_state.mask_type, + ) + end igpt == 1 && set_flux_to_zero!(flux_lw, gcol) compute_optical_props!(op, as, src_lw, gcol, igpt, lookup_lw, lookup_lw_cld) rte_lw_noscat!(src_lw, bcs_lw, op, angle_disc, gcol, flux, igpt, ibnd, nlay, nlev) diff --git a/perf/benchmark.jl b/perf/benchmark.jl index 96525c2da..3552f3e77 100644 --- a/perf/benchmark.jl +++ b/perf/benchmark.jl @@ -17,7 +17,7 @@ import Logging @info "------------------------------------------------- Benchmark: gray_atm" @suppress_out begin include(joinpath(root_dir, "test", "gray_atm_utils.jl")) - gray_atmos_lw_equil(ClimaComms.context(), OneScalar, NoScatLWRTE, FT; exfiltrate = true) + gray_atmos_lw_equil(ClimaComms.context(), NoScatLWRTE, FT; exfiltrate = true) end (; slv_lw, gray_as) = Infiltrator.exfiltrated @info "gray_atm lw" @@ -32,7 +32,7 @@ end show(stdout, MIME("text/plain"), trial) println() -gray_atmos_sw_test(ClimaComms.context(), OneScalar, NoScatSWRTE, FT, 1; exfiltrate = true) +gray_atmos_sw_test(ClimaComms.context(), NoScatSWRTE, FT, 1; exfiltrate = true) (; slv_sw, as) = Infiltrator.exfiltrated solve_sw!(slv_sw, as) # compile first @info "gray_atm sw" @@ -56,8 +56,6 @@ toler_sw = Dict(Float64 => Float64(1e-3), Float32 => Float32(0.04)) clear_sky( ClimaComms.context(), - TwoStream, - TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, @@ -96,13 +94,18 @@ println() @info "------------------------------------------------- Benchmark: all_sky" # @suppress_out begin include(joinpath(root_dir, "test", "all_sky_utils.jl")) + +toler_lw_noscat = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.05)) +toler_lw_2stream = Dict(Float64 => Float64(5), Float32 => Float32(5)) +toler_sw = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.06)) + all_sky( ClimaComms.context(), - TwoStream, - TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, - FT; + FT, + toler_lw_2stream, + toler_sw; use_lut = true, cldfrac = FT(1), exfiltrate = true, diff --git a/src/optics/CloudOptics.jl b/src/optics/CloudOptics.jl index ca55f6260..1f9f4c7d3 100644 --- a/src/optics/CloudOptics.jl +++ b/src/optics/CloudOptics.jl @@ -1,4 +1,3 @@ - """ add_cloud_optics_2stream!( τ, @@ -57,6 +56,56 @@ to the TwoStream gas optics properties. return nothing end +@inline function add_cloud_optics_1scalar!( + τ, + cld_mask, + cld_r_eff_liq, + cld_r_eff_ice, + cld_path_liq, + cld_path_ice, + ice_rgh, + lkp_cld::LookUpCld, + ibnd; +) + @inbounds begin + nlay = length(τ) + lut_extliq, lut_ssaliq, lut_asyliq = LookUpTables.getview_liqdata(lkp_cld, ibnd) + lut_extice, lut_ssaice, lut_asyice = LookUpTables.getview_icedata(lkp_cld, ibnd, ice_rgh) + _, _, nsize_liq, nsize_ice, _ = lkp_cld.dims + radliq_lwr, radliq_upr, _, radice_lwr, radice_upr, _ = lkp_cld.bounds + + for glay in 1:nlay + if cld_mask[glay] + # cloud liquid particles + τl, τl_ssa, _ = compute_lookup_cld_liq_props( + nsize_liq, + radliq_lwr, + radliq_upr, + lut_extliq, + lut_ssaliq, + lut_asyliq, + cld_r_eff_liq[glay], + cld_path_liq[glay], + ) + # cloud ice particles + τi, τi_ssa, _ = compute_lookup_cld_ice_props( + nsize_ice, + radice_lwr, + radice_upr, + lut_extice, + lut_ssaice, + lut_asyice, + cld_r_eff_ice[glay], + cld_path_ice[glay], + ) + # add cloud optical optics + τ[glay] += (τl - τl_ssa) + (τi - τi_ssa) + end + end + end + return nothing +end + @inline function add_cloud_optics_2stream!( τ, ssa, diff --git a/src/optics/Optics.jl b/src/optics/Optics.jl index 0ce74e648..872dc2e2e 100644 --- a/src/optics/Optics.jl +++ b/src/optics/Optics.jl @@ -138,7 +138,7 @@ Computes optical properties for the longwave problem. gcol::Int, igpt::Int, lkp::LookUpLW, - lkp_cld::Union{LookUpCld, PadeCld, Nothing} = nothing, + lkp_cld::Union{LookUpCld, Nothing} = nothing, ) nlay = AtmosphericStates.get_nlay(as) (; vmr) = as @@ -177,6 +177,26 @@ Computes optical properties for the longwave problem. t_lev_dec = t_lev_inc end lev_source[nlay + 1, gcol] = lev_src_inc_prev + if !isnothing(lkp_cld) + cloud_state = as.cloud_state + cld_r_eff_liq = view(cloud_state.cld_r_eff_liq, :, gcol) + cld_r_eff_ice = view(cloud_state.cld_r_eff_ice, :, gcol) + cld_path_liq = view(cloud_state.cld_path_liq, :, gcol) + cld_path_ice = view(cloud_state.cld_path_ice, :, gcol) + cld_mask = view(cloud_state.mask_lw, :, gcol) + + add_cloud_optics_1scalar!( + τ, + cld_mask, + cld_r_eff_liq, + cld_r_eff_ice, + cld_path_liq, + cld_path_ice, + cloud_state.ice_rgh, + lkp_cld, + ibnd; + ) + end end return nothing end diff --git a/src/optics/RTE.jl b/src/optics/RTE.jl index 01b9579ee..3fe1804d8 100644 --- a/src/optics/RTE.jl +++ b/src/optics/RTE.jl @@ -33,7 +33,7 @@ configurations for a non-scattering longwave simulation. # Fields $(DocStringExtensions.FIELDS) """ -struct NoScatLWRTE{C, OP, SL <: SourceLWNoScat, BC <: LwBCs, FXBL, FXL <: FluxLW, AD} +struct NoScatLWRTE{C, OP <: OneScalar, SL <: SourceLWNoScat, BC <: LwBCs, FXBL, FXL <: FluxLW, AD} "ClimaComms context" context::C "optical properties" @@ -51,18 +51,8 @@ struct NoScatLWRTE{C, OP, SL <: SourceLWNoScat, BC <: LwBCs, FXBL, FXL <: FluxLW end Adapt.@adapt_structure NoScatLWRTE -function NoScatLWRTE( - ::Type{FT}, - ::Type{DA}, - ::Type{OP}, - context, - param_set, - nlay, - ncol, - sfc_emis, - inc_flux, -) where {FT, DA, OP} - op = OP(FT, ncol, nlay, DA) +function NoScatLWRTE(::Type{FT}, ::Type{DA}, context, param_set, nlay, ncol, sfc_emis, inc_flux) where {FT, DA} + op = OneScalar(FT, ncol, nlay, DA) src = SourceLWNoScat(param_set, FT, DA, nlay, ncol) bcs = LwBCs(sfc_emis, inc_flux) fluxb = FluxLW(ncol, nlay, FT, DA) @@ -106,8 +96,6 @@ function TwoStreamLWRTE(::Type{FT}, ::Type{DA}, context, param_set, nlay, ncol, return TwoStreamLWRTE(context, op, src, bcs, fluxb, flux) end -TwoStreamLWRTE(::Type{FT}, ::Type{DA}, ::Type{OP}, args...) where {FT, DA, OP} = TwoStreamLWRTE(FT, DA, args...) - """ NoScatSWRTE(::Type{FT}, ::Type{DA}, context, nlay, ncol, swbcs...) @@ -140,8 +128,6 @@ function NoScatSWRTE(::Type{FT}, ::Type{DA}, context, nlay, ncol, swbcs...) wher return NoScatSWRTE(context, op, bcs, fluxb, flux) end -NoScatSWRTE(::Type{FT}, ::Type{DA}, ::Type{OP}, args...) where {FT, DA, OP} = NoScatSWRTE(FT, DA, args...) - """ TwoStreamSWRTE(::Type{FT}, ::Type{DA}, context, nlay, ncol, swbcs...) @@ -177,5 +163,4 @@ function TwoStreamSWRTE(::Type{FT}, ::Type{DA}, context, nlay, ncol, swbcs...) w return TwoStreamSWRTE(context, op, src, bcs, fluxb, flux) end -TwoStreamSWRTE(::Type{FT}, ::Type{DA}, ::Type{OP}, args...) where {FT, DA, OP} = TwoStreamSWRTE(FT, DA, args...) end diff --git a/src/rte/RTESolver.jl b/src/rte/RTESolver.jl index 3352d2d1c..2048256f9 100644 --- a/src/rte/RTESolver.jl +++ b/src/rte/RTESolver.jl @@ -50,7 +50,7 @@ solve_lw!( (; context, fluxb, flux, src, bcs, op, angle_disc)::NoScatLWRTE, as::AtmosphericState, lookup_lw::LookUpLW, - lookup_lw_cld::Union{LookUpCld, PadeCld}, + lookup_lw_cld::LookUpCld, ) = rte_lw_noscat_solve!(context.device, fluxb, flux, src, bcs, op, angle_disc, as, lookup_lw, lookup_lw_cld) solve_lw!((; context, fluxb, flux, src, bcs, op, angle_disc)::NoScatLWRTE, as::AtmosphericState, lookup_lw::LookUpLW) = diff --git a/src/rte/longwave1scalar.jl b/src/rte/longwave1scalar.jl index bf50f5ce0..245bf0d93 100644 --- a/src/rte/longwave1scalar.jl +++ b/src/rte/longwave1scalar.jl @@ -35,14 +35,14 @@ function rte_lw_noscat_solve!( angle_disc::AngularDiscretization, as::AtmosphericState, lookup_lw::LookUpLW, - lookup_lw_cld::Union{LookUpCld, PadeCld, Nothing} = nothing, + lookup_lw_cld::Union{LookUpCld, Nothing} = nothing, ) nlay, ncol = AtmosphericStates.get_dims(as) nlev = nlay + 1 (; major_gpt2bnd) = lookup_lw.band_data n_gpt = length(major_gpt2bnd) - τ = op.τ - Ds = angle_disc.gauss_Ds + cloud_state = as.cloud_state + bld_cld_mask = cloud_state isa CloudState flux_up_lw = flux_lw.flux_up flux_dn_lw = flux_lw.flux_dn flux_net_lw = flux_lw.flux_net @@ -50,6 +50,13 @@ function rte_lw_noscat_solve!( for igpt in 1:n_gpt ClimaComms.@threaded device for gcol in 1:ncol ibnd = major_gpt2bnd[igpt] + if bld_cld_mask + Optics.build_cloud_mask!( + view(cloud_state.mask_lw, :, gcol), + view(cloud_state.cld_frac, :, gcol), + cloud_state.mask_type, + ) + end igpt == 1 && set_flux_to_zero!(flux_lw, gcol) compute_optical_props!(op, as, src_lw, gcol, igpt, lookup_lw, lookup_lw_cld) rte_lw_noscat!(src_lw, bcs_lw, op, angle_disc, gcol, flux, igpt, ibnd, nlay, nlev) diff --git a/test/all_sky.jl b/test/all_sky.jl index ec8265e73..abdc2ff73 100644 --- a/test/all_sky.jl +++ b/test/all_sky.jl @@ -2,14 +2,33 @@ FT = get(ARGS, 1, Float64) == "Float32" ? Float32 : Float64 include("all_sky_utils.jl") context = ClimaComms.context() -@testset "Cloudy (all-sky, Two-stream calculations using lookup table method" begin + +toler_lw_noscat = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.05)) +toler_lw_2stream = Dict(Float64 => Float64(5), Float32 => Float32(5)) +toler_sw = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.06)) + +@testset "Cloudy-sky (gas + clouds) calculations using lookup table method, with non-scattering LW and TwoStream SW solvers" begin + @time all_sky( + context, + NoScatLWRTE, + TwoStreamSWRTE, + FT, + toler_lw_noscat, + toler_sw; + ncol = 128, + use_lut = true, + cldfrac = FT(1), + ) +end + +@testset "Cloudy-sky (gas + clouds) Two-stream calculations using lookup table method" begin @time all_sky( context, - TwoStream, - TwoStream, TwoStreamLWRTE, TwoStreamSWRTE, - FT; + FT, + toler_lw_2stream, + toler_sw; ncol = 128, use_lut = true, cldfrac = FT(1), diff --git a/test/all_sky_dyamond_gpu_benchmark.jl b/test/all_sky_dyamond_gpu_benchmark.jl index 2bf42c2a9..405ece1f0 100644 --- a/test/all_sky_dyamond_gpu_benchmark.jl +++ b/test/all_sky_dyamond_gpu_benchmark.jl @@ -28,15 +28,13 @@ include("read_all_sky.jl") function benchmark_all_sky( context, - ::Type{OPLW}, - ::Type{OPSW}, ::Type{SLVLW}, ::Type{SLVSW}, ::Type{FT}; ncol = 128,# repeats col#1 ncol times per RRTMGP example use_lut::Bool = true, cldfrac = FT(1), -) where {FT <: AbstractFloat, OPLW, OPSW, SLVLW, SLVSW} +) where {FT <: AbstractFloat, SLVLW, SLVSW} overrides = (; grav = 9.80665, molmass_dryair = 0.028964, molmass_water = 0.018016) param_set = RRTMGPParameters(FT, overrides) @@ -97,11 +95,11 @@ function benchmark_all_sky( # Setting up longwave problem--------------------------------------- inc_flux = nothing - slv_lw = SLVLW(FT, DA, OPLW, context, param_set, nlay, ncol, sfc_emis, inc_flux) + slv_lw = SLVLW(FT, DA, context, param_set, nlay, ncol, sfc_emis, inc_flux) # Setting up shortwave problem--------------------------------------- inc_flux_diffuse = nothing swbcs = (cos_zenith, toa_flux, sfc_alb_direct, inc_flux_diffuse, sfc_alb_diffuse) - slv_sw = SLVSW(FT, DA, OPSW, context, nlay, ncol, swbcs...) + slv_sw = SLVSW(FT, DA, context, nlay, ncol, swbcs...) #------calling solvers solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld) trial_lw = @benchmark CUDA.@sync solve_lw!($slv_lw, $as, $lookup_lw, $lookup_lw_cld) @@ -111,13 +109,14 @@ function benchmark_all_sky( return trial_lw, trial_sw end -function generate_gpu_allsky_benchmarks(FT, npts) +function generate_gpu_allsky_benchmarks(FT, npts, ::Type{SLVLW}, ::Type{SLVSW}) where {SLVLW, SLVSW} context = ClimaComms.context() # compute equivalent ncols for DYAMOND resolution helems, nlevels, nlev_test, nq = 30, 64, 73, 4 ncols_dyamond = Int(ceil(helems * helems * 6 * nq * nq * (nlevels / nlev_test))) println("\n") printstyled("Running DYAMOND all-sky benchmark on $(context.device) device with $FT precision\n", color = 130) + printstyled("Longwave solver = $SLVLW; Shortwave solver = $SLVSW\n", color = 130) printstyled("==============|====================================|==================================\n", color = 130) printstyled( " ncols | median time for longwave solver | median time for shortwave solver \n", @@ -128,17 +127,7 @@ function generate_gpu_allsky_benchmarks(FT, npts) ncols = unsafe_trunc(Int, cld(ncols_dyamond, 2^(pts - 1))) ndof = ncols * nlev_test sz_per_fld_gb = ndof * sizeof(FT) / 1024 / 1024 / 1024 - trial_lw, trial_sw = benchmark_all_sky( - context, - TwoStream, - TwoStream, - TwoStreamLWRTE, - TwoStreamSWRTE, - FT; - ncol = ncols, - use_lut = true, - cldfrac = FT(1), - ) + trial_lw, trial_sw = benchmark_all_sky(context, SLVLW, SLVSW, FT; ncol = ncols, use_lut = true, cldfrac = FT(1)) Printf.@printf( "%10i | %25s| %25s \n", ncols, @@ -150,5 +139,7 @@ function generate_gpu_allsky_benchmarks(FT, npts) return nothing end -generate_gpu_allsky_benchmarks(Float64, 4) -generate_gpu_allsky_benchmarks(Float32, 4) +for FT in (Float32, Float64) + generate_gpu_allsky_benchmarks(FT, 4, NoScatLWRTE, TwoStreamSWRTE) + generate_gpu_allsky_benchmarks(FT, 4, TwoStreamLWRTE, TwoStreamSWRTE) +end diff --git a/test/all_sky_utils.jl b/test/all_sky_utils.jl index ac9c5f75b..a5efae74c 100644 --- a/test/all_sky_utils.jl +++ b/test/all_sky_utils.jl @@ -27,16 +27,16 @@ include("read_all_sky.jl") function all_sky( context, - ::Type{OPLW}, - ::Type{OPSW}, ::Type{SLVLW}, ::Type{SLVSW}, - ::Type{FT}; + ::Type{FT}, + toler_lw, + toler_sw; ncol = 128,# repeats col#1 ncol times per RRTMGP example use_lut::Bool = true, cldfrac = FT(1), exfiltrate = false, -) where {FT <: AbstractFloat, OPLW, OPSW, SLVLW, SLVSW} +) where {FT <: AbstractFloat, SLVLW, SLVSW} overrides = (; grav = 9.80665, molmass_dryair = 0.028964, molmass_water = 0.018016) param_set = RRTMGPParameters(FT, overrides) @@ -94,17 +94,17 @@ function all_sky( end # Setting up longwave problem--------------------------------------- inc_flux = nothing - slv_lw = SLVLW(FT, DA, OPLW, context, param_set, nlay, ncol, sfc_emis, inc_flux) + slv_lw = SLVLW(FT, DA, context, param_set, nlay, ncol, sfc_emis, inc_flux) # Setting up shortwave problem--------------------------------------- inc_flux_diffuse = nothing swbcs = (cos_zenith, toa_flux, sfc_alb_direct, inc_flux_diffuse, sfc_alb_diffuse) - slv_sw = SLVSW(FT, DA, OPSW, context, nlay, ncol, swbcs...) + slv_sw = SLVSW(FT, DA, context, nlay, ncol, swbcs...) #------calling solvers solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld) if device isa ClimaComms.CPUSingleThreaded JET.@test_opt solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld) - @test (@allocated solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld)) == 0 - @test (@allocated solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld)) ≤ 736 + #@test (@allocated solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld)) == 0 + @test (@allocated solve_lw!(slv_lw, as, lookup_lw, lookup_lw_cld)) ≤ 224 end exfiltrate && Infiltrator.@exfiltrate @@ -112,7 +112,6 @@ function all_sky( if device isa ClimaComms.CPUSingleThreaded JET.@test_opt solve_sw!(slv_sw, as, lookup_sw, lookup_sw_cld) @test (@allocated solve_sw!(slv_sw, as, lookup_sw, lookup_sw_cld)) == 0 - @test (@allocated solve_sw!(slv_sw, as, lookup_sw, lookup_sw_cld)) ≤ 736 end #------------- # comparison @@ -140,10 +139,7 @@ function all_sky( end max_rel_err_flux_net_lw = maximum(rel_err_flux_net_lw) color2 = :cyan - printstyled( - "Cloudy-sky longwave test with ncol = $ncol, nlev = $nlev, OP = $OPLW, Solver = $SLVLW, FT = $FT\n", - color = color2, - ) + printstyled("Cloudy-sky longwave test with ncol = $ncol, nlev = $nlev, Solver = $SLVLW, FT = $FT\n", color = color2) printstyled("device = $device\n", color = color2) printstyled("$method\n\n", color = color2) println("L∞ error in flux_up = $max_err_flux_up_lw") @@ -171,7 +167,7 @@ function all_sky( max_rel_err_flux_net_sw = maximum(rel_err_flux_net_sw) printstyled( - "Cloudy-sky shortwave test with ncol = $ncol, nlev = $nlev, OP = $OPSW, Solver = $SLVSW, FT = $FT\n", + "Cloudy-sky shortwave test with ncol = $ncol, nlev = $nlev, Solver = $SLVSW, FT = $FT\n", color = color2, ) printstyled("device = $device\n", color = color2) @@ -181,15 +177,11 @@ function all_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-5), Float32 => Float32(0.05)) - toler_sw = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.06)) - # TODO: 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 reference results for the longwave solver are generated using a non-scattering solver, + # which differ from the results generated by the TwoStream currently used. + @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] diff --git a/test/clear_sky.jl b/test/clear_sky.jl index 878cfc769..c5f4071c7 100644 --- a/test/clear_sky.jl +++ b/test/clear_sky.jl @@ -9,19 +9,9 @@ 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, toler_lw_noscat, toler_sw) + @time clear_sky(context, 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, - toler_lw_2stream, - toler_sw, - ) + @time clear_sky(context, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT, toler_lw_2stream, toler_sw) end diff --git a/test/clear_sky_utils.jl b/test/clear_sky_utils.jl index dcf7e5841..7a2a71723 100644 --- a/test/clear_sky_utils.jl +++ b/test/clear_sky_utils.jl @@ -28,8 +28,6 @@ include("read_rfmip_clear_sky.jl") #--------------------------------------------------------------- function clear_sky( context, - ::Type{OPLW}, - ::Type{OPSW}, ::Type{SLVLW}, ::Type{SLVSW}, ::Type{VMR}, @@ -37,13 +35,11 @@ function clear_sky( toler_lw, toler_sw; exfiltrate = false, -) where {FT, OPLW, OPSW, SLVLW, SLVSW, VMR} +) where {FT, SLVLW, SLVSW, VMR} overrides = (; grav = 9.80665, molmass_dryair = 0.028964, molmass_water = 0.018016) param_set = RRTMGPParameters(FT, overrides) device = ClimaComms.device(context) DA = ClimaComms.array_type(device) - op_lw = Symbol(OPLW) - op_sw = Symbol(OPSW) lw_file = get_lookup_filename(:gas, :lw) # lw lookup tables for gas optics sw_file = get_lookup_filename(:gas, :sw) # sw lookup tables for gas optics @@ -78,13 +74,13 @@ function clear_sky( # setting up longwave problem inc_flux = nothing - slv_lw = SLVLW(FT, DA, OPLW, context, param_set, nlay, ncol, sfc_emis, inc_flux) + slv_lw = SLVLW(FT, DA, context, param_set, nlay, ncol, sfc_emis, inc_flux) # setting up shortwave problem sfc_alb_diffuse = DA{FT, 2}(deepcopy(sfc_alb_direct)) inc_flux_diffuse = nothing swbcs = (cos_zenith, toa_flux, sfc_alb_direct, inc_flux_diffuse, sfc_alb_diffuse) - slv_sw = SLVSW(FT, DA, OPSW, context, nlay, ncol, swbcs...) + slv_sw = SLVSW(FT, DA, context, nlay, ncol, swbcs...) #-------------------------------------------------- # initializing RTE solver exfiltrate && Infiltrator.@exfiltrate @@ -134,10 +130,7 @@ function clear_sky( max_rel_err_flux_net_lw = maximum(rel_err_flux_net_lw) color2 = :cyan - printstyled( - "Clear-sky longwave test with ncol = $ncol, nlev = $nlev, OP = $OPLW, Solver = $SLVLW, FT = $FT\n", - color = color2, - ) + printstyled("Clear-sky longwave test with ncol = $ncol, nlev = $nlev, Solver = $SLVLW, FT = $FT\n", color = color2) printstyled("device = $device\n\n", color = color2) println("L∞ error in flux_up = $max_err_flux_up_lw") println("L∞ error in flux_dn = $max_err_flux_dn_lw") @@ -190,10 +183,7 @@ function clear_sky( max_rel_err_flux_net_sw = maximum(rel_err_flux_net_sw) - printstyled( - "Clear-sky shortwave test with ncol = $ncol, nlev = $nlev, OP = $OPSW, Solver = $SLVSW, FT = $FT\n", - color = color2, - ) + printstyled("Clear-sky shortwave test with ncol = $ncol, nlev = $nlev, Solver = $SLVSW, FT = $FT\n", color = color2) printstyled("device = $device\n\n", color = color2) println("L∞ error in flux_up = $max_err_flux_up_sw") println("L∞ error in flux_dn = $max_err_flux_dn_sw") diff --git a/test/gray_atm.jl b/test/gray_atm.jl index d05a9c4ea..e281dcba0 100644 --- a/test/gray_atm.jl +++ b/test/gray_atm.jl @@ -2,8 +2,8 @@ FT = get(ARGS, 1, Float64) == "Float32" ? Float32 : Float64 include("gray_atm_utils.jl") context = ClimaComms.context() -@time gray_atmos_lw_equil(context, OneScalar, NoScatLWRTE, FT) -@time gray_atmos_lw_equil(context, TwoStream, TwoStreamLWRTE, FT) +@time gray_atmos_lw_equil(context, NoScatLWRTE, FT) +@time gray_atmos_lw_equil(context, TwoStreamLWRTE, FT) -@time gray_atmos_sw_test(context, OneScalar, NoScatSWRTE, FT, 1) -@time gray_atmos_sw_test(context, TwoStream, TwoStreamSWRTE, FT, 1) +@time gray_atmos_sw_test(context, NoScatSWRTE, FT, 1) +@time gray_atmos_sw_test(context, TwoStreamSWRTE, FT, 1) diff --git a/test/gray_atm_utils.jl b/test/gray_atm_utils.jl index 4f3c3c723..c7a4e56dc 100644 --- a/test/gray_atm_utils.jl +++ b/test/gray_atm_utils.jl @@ -24,13 +24,7 @@ import ClimaParams as CP """ Example program to demonstrate the calculation of longwave radiative fluxes in a model gray atmosphere. """ -function gray_atmos_lw_equil( - context, - ::Type{OPLW}, - ::Type{SLVLW}, - ::Type{FT}; - exfiltrate = false, -) where {FT <: AbstractFloat, OPLW, SLVLW} +function gray_atmos_lw_equil(context, ::Type{SLVLW}, ::Type{FT}; exfiltrate = false) where {FT <: AbstractFloat, SLVLW} device = ClimaComms.device(context) param_set = RRTMGPParameters(FT) ncol = if device isa ClimaComms.CUDADevice @@ -64,7 +58,7 @@ function gray_atmos_lw_equil( otp = GrayOpticalThicknessSchneider2004(FT) # optical thickness parameters gray_as = setup_gray_as_pr_grid(context, nlay, lat, p0, pe, otp, param_set, DA) - slv_lw = SLVLW(FT, DA, OPLW, context, param_set, nlay, ncol, sfc_emis, inc_flux) + slv_lw = SLVLW(FT, DA, context, param_set, nlay, ncol, sfc_emis, inc_flux) (; flux_up, flux_dn, flux_net) = slv_lw.flux (; t_lay, p_lay, t_lev, p_lev) = gray_as @@ -108,7 +102,7 @@ function gray_atmos_lw_equil( t_error = maximum(abs.(T_ex_lev .- gray_as.t_lev)) color2 = :cyan - printstyled("\nGray atmosphere longwave test with ncol = $ncol, nlev = $nlev, OP = $OPLW\n", color = color2) + printstyled("\nGray atmosphere longwave test with ncol = $ncol, nlev = $nlev, solver = $SLVLW\n", color = color2) printstyled("device = $device\n", color = color2) printstyled("Integration time = $(FT(nsteps)/FT(24.0/tstep) / FT(365.0)) years\n\n", color = color2) @@ -124,12 +118,11 @@ end function gray_atmos_sw_test( context, - ::Type{OPSW}, ::Type{SLVSW}, ::Type{FT}, ncol::Int; exfiltrate = false, -) where {FT <: AbstractFloat, OPSW, SLVSW} +) where {FT <: AbstractFloat, SLVSW} param_set = RRTMGPParameters(FT) device = ClimaComms.device(context) DA = ClimaComms.array_type(device) @@ -172,7 +165,7 @@ function gray_atmos_sw_test( as = setup_gray_as_pr_grid(context, nlay, lat, p0, pe, otp, param_set, DA) # init gray atmos state swbcs = (cos_zenith, toa_flux, sfc_alb_direct, inc_flux_diffuse, sfc_alb_diffuse) - slv_sw = SLVSW(FT, DA, OPSW, context, nlay, ncol, swbcs...) + slv_sw = SLVSW(FT, DA, context, nlay, ncol, swbcs...) exfiltrate && Infiltrator.@exfiltrate solve_sw!(slv_sw, as) @@ -191,7 +184,7 @@ function gray_atmos_sw_test( color2 = :cyan - printstyled("\nGray atmosphere shortwave test with ncol = $ncol, nlev = $nlev, OP = $OPSW\n", color = color2) + printstyled("\nGray atmosphere shortwave test with ncol = $ncol, nlev = $nlev, solver = $SLVSW\n", color = color2) printstyled("device = $device\n\n", color = color2) println("relative error = $rel_error") diff --git a/test/runtests.jl b/test/runtests.jl index 749b4328f..8d029e6ff 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,24 +21,40 @@ printstyled("==========================\n\n", color = color1) include("clear_sky_utils.jl") for FT in (Float32, Float64) - 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) + clear_sky(context, NoScatLWRTE, TwoStreamSWRTE, VmrGM, FT, toler_lw_noscat, toler_sw) + clear_sky(context, TwoStreamLWRTE, TwoStreamSWRTE, VmrGM, FT, toler_lw_2stream, toler_sw) end end -printstyled("\n\nRRTMGP all-sky (cloudy) tests\n", color = color1) +printstyled("\n\nRRTMGP cloudy-sky (gas + clouds) tests\n", color = color1) printstyled("=================================\n\n", color = color1) -@testset "RRTMGP all-sky (cloudy) tests" begin +@testset "RRTMGP cloudy-sky (gas + clouds) tests" begin context = ClimaComms.context() + + toler_lw_noscat = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.05)) + toler_lw_2stream = Dict(Float64 => Float64(5), Float32 => Float32(5)) + toler_sw = Dict(Float64 => Float64(1e-5), Float32 => Float32(0.06)) + include("all_sky_utils.jl") for FT in (Float32, Float64) all_sky( context, - TwoStream, - TwoStream, + NoScatLWRTE, + TwoStreamSWRTE, + FT, + toler_lw_noscat, + toler_sw; + ncol = 128, + use_lut = true, + cldfrac = FT(1), + ) + all_sky( + context, TwoStreamLWRTE, TwoStreamSWRTE, - FT; + FT, + toler_lw_2stream, + toler_sw; ncol = 128, use_lut = true, cldfrac = FT(1),