From ec5300713a3e0aac2804a32038c80e15e6925033 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Wed, 13 Nov 2024 22:28:37 -0500 Subject: [PATCH] Implement and test vertical mass borrowing limiters --- NEWS.md | 3 + docs/refs.bib | 11 ++ docs/src/api.md | 3 +- src/Limiters/Limiters.jl | 1 + .../vertical_mass_borrowing_limiter.jl | 119 ++++++++++++++ .../vertical_mass_borrowing_limiter.jl | 91 +++++++++++ ...rtical_mass_borrowing_limiter_advection.jl | 145 ++++++++++++++++++ 7 files changed, 372 insertions(+), 1 deletion(-) create mode 100644 src/Limiters/vertical_mass_borrowing_limiter.jl create mode 100644 test/Limiters/vertical_mass_borrowing_limiter.jl create mode 100644 test/Limiters/vertical_mass_borrowing_limiter_advection.jl diff --git a/NEWS.md b/NEWS.md index b4eb35a610..782153802e 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) + v0.14.20 ------- diff --git a/docs/refs.bib b/docs/refs.bib index 8d3111f485..a2b55579f4 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -257,3 +257,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 cfb9fd029b..02e80acbcf 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -332,6 +332,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). @@ -339,7 +341,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 55b8deb36c..2a3fd07f3e 100644 --- a/src/Limiters/Limiters.jl +++ b/src/Limiters/Limiters.jl @@ -19,5 +19,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..1a60d6f9ef --- /dev/null +++ b/src/Limiters/vertical_mass_borrowing_limiter.jl @@ -0,0 +1,119 @@ +import .DataLayouts as DL + +""" + VerticalMassBorrowingLimiter(f::Fields.Field, q_min) + +A vertical-only mass borrowing limier. + +The mass borrower borrows tracer mass from an adjacent layer. +It conserves the 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, lim::VerticalMassBorrowingLimiter) + +Apply the vertical mass borrowing +limiter `lim` to field `q`. +""" +apply_limiter!(q::Fields.Field, lim::VerticalMassBorrowingLimiter) = + apply_limiter!(q, axes(q), lim, ClimaComms.device(axes(q))) + +function apply_limiter!( + q::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, + 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], 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, cache) # column fields + (; bmass, ic, q_min) = cache + + pdel = Fields.field_values(Fields.Δz_field(q)) + # 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 + pdel_v = DL.getindex_field(pdel, CI) + nmass = DL.getindex_field(q_vals, CI) + bmass[] / pdel_v + 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]) * pdel_v + 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 + pdel_v = DL.getindex_field(pdel, CI) + # new mass in the current layer + nmass = DL.getindex_field(q_vals, CI) + bmass[] / pdel_v + 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]) * pdel_v + 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..fa7904ac8e --- /dev/null +++ b/test/Limiters/vertical_mass_borrowing_limiter.jl @@ -0,0 +1,91 @@ +#= +julia --project +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.Grids +using ClimaCore.CommonGrids +using Test +using Random + +@testset "Vertical mass borrowing limiter - column" begin + Random.seed!(1234) + FT = Float64 + z_elem = 10 + z_min = 0 + z_max = 1 + device = ClimaComms.device() + grid = ColumnGrid(; z_elem, z_min, z_max, device) + cspace = Spaces.FiniteDifferenceSpace(grid, Grids.CellCenter()) + context = ClimaComms.context(device) + ArrayType = ClimaComms.array_type(device) + tol = 0.01 + perturb = 0.2 + + # Initialize fields + coords = Fields.coordinate_field(cspace) + q = map(coord -> 0.1, coords) + (; z) = coords + rand_data = ArrayType(rand(size(parent(q))...)) # [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 # rand_data now in [-0.2:0.2] + parent(q) .= parent(q) .+ rand_data # q in [0.1 ± 0.2] + sum_q_init = sum(q) + # 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 ≤ minimum(q) ≤ -tol + limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,)) + Limiters.apply_limiter!(q, limiter) + @test 0 ≤ minimum(q) + @test isapprox(sum(q), sum_q_init; atol = 0.00000000001) + @test isapprox(sum(q), sum_q_init; rtol = 0.01) +end + +@testset "Vertical mass borrowing limiter" begin + FT = Float64 + z_elem = 10 + z_min = 0 + z_max = 1 + radius = 10 + h_elem = 10 + n_quad_points = 4 + grid = ExtrudedCubedSphereGrid(; + z_elem, + z_min, + z_max, + radius, + h_elem, + n_quad_points, + ) + cspace = Spaces.ExtrudedFiniteDifferenceSpace(grid, Grids.CellCenter()) + context = ClimaComms.context(device) + ArrayType = ClimaComms.array_type(device) + tol = 0.01 + perturb = 0.2 + + # Initialize fields + coords = Fields.coordinate_field(cspace) + q = map(coord -> 0.1, coords) + + rand_data = ArrayType(rand(size(parent(q))...)) # [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 # rand_data now in [-0.2:0.2] + parent(q) .= parent(q) .+ rand_data # q in [0.1 ± 0.2] + sum_q_init = sum(q) + + # 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 -0.11 ≤ minimum(q) ≤ -0.08 + limiter = Limiters.VerticalMassBorrowingLimiter(q, (0.0,)) + Limiters.apply_limiter!(q, limiter) + @test 0 ≤ minimum(q) + @test isapprox(sum(q), sum_q_init; atol = 0.013) + @test isapprox(sum(q), sum_q_init; rtol = 0.01) +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..5068612d07 --- /dev/null +++ b/test/Limiters/vertical_mass_borrowing_limiter_advection.jl @@ -0,0 +1,145 @@ +#= +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 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" +using ClimaCorePlots, 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, limiter) + 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(w * If(y.q)) + -divf2c(upwind1(w, y.q) * If(y.q)) + # -divf2c(upwind3(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₀ + +n = 2 .^ 6 + +domain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(-π), + Geometry.ZPoint{FT}(π); + boundary_names = (:bottom, :top), +) + +# stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(FT(7.0))) +stretch_fns = (Meshes.Uniform(),) +plot_string = ["uniform", "stretched"] + +@testset "VerticalMassBorrowingLimiter on advection" begin + for (i, stretch_fn) in enumerate(stretch_fns) + mesh = Meshes.IntervalMesh(domain, stretch_fn; nelems = n) + cent_space = Spaces.CenterFiniteDifferenceSpace(mesh) + face_space = Spaces.FaceFiniteDifferenceSpace(cent_space) + z = Fields.coordinate_field(cent_space).z + O = ones(FT, face_space) + + # Initial condition + q_init = pulse.(z, 0.0, z₀, zₕ, z₁) + q = q_init + 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)) + + + p = 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) + @test err ≤ 0.25 + @test rel_mass_err ≤ 10eps() + @test minimum(q_final) ≥ 0 + end +end