diff --git a/NEWS.md b/NEWS.md index fd10c026f4..8a11fa1a90 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,9 @@ ClimaCore.jl Release Notes main ------- + - A new limiter, `VerticalMassBorrowingLimiter`, was added. Which redistributes all negative values from a given field while preserving mass. + PR [2084](https://github.com/CliMA/ClimaCore.jl/pull/2084) + ### ![][badge-✨feature/enhancement] Various improvements to `Remapper` [2060](https://github.com/CliMA/ClimaCore.jl/pull/2060) The `ClimaCore.Remapping` module received two improvements. First, `Remapper` is diff --git a/docs/refs.bib b/docs/refs.bib index fb5bab7d2c..3700198cd9 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -271,3 +271,14 @@ @article{Wiin1967 year = {1967}, publisher = {Wiley Online Library} } + +@article{zhang2018impact, + title={Impact of numerical choices on water conservation in the E3SM Atmosphere Model version 1 (EAMv1)}, + author={Zhang, Kai and Rasch, Philip J and Taylor, Mark A and Wan, Hui and Leung, Ruby and Ma, Po-Lun and Golaz, Jean-Christophe and Wolfe, Jon and Lin, Wuyin and Singh, Balwinder and others}, + journal={Geoscientific Model Development}, + volume={11}, + number={5}, + pages={1971--1988}, + year={2018}, + publisher={Copernicus Publications G{\"o}ttingen, Germany} +} diff --git a/docs/src/api.md b/docs/src/api.md index c20022a44a..4256ab9b26 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -356,6 +356,8 @@ Hypsography.ref_z_to_physical_z The limiters supertype is ```@docs Limiters.AbstractLimiter +Limiters.QuasiMonotoneLimiter +Limiters.VerticalMassBorrowingLimiter ``` This class of flux-limiters is applied only in the horizontal direction (on spectral advection operators). @@ -363,7 +365,6 @@ This class of flux-limiters is applied only in the horizontal direction (on spec ### Interfaces ```@docs -Limiters.QuasiMonotoneLimiter Limiters.compute_bounds! Limiters.apply_limiter! ``` diff --git a/src/Limiters/Limiters.jl b/src/Limiters/Limiters.jl index bd9116cf0d..a731f38fcc 100644 --- a/src/Limiters/Limiters.jl +++ b/src/Limiters/Limiters.jl @@ -20,5 +20,6 @@ abstract type AbstractLimiter end # implementations include("quasimonotone.jl") +include("vertical_mass_borrowing_limiter.jl") end # end module diff --git a/src/Limiters/vertical_mass_borrowing_limiter.jl b/src/Limiters/vertical_mass_borrowing_limiter.jl new file mode 100644 index 0000000000..c33d79d148 --- /dev/null +++ b/src/Limiters/vertical_mass_borrowing_limiter.jl @@ -0,0 +1,136 @@ +import .DataLayouts as DL + +""" + VerticalMassBorrowingLimiter(f::Fields.Field, q_min) + +A vertical-only mass borrowing limiter. + +The mass borrower borrows tracer mass from an adjacent, lower layer. +It conserves the total tracer mass and can avoid negative tracers. + +At level k, it will first borrow the mass from the layer k+1 (lower level). +If the mass is not sufficient in layer k+1, it will borrow mass from +layer k+2. The borrower will proceed this process until the bottom layer. +If the tracer mass in the bottom layer goes negative, it will repeat the +process from the bottom to the top. In this way, the borrower works for +any shape of mass profiles. + +This code was adapted from [E3SM](https://github.com/E3SM-Project/E3SM/blob/2c377c5ec9a5585170524b366ad85074ab1b1a5c/components/eam/src/physics/cam/massborrow.F90) + +References: + - [zhang2018impact](@cite) +""" +struct VerticalMassBorrowingLimiter{F, T} + bmass::F + ic::F + q_min::T +end +function VerticalMassBorrowingLimiter(f::Fields.Field, q_min) + bmass = similar(Spaces.level(f, 1)) + ic = similar(Spaces.level(f, 1)) + return VerticalMassBorrowingLimiter(bmass, ic, q_min) +end + + +""" + apply_limiter!(q::Fields.Field, ρ::Fields.Field, lim::VerticalMassBorrowingLimiter) + +Apply the vertical mass borrowing +limiter `lim` to field `q`, given +density `ρ`. +""" +apply_limiter!( + q::Fields.Field, + ρ::Fields.Field, + lim::VerticalMassBorrowingLimiter, +) = apply_limiter!(q, ρ, axes(q), lim, ClimaComms.device(axes(q))) + +function apply_limiter!( + q::Fields.Field, + ρ::Fields.Field, + space::Spaces.FiniteDifferenceSpace, + lim::VerticalMassBorrowingLimiter, + device::ClimaComms.AbstractCPUDevice, +) + cache = (; bmass = lim.bmass, ic = lim.ic, q_min = lim.q_min) + columnwise_massborrow_cpu(q, ρ, cache) +end + +function apply_limiter!( + q::Fields.Field, + ρ::Fields.Field, + space::Spaces.ExtrudedFiniteDifferenceSpace, + lim::VerticalMassBorrowingLimiter, + device::ClimaComms.AbstractCPUDevice, +) + Fields.bycolumn(axes(q)) do colidx + cache = (; + bmass = lim.bmass[colidx], + ic = lim.ic[colidx], + q_min = lim.q_min, + ) + columnwise_massborrow_cpu(q[colidx], ρ[colidx], cache) + end +end + +# TODO: can we improve the performance? +# `bycolumn` on the CPU may be better here since we could multithread it. +function columnwise_massborrow_cpu(q::Fields.Field, ρ::Fields.Field, cache) # column fields + (; bmass, ic, q_min) = cache + + Δz = Fields.Δz_field(q) + Δz_vals = Fields.field_values(Δz) + (; J) = Fields.local_geometry_field(ρ) + # ΔV_vals = Fields.field_values(J) + ΔV_vals = Δz_vals + ρ_vals = Fields.field_values(ρ) + # loop over tracers + nlevels = Spaces.nlevels(axes(q)) + @. ic = 0 + @. bmass = 0 + q_vals = Fields.field_values(q) + # top to bottom + for f in 1:DataLayouts.ncomponents(q_vals) + for v in 1:nlevels + CI = CartesianIndex(1, 1, f, v, 1) + # new mass in the current layer + ρΔV_lev = + DL.getindex_field(ΔV_vals, CI) * DL.getindex_field(ρ_vals, CI) + nmass = DL.getindex_field(q_vals, CI) + bmass[] / ρΔV_lev + if nmass > q_min[f] + # if new mass in the current layer is positive, don't borrow mass any more + DL.setindex_field!(q_vals, nmass, CI) + bmass[] = 0 + else + # set mass to q_min in the current layer, and save bmass + bmass[] = (nmass - q_min[f]) * ρΔV_lev + DL.setindex_field!(q_vals, q_min[f], CI) + ic[] = ic[] + 1 + end + end + + # bottom to top + for v in nlevels:-1:1 + CI = CartesianIndex(1, 1, f, v, 1) + # if the surface layer still needs to borrow mass + if bmass[] < 0 + ρΔV_lev = + DL.getindex_field(ΔV_vals, CI) * + DL.getindex_field(ρ_vals, CI) + # new mass in the current layer + nmass = DL.getindex_field(q_vals, CI) + bmass[] / ρΔV_lev + if nmass > q_min[f] + # if new mass in the current layer is positive, don't borrow mass any more + DL.setindex_field!(q_vals, nmass, CI) + bmass[] = 0 + else + # if new mass in the current layer is negative, continue to borrow mass + bmass[] = (nmass - q_min[f]) * ρΔV_lev + DL.setindex_field!(q_vals, q_min[f], CI) + end + end + end + end + + return nothing +end diff --git a/test/Limiters/vertical_mass_borrowing_limiter.jl b/test/Limiters/vertical_mass_borrowing_limiter.jl new file mode 100644 index 0000000000..e8c9702ce4 --- /dev/null +++ b/test/Limiters/vertical_mass_borrowing_limiter.jl @@ -0,0 +1,168 @@ +#= +julia --project=.buildkite +using Revise; include(joinpath("test", "Limiters", "vertical_mass_borrowing_limiter.jl")) +=# +using ClimaComms +ClimaComms.@import_required_backends +using ClimaCore: Fields, Spaces, Limiters +using ClimaCore.RecursiveApply +using ClimaCore.Geometry +using ClimaCore.Grids +using ClimaCore.CommonGrids +using Test +using Random + +#= +import Plots +import ClimaCorePlots +function plot_results(f, f₀) + col = Fields.ColumnIndex((1, 1), 1) + fcol = f[col] + f₀col = f₀[col] + p = Plots.plot() + Plots.plot(fcol; label = "field") + Plots.plot!(f₀col; label = "initial") +end +plot_results(ρq, ρq_init) +=# + +function perturb_field!(f::Fields.Field; perturb_radius) + device = ClimaComms.device(f) + ArrayType = ClimaComms.array_type(device) + rand_data = ArrayType(rand(size(parent(f))...)) # [0-1] + rand_data = rand_data .- sum(rand_data) / length(rand_data) # make centered about 0 ([-0.5:0.5]) + rand_data = (rand_data ./ maximum(rand_data)) .* perturb_radius # rand_data now in [-perturb_radius:perturb_radius] + parent(f) .= parent(f) .+ rand_data # f in [f ± perturb_radius] + return nothing +end + +@testset "Vertical mass borrowing limiter - column" begin + FT = Float64 + Random.seed!(1234) + z_elem = 10 + z_min = 0 + z_max = 1 + device = ClimaComms.device() + grid = ColumnGrid(FT; z_elem, z_min, z_max, device) + cspace = Spaces.FiniteDifferenceSpace(grid, Grids.CellCenter()) + tol = 0.01 + perturb_q = 0.3 + perturb_ρ = 0.2 + + # Initialize fields + coords = Fields.coordinate_field(cspace) + ρ = map(coord -> 1.0, coords) + q = map(coord -> 0.1, coords) + (; z) = coords + perturb_field!(q; perturb_radius = perturb_q) + perturb_field!(ρ; perturb_radius = perturb_ρ) + ρq_init = ρ .⊠ q + sum_ρq_init = sum(ρq_init) + + # Test that the minimum is below 0 + @test length(parent(q)) == z_elem == 10 + @test 0.3 ≤ count(x -> x < 0, parent(q)) / 10 ≤ 0.5 # ensure reasonable percentage of points are negative + + @test -2 * perturb_q ≤ minimum(q) ≤ -tol + limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,)) + Limiters.apply_limiter!(q, ρ, limiter) + @test 0 ≤ minimum(q) + ρq = ρ .⊠ q + @test isapprox(sum(ρq), sum_ρq_init; atol = 1e-15) + @test isapprox(sum(ρq), sum_ρq_init; rtol = 1e-10) + # @show sum(ρq) # 0.07388931313511024 + # @show sum_ρq_init # 0.07388931313511025 +end + +@testset "Vertical mass borrowing limiter - sphere" begin + FT = Float64 + z_elem = 10 + z_min = 0 + z_max = 1 + radius = 10 + h_elem = 10 + n_quad_points = 4 + grid = ExtrudedCubedSphereGrid(FT; + z_elem, + z_min, + z_max, + radius, + h_elem, + n_quad_points, + ) + cspace = Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter()) + tol = 0.01 + perturb_q = 0.2 + perturb_ρ = 0.1 + + # Initialize fields + coords = Fields.coordinate_field(cspace) + ρ = map(coord -> 1.0, coords) + q = map(coord -> 0.1, coords) + + perturb_field!(q; perturb_radius = perturb_q) + perturb_field!(ρ; perturb_radius = perturb_ρ) + ρq_init = ρ .⊠ q + sum_ρq_init = sum(ρq_init) + + # Test that the minimum is below 0 + @test length(parent(q)) == z_elem * h_elem^2 * 6 * n_quad_points^2 == 96000 + @test 0.10 ≤ count(x -> x < 0.0001, parent(q)) / 96000 ≤ 1 # ensure 10% of points are negative + + @test -2 * perturb_q ≤ minimum(q) ≤ -tol + limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,)) + Limiters.apply_limiter!(q, ρ, limiter) + @test 0 ≤ minimum(q) + ρq = ρ .⊠ q + @test isapprox(sum(ρq), sum_ρq_init; atol = 0.07) + @test isapprox(sum(ρq), sum_ρq_init; rtol = 0.07) + # @show sum(ρq) # 125.5483524237572 + # @show sum_ρq_init # 125.52052986152977 +end + +@testset "Vertical mass borrowing limiter - deep atmosphere" begin + FT = Float64 + z_elem = 10 + z_min = 0 + z_max = 1 + radius = 10 + h_elem = 10 + n_quad_points = 4 + grid = ExtrudedCubedSphereGrid(FT; + z_elem, + z_min, + z_max, + radius, + global_geometry = Geometry.DeepSphericalGlobalGeometry{FT}(radius), + h_elem, + n_quad_points, + ) + cspace = Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter()) + tol = 0.01 + perturb_q = 0.2 + perturb_ρ = 0.1 + + # Initialize fields + coords = Fields.coordinate_field(cspace) + ρ = map(coord -> 1.0, coords) + q = map(coord -> 0.1, coords) + + perturb_field!(q; perturb_radius = perturb_q) + perturb_field!(ρ; perturb_radius = perturb_ρ) + ρq_init = ρ .⊠ q + sum_ρq_init = sum(ρq_init) + + # Test that the minimum is below 0 + @test length(parent(q)) == z_elem * h_elem^2 * 6 * n_quad_points^2 == 96000 + @test 0.10 ≤ count(x -> x < 0.0001, parent(q)) / 96000 ≤ 1 # ensure 10% of points are negative + + @test -2 * perturb_q ≤ minimum(q) ≤ -tol + limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,)) + Limiters.apply_limiter!(q, ρ, limiter) + @test 0 ≤ minimum(q) + ρq = ρ .⊠ q + @test isapprox(sum(ρq), sum_ρq_init; atol = 0.269) + @test isapprox(sum(ρq), sum_ρq_init; rtol = 0.00199) + # @show sum(ρq) # 138.90494931641584 + # @show sum_ρq_init # 139.1735731377394 +end diff --git a/test/Limiters/vertical_mass_borrowing_limiter_advection.jl b/test/Limiters/vertical_mass_borrowing_limiter_advection.jl new file mode 100644 index 0000000000..884d1f2076 --- /dev/null +++ b/test/Limiters/vertical_mass_borrowing_limiter_advection.jl @@ -0,0 +1,207 @@ +#= +julia --project=.buildkite +using Revise; include(joinpath("test", "Limiters", "vertical_mass_borrowing_limiter_advection.jl")) +=# +using Test +using LinearAlgebra +import ClimaComms +ClimaComms.@import_required_backends +using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 +using ClimaCore.CommonGrids +using ClimaCore.Grids +using ClimaTimeSteppers + +import ClimaCore: + Fields, + Domains, + Limiters, + Topologies, + Meshes, + DataLayouts, + Operators, + Geometry, + Spaces + + +# Advection Equation, with constant advective velocity (so advection form = flux form) +# ∂_t y + w ∂_z y = 0 +# the solution translates to the right at speed w, +# so at time t, the solution is y(z - w * t) + +# visualization artifacts +ENV["GKSwstype"] = "nul" +import ClimaCorePlots +import Plots +Plots.GRBackend() +dir = "vert_mass_borrow_advection" +path = joinpath(@__DIR__, "output", dir) +mkpath(path) + +function lim!(y, parameters, t, y_ref) + (; w, Δt, limiter) = parameters + Limiters.apply_limiter!(y.q, y.ρ, limiter) + return nothing +end + +function perturb_field!(f::Fields.Field; perturb_radius) + device = ClimaComms.device(f) + ArrayType = ClimaComms.array_type(device) + rand_data = ArrayType(rand(size(parent(f))...)) # [0-1] + rand_data = rand_data .- sum(rand_data) / length(rand_data) # make centered about 0 ([-0.5:0.5]) + rand_data = (rand_data ./ maximum(rand_data)) .* perturb_radius # rand_data now in [-perturb_radius:perturb_radius] + parent(f) .= parent(f) .+ rand_data # f in [f ± perturb_radius] + return nothing +end + +function tendency!(yₜ, y, parameters, t) + (; w, Δt, limiter) = parameters + FT = Spaces.undertype(axes(y.q)) + bcvel = pulse(-π, t, z₀, zₕ, z₁) + divf2c = Operators.DivergenceF2C( + bottom = Operators.SetValue(Geometry.WVector(FT(bcvel))), + top = Operators.SetValue(Geometry.WVector(FT(0))), + ) + upwind1 = Operators.UpwindBiasedProductC2F( + bottom = Operators.Extrapolate(), + top = Operators.Extrapolate(), + ) + upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( + bottom = Operators.ThirdOrderOneSided(), + top = Operators.ThirdOrderOneSided(), + ) + If = Operators.InterpolateC2F() + @. yₜ.q = -divf2c(upwind1(w, y.q) * If(y.q)) + return nothing +end + +# Define a pulse wave or square wave + +FT = Float64 +t₀ = FT(0.0) +Δt = 0.0001 * 25 +t₁ = 200Δt +z₀ = FT(0.0) +zₕ = FT(1.0) +z₁ = FT(1.0) +speed = FT(-1.0) +pulse(z, t, z₀, zₕ, z₁) = abs(z - speed * t) ≤ zₕ ? z₁ : z₀ + +stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(FT(7.0))) +plot_string = ["uniform", "stretched"] + +@testset "VerticalMassBorrowingLimiter on advection" begin + for (i, stretch) in enumerate(stretch_fns) + i = 1 + stretch = Meshes.Uniform() + + z_elem = 2^6 + z_min = -π + z_max = π + device = ClimaComms.device() + + # use_column = true + use_column = false + if use_column + grid = ColumnGrid(; z_elem, z_min, z_max, stretch, device) + cspace = Spaces.FiniteDifferenceSpace(grid, Grids.CellCenter()) + fspace = Spaces.FaceFiniteDifferenceSpace(cspace) + else + grid = ExtrudedCubedSphereGrid(; + z_elem, + z_min, + z_max, + stretch, + radius = 10, + h_elem = 10, + n_quad_points = 4, + device, + ) + cspace = + Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter()) + fspace = Spaces.FaceExtrudedFiniteDifferenceSpace(cspace) + end + + z = Fields.coordinate_field(cspace).z + O = ones(FT, fspace) + + # Initial condition + q_init = pulse.(z, 0.0, z₀, zₕ, z₁) + q = q_init + coords = Fields.coordinate_field(q) + ρ = map(coord -> 1.0, coords) + perturb_field!(ρ; perturb_radius = 0.1) + y = Fields.FieldVector(; q, ρ) + limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,)) + + # Unitary, constant advective velocity + w = Geometry.WVector.(speed .* O) + + # Solve the ODE + parameters = (; w, Δt, limiter) + prob = ODEProblem( + ClimaODEFunction(; lim!, T_lim! = tendency!), + y, + (t₀, t₁), + parameters, + ) + sol = solve( + prob, + ExplicitAlgorithm(SSP33ShuOsher()), + dt = Δt, + saveat = Δt, + ) + + q_init = sol.u[1].q + q_final = sol.u[end].q + q_analytic = pulse.(z, t₁, z₀, zₕ, z₁) + err = norm(q_final .- q_analytic) + rel_mass_err = norm((sum(q_final) - sum(q_init)) / sum(q_init)) + + + if use_column + p = Plots.plot() + Plots.plot!(q_init, label = "init") + Plots.plot!(q_final, label = "computed") + Plots.plot!(q_analytic, label = "analytic") + Plots.plot!(; legend = :topright) + Plots.plot!(; xlabel = "q", title = "VerticalMassBorrowingLimiter") + f = joinpath( + path, + "VerticalMassBorrowingLimiter_comparison_$(plot_string[i]).png", + ) + Plots.png(p, f) + else + colidx = Fields.ColumnIndex((1, 1), 1) + p = Plots.plot() + Plots.plot!( + vec(parent(z[colidx])), + vec(parent(q_init[colidx])), + label = "init", + ) + Plots.plot!( + vec(parent(z[colidx])), + vec(parent(q_final[colidx])), + label = "computed", + ) + Plots.plot!( + vec(parent(z[colidx])), + vec(parent(q_analytic[colidx])), + label = "analytic", + ) + Plots.plot!(; legend = :topright) + Plots.plot!(; + xlabel = "q", + ylabel = "z", + title = "VerticalMassBorrowingLimiter", + ) + f = joinpath( + path, + "VerticalMassBorrowingLimiter_comparison_$(plot_string[i]).png", + ) + Plots.png(p, f) + end + @test err ≤ 0.25 + @test rel_mass_err ≤ 10eps() + @test minimum(q_final) ≥ 0 + end +end