diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7fdc1b51f8..506f3141a4 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -773,17 +773,6 @@ steps: agents: slurm_gpus: 1 - - label: "Unit: FD operator (shmem)" - key: unit_fd_operator_shmem - retry: *retry_policy - command: - - "julia --color=yes --check-bounds=yes --project=.buildkite test/Operators/finitedifference/unit_fd_ops_shared_memory.jl" - - "julia --color=yes --project=.buildkite test/Operators/finitedifference/benchmark_fd_ops_shared_memory.jl" - env: - CLIMACOMMS_DEVICE: "CUDA" - agents: - slurm_gpus: 1 - - label: "Unit: gpu columnwise" key: unit_gpu_columnwise retry: *retry_policy @@ -1848,14 +1837,6 @@ steps: artifact_paths: - "examples/column/output/heat/*" - - label: ":computer: Column Advection Step Eq" - key: "cpu_column_step_advect" - retry: *retry_policy - command: - - "julia --color=yes --project=.buildkite examples/column/step.jl" - artifact_paths: - - "examples/column/output/advect_step_function/*" - - label: ":computer: Column Advection Eq" key: "cpu_column_advect" retry: *retry_policy diff --git a/NEWS.md b/NEWS.md index e3c247e251..7d28f80559 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,11 @@ ClimaCore.jl Release Notes ======================== +- Removed unused finite difference operators and boundary condition [2521](https://github.com/CliMA/ClimaCore.jl/pull/2521) + - Removed `SetValue` from `GradientC2F`, `DivergenceC2F`, `CurlC2F`, `UpwindBiasedProductC2F` + - Removed `SetGradient` from `InterpolateC2F` and `WeightedInterpolateC2F` + - Removed `AdvectionC2C`, `AdvectionF2F`, `FluxCorrectionF2F`, and `FluxCorrectionC2C` operators + main ------- diff --git a/docs/Project.toml b/docs/Project.toml index 38e051b3df..1e1a333d07 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -11,6 +11,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" +LazyBroadcast = "9dccce8e-a116-406d-9fcc-a88ed4f510c8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" diff --git a/docs/src/examples.md b/docs/src/examples.md index a00b801f1a..ebdba0bdd5 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -84,9 +84,9 @@ The example code solves the equation for 4 different tendencies with the followi where $\textrm{fcc}(v, \theta)$ is the [`center-to-center flux correction`](https://github.com/CliMA/ClimaCore.jl/blob/main/src/Operators/finitedifference.jl#L2617) operator. - Tendency 3: - $$D = A,$$ + $$D = v \cdot G(T)$$ - where $A$ is the [`center-to-center vertical advection`](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.AdvectionC2C) operator. + where``G`` is the [center-to-face gradient](https://clima.github.io/ClimaCore.jl/dev/operators/#ClimaCore.Operators.GradientC2F) operator, called `gradc2f` in the example code - Tendency 4: $$D = A + \textrm{fcc}(v, \theta),$$ diff --git a/docs/src/operators.md b/docs/src/operators.md index b3b6a25d8b..c4066e7176 100644 --- a/docs/src/operators.md +++ b/docs/src/operators.md @@ -74,8 +74,6 @@ AbstractTVDSlopeLimiter ```@docs GradientF2C GradientC2F -AdvectionF2F -AdvectionC2C DivergenceF2C DivergenceC2F CurlC2F diff --git a/docs/tutorials/introduction.jl b/docs/tutorials/introduction.jl index 70dfdaa6b3..e684782c26 100644 --- a/docs/tutorials/introduction.jl +++ b/docs/tutorials/introduction.jl @@ -11,7 +11,7 @@ # - currently under development using ClimaComms, ClimaCore, ClimaCorePlots, LinearAlgebra, IntervalSets, Plots - +import LazyBroadcast: lazy #---------------------------------------------------------------------------- # ## 1. Constructing a discretization @@ -297,13 +297,8 @@ gradc2f = ClimaCore.Operators.GradientC2F() # ???? # ``` # -# To handle boundaries we need to *modify the stencil*. Two options: -# - provide the _value_ $\theta^*$ of $\theta$ at the boundary: -# ```math -# \nabla\theta[\tfrac{1}{2}] = \frac{\theta[1] - \theta^*}{\Delta z /2} -# ``` +# To handle boundaries we need to *modify the stencil*. This is done by providing the *gradient* $\nabla\theta^*$ of $\theta$ at the boundary: # -# - provide the *gradient* $\nabla\theta^*$ of $\theta$ at the boundary: # ```math # \nabla\theta[\tfrac{1}{2}] = \nabla\theta^* # ``` @@ -312,7 +307,9 @@ gradc2f = ClimaCore.Operators.GradientC2F() sinz = sin.(column_center_coords.z) gradc2f = ClimaCore.Operators.GradientC2F( - bottom = ClimaCore.Operators.SetValue(sin(0.0)), + bottom = ClimaCore.Operators.SetGradient( + ClimaCore.Geometry.WVector(cos(0.0)), + ), top = ClimaCore.Operators.SetGradient( ClimaCore.Geometry.WVector(cos(10.0)), ), @@ -370,7 +367,8 @@ import OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 # \frac{\partial y}{\partial t} = \alpha \nabla \cdot \nabla y # ``` # -# At the bottom we will use a Dirichlet condition ``y(0) = 1`` at the bottom: since we don't actually have a value located at the bottom, we will use a `SetValue` boundary modifier on the inner gradient. +# At the bottom we will use a Dirichlet condition ``y(0) = 1`` at the bottom: since we don't actually have a value located at the bottom, we will use a `SetGradient` boundary modifier on the inner gradient. +# If ``y(0) = 1``, then the gradient at the first face is ``\frac{\partial y}{\partial z}[\tfrac{1}{2}] = 2(y[\tfrac{1}{2}] - 1)``. # # At the top we will use a Neumann condition ``\frac{\partial y}{\partial z}(10) = 0``. We can do this two equivalent ways: # - a `SetGradient` on the gradient operator @@ -382,9 +380,11 @@ y0 = zeros(column_center_space) ## define the tendency function function heat_fd_tendency!(dydt, y, α, t) + bottom_level_y = ClimaCore.Fields.level(y, 1) + bottom_grad = @. lazy(ClimaCore.Geometry.Covariant3Vector(2.0 * (bottom_level_y - 1.0))) gradc2f = ClimaCore.Operators.GradientC2F( - bottom = ClimaCore.Operators.SetValue(1.0), - top = ClimaCore.Operators.SetGradient(ClimaCore.Geometry.WVector(0.0)), + bottom = ClimaCore.Operators.SetGradient(bottom_grad), + top = ClimaCore.Operators.SetGradient(ClimaCore.Geometry.Covariant3Vector(0.0)), ) divf2c = ClimaCore.Operators.DivergenceF2C() ## the @. macro "dots" the whole expression diff --git a/examples/column/advect.jl b/examples/column/advect.jl index 5d8aee9a46..aeee2e122f 100644 --- a/examples/column/advect.jl +++ b/examples/column/advect.jl @@ -8,9 +8,13 @@ import ClimaCore: DataLayouts, Operators, Geometry, - Spaces + Spaces, + Utilities +import LazyBroadcast: lazy using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 +using ClimaCorePlots +using Plots import Logging import TerminalLoggers @@ -27,187 +31,268 @@ domain = Domains.IntervalDomain( Geometry.ZPoint{FT}(b), boundary_names = (:left, :right), ) -mesh = Meshes.IntervalMesh(domain, nelems = n) +mesh_sin = Meshes.IntervalMesh(domain, nelems = n) + +a = FT(-20.0) +b = FT(20.0) +n = 64 +α = FT(0.1) +domain = Domains.IntervalDomain( + Geometry.ZPoint(a), + Geometry.ZPoint(b), + boundary_names = (:left, :right), +) +mesh_step = Meshes.IntervalMesh(domain, nelems = n) device = ClimaComms.device() -cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) -fs = Spaces.FaceFiniteDifferenceSpace(cs) +for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step)) -V = Geometry.WVector.(ones(FT, fs)) -θ = sin.(Fields.coordinate_field(cs).z) + cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) + fs = Spaces.FaceFiniteDifferenceSpace(cs) -# Solve advection Equation: ∂θ/dt = -∂(vθ) + V = Geometry.WVector.(ones(FT, fs)) + if fn == "sin" + θ = sin.(Fields.coordinate_field(cs).z) + elseif fn == "step" + function heaviside(pt) + 0.5 * (sign(pt.z) + 1) + end + θ = heaviside.(Fields.coordinate_field(cs)) + end -# upwinding -function tendency1!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a - t)), - right = Operators.SetValue(sin(b - t)), - ) - ∂ = Operators.DivergenceF2C() + # Solve advection Equation: ∂θ/dt = -∂(vθ) - return @. dθ = -∂(UB(V, θ)) -end -function tendency2!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a - t)), - right = Operators.SetValue(sin(b - t)), - ) - ∂ = Operators.DivergenceF2C() - return @. dθ = -∂(UB(V, θ)) + fcc(V, θ) -end -# use the advection operator -function tendency3!(dθ, θ, _, t) + # upwinding + function tendency1!(dθ, θ, _, t) + lg_field = Fields.local_geometry_field(fs) + lg_left = Fields.level(lg_field, Utilities.PlusHalf(0)) + lg_right = Fields.level(lg_field, Fields.nlevels(lg_field) - Utilities.PlusHalf(0)) + v_left = Fields.field_values( + Geometry.contravariant3.(Fields.level(V, Utilities.PlusHalf(0)), lg_left), + )[] + aᴸᴮ = sin(a - t) + aᴸ = Fields.field_values(Fields.level(θ, 1))[] + left_bc = Operators.SetValue( + Geometry.Contravariant3Vector(Operators.upwind_biased_product(v_left, aᴸᴮ, aᴸ)), + ) + v_right = Fields.field_values( + Geometry.contravariant3.( + Fields.level(V, Fields.nlevels(V) - Utilities.PlusHalf(0)), + lg_right, + ), + )[] + aᴿᴮ = sin(b - t) + aᴿ = Fields.field_values(Fields.level(θ, Fields.nlevels(θ)))[] + right_bc = Operators.SetValue( + Geometry.Contravariant3Vector( + Operators.upwind_biased_product(v_right, aᴿ, aᴿᴮ), + ), + ) + set_bcs = Operators.SetBoundaryOperator(; left = left_bc, right = right_bc) + UB = Operators.UpwindBiasedProductC2F() + ∂ = Operators.DivergenceF2C() - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), + return @. dθ = -∂(set_bcs(UB(V, θ))) + end + function tendency2!(dθ, θ, _, t) + lg_field = Fields.local_geometry_field(fs) + lg_left = Fields.level(lg_field, Utilities.PlusHalf(0)) + lg_right = Fields.level(lg_field, Fields.nlevels(lg_field) - Utilities.PlusHalf(0)) + v_left = Fields.field_values( + Geometry.contravariant3.(Fields.level(V, Utilities.PlusHalf(0)), lg_left), + )[] + aᴸᴮ = sin(a - t) + aᴸ = Fields.field_values(Fields.level(θ, 1))[] + left_bc = Operators.SetValue( + Geometry.Contravariant3Vector(Operators.upwind_biased_product(v_left, aᴸᴮ, aᴸ)), + ) + v_right = Fields.field_values( + Geometry.contravariant3.( + Fields.level(V, Fields.nlevels(V) - Utilities.PlusHalf(0)), + lg_right, + ), + )[] + aᴿᴮ = sin(b - t) + aᴿ = Fields.field_values(Fields.level(θ, Fields.nlevels(θ)))[] + right_bc = Operators.SetValue( + Geometry.Contravariant3Vector( + Operators.upwind_biased_product(v_right, aᴿ, aᴿᴮ), + ), + ) + set_bcs = Operators.SetBoundaryOperator(; left = left_bc, right = right_bc) + UB = Operators.UpwindBiasedProductC2F() + ∂ = Operators.DivergenceF2C() + left_center = Fields.level(θ, 1) + θ_top = Fields.level(θ, Fields.nlevels(θ)) + θ_top_m1 = Fields.level(θ, Fields.nlevels(θ) - 1) + right_center_left_biased_grad = + @. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1)) + right_gradient_extrapolate = Operators.SetGradient(right_center_left_biased_grad) + θ_2 = Fields.level(θ, 2) + _left_lazy_grad2 = @. lazy(Geometry.Covariant3Vector(θ_2 - left_center)) + left_gradient_extrapolate = Operators.SetGradient(_left_lazy_grad2) + gradc2f_fcc = Operators.GradientC2F( + left = left_gradient_extrapolate, + right = right_gradient_extrapolate, + ) + gradf2c = Operators.GradientF2C() + return @. dθ = + -∂(set_bcs(UB(V, θ))) + + parent( + gradf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f_fcc(θ))), + ).data.:1 + end + # use the advection operator + function tendency3!(dθ, θ, _, t) + left_center = Fields.level(θ, 1) + θ_top = Fields.level(θ, Fields.nlevels(θ)) + θ_top_m1 = Fields.level(θ, Fields.nlevels(θ) - 1) + right_center_left_biased_grad = + @. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1)) + left_gradient = + Operators.SetGradient( + @. lazy(Geometry.Covariant3Vector(2 * (left_center - sin(-t)))) + ) + right_gradient = Operators.SetGradient(right_center_left_biased_grad) + gradc2f = Operators.GradientC2F(left = left_gradient, right = right_gradient) + interpf2c = Operators.InterpolateF2C() + return @. dθ = + -1 * interpf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f(θ))) + end + # use the advection operator + function tendency4!(dθ, θ, _, t) + left_center = Fields.level(θ, 1) + θ_top = Fields.level(θ, Fields.nlevels(θ)) + θ_top_m1 = Fields.level(θ, Fields.nlevels(θ) - 1) + right_center_left_biased_grad = + @. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1)) + left_gradient = + Operators.SetGradient( + @. lazy(Geometry.Covariant3Vector(2 * (left_center - sin(-t)))) + ) + right_gradient_extrapolate = Operators.SetGradient(right_center_left_biased_grad) + gradc2f = + Operators.GradientC2F(left = left_gradient, right = right_gradient_extrapolate) + interpf2c = Operators.InterpolateF2C() + θ_2 = Fields.level(θ, 2) + _left_lazy_grad4 = @. lazy(Geometry.Covariant3Vector(θ_2 - left_center)) + left_gradient_extrapolate = Operators.SetGradient(_left_lazy_grad4) + gradc2f_fcc = Operators.GradientC2F( + left = left_gradient_extrapolate, + right = right_gradient_extrapolate, + ) + gradf2c = Operators.GradientF2C() + return @. dθ = + -1 * interpf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f(θ))) + + parent( + gradf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f_fcc(θ))), + ).data.:1 + end + + # use the advection operator + + @show tendency1!(similar(θ), θ, nothing, 0.0) + # Solve the ODE operator + Δt = 0.001 + t_end = fn == "sin" ? 10.0 : 5.0 + prob1 = ODEProblem(tendency1!, θ, (0.0, t_end)) + prob2 = ODEProblem(tendency2!, θ, (0.0, t_end)) + prob3 = ODEProblem(tendency3!, θ, (0.0, t_end)) + prob4 = ODEProblem(tendency4!, θ, (0.0, t_end)) + sol1 = solve( + prob1, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, ) - A = Operators.AdvectionC2C( - left = Operators.SetValue(sin(-t)), - right = Operators.Extrapolate(), + sol2 = solve( + prob2, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, ) - return @. dθ = -A(V, θ) -end -# use the advection operator -function tendency4!(dθ, θ, _, t) - - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), + sol3 = solve( + prob3, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), + sol4 = solve( + prob4, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, ) - A = Operators.AdvectionC2C( - left = Operators.SetValue(sin(-t)), - right = Operators.Extrapolate(), + + ENV["GKSwstype"] = "nul" + + Plots.GRBackend() + + sim_type = fn == "sin" ? "advect" : "advect_step_function" + dir = sim_type + path = joinpath(@__DIR__, "output", dir) + mkpath(path) + + anim = Plots.@animate for u in sol1.u + Plots.plot(u, xlim = (-1, 1)) + end + Plots.mp4(anim, joinpath(path, "UBP_$(sim_type).mp4"), fps = 10) + Plots.png( + Plots.plot(sol1.u[end], xlim = (-1, 1)), + joinpath(path, "sol1_$(sim_type)_end.png"), ) - return @. dθ = -A(V, θ) + fcc(V, θ) -end -# use the advection operator - -@show tendency1!(similar(θ), θ, nothing, 0.0) -# Solve the ODE operator -Δt = 0.001 -prob1 = ODEProblem(tendency1!, θ, (0.0, 10.0)) -prob2 = ODEProblem(tendency2!, θ, (0.0, 10.0)) -prob3 = ODEProblem(tendency3!, θ, (0.0, 10.0)) -prob4 = ODEProblem(tendency4!, θ, (0.0, 10.0)) -sol1 = solve( - prob1, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):10.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); -sol2 = solve( - prob2, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):10.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); -sol3 = solve( - prob3, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):10.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); -sol4 = solve( - prob4, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):10.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); - -ENV["GKSwstype"] = "nul" -using ClimaCorePlots, Plots -Plots.GRBackend() - -dir = "advect" -path = joinpath(@__DIR__, "output", dir) -mkpath(path) - -anim = Plots.@animate for u in sol1.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "UBP_advect.mp4"), fps = 10) -Plots.png( - Plots.plot(sol1.u[end], xlim = (-1, 1)), - joinpath(path, "sol1_advect_end.png"), -) + anim = Plots.@animate for u in sol2.u + Plots.plot(u, xlim = (-1, 1)) + end + Plots.mp4(anim, joinpath(path, "UBP_$(sim_type)_fc.mp4"), fps = 10) + Plots.png( + Plots.plot(sol2.u[end], xlim = (-1, 1)), + joinpath(path, "sol2_$(sim_type)_end.png"), + ) -anim = Plots.@animate for u in sol2.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "UBP_advect_fc.mp4"), fps = 10) -Plots.png( - Plots.plot(sol2.u[end], xlim = (-1, 1)), - joinpath(path, "sol2_advect_end.png"), -) + anim = Plots.@animate for u in sol3.u + Plots.plot(u, xlim = (-1, 1)) + end + Plots.mp4(anim, joinpath(path, "C2C_$(sim_type).mp4"), fps = 10) + Plots.png( + Plots.plot(sol3.u[end], xlim = (-1, 1)), + joinpath(path, "sol3_$(sim_type)_end.png"), + ) -anim = Plots.@animate for u in sol3.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "C2C_advect.mp4"), fps = 10) -Plots.png( - Plots.plot(sol3.u[end], xlim = (-1, 1)), - joinpath(path, "sol3_advect_end.png"), -) + anim = Plots.@animate for u in sol4.u + Plots.plot(u, xlim = (-1, 1)) + end + Plots.mp4(anim, joinpath(path, "C2C_$(sim_type)_fc.mp4"), fps = 10) + Plots.png( + Plots.plot(sol4.u[end], xlim = (-1, 1)), + joinpath(path, "sol4_$(sim_type)_end.png"), + ) -anim = Plots.@animate for u in sol4.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "C2C_advect_fc.mp4"), fps = 10) -Plots.png( - Plots.plot(sol4.u[end], xlim = (-1, 1)), - joinpath(path, "sol4_advect_end.png"), -) + p = Plots.plot(sol1.u[end], xlim = (-1, 1), ls = :dash, label = "UBP") + p = Plots.plot!(sol2.u[end], xlim = (-1, 1), ls = :dot, label = "UBP_FC") + p = Plots.plot!(sol3.u[end], xlim = (-1, 1), ls = :solid, label = "C2C") + p = Plots.plot!(sol4.u[end], xlim = (-1, 1), ls = :dashdot, label = "C2C_FC") + Plots.png(p, joinpath(path, "all_$(sim_type)_end.png")) -p = Plots.plot(sol1.u[end], xlim = (-1, 1), ls = :dash, label = "UBP") -p = Plots.plot!(sol2.u[end], xlim = (-1, 1), ls = :dot, label = "UBP_FC") -p = Plots.plot!(sol3.u[end], xlim = (-1, 1), ls = :solid, label = "C2C") -p = Plots.plot!(sol4.u[end], xlim = (-1, 1), ls = :dashdot, label = "C2C_FC") -Plots.png(p, joinpath(path, "all_advect_end.png")) - -function linkfig(figpath, alt = "") - # buildkite-agent upload figpath - # link figure in logs if we are running on CI - if get(ENV, "BUILDKITE", "") == "true" - artifact_url = "artifact://$figpath" - print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n") + function linkfig(figpath, alt = "") + # buildkite-agent upload figpath + # link figure in logs if we are running on CI + if get(ENV, "BUILDKITE", "") == "true" + artifact_url = "artifact://$figpath" + print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n") + end end -end -linkfig( - relpath(joinpath(path, "advect_end.png"), joinpath(@__DIR__, "../..")), - "Advect End Simulation", -) + linkfig( + relpath(joinpath(path, "$(sim_type)_end.png"), joinpath(@__DIR__, "../..")), + "Advect End Simulation", + ) +end diff --git a/examples/column/advect_diffusion.jl b/examples/column/advect_diffusion.jl index e2bf74f6d6..9f4d61d6fc 100644 --- a/examples/column/advect_diffusion.jl +++ b/examples/column/advect_diffusion.jl @@ -10,6 +10,7 @@ import ClimaCore: Geometry, Spaces +import LazyBroadcast: lazy using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 import Logging @@ -52,37 +53,30 @@ end T = gaussian.(zc.z, -0; μ = μ, δ = δ, ν = ν, 𝓌 = 𝓌) V = Geometry.WVector.(ones(FT, fs)) -# Solve Adv-Diff Equation: ∂_t T = α ∇²T +# Solve Adv-Diff Equation: ∂_t T = α ∇²T - v⋅∇T z₀ = FT(0) z₁ = FT(10) function ∑tendencies!(dT, T, z, t) - - ic2f = Operators.InterpolateC2F() - bc_vb = Operators.SetValue(FT(gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))) - bc_vt = Operators.SetValue(FT(gaussian(z₁, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))) bc_gb = Operators.SetGradient( Geometry.WVector(FT(∇gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))), ) bc_gt = Operators.SetGradient( Geometry.WVector(FT(∇gaussian(z₁, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))), ) - - # Upwind Biased Product - # UB = Operators.UpwindBiasedProductC2F( - # bottom = Operators.Extrapolate(), - # top = bc_vt, - # ) - # ∂ = Operators.GradientF2C() - # return @. dT = -∂(UB(V, ic2f(T))) - - A = Operators.AdvectionC2C(bottom = bc_vb, top = Operators.Extrapolate()) - - - gradc2f = Operators.GradientC2F(bottom = bc_vb, top = bc_gt) + T_top = Fields.level(T, Fields.nlevels(T)) + T_top_m1 = Fields.level(T, Fields.nlevels(T) - 1) + top_center_left_biased_grad = @. lazy(Geometry.Covariant3Vector(T_top - T_top_m1)) + + bc_gt_lb = Operators.SetGradient(top_center_left_biased_grad) + gradc2f = Operators.GradientC2F(bottom = bc_gb, top = bc_gt) + gradc2f_advect = Operators.GradientC2F(bottom = bc_gb, top = bc_gt_lb) + interpf2c = Operators.InterpolateF2C() divf2c = Operators.DivergenceF2C() - return @. dT = divf2c(ν * gradc2f(T)) - A(V, T) + return @. dT = + divf2c(ν * gradc2f(T)) - + interpf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f_advect(T))) end @show ∑tendencies!(similar(T), T, nothing, 0.0); diff --git a/examples/column/ekman.jl b/examples/column/ekman.jl index 766ada9cb9..bf1367c72b 100644 --- a/examples/column/ekman.jl +++ b/examples/column/ekman.jl @@ -13,6 +13,7 @@ import ClimaCore: using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 +import LazyBroadcast: lazy import Logging import TerminalLoggers Logging.global_logger(TerminalLoggers.TerminalLogger()) @@ -85,24 +86,42 @@ function tendency!(dY, Y, _, t) u_1 = parent(u)[1] v_1 = parent(v)[1] u_wind = sqrt(u_1^2 + v_1^2) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(0.0), - top = Operators.SetValue(0.0), - ) - # u-momentum bcs_bottom = Operators.SetValue(Geometry.WVector(Cd * u_wind * u_1)) # Eq. 4.16 - bcs_top = Operators.SetValue(FT(ug)) # Eq. 4.18 + u_top = Fields.level(u, Fields.nlevels(u)) + u_bottom = Fields.level(u, 1) + bcs_top = + Operators.SetGradient(@. lazy(Geometry.Covariant3Vector(2 * (FT(ug) - u_top)))) # Eq. 4.18 gradc2f = Operators.GradientC2F(top = bcs_top) divf2c = Operators.DivergenceF2C(bottom = bcs_bottom) - @. du = divf2c(ν * gradc2f(u)) + f * (v - vg) - A(w, u) # Eq. 4.8 + bcs_bottom_advection = + Operators.SetGradient(@. lazy(Geometry.Covariant3Vector(2 * u_bottom))) + bcs_top_advection = + Operators.SetGradient(@. lazy(Geometry.Covariant3Vector(-2 * u_top))) + interpf2c = Operators.InterpolateF2C() + gradc2f_advect = + Operators.GradientC2F(top = bcs_top_advection, bottom = bcs_bottom_advection) + @. du = + divf2c(ν * gradc2f(u)) + f * (v - vg) - + interpf2c(Geometry.dot(Geometry.Contravariant3Vector(w), gradc2f_advect(u))) # Eq. 4.8 # v-momentum bcs_bottom = Operators.SetValue(Geometry.WVector(Cd * u_wind * v_1)) # Eq. 4.17 - bcs_top = Operators.SetValue(FT(vg)) # Eq. 4.19 + v_top = Fields.level(v, Fields.nlevels(v)) + v_bottom = Fields.level(v, 1) + bcs_top = + Operators.SetGradient(@. lazy(Geometry.Covariant3Vector(2 * (FT(vg) - v_top)))) # Eq. 4.19 gradc2f = Operators.GradientC2F(top = bcs_top) divf2c = Operators.DivergenceF2C(bottom = bcs_bottom) - @. dv = divf2c(ν * gradc2f(v)) - f * (u - ug) - A(w, v) # Eq. 4.9 + bcs_bottom_advection = + Operators.SetGradient(@. lazy(Geometry.Covariant3Vector(2 * v_bottom))) + bcs_top_advection = + Operators.SetGradient(@. lazy(Geometry.Covariant3Vector(-2 * v_top))) + gradc2f_advect = + Operators.GradientC2F(top = bcs_top_advection, bottom = bcs_bottom_advection) + @. dv = + divf2c(ν * gradc2f(v)) - f * (u - ug) - + interpf2c(Geometry.dot(Geometry.Contravariant3Vector(w), gradc2f_advect(v))) # Eq. 4.9 return dY end diff --git a/examples/column/heat.jl b/examples/column/heat.jl index fb252db960..6dc5ed16b2 100644 --- a/examples/column/heat.jl +++ b/examples/column/heat.jl @@ -9,7 +9,7 @@ import ClimaCore: Operators, Geometry, Spaces - +import LazyBroadcast: lazy using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 import Logging @@ -38,6 +38,9 @@ T = Fields.zeros(FT, cs) function ∑tendencies!(dT, T, _, t) bcs_bottom = Operators.SetValue(FT(0.0)) + bottom_level_T = Fields.level(T, 1) + bottom_grad = @. lazy(2.0 * Geometry.Covariant3Vector(bottom_level_T)) + bcs_bottom = Operators.SetGradient(bottom_grad) bcs_top = Operators.SetGradient(Geometry.WVector(FT(1.0))) gradc2f = Operators.GradientC2F(bottom = bcs_bottom, top = bcs_top) diff --git a/examples/column/hydrostatic_ekman.jl b/examples/column/hydrostatic_ekman.jl index b73fa057f3..5449a7bb34 100644 --- a/examples/column/hydrostatic_ekman.jl +++ b/examples/column/hydrostatic_ekman.jl @@ -136,18 +136,26 @@ function tendency!(dY, Y, _, t) uv_1 = Operators.getidx(axes(uv), uv, 1) u_wind = LinearAlgebra.norm(uv_1) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(Geometry.UVVector(0.0, 0.0)), - top = Operators.SetValue(Geometry.UVVector(0.0, 0.0)), - ) # uv bcs_bottom = Operators.SetValue(Geometry.WVector(Cd * u_wind) ⊗ uv_1) - bcs_top = Operators.SetValue(uvg) + uv_top = Fields.level(uv, Fields.nlevels(uv)) + uv_top_val = Fields.field_values(uv_top)[] + uv_bottom = Fields.level(uv, Fields.nlevels(uv)) + uv_bottom_val = Fields.field_values(uv_bottom)[] + bcs_top = Operators.SetGradient(Geometry.Covariant3Vector(1) ⊗ ((uvg - uv_top_val))) ∂c = Operators.DivergenceF2C(bottom = bcs_bottom) ∂f = Operators.GradientC2F(top = bcs_top) + bcs_bottom_advection = + Operators.SetGradient(Geometry.Covariant3Vector(2) ⊗ uv_bottom_val) + bcs_top_advection = Operators.SetGradient(Geometry.Covariant3Vector(-2) ⊗ uv_top_val) + interpf2c = Operators.InterpolateF2C() + gradc2f_advect = + Operators.GradientC2F(top = bcs_top_advection, bottom = bcs_bottom_advection) duv .= (uv .- Ref(uvg)) .× Ref(Geometry.WVector(f)) - @. duv += ∂c(ν * ∂f(uv)) - A(w, uv) + @. duv += + ∂c(ν * ∂f(uv)) - + adjoint(interpf2c(adjoint(Geometry.Contravariant3Vector(w)) * gradc2f_advect(uv))) # w If = Operators.InterpolateC2F( @@ -156,7 +164,8 @@ function tendency!(dY, Y, _, t) ) ∂f = Operators.GradientC2F() ∂c = Operators.GradientF2C() - Af = Operators.AdvectionF2F() + gradf2c = Operators.GradientF2C() + interpc2f = Operators.InterpolateC2F() divf = Operators.DivergenceC2F() B = Operators.SetBoundaryOperator( bottom = Operators.SetValue(Geometry.WVector(zero(FT))), @@ -164,7 +173,8 @@ function tendency!(dY, Y, _, t) ) @. dw = B( Geometry.WVector(-(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z))) + - divf(ν * ∂c(w)) - Af(w, w), + divf(ν * ∂c(w)) - + adjoint(interpc2f(adjoint(Geometry.Contravariant3Vector(uv)) * gradf2c(w))), ) return dY @@ -175,12 +185,12 @@ dY = tendency!(similar(Y), Y, nothing, 0.0) Δt = 1.0 / 100.0 # Solve the ODE operator -prob = ODEProblem(tendency!, Y, (0.0, 60 * 60)) +prob = ODEProblem(tendency!, Y, (0.0, 60 * 60 * 10)) sol = solve( prob, SSPRK33(), dt = Δt, - saveat = collect(0.0:600:(60 * 60)), # save every hour + saveat = collect(0.0:600:(60 * 60 * 10)), # save every hour progress = true, progress_message = (dt, u, p, t) -> t, ); diff --git a/examples/column/step.jl b/examples/column/step.jl deleted file mode 100644 index 12e68a743f..0000000000 --- a/examples/column/step.jl +++ /dev/null @@ -1,223 +0,0 @@ -import ClimaComms -ClimaComms.@import_required_backends -using LinearAlgebra -import ClimaCore: - Fields, - Domains, - Topologies, - Meshes, - DataLayouts, - Operators, - Geometry, - Spaces - -using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 - -import Logging -import TerminalLoggers -Logging.global_logger(TerminalLoggers.TerminalLogger()) - -const FT = Float64 - -a = Geometry.ZPoint(FT(-20.0)) -b = Geometry.ZPoint(FT(20.0)) -n = 64 -α = FT(0.1) - -function heaviside(pt) - 0.5 * (sign(pt.z) + 1) -end - -domain = Domains.IntervalDomain(a, b, boundary_names = (:left, :right)) -mesh = Meshes.IntervalMesh(domain, nelems = n) -device = ClimaComms.device() -cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) -fs = Spaces.FaceFiniteDifferenceSpace(cs) - -V = Geometry.WVector.(ones(FT, fs)) -θ = heaviside.(Fields.coordinate_field(cs)) - -# Solve advection Equation: ∂θ/dt = -∂(vθ) - -# upwinding -function tendency1!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a.z - t)), - right = Operators.SetValue(sin(b.z - t)), - ) - ∂ = Operators.DivergenceF2C() - - return @. dθ = -∂(UB(V, θ)) -end -function tendency2!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a.z - t)), - right = Operators.SetValue(sin(b.z - t)), - ) - ∂ = Operators.DivergenceF2C() - return @. dθ = -∂(UB(V, θ)) + fcc(V, θ) -end -# use the advection operator -function tendency3!(dθ, θ, _, t) - - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - A = Operators.AdvectionC2C( - left = Operators.SetValue(sin(-t)), - right = Operators.Extrapolate(), - ) - return @. dθ = -A(V, θ) -end -# use the advection operator -function tendency4!(dθ, θ, _, t) - - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - A = Operators.AdvectionC2C( - left = Operators.SetValue(sin(-t)), - right = Operators.Extrapolate(), - ) - return @. dθ = -A(V, θ) + fcc(V, θ) -end - -# use the advection operator - -@show tendency1!(similar(θ), θ, nothing, 0.0) -# Solve the ODE operator -Δt = 0.001 -prob1 = ODEProblem(tendency1!, θ, (0.0, 5.0)) -prob2 = ODEProblem(tendency2!, θ, (0.0, 5.0)) -prob3 = ODEProblem(tendency3!, θ, (0.0, 5.0)) -prob4 = ODEProblem(tendency4!, θ, (0.0, 5.0)) -sol1 = solve( - prob1, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):5.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); -sol2 = solve( - prob2, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):5.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); -sol3 = solve( - prob3, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):5.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); -sol4 = solve( - prob4, - SSPRK33(), - dt = Δt, - saveat = collect(0.0:(10 * Δt):5.0), - progress = true, - progress_message = (dt, u, p, t) -> t, -); - -ENV["GKSwstype"] = "nul" -import ClimaCorePlots, Plots -Plots.GRBackend() - -dir = "advect_step_function" -path = joinpath(@__DIR__, "output", dir) -mkpath(path) - -anim = Plots.@animate for u in sol1.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "UBP_advect_step_function.mp4"), fps = 10) -Plots.png( - Plots.plot(sol1.u[end], xlim = (-1, 1)), - joinpath(path, "sol1_advect_step_function_end.png"), -) - -anim = Plots.@animate for u in sol2.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "UBP_advect_step_function_fc.mp4"), fps = 10) -Plots.png( - Plots.plot(sol2.u[end], xlim = (-1, 1)), - joinpath(path, "sol2_advect_step_function_end.png"), -) - -anim = Plots.@animate for u in sol3.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "C2C_advect_step_function.mp4"), fps = 10) -Plots.png( - Plots.plot(sol3.u[end], xlim = (-1, 1)), - joinpath(path, "sol3_advect_step_function_end.png"), -) - -anim = Plots.@animate for u in sol4.u - Plots.plot(u, xlim = (-1, 1)) -end -Plots.mp4(anim, joinpath(path, "C2C_advect_step_function_fc.mp4"), fps = 10) -Plots.png( - Plots.plot(sol4.u[end], xlim = (-1, 1)), - joinpath(path, "sol4_advect_step_function_end.png"), -) - -p = Plots.plot(sol1.u[end], xlim = (-1, 1), label = "UBP") -p = Plots.plot!(sol2.u[end], xlim = (-1, 1), label = "UBP_FC") -p = Plots.plot!(sol3.u[end], xlim = (-1, 1), label = "C2C") -p = Plots.plot!(sol4.u[end], xlim = (-1, 1), label = "C2C_FC") -Plots.png(p, joinpath(path, "all_advect_step_function_end.png")) - -function linkfig(figpath, alt = "") - # buildkite-agent upload figpath - # link figure in logs if we are running on CI - if get(ENV, "BUILDKITE", "") == "true" - artifact_url = "artifact://$figpath" - print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n") - end -end - -linkfig( - "examples/column/output/$(dir)/advect_step_function_end.png", - "Advect End Simulation", -) - -linkfig( - relpath( - joinpath(path, "advect_step_function_end.png"), - joinpath(@__DIR__, "../.."), - ), - "Advect End Simulation", -) diff --git a/examples/column/wave.jl b/examples/column/wave.jl index cadd4358dc..5b81d13f62 100644 --- a/examples/column/wave.jl +++ b/examples/column/wave.jl @@ -40,10 +40,11 @@ function tendency!(dY, Y, _, t) du = dY.u dp = dY.p - + u_left = Fields.field_values(Fields.level(u, 1))[] + u_right = Fields.field_values(Fields.level(u, Fields.nlevels(u)))[] ∂f = Operators.GradientC2F( - left = Operators.SetValue(0.0), - right = Operators.SetValue(0.0), + left = Operators.SetGradient(Geometry.Covariant3Vector(2 * u_left)), + right = Operators.SetGradient(Geometry.Covariant3Vector(-2 * u_right)), ) ∂c = Operators.DivergenceF2C() diff --git a/examples/hybrid/box/bubble_3d_invariant_rhoe.jl b/examples/hybrid/box/bubble_3d_invariant_rhoe.jl index fbd29e0484..32089e40ae 100644 --- a/examples/hybrid/box/bubble_3d_invariant_rhoe.jl +++ b/examples/hybrid/box/bubble_3d_invariant_rhoe.jl @@ -36,6 +36,7 @@ using Logging using CUDA CUDA.allowscalar(false) +include("../flux_correction_utils.jl") """ SimulationParameters{FT} @@ -330,19 +331,8 @@ function rhs_invariant!(dY, Y, ghost_buffer, t) @. dρe -= hdiv(cuvw * (cρe + cp)) @. dρe -= vdivf2c(fw * Ic2f(cρe + cp)) @. dρe -= vdivf2c(Ic2f(cuₕ * (cρe + cp))) - - fcc = Operators.FluxCorrectionC2C( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - - @. dρ += fcc(fw, cρ) - @. dρe += fcc(fw, cρe) - # dYc.ρuₕ += fcc(w, Yc.ρuₕ) + add_flux_correction_c2c(dρe, fw, cρe) + add_flux_correction_c2c(dρ, fw, cρ) Spaces.weighted_dss_start!(dY.Yc, ghost_buffer.Yc) Spaces.weighted_dss_start!(dY.uₕ, ghost_buffer.uₕ) diff --git a/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl b/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl index d5486fed7c..34e6c48852 100644 --- a/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl +++ b/examples/hybrid/box/bubble_3d_invariant_rhotheta.jl @@ -23,6 +23,7 @@ import TerminalLoggers Logging.global_logger(TerminalLoggers.TerminalLogger()) const context = ClimaComms.SingletonCommsContext() +include("../flux_correction_utils.jl") function hvspace_3D( xlim = (-π, π), @@ -279,18 +280,8 @@ function rhs_invariant!(dY, Y, _, t) @. dρθ -= vdivf2c(fw * Ic2f(cρθ)) @. dρθ -= vdivf2c(Ic2f(cuₕ * cρθ)) - fcc = Operators.FluxCorrectionC2C( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - - @. dρ += fcc(fw, cρ) - @. dρθ += fcc(fw, cρθ) - # dYc.ρuₕ += fcc(w, Yc.ρuₕ) + add_flux_correction_c2c(dρ, fw, cρ) + add_flux_correction_c2c(dρθ, fw, cρθ) Spaces.weighted_dss!(dY.Yc) Spaces.weighted_dss!(dY.uₕ) diff --git a/examples/hybrid/box/bubble_3d_rhotheta.jl b/examples/hybrid/box/bubble_3d_rhotheta.jl index eaedf2184c..ffd1fa254d 100644 --- a/examples/hybrid/box/bubble_3d_rhotheta.jl +++ b/examples/hybrid/box/bubble_3d_rhotheta.jl @@ -22,6 +22,7 @@ import Logging import TerminalLoggers Logging.global_logger(TerminalLoggers.TerminalLogger()) const context = ClimaComms.SingletonCommsContext() +include("../flux_correction_utils.jl") # set up function space function hvspace_3D( @@ -213,14 +214,6 @@ function rhs!(dY, Y, _, t) top = Operators.SetValue(Geometry.WVector(0.0)), ) - fcc = Operators.FluxCorrectionC2C( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) uₕ = @. ρuₕ / ρ w = @. ρw / If(ρ) @@ -272,10 +265,11 @@ function rhs!(dY, Y, _, t) ### UPWIND FLUX CORRECTION upwind_correction = true if upwind_correction - @. dρ += fcc(w, ρ) - @. dρθ += fcc(w, ρθ) - @. dρuₕ += fcc(w, ρuₕ) - @. dρw += fcf(wc, ρw) + add_flux_correction_c2c(dρ, w, ρ) + add_flux_correction_c2c(dρθ, w, ρθ) + add_flux_correction_c2c(dρuₕ, w, ρuₕ) + add_flux_correction_f2f(dρw, wc, ρw) + add_flux_correction_f2f(dρw, wc, ρw) end ### DIFFUSION diff --git a/examples/hybrid/flux_correction_utils.jl b/examples/hybrid/flux_correction_utils.jl new file mode 100644 index 0000000000..23aa0e2a99 --- /dev/null +++ b/examples/hybrid/flux_correction_utils.jl @@ -0,0 +1,100 @@ +using LinearAlgebra +import ClimaCore: + ClimaCore, + DataLayouts, + Spaces, + Domains, + Meshes, + Geometry, + Topologies, + Spaces, + Quadratures, + Fields, + Operators +using ClimaCore.Geometry +import ClimaCore.Geometry: ⊗ + +function add_flux_correction_c2c(d_, velocity, quantity) + n_levels = Fields.nlevels(quantity) + top_level_space = axes(Fields.level(quantity, n_levels)) + bottom_level_space = axes(Fields.level(quantity, 1)) + quantity_top_gradient_extrapolate = Operators.SetGradient( + Geometry.outer.( + (Geometry.Covariant3Vector(1),), + Fields.Field( + Fields.field_values(Fields.level(quantity, n_levels)) .- + Fields.field_values(Fields.level(quantity, n_levels - 1)), + top_level_space, + ), + ), + ) + quantity_bottom_gradient_extrapolate = Operators.SetGradient( + Geometry.outer.( + (Geometry.Covariant3Vector(1),), + Fields.Field( + ( + Fields.field_values(Fields.level(quantity, 2)) .- + Fields.field_values(Fields.level(quantity, 1)) + ), + bottom_level_space, + ), + )) + quantity_gradc2f = Operators.GradientC2F( + bottom = quantity_bottom_gradient_extrapolate, + top = quantity_top_gradient_extrapolate, + ) + lg_field = Fields.local_geometry_field(axes(velocity)) + gradf2c = Operators.GradientF2C() + @. d_ += + adjoint( + gradf2c( + adjoint(quantity_gradc2f(quantity)) * Geometry.Contravariant3Vector( + abs(Geometry.contravariant3(velocity, lg_field)), + ), + ), + ) * Geometry.Contravariant3Vector(1) + return +end + +function add_flux_correction_f2f(d_, velocity, quantity) + gradf2c = Operators.GradientF2C() + lg_field = Fields.local_geometry_field(axes(velocity)) + inner_grad = @. adjoint(gradf2c(quantity)) * Geometry.Contravariant3Vector( + abs(Geometry.contravariant3(velocity, lg_field)), + ) + n_levels = Fields.nlevels(inner_grad) + top_level_space = axes(Fields.level(inner_grad, n_levels)) + bottom_level_space = axes(Fields.level(inner_grad, 1)) + top_gradient_extrapolate = Operators.SetGradient( + Geometry.outer.( + (Geometry.Covariant3Vector(1),), + Fields.Field( + Fields.field_values(Fields.level(inner_grad, n_levels)) .- + Fields.field_values(Fields.level(inner_grad, n_levels - 1)), + top_level_space, + ), + ), + ) + bottom_gradient_extrapolate = Operators.SetGradient( + Geometry.outer.( + (Geometry.Covariant3Vector(1),), + Fields.Field( + ( + Fields.field_values(Fields.level(inner_grad, 1)) .- + Fields.field_values(Fields.level(inner_grad, 2)) + ), + bottom_level_space, + ), + )) + gradc2f = Operators.GradientC2F( + bottom = bottom_gradient_extrapolate, + top = top_gradient_extrapolate, + ) + lg_field = Fields.local_geometry_field(axes(velocity)) + @. d_ += + adjoint( + gradc2f(inner_grad + ), + ) * Geometry.Contravariant3Vector(1) + return +end diff --git a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl index b2f4da689f..6e2cc8c998 100644 --- a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl +++ b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl @@ -228,18 +228,6 @@ function rhs_invariant!(dY, Y, _, t) @. dρe -= vdivf2c((Ic2f(cρ) * third_order_upwind_c2f(fw, (cρe + cp) / cρ))) @. dρe -= vdivf2c(Ic2f(cuₕ * (cρe + cp))) - fcc = Operators.FluxCorrectionC2C( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - - # Flux correction (Upwind Correction to Central scheme) - # @. dρ += fcc(fw, cρ) - # @. dρe += fcc(fw, cρe) Spaces.weighted_dss!(dY.Yc) Spaces.weighted_dss!(dY.uₕ) diff --git a/examples/hybrid/plane/density_current_2d_flux_form.jl b/examples/hybrid/plane/density_current_2d_flux_form.jl index faf1f18630..81dff72203 100644 --- a/examples/hybrid/plane/density_current_2d_flux_form.jl +++ b/examples/hybrid/plane/density_current_2d_flux_form.jl @@ -20,6 +20,7 @@ import ClimaCore.Geometry: ⊗ import Logging import TerminalLoggers Logging.global_logger(TerminalLoggers.TerminalLogger()) +include("../flux_correction_utils.jl") # set up function space function hvspace_2D( @@ -201,15 +202,6 @@ function rhs!(dY, Y, _, t) top = Operators.SetValue(Geometry.WVector(0.0)), ) - fcc = Operators.FluxCorrectionC2C( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) - uₕ = @. ρuₕ / ρ w = @. ρw / If(ρ) wc = @. Ic(ρw) / ρ @@ -255,15 +247,6 @@ function rhs!(dY, Y, _, t) uₕf = @. If(ρuₕ / ρ) # requires boundary conditions @. dρw -= hdiv(uₕf ⊗ ρw) - ### UPWIND FLUX CORRECTION - upwind_correction = false - if upwind_correction - @. dρ += fcc(w, ρ) - @. dρθ += fcc(w, ρθ) - @. dρuₕ += fcc(w, ρuₕ) - @. dρw += fcf(wc, ρw) - end - ### DIFFUSION κ₂ = 75.0 # m^2/s # 1a) horizontal div of horizontal grad of horiz momentun diff --git a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl index 20df8bb6a8..7fc4e18c4f 100644 --- a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl +++ b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl @@ -16,8 +16,10 @@ import ClimaCore: Spaces, Quadratures, Fields, - Operators + Operators, + Utilities using ClimaCore.Geometry +import LazyBroadcast: lazy using Logging: global_logger using TerminalLoggers: TerminalLogger @@ -177,10 +179,7 @@ function rhs_invariant!(dY, Y, _, t) top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), ) - vdivc2f = Operators.DivergenceC2F( - top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - ) + vdivc2f = Operators.DivergenceC2F() # we want the total u³ at the boundary to be zero: we can either constrain # both to be zero, or allow one to be non-zero and set the other to be its # negation @@ -246,7 +245,32 @@ function rhs_invariant!(dY, Y, _, t) hκ₂∇²uₕ = @. hwdiv(κ₂ * ᶜ∇ₕuₕ) vκ₂∇²uₕ = @. vdivf2c(κ₂ * ᶠ∇ᵥuₕ) hκ₂∇²w = @. hwdiv(κ₂ * ᶠ∇ₕw) - vκ₂∇²w = @. vdivc2f(κ₂ * ᶜ∇ᵥw) + + lg_field_faces = Fields.local_geometry_field(axes(fw)) + lg_field_centers = Fields.local_geometry_field(axes(cρ)) + lg_bottom_face = Operators.Fields.level(Operators.RightBiasedF2C().(lg_field_faces), 1) + lg_top_face = Fields.level( + Operators.LeftBiasedF2C().(lg_field_faces), + Fields.nlevels(lg_field_centers), + ) + lg_bottom_center = Fields.level(lg_field_centers, 1) + lg_top_center = Fields.level(lg_field_centers, Fields.nlevels(lg_field_centers)) + ᶜ∇ᵥw_bottom = Fields.level(ᶜ∇ᵥw, 1) + ᶜ∇ᵥw_top = Fields.level(ᶜ∇ᵥw, Fields.nlevels(ᶜ∇ᵥw)) + bottom_divergence = @. lazy( + Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * + (2 * inv(lg_bottom_face.J)), + ) + top_divergence = + @. lazy( + Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * inv(lg_top_face.J)), + ) + set_bcs = Operators.SetBoundaryOperator( + bottom = Operators.SetValue(bottom_divergence), + top = Operators.SetValue(top_divergence), + ) + + vκ₂∇²w = @. set_bcs(vdivc2f(κ₂ * ᶜ∇ᵥw)) hκ₂∇²h_tot = @. hwdiv(cρ * κ₂ * ᶜ∇ₕh_tot) vκ₂∇²h_tot = @. vdivf2c(fρ * κ₂ * ᶠ∇ᵥh_tot) diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl index e4396f1960..f254ebda84 100644 --- a/examples/hybrid/plane/topo_agnesi_nh.jl +++ b/examples/hybrid/plane/topo_agnesi_nh.jl @@ -16,6 +16,7 @@ import ClimaCore: Operators, Hypsography using ClimaCore.Geometry +import LazyBroadcast: lazy using DiffEqCallbacks @@ -233,10 +234,7 @@ function rhs_invariant!(dY, Y, _, t) top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), ) - vdivc2f = Operators.DivergenceC2F( - top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - ) + vdivc2f = Operators.DivergenceC2F() # we want the total u³ at the boundary to be zero: we can either constrain # both to be zero, or allow one to be non-zero and set the other to be its # negation @@ -303,7 +301,30 @@ function rhs_invariant!(dY, Y, _, t) hκ₂∇²uₕ = @. hwdiv(κ₂ * ᶜ∇ₕuₕ) vκ₂∇²uₕ = @. vdivf2c(κ₂ * ᶠ∇ᵥuₕ) hκ₂∇²w = @. hwdiv(κ₂ * ᶠ∇ₕw) - vκ₂∇²w = @. vdivc2f(κ₂ * ᶜ∇ᵥw) + lg_field_faces = Fields.local_geometry_field(axes(fw)) + lg_field_centers = Fields.local_geometry_field(axes(cρ)) + lg_bottom_face = Operators.Fields.level(Operators.RightBiasedF2C().(lg_field_faces), 1) + lg_top_face = Fields.level( + Operators.LeftBiasedF2C().(lg_field_faces), + Fields.nlevels(lg_field_centers), + ) + lg_bottom_center = Fields.level(lg_field_centers, 1) + lg_top_center = Fields.level(lg_field_centers, Fields.nlevels(lg_field_centers)) + ᶜ∇ᵥw_bottom = Fields.level(ᶜ∇ᵥw, 1) + ᶜ∇ᵥw_top = Fields.level(ᶜ∇ᵥw, Fields.nlevels(ᶜ∇ᵥw)) + bottom_divergence = @. lazy( + Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * + (2 * inv(lg_bottom_face.J)), + ) + top_divergence = + @. lazy( + Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * inv(lg_top_face.J)), + ) + set_bcs = Operators.SetBoundaryOperator( + bottom = Operators.SetValue(bottom_divergence), + top = Operators.SetValue(top_divergence), + ) + vκ₂∇²w = @. set_bcs(vdivc2f(κ₂ * ᶜ∇ᵥw)) hκ₂∇²h_tot = @. hwdiv(cρ * κ₂ * ᶜ∇ₕh_tot) vκ₂∇²h_tot = @. vdivf2c(fρ * κ₂ * ᶠ∇ᵥh_tot) diff --git a/examples/hybrid/plane/topo_schar_nh.jl b/examples/hybrid/plane/topo_schar_nh.jl index b1efe08ce2..c83b5ec7db 100644 --- a/examples/hybrid/plane/topo_schar_nh.jl +++ b/examples/hybrid/plane/topo_schar_nh.jl @@ -308,10 +308,6 @@ function rhs_invariant!(dY, Y, _, t) top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), ) - vdivc2f = Operators.DivergenceC2F( - top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - ) # we want the total u³ at the boundary to be zero: we can either constrain # both to be zero, or allow one to be non-zero and set the other to be its # negation @@ -475,10 +471,6 @@ if tendency_plots top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), ) - vdivc2f = Operators.DivergenceC2F( - top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - ) cw = @. If2c(fw) fuₕ = @. Ic2f(cuₕ) @@ -593,10 +585,6 @@ if tendency_plots top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), ) - vdivc2f = Operators.DivergenceC2F( - top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)), - ) cw = @. If2c(fw) fuₕ = @. Ic2f(cuₕ) diff --git a/examples/hybrid/staggered_nonhydrostatic_model.jl b/examples/hybrid/staggered_nonhydrostatic_model.jl index 52ed153cbe..18bfe397ea 100644 --- a/examples/hybrid/staggered_nonhydrostatic_model.jl +++ b/examples/hybrid/staggered_nonhydrostatic_model.jl @@ -53,10 +53,6 @@ const ᶠcurlᵥ = Operators.CurlC2F( bottom = Operators.SetCurl(CT12(FT(0), FT(0))), top = Operators.SetCurl(CT12(FT(0), FT(0))), ) -const ᶜFC = Operators.FluxCorrectionC2C( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), -) const ᶠupwind_product1 = Operators.UpwindBiasedProductC2F() const ᶠupwind_product3 = Operators.Upwind3rdOrderBiasedProductC2F( bottom = Operators.ThirdOrderOneSided(), @@ -190,16 +186,6 @@ function implicit_tendency!(Yₜ, Y, p, t) @. Yₜ.f.w = -(ᶠgradᵥ(ᶜp) / ᶠinterp(ᶜρ) + ᶠgradᵥ(ᶜK + ᶜΦ)) - # TODO: Add flux correction to the Jacobian - # @. Yₜ.c.ρ += ᶜFC(ᶠw, ᶜρ) - # if :ρθ in propertynames(Y.c) - # @. Yₜ.c.ρθ += ᶜFC(ᶠw, ᶜρθ) - # elseif :ρe in propertynames(Y.c) - # @. Yₜ.c.ρe += ᶜFC(ᶠw, ᶜρe) - # elseif :ρe_int in propertynames(Y.c) - # @. Yₜ.c.ρe_int += ᶜFC(ᶠw, ᶜρe_int) - # end - return Yₜ end diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index f84d1ee0b7..e445e2cdcf 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -19,13 +19,9 @@ const TwoArgFDOperatorWithCenterInput = Union{ Operators.WeightedInterpolateC2F, Operators.UpwindBiasedProductC2F, Operators.Upwind3rdOrderBiasedProductC2F, - Operators.AdvectionC2C, - Operators.FluxCorrectionC2C, } const TwoArgFDOperatorWithFaceInput = Union{ Operators.WeightedInterpolateF2C, - Operators.AdvectionF2F, - Operators.FluxCorrectionF2F, } const ErroneousFDOperator = Union{ Operators.LeftBiased3rdOrderC2F, @@ -33,7 +29,8 @@ const ErroneousFDOperator = Union{ Operators.RightBiased3rdOrderC2F, Operators.RightBiased3rdOrderF2C, } -const NonlinearFDOperator = Union{Operators.FCTBorisBook, Operators.FCTZalesak} +const NonlinearFDOperator = + Union{Operators.FCTBorisBook, Operators.FCTZalesak, Operators.LinVanLeerC2F} const OneArgFDOperator = Union{OneArgFDOperatorWithCenterInput, OneArgFDOperatorWithFaceInput} @@ -390,16 +387,6 @@ op_matrix_last_row( ::Operators.SetValue, ::Type{FT}, ) where {FT} = LowerDiagonalMatrixRow(FT(0)) -op_matrix_first_row( - ::Operators.InterpolateC2F, - ::Operators.SetGradient, - ::Type{FT}, -) where {FT} = UpperDiagonalMatrixRow(FT(1)) -op_matrix_last_row( - ::Operators.InterpolateC2F, - ::Operators.SetGradient, - ::Type{FT}, -) where {FT} = LowerDiagonalMatrixRow(FT(1)) op_matrix_first_row( ::Operators.InterpolateC2F, ::Operators.Extrapolate, @@ -468,16 +455,6 @@ op_matrix_last_row( ::Operators.SetValue, ::Type{FT}, ) where {FT} = LowerDiagonalMatrixRow(FT(0)) -op_matrix_first_row( - ::Operators.WeightedInterpolateC2F, - ::Operators.SetGradient, - ::Type{FT}, -) where {FT} = UpperDiagonalMatrixRow(FT(1)) -op_matrix_last_row( - ::Operators.WeightedInterpolateC2F, - ::Operators.SetGradient, - ::Type{FT}, -) where {FT} = LowerDiagonalMatrixRow(FT(1)) op_matrix_first_row( ::Operators.WeightedInterpolateC2F, ::Operators.Extrapolate, @@ -507,30 +484,6 @@ Base.@propagate_inbounds function op_matrix_interior_row( av³ = CT3(abs(v³.u³)) return BidiagonalMatrixRow(v³ + av³, v³ - av³) / 2 end -Base.@propagate_inbounds function op_matrix_first_row( - ::Operators.UpwindBiasedProductC2F, - ::Operators.SetValue, - space, - idx, - hidx, - velocity, -) - v³ = CT3(ct3_data(velocity, space, idx, hidx)) - av³ = CT3(abs(v³.u³)) - return UpperDiagonalMatrixRow(v³ - av³) / 2 -end -Base.@propagate_inbounds function op_matrix_last_row( - ::Operators.UpwindBiasedProductC2F, - ::Operators.SetValue, - space, - idx, - hidx, - velocity, -) - v³ = CT3(ct3_data(velocity, space, idx, hidx)) - av³ = CT3(abs(v³.u³)) - return LowerDiagonalMatrixRow(v³ + av³) / 2 -end Base.@propagate_inbounds function op_matrix_first_row( ::Operators.UpwindBiasedProductC2F, ::Operators.Extrapolate, @@ -622,107 +575,6 @@ end op_matrix_row_type(::Operators.AdvectionOperator, ::Type{FT}, _) where {FT} = TridiagonalMatrixRow{FT} -Base.@propagate_inbounds function op_matrix_interior_row( - ::Operators.AdvectionC2C, - space, - idx, - hidx, - velocity, -) - v³⁻_data = ct3_data(velocity, space, idx - half, hidx) - v³⁺_data = ct3_data(velocity, space, idx + half, hidx) - return TridiagonalMatrixRow(-v³⁻_data, v³⁻_data - v³⁺_data, v³⁺_data) / 2 -end -Base.@propagate_inbounds function op_matrix_first_row( - ::Operators.AdvectionC2C, - ::Operators.SetValue, - space, - idx, - hidx, - velocity, -) - v³⁻_data = ct3_data(velocity, space, idx - half, hidx) - v³⁺_data = ct3_data(velocity, space, idx + half, hidx) - return UpperBidiagonalSquareMatrixRow(2v³⁻_data - v³⁺_data, v³⁺_data) / 2 -end -Base.@propagate_inbounds function op_matrix_last_row( - ::Operators.AdvectionC2C, - ::Operators.SetValue, - space, - idx, - hidx, - velocity, -) - v³⁻_data = ct3_data(velocity, space, idx - half, hidx) - v³⁺_data = ct3_data(velocity, space, idx + half, hidx) - return LowerBidiagonalSquareMatrixRow(-v³⁻_data, v³⁻_data - 2v³⁺_data) / 2 -end -Base.@propagate_inbounds function op_matrix_first_row( - ::Operators.AdvectionC2C, - ::Operators.Extrapolate, - space, - idx, - hidx, - velocity, -) - v³⁺_data = ct3_data(velocity, space, idx + half, hidx) - return UpperBidiagonalSquareMatrixRow(-v³⁺_data, v³⁺_data) -end -Base.@propagate_inbounds function op_matrix_last_row( - ::Operators.AdvectionC2C, - ::Operators.Extrapolate, - space, - idx, - hidx, - velocity, -) - v³⁻_data = ct3_data(velocity, space, idx - half, hidx) - return LowerBidiagonalSquareMatrixRow(-v³⁻_data, v³⁻_data) -end -Base.@propagate_inbounds function op_matrix_interior_row( - ::Operators.AdvectionF2F, - space, - idx, - hidx, - velocity, -) - FT = Spaces.undertype(space) - v³_data = ct3_data(velocity, space, idx, hidx) - return TridiagonalMatrixRow(-v³_data, FT(0), v³_data) / 2 -end -Base.@propagate_inbounds function op_matrix_interior_row( - ::Union{Operators.FluxCorrectionC2C, Operators.FluxCorrectionF2F}, - space, - idx, - hidx, - velocity, -) - av³⁻_data = abs(ct3_data(velocity, space, idx - half, hidx)) - av³⁺_data = abs(ct3_data(velocity, space, idx + half, hidx)) - return TridiagonalMatrixRow(av³⁻_data, -av³⁻_data - av³⁺_data, av³⁺_data) -end -Base.@propagate_inbounds function op_matrix_first_row( - ::Union{Operators.FluxCorrectionC2C, Operators.FluxCorrectionF2F}, - ::Operators.Extrapolate, - space, - idx, - hidx, - velocity, -) - av³⁺_data = abs(ct3_data(velocity, space, idx + half, hidx)) - return UpperBidiagonalSquareMatrixRow(-av³⁺_data, av³⁺_data) -end -Base.@propagate_inbounds function op_matrix_last_row( - ::Union{Operators.FluxCorrectionC2C, Operators.FluxCorrectionF2F}, - ::Operators.Extrapolate, - space, - idx, - hidx, - velocity, -) - av³⁻_data = abs(ct3_data(velocity, space, idx - half, hidx)) - return LowerBidiagonalSquareMatrixRow(av³⁻_data, -av³⁻_data) -end op_matrix_interior_row(::Operators.SetBoundaryOperator, ::Type{FT}) where {FT} = DiagonalMatrixRow(FT(1)) @@ -742,16 +594,6 @@ op_matrix_row_type(op::Operators.GradientOperator, ::Type{FT}) where {FT} = BidiagonalMatrixRow{C3{FT}} op_matrix_interior_row(::Operators.GradientOperator, ::Type{FT}) where {FT} = BidiagonalMatrixRow(-C3(FT(1)), C3(FT(1))) -op_matrix_first_row( - ::Operators.GradientC2F, - ::Operators.SetValue, - ::Type{FT}, -) where {FT} = UpperDiagonalMatrixRow(C3(FT(2))) -op_matrix_last_row( - ::Operators.GradientC2F, - ::Operators.SetValue, - ::Type{FT}, -) where {FT} = LowerDiagonalMatrixRow(-C3(FT(2))) op_matrix_first_row( ::Operators.GradientC2F, ::Operators.SetGradient, @@ -797,28 +639,6 @@ Base.@propagate_inbounds function op_matrix_interior_row( J⁺ = Geometry.LocalGeometry(space, idx + half, hidx).J return BidiagonalMatrixRow(-C3(J⁻)', C3(J⁺)') * invJ end -Base.@propagate_inbounds function op_matrix_first_row( - ::Operators.DivergenceC2F, - ::Operators.SetValue, - space, - idx, - hidx, -) - invJ = Geometry.LocalGeometry(space, idx, hidx).invJ - J⁺ = Geometry.LocalGeometry(space, idx + half, hidx).J - return UpperDiagonalMatrixRow(C3(J⁺)') * 2invJ -end -Base.@propagate_inbounds function op_matrix_last_row( - ::Operators.DivergenceC2F, - ::Operators.SetValue, - space, - idx, - hidx, -) - invJ = Geometry.LocalGeometry(space, idx, hidx).invJ - J⁻ = Geometry.LocalGeometry(space, idx - half, hidx).J - return LowerDiagonalMatrixRow(-C3(J⁻)') * 2invJ -end op_matrix_first_row( ::Operators.DivergenceC2F, ::Operators.SetDivergence, @@ -903,26 +723,6 @@ Base.@propagate_inbounds function op_matrix_interior_row( invJ = Geometry.LocalGeometry(space, idx, hidx).invJ return BidiagonalMatrixRow(-εⁱʲ, εⁱʲ) * invJ end -Base.@propagate_inbounds function op_matrix_first_row( - ::Operators.CurlC2F, - ::Operators.SetValue, - space, - idx, - hidx, -) - invJ = Geometry.LocalGeometry(space, idx, hidx).invJ - return UpperDiagonalMatrixRow(εⁱʲ) * 2invJ -end -Base.@propagate_inbounds function op_matrix_last_row( - ::Operators.CurlC2F, - ::Operators.SetValue, - space, - idx, - hidx, -) - invJ = Geometry.LocalGeometry(space, idx, hidx).invJ - return LowerDiagonalMatrixRow(-εⁱʲ) * 2invJ -end op_matrix_first_row( ::Operators.CurlC2F, ::Operators.SetCurl, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 5896625d6b..1ba1bd2645 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -422,11 +422,6 @@ Supported boundary conditions are: ```math I(x)[\\tfrac{1}{2}] = x₀ ``` -- [`SetGradient(v)`](@ref): set the value at the boundary such that the gradient - is `v`. At the left boundary the stencil is -```math -I(x)[\\tfrac{1}{2}] = x[1] - \\frac{1}{2} v³ -``` - [`Extrapolate`](@ref): use the closest interior point as the boundary value. At the left boundary the stencil is ```math @@ -439,7 +434,7 @@ struct InterpolateC2F{BCS} <: InterpolationOperator assert_valid_bcs( "InterpolateC2F", kwargs, - (SetValue, SetGradient, Extrapolate), + (SetValue, Extrapolate), ) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end @@ -486,39 +481,6 @@ Base.@propagate_inbounds function stencil_right_boundary( getidx(space, bc.val, nothing, hidx) end -Base.@propagate_inbounds function stencil_left_boundary( - ::InterpolateC2F, - bc::SetGradient, - space, - idx, - hidx, - arg, -) - @assert idx == left_face_boundary_idx(space) - a⁺ = getidx(space, arg, idx + half, hidx) - v₃ = Geometry.covariant3( - getidx(space, bc.val, nothing, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁺ - v₃ / 2 -end -Base.@propagate_inbounds function stencil_right_boundary( - ::InterpolateC2F, - bc::SetGradient, - space, - idx, - hidx, - arg, -) - @assert idx == right_face_boundary_idx(space) - a⁻ = getidx(space, arg, idx - half, hidx) - v₃ = Geometry.covariant3( - getidx(space, bc.val, nothing, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁻ + v₃ / 2 -end - Base.@propagate_inbounds function stencil_left_boundary( ::InterpolateC2F, bc::Extrapolate, @@ -1108,7 +1070,6 @@ WI(w, x)[i] = \\frac{ Supported boundary conditions are: - [`SetValue(val)`](@ref): set the value at the boundary face to be `val`. -- [`SetGradient`](@ref): set the value at the boundary such that the gradient is `val`. - [`Extrapolate`](@ref): use the closest interior point as the boundary value. These have the same stencil as in [`InterpolateC2F`](@ref). @@ -1119,7 +1080,7 @@ struct WeightedInterpolateC2F{BCS} <: WeightedInterpolationOperator assert_valid_bcs( "WeightedInterpolateC2F", kwargs, - (SetValue, SetGradient, Extrapolate), + (SetValue, Extrapolate), ) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end @@ -1175,41 +1136,6 @@ Base.@propagate_inbounds function stencil_right_boundary( getidx(space, bc.val, nothing, hidx) end -Base.@propagate_inbounds function stencil_left_boundary( - ::WeightedInterpolateC2F, - bc::SetGradient, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == left_face_boundary_idx(space) - a⁺ = getidx(space, arg, idx + half, hidx) - v₃ = Geometry.covariant3( - getidx(space, bc.val, nothing, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁺ - v₃ / 2 -end -Base.@propagate_inbounds function stencil_right_boundary( - ::WeightedInterpolateC2F, - bc::SetGradient, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == right_face_boundary_idx(space) - a⁻ = getidx(space, arg, idx - half, hidx) - v₃ = Geometry.covariant3( - getidx(space, bc.val, nothing, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁻ + v₃ / 2 -end - Base.@propagate_inbounds function stencil_left_boundary( ::WeightedInterpolateC2F, bc::Extrapolate, @@ -1259,15 +1185,6 @@ U(\\boldsymbol{v},x)[i] = \\begin{cases} where ``\\boldsymbol{e}_3`` is the 3rd covariant basis vector. Supported boundary conditions are: -- [`SetValue(x₀)`](@ref): set the value of `x` to be `x₀` in a hypothetical - ghost cell on the other side of the boundary. On the left boundary the stencil - is - ```math - U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = \\begin{cases} - v^3[\\tfrac{1}{2}] x_0 \\boldsymbol{e}_3 \\textrm{, if } v^3[\\tfrac{1}{2}] > 0 \\\\ - v^3[\\tfrac{1}{2}] x[1] \\boldsymbol{e}_3 \\textrm{, if } v^3[\\tfrac{1}{2}] < 0 - \\end{cases} - ``` - [`Extrapolate()`](@ref): set the value of `x` to be the same as the closest interior point. On the left boundary, the stencil is ```math @@ -1280,7 +1197,7 @@ struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator assert_valid_bcs( "UpwindBiasedProductC2F", kwargs, - (SetValue, Extrapolate), + (Extrapolate,), ) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end @@ -1320,44 +1237,6 @@ end boundary_width(::UpwindBiasedProductC2F, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::UpwindBiasedProductC2F, - bc::SetValue, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_face_boundary_idx(space) - aᴸᴮ = getidx(space, bc.val, nothing, hidx) - a⁺ = stencil_interior(RightBiasedC2F(), space, idx, hidx, arg) - vᶠ = Geometry.contravariant3( - getidx(space, velocity, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, aᴸᴮ, a⁺)) -end - -Base.@propagate_inbounds function stencil_right_boundary( - ::UpwindBiasedProductC2F, - bc::SetValue, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_face_boundary_idx(space) - a⁻ = stencil_interior(LeftBiasedC2F(), space, idx, hidx, arg) - aᴿᴮ = getidx(space, bc.val, nothing, hidx) - vᶠ = Geometry.contravariant3( - getidx(space, velocity, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, a⁻, aᴿᴮ)) -end - Base.@propagate_inbounds function stencil_left_boundary( op::UpwindBiasedProductC2F, ::Extrapolate, @@ -2236,409 +2115,6 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) end -""" - A = AdvectionF2F(;boundaries) - A.(v, θ) - -Vertical advection operator at cell faces, for a face-valued velocity field `v` and face-valued -variables `θ`, approximating ``v^3 \\partial_3 \\theta``. - -It uses the following stencil -```math -A(v,θ)[i] = \\frac{1}{2} (θ[i+1] - θ[i-1]) v³[i] -``` - -No boundary conditions are currently supported. -""" -struct AdvectionF2F{BCS <: @NamedTuple{}} <: AdvectionOperator - bcs::BCS -end - -function AdvectionF2F(; kwargs...) - assert_no_bcs("AdvectionF2F", kwargs) - AdvectionF2F(NamedTuple(kwargs)) -end - -return_space( - ::AdvectionF2F, - velocity_space::AllFaceFiniteDifferenceSpace, - arg_space::AllFaceFiniteDifferenceSpace, -) = arg_space - -stencil_interior_width(::AdvectionF2F, velocity, arg) = ((0, 0), (-1, 1)) -Base.@propagate_inbounds function stencil_interior( - ::AdvectionF2F, - space, - idx, - hidx, - velocity, - arg, -) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³ = Geometry.contravariant3( - getidx(space, velocity, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - ∂θ₃ = (θ⁺ - θ⁻) / 2 - return w³ * ∂θ₃ -end -boundary_width(::AdvectionF2F, ::AbstractBoundaryCondition) = 1 - -""" - A = AdvectionC2C(;boundaries) - A.(v, θ) - -Vertical advection operator at cell centers, for cell face velocity field `v` -cell center variables `θ`, approximating ``v^3 \\partial_3 \\theta``. - -It uses the following stencil -```math -A(v,θ)[i] = \\frac{1}{2} \\{ (θ[i+1] - θ[i]) v³[i+\\tfrac{1}{2}] + (θ[i] - θ[i-1])v³[i-\\tfrac{1}{2}]\\} -``` - -Supported boundary conditions: - -- [`SetValue(θ₀)`](@ref): set the value of `θ` at the boundary face to be `θ₀`. - At the lower boundary, this is: -```math -A(v,θ)[1] = \\frac{1}{2} \\{ (θ[2] - θ[1]) v³[1 + \\tfrac{1}{2}] + (θ[1] - θ₀)v³[\\tfrac{1}{2}]\\} -``` -- [`Extrapolate`](@ref): use the closest interior point as the boundary value. - At the lower boundary, this is: -```math -A(v,θ)[1] = (θ[2] - θ[1]) v³[1 + \\tfrac{1}{2}] \\} -``` -""" -struct AdvectionC2C{BCS} <: AdvectionOperator - bcs::BCS - function AdvectionC2C(; kwargs...) - assert_valid_bcs("AdvectionC2C", kwargs, (SetValue, Extrapolate)) - new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) - end - AdvectionC2C(bcs) = AdvectionC2C(; bcs...) -end - -return_space( - ::AdvectionC2C, - velocity_space::AllFaceFiniteDifferenceSpace, - arg_space::AllCenterFiniteDifferenceSpace, -) = arg_space - -stencil_interior_width(::AdvectionC2C, velocity, arg) = - ((-half, +half), (-1, 1)) -Base.@propagate_inbounds function stencil_interior( - ::AdvectionC2C, - space, - idx, - hidx, - velocity, - arg, -) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - ∂θ₃⁻ = θ - θ⁻ - return (w³⁺ * ∂θ₃⁺ + w³⁻ * ∂θ₃⁻) / 2 -end - -boundary_width(::AdvectionC2C, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::AdvectionC2C, - bc::SetValue, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_center_boundary_idx(space) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, bc.val, nothing, hidx) # defined at face, not the center - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - ∂θ₃⁻ = 2 * (θ - θ⁻) - return (w³⁺ * ∂θ₃⁺ + w³⁻ * ∂θ₃⁻) / 2 -end -Base.@propagate_inbounds function stencil_right_boundary( - ::AdvectionC2C, - bc::SetValue, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_center_boundary_idx(space) - θ⁺ = getidx(space, bc.val, nothing, hidx) # value at the face - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁺ = 2 * (θ⁺ - θ) - ∂θ₃⁻ = θ - θ⁻ - return (w³⁺ * ∂θ₃⁺ + w³⁻ * ∂θ₃⁻) / 2 -end - -Base.@propagate_inbounds function stencil_left_boundary( - ::AdvectionC2C, - ::Extrapolate, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_center_boundary_idx(space) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - return (w³⁺ * ∂θ₃⁺) -end -Base.@propagate_inbounds function stencil_right_boundary( - ::AdvectionC2C, - ::Extrapolate, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_center_boundary_idx(space) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁻ = θ - θ⁻ - return (w³⁻ * ∂θ₃⁻) -end - -""" - A = FluxCorrectionC2C(;boundaries) - A.(v, θ) - -Vertical advection operator at cell centers, for cell center velocity field `v` -cell center variables `θ`, approximating ``v^3 \\partial_3 \\theta``. - -It uses the following stencil (TODO) - -```math -``` - -Supported boundary conditions: - -- [`Extrapolate`](@ref): use the closest interior point as the boundary value. - At the lower boundary. -""" -struct FluxCorrectionC2C{BCS} <: AdvectionOperator - bcs::BCS - function FluxCorrectionC2C(; kwargs...) - assert_valid_bcs("FluxCorrectionC2C", kwargs, (Extrapolate,)) - new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) - end - FluxCorrectionC2C(bcs) = FluxCorrectionC2C(; bcs...) -end - -return_space( - ::FluxCorrectionC2C, - velocity_space::AllFaceFiniteDifferenceSpace, - arg_space::AllCenterFiniteDifferenceSpace, -) = arg_space - -stencil_interior_width(::FluxCorrectionC2C, velocity, arg) = - ((-half, +half), (-1, 1)) -Base.@propagate_inbounds function stencil_interior( - ::FluxCorrectionC2C, - space, - idx, - hidx, - velocity, - arg, -) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - ∂θ₃⁻ = θ - θ⁻ - return abs(w³⁺) * ∂θ₃⁺ - abs(w³⁻) * ∂θ₃⁻ -end - -boundary_width(::FluxCorrectionC2C, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::FluxCorrectionC2C, - ::Extrapolate, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_center_boundary_idx(space) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - return abs(w³⁺) * ∂θ₃⁺ -end -Base.@propagate_inbounds function stencil_right_boundary( - ::FluxCorrectionC2C, - ::Extrapolate, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_center_boundary_idx(space) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁻ = θ - θ⁻ - return -abs(w³⁻) * ∂θ₃⁻ -end - -""" - A = FluxCorrectionF2F(;boundaries) - A.(v, θ) - -Vertical advection operator at cell faces, for cell face velocity field `v` -cell face variables `θ`, approximating ``v^3 \\partial_3 \\theta``. - -It uses the following stencil (TODO) - -```math -``` - -Supported boundary conditions: - -- [`Extrapolate`](@ref): use the closest interior point as the boundary value. - At the lower boundary. -""" -struct FluxCorrectionF2F{BCS} <: AdvectionOperator - bcs::BCS - function FluxCorrectionF2F(; kwargs...) - assert_valid_bcs("FluxCorrectionF2F", kwargs, (Extrapolate,)) - new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) - end - FluxCorrectionF2F(bcs) = FluxCorrectionF2F(; bcs...) -end - -return_space( - ::FluxCorrectionF2F, - velocity_space::AllCenterFiniteDifferenceSpace, - arg_space::AllFaceFiniteDifferenceSpace, -) = arg_space - -stencil_interior_width(::FluxCorrectionF2F, velocity, arg) = - ((-half, +half), (-1, 1)) -Base.@propagate_inbounds function stencil_interior( - ::FluxCorrectionF2F, - space, - idx, - hidx, - velocity, - arg, -) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - ∂θ₃⁻ = θ - θ⁻ - return abs(w³⁺) * ∂θ₃⁺ - abs(w³⁻) * ∂θ₃⁻ -end - -boundary_width(::FluxCorrectionF2F, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::FluxCorrectionF2F, - ::Extrapolate, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_face_boundary_idx(space) - θ⁺ = getidx(space, arg, idx + 1, hidx) - θ = getidx(space, arg, idx, hidx) - w³⁺ = Geometry.contravariant3( - getidx(space, velocity, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - ∂θ₃⁺ = θ⁺ - θ - return abs(w³⁺) * ∂θ₃⁺ -end -Base.@propagate_inbounds function stencil_right_boundary( - ::FluxCorrectionF2F, - ::Extrapolate, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_face_boundary_idx(space) - θ = getidx(space, arg, idx, hidx) - θ⁻ = getidx(space, arg, idx - 1, hidx) - w³⁻ = Geometry.contravariant3( - getidx(space, velocity, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - ∂θ₃⁻ = θ - θ⁻ - return -abs(w³⁻) * ∂θ₃⁻ -end - - abstract type BoundaryOperator <: FiniteDifferenceOperator end """ @@ -2829,11 +2305,6 @@ G(x)[i]^3 = x[i+\\tfrac{1}{2}] - x[i-\\tfrac{1}{2}] ``` The following boundary conditions are supported: -- [`SetValue(x₀)`](@ref): calculate the gradient assuming the value at the - boundary is `x₀`. For the left boundary, this becomes: - ```math - G(x)[\\tfrac{1}{2}]³ = 2 (x[1] - x₀) - ``` - [`SetGradient(v₀)`](@ref): set the value of the gradient at the boundary to be `v₀`. For the left boundary, this becomes: ```math @@ -2843,7 +2314,7 @@ The following boundary conditions are supported: struct GradientC2F{BC} <: GradientOperator bcs::BC function GradientC2F(; kwargs...) - assert_valid_bcs("GradientC2F", kwargs, (SetValue, SetGradient)) + assert_valid_bcs("GradientC2F", kwargs, (SetGradient,)) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end GradientC2F(bcs) = GradientC2F(; bcs...) @@ -2867,36 +2338,6 @@ Base.@propagate_inbounds function stencil_interior( end boundary_width(::GradientC2F, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::GradientC2F, - bc::SetValue, - space, - idx, - hidx, - arg, -) - @assert idx == left_face_boundary_idx(space) - # ∂x[i] = 2(∂x[i + half] - val) - Geometry.Covariant3Vector(2) ⊗ ( - getidx(space, arg, idx + half, hidx) - - getidx(space, bc.val, nothing, hidx) - ) -end -Base.@propagate_inbounds function stencil_right_boundary( - ::GradientC2F, - bc::SetValue, - space, - idx, - hidx, - arg, -) - @assert idx == right_face_boundary_idx(space) - Geometry.Covariant3Vector(2) ⊗ ( - getidx(space, bc.val, nothing, hidx) - - getidx(space, arg, idx - half, hidx) - ) -end - # left / right SetGradient boundary conditions Base.@propagate_inbounds function stencil_left_boundary( @@ -3230,11 +2671,6 @@ where `Jv³` is the Jacobian multiplied by the third contravariant component of `v`. The following boundary conditions are supported: -- [`SetValue(v₀)`](@ref): calculate the divergence assuming the value at the - boundary is `v₀`. For the left boundary, this becomes: - ```math - D(v)[\\tfrac{1}{2}] = \\frac{1}{2} (Jv³[1] - Jv³₀) / J[i] - ``` - [`SetDivergence(x)`](@ref): set the value of the divergence at the boundary to be `x`. ```math D(v)[\\tfrac{1}{2}] = x @@ -3243,7 +2679,7 @@ The following boundary conditions are supported: struct DivergenceC2F{BC} <: DivergenceOperator bcs::BC function DivergenceC2F(; kwargs...) - assert_valid_bcs("DivergenceC2F", kwargs, (SetValue, SetDivergence)) + assert_valid_bcs("DivergenceC2F", kwargs, (SetDivergence,)) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end DivergenceC2F(bcs) = DivergenceC2F(; bcs...) @@ -3273,48 +2709,6 @@ Base.@propagate_inbounds function stencil_interior( end boundary_width(::DivergenceC2F, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::DivergenceC2F, - bc::SetValue, - space, - idx, - hidx, - arg, -) - @assert idx == left_face_boundary_idx(space) - # ∂x[i] = 2(∂x[i + half] - val) - local_geometry = Geometry.LocalGeometry(space, idx, hidx) - Ju³₊ = Geometry.Jcontravariant3( - getidx(space, arg, idx + half, hidx), - Geometry.LocalGeometry(space, idx + half, hidx), - ) - Ju³ = Geometry.Jcontravariant3( - getidx(space, bc.val, nothing, hidx), - local_geometry, - ) - (Ju³₊ - Ju³) * (2 * local_geometry.invJ) -end -Base.@propagate_inbounds function stencil_right_boundary( - ::DivergenceC2F, - bc::SetValue, - space, - idx, - hidx, - arg, -) - @assert idx == right_face_boundary_idx(space) - local_geometry = Geometry.LocalGeometry(space, idx, hidx) - Ju³ = Geometry.Jcontravariant3( - getidx(space, bc.val, nothing, hidx), - local_geometry, - ) - Ju³₋ = Geometry.Jcontravariant3( - getidx(space, arg, idx - half, hidx), - Geometry.LocalGeometry(space, idx - half, hidx), - ) - (Ju³ - Ju³₋) * (2 * local_geometry.invJ) -end - # left / right SetDivergence boundary conditions Base.@propagate_inbounds function stencil_left_boundary( ::DivergenceC2F, @@ -3374,19 +2768,13 @@ where ``v₁`` and ``v₂`` are the 1st and 2nd covariant components of ``v``, a The following boundary conditions are supported: -- [`SetValue(v₀)`](@ref): calculate the curl assuming the value of ``v`` at the - boundary is `v₀`. For the left boundary, this becomes: - ```math - C(v)[\\tfrac{1}{2}]^1 = -\\frac{2}{J[i]} (v_2[1] - (v₀)_2) - C(v)[\\tfrac{1}{2}]^2 = \\frac{2}{J[i]} (v_1[1] - (v₀)_1) - ``` - [`SetCurl(v⁰)`](@ref): enforce the curl operator output at the boundary to be the contravariant vector `v⁰`. """ struct CurlC2F{BC} <: CurlFiniteDifferenceOperator bcs::BC function CurlC2F(; kwargs...) - assert_valid_bcs("CurlC2F", kwargs, (SetValue, SetCurl)) + assert_valid_bcs("CurlC2F", kwargs, (SetCurl,)) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end CurlC2F(bcs) = CurlC2F(; bcs...) @@ -3418,32 +2806,6 @@ Base.@propagate_inbounds function stencil_interior( end boundary_width(::CurlC2F, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::CurlC2F, - bc::SetValue, - space, - idx, - hidx, - arg, -) - u₊ = getidx(space, arg, idx + half, hidx) - u = getidx(space, bc.val, nothing, hidx) - local_geometry = Geometry.LocalGeometry(space, idx, hidx) - return fd3_curl(u₊, u, local_geometry.invJ * 2) -end -Base.@propagate_inbounds function stencil_right_boundary( - ::CurlC2F, - bc::SetValue, - space, - idx, - hidx, - arg, -) - u = getidx(space, bc.val, nothing, hidx) - u₋ = getidx(space, arg, idx - half, hidx) - local_geometry = Geometry.LocalGeometry(space, idx, hidx) - return fd3_curl(u, u₋, local_geometry.invJ * 2) -end # Project the user-supplied curl value onto the full Contravariant123 axis so # the boundary output matches `return_eltype` (uniformly Contravariant123Vector). diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index cb6deadee9..bd426ce662 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -27,10 +27,6 @@ import ClimaCore.Operators: Upwind3rdOrderBiasedProductC2F, FCTBorisBook, FCTZalesak, - AdvectionC2C, - AdvectionF2F, - FluxCorrectionC2C, - FluxCorrectionF2F, SetBoundaryOperator, GradientC2F, GradientF2C, @@ -126,7 +122,6 @@ end # Note: The Curl operator currently only works with C12, C1, or C2 inputs. test_op_matrix(InterpolateC2F, Nothing, (ᶜnested,), true) test_op_matrix(InterpolateC2F, SetValue, (ᶜnested,)) - test_op_matrix(InterpolateC2F, SetGradient, (ᶜnested,)) test_op_matrix(InterpolateC2F, Extrapolate, (ᶜnested,)) test_op_matrix(InterpolateF2C, Nothing, (ᶠnested,)) test_op_matrix(LeftBiasedC2F, Nothing, (ᶜnested,), true) @@ -139,11 +134,9 @@ end test_op_matrix(RightBiasedF2C, SetValue, (ᶠnested,)) test_op_matrix(WeightedInterpolateC2F, Nothing, (ᶜscalar, ᶜnested), true) test_op_matrix(WeightedInterpolateC2F, SetValue, (ᶜscalar, ᶜnested)) - test_op_matrix(WeightedInterpolateC2F, SetGradient, (ᶜscalar, ᶜnested)) test_op_matrix(WeightedInterpolateC2F, Extrapolate, (ᶜscalar, ᶜnested)) test_op_matrix(WeightedInterpolateF2C, Nothing, (ᶠscalar, ᶠnested)) test_op_matrix(UpwindBiasedProductC2F, Nothing, (ᶠuvw, ᶜscalar), true) - test_op_matrix(UpwindBiasedProductC2F, SetValue, (ᶠuvw, ᶜscalar)) test_op_matrix(UpwindBiasedProductC2F, Extrapolate, (ᶠuvw, ᶜscalar)) test_op_matrix( Upwind3rdOrderBiasedProductC2F, @@ -157,27 +150,19 @@ end (ᶠuvw, ᶜscalar), true, ) - test_op_matrix(AdvectionC2C, SetValue, (ᶠuvw, ᶜnested)) - test_op_matrix(AdvectionC2C, Extrapolate, (ᶠuvw, ᶜnested)) - test_op_matrix(AdvectionF2F, Nothing, (ᶠuvw, ᶠnested), true) - test_op_matrix(FluxCorrectionC2C, Extrapolate, (ᶠuvw, ᶜnested)) - test_op_matrix(FluxCorrectionF2F, Extrapolate, (ᶜuvw, ᶠnested)) test_op_matrix(SetBoundaryOperator, SetValue, (ᶠnested,)) test_op_matrix(GradientC2F, Nothing, (ᶜscalar,), true) - test_op_matrix(GradientC2F, SetValue, (ᶜscalar,)) test_op_matrix(GradientC2F, SetGradient, (ᶜscalar,)) test_op_matrix(GradientF2C, Nothing, (ᶠscalar,)) test_op_matrix(GradientF2C, SetValue, (ᶠscalar,)) test_op_matrix(GradientF2C, Extrapolate, (ᶠscalar,)) test_op_matrix(DivergenceC2F, Nothing, (ᶜuvw,), true) - test_op_matrix(DivergenceC2F, SetValue, (ᶜuvw,)) test_op_matrix(DivergenceC2F, SetDivergence, (ᶜuvw,)) test_op_matrix(DivergenceF2C, Nothing, (ᶠuvw,)) test_op_matrix(DivergenceF2C, SetValue, (ᶠuvw,)) test_op_matrix(DivergenceF2C, SetDivergence, (ᶠuvw,)) test_op_matrix(DivergenceF2C, Extrapolate, (ᶠuvw,)) test_op_matrix(CurlC2F, Nothing, (ᶜc12,), true) - test_op_matrix(CurlC2F, SetValue, (ᶜc12,)) test_op_matrix(CurlC2F, SetCurl, (ᶜc12,)) @test_throws "nonlinear" MatrixFields.operator_matrix(FCTBorisBook()) @@ -202,29 +187,18 @@ end set_nested_values = (; bottom = SetValue(nested_zero), top = SetValue(nested_zero)) c12_zero = zero(Geometry.Covariant12Vector{FT}) - set_c12_values = (; bottom = SetValue(c12_zero), top = SetValue(c12_zero)) extrapolate = (; bottom = Extrapolate(), top = Extrapolate()) ᶠinterp = InterpolateC2F(; set_nested_values...) ᶜlbias = LeftBiasedF2C() ᶠrbias = RightBiasedC2F(; set_nested_values.top) ᶜwinterp = WeightedInterpolateF2C() - ᶠupwind = UpwindBiasedProductC2F(; set_scalar_values...) - ᶜadvect = AdvectionC2C(; extrapolate...) - ᶜflux_correct = FluxCorrectionC2C(; extrapolate...) - ᶠgrad = GradientC2F(; set_scalar_values...) ᶜdiv = DivergenceF2C() - ᶠcurl = CurlC2F(; set_c12_values...) ᶠinterp_matrix = MatrixFields.operator_matrix(ᶠinterp) ᶜlbias_matrix = MatrixFields.operator_matrix(ᶜlbias) ᶠrbias_matrix = MatrixFields.operator_matrix(ᶠrbias) ᶜwinterp_matrix = MatrixFields.operator_matrix(ᶜwinterp) - ᶠupwind_matrix = MatrixFields.operator_matrix(ᶠupwind) - ᶜadvect_matrix = MatrixFields.operator_matrix(ᶜadvect) - ᶜflux_correct_matrix = MatrixFields.operator_matrix(ᶜflux_correct) - ᶠgrad_matrix = MatrixFields.operator_matrix(ᶠgrad) ᶜdiv_matrix = MatrixFields.operator_matrix(ᶜdiv) - ᶠcurl_matrix = MatrixFields.operator_matrix(ᶠcurl) @test_throws "does not contain any Fields" @. ᶜlbias_matrix() * ᶠinterp_matrix() @@ -249,17 +223,13 @@ end test_field_broadcast(; test_name = "product of six operator matrices", get_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() ), set_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() @@ -270,31 +240,21 @@ end test_name = "applying six operators to a nested field using operator \ matrices", get_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() * ᶜnested ), set_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() * ᶜnested ), ref_set_result = @lazy( - @. ᶜflux_correct( - ᶠuvw, - ᶜadvect( - ᶠuvw, - ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))), - ), - ) + @. ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))) ), ) # this test is will fail because of incorrect results, not InvalidIRError @@ -302,151 +262,30 @@ end test_name = "applying six operators to a nested field using operator \ matrices, but with forced right associativity", get_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * ( - ᶜadvect_matrix(ᶠuvw) * ( - ᶜwinterp_matrix(ᶠscalar) * ( - ᶠrbias_matrix() * - (ᶜlbias_matrix() * (ᶠinterp_matrix() * ᶜnested)) - ) + @. ( + ( + ᶜwinterp_matrix(ᶠscalar) * ( + ᶠrbias_matrix() * + (ᶜlbias_matrix() * (ᶠinterp_matrix() * ᶜnested)) ) ) + ) ), set_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * ( - ᶜadvect_matrix(ᶠuvw) * ( - ᶜwinterp_matrix(ᶠscalar) * ( - ᶠrbias_matrix() * - (ᶜlbias_matrix() * (ᶠinterp_matrix() * ᶜnested)) - ) + @. ( + ( + ᶜwinterp_matrix(ᶠscalar) * ( + ᶠrbias_matrix() * + (ᶜlbias_matrix() * (ᶠinterp_matrix() * ᶜnested)) ) ) + ) ), ref_set_result = @lazy( - @. ᶜflux_correct( - ᶠuvw, - ᶜadvect( - ᶠuvw, - ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))), - ), - ) + @. ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))) ), time_ratio_limit = 30, # This case's ref function is fast on Buildkite. test_broken_with_cuda = true, # TODO: Fix this. ) - # TODO: For some reason, we need to compile and run @test_opt on several - # simpler broadcast expressions before we can run the remaining two test - # cases. As of Julia 1.8.5, the tests fail if we skip this step. Is this a - # false positive, a compiler issue, or a sign that the code can be improved? - for get_result in ( - @lazy( - @. (c12_b',) * - ᶜwinterp_matrix(ᶠscalar) * - ᶠcurl_matrix() * - (c12_a,) + - (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5 - ), - @lazy( - @. ᶜdiv_matrix() * - DiagonalMatrixRow(ᶠscalar) * - ᶠgrad_matrix() * - ( - (c12_b',) * - ᶜwinterp_matrix(ᶠscalar) * - ᶠcurl_matrix() * - (c12_a,) + - (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5 - ) - ), - ) - materialize(get_result) - @test_opt ignored_modules = CUDA_FRAMES materialize(get_result) - end - - test_field_broadcast(; - test_name = "non-trivial combination of operator matrices and other \ - matrix fields", - get_result = @lazy( - @. ᶠupwind_matrix(ᶠuvw) * ( - ᶜdiv_matrix() * - DiagonalMatrixRow(ᶠscalar) * - ᶠgrad_matrix() * - ( - (c12_b',) * - ᶜwinterp_matrix(ᶠscalar) * - ᶠcurl_matrix() * - (c12_a,) + - (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5 - ) - (2I,) - ) - ), - set_result = @lazy( - @. ᶠupwind_matrix(ᶠuvw) * ( - ᶜdiv_matrix() * - DiagonalMatrixRow(ᶠscalar) * - ᶠgrad_matrix() * - ( - (c12_b',) * - ᶜwinterp_matrix(ᶠscalar) * - ᶠcurl_matrix() * - (c12_a,) + - (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / 5 - ) - (2I,) - ) - ), - ) - - # TODO: This case's reference function takes too long to compile on both - # CPUs and GPUs (more than half an hour), as of Julia 1.9. This might be - # happening because of excessive inlining---aside from *, all other finite - # difference operators use @propagate_inbounds. So, the reference function - # is currently disabled, although the test does pass when it is enabled. - test_field_broadcast(; - test_name = "applying a non-trivial sequence of operations to a scalar \ - field using operator matrices and other matrix fields", - get_result = @lazy( - @. ᶠupwind_matrix(ᶠuvw) * - ( - ᶜdiv_matrix() * - DiagonalMatrixRow(ᶠscalar) * - ᶠgrad_matrix() * - ( - (c12_b',) * - ᶜwinterp_matrix(ᶠscalar) * - ᶠcurl_matrix() * - (c12_a,) + - (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / - 5 - ) - (2I,) - ) * - ᶜscalar - ), - set_result = @lazy( - @. ᶠupwind_matrix(ᶠuvw) * - ( - ᶜdiv_matrix() * - DiagonalMatrixRow(ᶠscalar) * - ᶠgrad_matrix() * - ( - (c12_b',) * - ᶜwinterp_matrix(ᶠscalar) * - ᶠcurl_matrix() * - (c12_a,) + - (DiagonalMatrixRow(ᶜdiv(ᶠuvw)) - ᶜadvect_matrix(ᶠuvw)) / - 5 - ) - (2I,) - ) * - ᶜscalar - ), - # ref_set_result = @lazy(@. ᶠupwind( - # ᶠuvw, - # ᶜdiv( - # ᶠscalar * ᶠgrad( - # (c12_b',) * ᶜwinterp(ᶠscalar, ᶠcurl((c12_a,) * ᶜscalar)) + - # (ᶜdiv(ᶠuvw) * ᶜscalar - ᶜadvect(ᶠuvw, ᶜscalar)) / 5, - # ), - # ) - 2 * ᶜscalar, - # )), - # max_eps_error_limit = 20, # This case's roundoff error is large. - ) end diff --git a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl index 78668a8549..ff74e0676e 100644 --- a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl +++ b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl @@ -82,13 +82,6 @@ function op_Upwind3rdOrderBiasedProductC2F!(c, f, bcs = (;)) end n_reads_writes(::Type{typeof(op_Upwind3rdOrderBiasedProductC2F!)}) = -1 # todo #### Simple composed (non-exhaustive due to combinatorial explosion) -function op_divgrad_CC!(c, f, bcs) - grad = Operators.GradientC2F(bcs.inner) - div = Operators.DivergenceF2C(bcs.outer) - @. c.y = div(grad(c.x)) - return nothing -end -n_reads_writes(::Type{typeof(op_divgrad_CC!)}) = 3 # 1 write, 2 reads (1 metric term) function op_divgrad_FF!(c, f, bcs) grad = Operators.GradientF2C(bcs.inner) div = Operators.DivergenceC2F(bcs.outer) @@ -194,9 +187,5 @@ LeftBiased3rdOrderC2F # unused LeftBiased3rdOrderF2C # unused RightBiased3rdOrderC2F # unused RightBiased3rdOrderF2C # unused -AdvectionF2F # only used in ClimaAtmos src/ -AdvectionC2C # only used in ClimaAtmos src/ -FluxCorrectionC2C # unused at the moment -FluxCorrectionF2F # only used in ClimaAtmos src/ ``` =# diff --git a/test/Operators/finitedifference/benchmark_stencils_utils.jl b/test/Operators/finitedifference/benchmark_stencils_utils.jl index fc91914d59..7a925490dd 100755 --- a/test/Operators/finitedifference/benchmark_stencils_utils.jl +++ b/test/Operators/finitedifference/benchmark_stencils_utils.jl @@ -164,13 +164,6 @@ function set_curl_bcs(c) top = Operators.SetCurl(contra12(FT(0), FT(0)))) end -function set_curl_value_bcs(c) - FT = Spaces.undertype(axes(c)) - cov12 = Geometry.Covariant12Vector - return (;bottom = Operators.SetValue(cov12(FT(0), FT(0))), - top = Operators.SetValue(cov12(FT(0), FT(0)))) -end - function bc_name(bcs::NamedTuple) if haskey(bcs, :inner) && haskey(bcs, :outer) return (bc_name_base(bcs.inner)..., bc_name_base(bcs.outer)...) @@ -200,7 +193,7 @@ bcs_tested(c, ::typeof(op_broadcast_example1!)) = ((;), ) bcs_tested(c, ::typeof(op_broadcast_example2!)) = ((;), ) bcs_tested(c, ::typeof(op_GradientF2C!)) = ((;), set_value_bcs(c)) -bcs_tested(c, ::typeof(op_GradientC2F!)) = (set_gradient_value_bcs(c), set_value_bcs(c)) +bcs_tested(c, ::typeof(op_GradientC2F!)) = (set_gradient_value_bcs(c),) bcs_tested(c, ::typeof(op_DivergenceF2C!)) = ((;), extrapolate_bcs(c)) bcs_tested(c, ::typeof(op_DivergenceC2F!)) = (set_divergence_bcs(c), ) bcs_tested(c, ::typeof(op_InterpolateF2C!)) = ((;), ) @@ -209,21 +202,19 @@ bcs_tested(c, ::typeof(op_LeftBiasedC2F!)) = (set_bot_value_bc(c),) bcs_tested(c, ::typeof(op_LeftBiasedF2C!)) = ((;), set_bot_value_bc(c)) bcs_tested(c, ::typeof(op_RightBiasedC2F!)) = (set_top_value_bc(c),) bcs_tested(c, ::typeof(op_RightBiasedF2C!)) = ((;), set_top_value_bc(c)) -bcs_tested(c, ::typeof(op_CurlC2F!)) = (set_curl_bcs(c), set_curl_value_bcs(c)) -bcs_tested(c, ::typeof(op_UpwindBiasedProductC2F!)) = (set_value_bcs(c), extrapolate_bcs(c)) +bcs_tested(c, ::typeof(op_CurlC2F!)) = (set_curl_bcs(c),) +bcs_tested(c, ::typeof(op_UpwindBiasedProductC2F!)) = (extrapolate_bcs(c),) bcs_tested(c, ::typeof(op_Upwind3rdOrderBiasedProductC2F!)) = (set_upwind_biased_3_bcs(c), extrapolate_bcs(c)) # Composed operators (bcs handled case-by-case) bcs_tested(c, ::typeof(op_divUpwind3rdOrderBiasedProductC2F!)) = ((; inner = set_upwind_biased_3_bcs(c), outer = set_value_contra3_bcs(c)), ) -bcs_tested(c, ::typeof(op_divgrad_CC!)) = - ((; inner = set_value_bcs(c), outer = (;)), ) bcs_tested(c, ::typeof(op_divgrad_FF!)) = ((; inner = (;), outer = set_divergence_bcs(c)), ) bcs_tested(c, ::typeof(op_div_interp_CC!)) = ((; inner = set_value_contra3_bcs(c), outer = (;)), ) bcs_tested(c, ::typeof(op_div_interp_FF!)) = - ((; inner = (;), outer = set_value_contra3_bcs(c)), ) + ((; inner = (;), outer = set_divergence_bcs(c)), ) bcs_tested(c, ::typeof(op_divgrad_uₕ!)) = ( (; inner = (;), outer = set_value_divgrad_uₕ_bcs(c)), @@ -406,7 +397,6 @@ function benchmark_operators_base(bm, trials, t_min, cfield, ffield, name; compi # op_Upwind3rdOrderBiasedProductC2F!, # TODO: do we need to test this for different w values? #### Composed op_divUpwind3rdOrderBiasedProductC2F!, - op_divgrad_CC!, op_divgrad_FF!, op_div_interp_CC!, op_div_interp_FF!, @@ -435,7 +425,6 @@ function test_results_column(t_min) [(op_GradientF2C!, :none), 253.100*ns*buffer], [(op_GradientF2C!, :SetValue, :SetValue), 270.448*ns*buffer], [(op_GradientC2F!, :SetGradient, :SetGradient), 242.053*ns*buffer], - [(op_GradientC2F!, :SetValue, :SetValue), 241.647*ns*buffer], [(op_DivergenceF2C!, :none), 1.005*μs*buffer], [(op_DivergenceF2C!, :Extrapolate, :Extrapolate), 1.076*μs*buffer], [(op_DivergenceC2F!, :SetDivergence, :SetDivergence), 878.028*ns*buffer], @@ -450,14 +439,11 @@ function test_results_column(t_min) [(op_RightBiasedF2C!, :none), 142.120*ns*buffer], [(op_RightBiasedF2C!, :SetValue), 141.446*ns*buffer], [(op_CurlC2F!, :SetCurl, :SetCurl), 1.692*μs*buffer], - [(op_CurlC2F!, :SetValue, :SetValue), 1.616*μs*buffer], - [(op_UpwindBiasedProductC2F!, :SetValue, :SetValue), 754.856*ns*buffer], [(op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate), 765.401*ns*buffer], [(op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue), 2.540*μs*buffer], - [(op_divgrad_CC!, :SetValue, :SetValue, :none), 924.147*ns*buffer], [(op_divgrad_FF!, :none, :SetDivergence, :SetDivergence), 876.510*ns*buffer], [(op_div_interp_CC!, :SetValue, :SetValue, :none), 721.119*ns*buffer], - [(op_div_interp_FF!, :none, :SetValue, :SetValue), 686.581*ns*buffer], + [(op_div_interp_FF!, :none, :SetDivergence, :SetDivergence), 686.581*ns*buffer], [(op_divgrad_uₕ!, :none, :SetValue, :Extrapolate), 4.960*μs*buffer], [(op_divgrad_uₕ!, :none, :SetValue, :SetValue), 5.047*μs*buffer], ] @@ -479,7 +465,6 @@ function test_results_sphere(t_min) [(op_GradientF2C!, :none), 1.746*ms*buffer], [(op_GradientF2C!, :SetValue, :SetValue), 1.754*ms*buffer], [(op_GradientC2F!, :SetGradient, :SetGradient), 1.899*ms*buffer], - [(op_GradientC2F!, :SetValue, :SetValue), 1.782*ms*buffer], [(op_DivergenceF2C!, :none), 6.792*ms*buffer], [(op_DivergenceF2C!, :Extrapolate, :Extrapolate), 6.776*ms*buffer], [(op_DivergenceC2F!, :SetDivergence, :SetDivergence), 6.720*ms*buffer], @@ -496,14 +481,11 @@ function test_results_sphere(t_min) [(op_RightBiasedF2C!, :none), 1.582*ms*buffer], [(op_RightBiasedF2C!, :SetValue), 1.551*ms*buffer], [(op_CurlC2F!, :SetCurl, :SetCurl), 4.669*ms*buffer], - [(op_CurlC2F!, :SetValue, :SetValue), 4.568*ms*buffer], - [(op_UpwindBiasedProductC2F!, :SetValue, :SetValue), 3.444*ms*buffer], [(op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate), 3.432*ms*buffer], [(op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue), 5.650*ms*buffer], - [(op_divgrad_CC!, :SetValue, :SetValue, :none), 4.474*ms*buffer], [(op_divgrad_FF!, :none, :SetDivergence, :SetDivergence), 4.470*ms*buffer], [(op_div_interp_CC!, :SetValue, :SetValue, :none), 3.566*ms*buffer], - [(op_div_interp_FF!, :none, :SetValue, :SetValue), 3.663*ms*buffer], + [(op_div_interp_FF!, :none, :SetDivergence, :SetDivergence), 3.663*ms*buffer], [(op_divgrad_uₕ!, :none, :SetValue, :Extrapolate), 7.470*ms*buffer], [(op_divgrad_uₕ!, :none, :SetValue, :SetValue), 7.251*ms*buffer], ] diff --git a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl index 1ca603b8d2..510e1418e8 100644 --- a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl +++ b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl @@ -61,21 +61,26 @@ convergence_rate(err, Δh) = V = Geometry.WVector.(ones(FT, fs)) function ∑tendencies!(dT, T, z, t) - bc_vb = Operators.SetValue( - FT(gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ)), + bc_gb = Operators.SetGradient( + Geometry.WVector(FT(∇gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))), ) bc_gt = Operators.SetGradient( - Geometry.WVector( - FT(∇gaussian(z₁, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ)), - ), + Geometry.WVector(FT(∇gaussian(z₁, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))), ) - A = Operators.AdvectionC2C( - bottom = bc_vb, - top = Operators.Extrapolate(), - ) - gradc2f = Operators.GradientC2F(; bottom = bc_vb, top = bc_gt) + top_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(T, Fields.nlevels(T)) .- + Fields.level(T, Fields.nlevels(T) - 1) + ) + + bc_gt_lb = Operators.SetGradient(top_center_left_biased_grad) + gradc2f = Operators.GradientC2F(bottom = bc_gb, top = bc_gt) + gradc2f_advect = Operators.GradientC2F(bottom = bc_gb, top = bc_gt_lb) + interpf2c = Operators.InterpolateF2C() divf2c = Operators.DivergenceF2C() - return @. dT = divf2c(ν * gradc2f(T)) - A(V, T) + return @. dT = + divf2c(ν * gradc2f(T)) - + interpf2c(Geometry.dot(Geometry.Contravariant3Vector(V), gradc2f_advect(T))) end # Solve the ODE operator diff --git a/test/Operators/finitedifference/convergence_column.jl b/test/Operators/finitedifference/convergence_column.jl index 0d2c4ff65c..8e642d4f32 100644 --- a/test/Operators/finitedifference/convergence_column.jl +++ b/test/Operators/finitedifference/convergence_column.jl @@ -7,6 +7,7 @@ using Random, StaticArrays, IntervalSets, LinearAlgebra using ClimaComms ClimaComms.@import_required_backends +import ClimaCore import ClimaCore: slab, Domains, Meshes, Topologies, Spaces, Fields, Operators import ClimaCore.Domains: Geometry import ClimaCore.DataLayouts: vindex @@ -249,9 +250,7 @@ end err_grad_sin_c = zeros(FT, length(n_elems_seq)) err_div_sin_c = zeros(FT, length(n_elems_seq)) err_grad_z_f = zeros(FT, length(n_elems_seq)) - err_grad_cos_f1 = zeros(FT, length(n_elems_seq)) err_grad_cos_f2 = zeros(FT, length(n_elems_seq)) - err_div_sin_f = zeros(FT, length(n_elems_seq)) err_div_cos_f = zeros(FT, length(n_elems_seq)) err_curl_sin_f = zeros(FT, length(n_elems_seq)) Δh = zeros(FT, length(n_elems_seq)) @@ -282,23 +281,6 @@ end divᶜ = Operators.DivergenceF2C() divsinᶜ = divᶜ.(Geometry.WVector.(sin.(faces))) - # Center -> Face operators: - # GradientC2F, SetValue - # f(z) = z - ∇ᶠ⁰ = Operators.GradientC2F( - left = Operators.SetValue(FT(0)), - right = Operators.SetValue(FT(pi)), - ) - ∂zᶠ = Geometry.WVector.(∇ᶠ⁰.(centers)) - - # GradientC2F, SetValue - # f(z) = cos(z) - ∇ᶠ¹ = Operators.GradientC2F( - left = Operators.SetValue(FT(1)), - right = Operators.SetValue(FT(-1)), - ) - gradcosᶠ¹ = Geometry.WVector.(∇ᶠ¹.(cos.(centers))) - # GradientC2F, SetGradient # f(z) = cos(z) ∇ᶠ² = Operators.GradientC2F( @@ -307,14 +289,6 @@ end ) gradcosᶠ² = Geometry.WVector.(∇ᶠ².(cos.(centers))) - # DivergenceC2F, SetValue - # f(z) = sin(z) - divᶠ⁰ = Operators.DivergenceC2F( - left = Operators.SetValue(Geometry.WVector(zero(FT))), - right = Operators.SetValue(Geometry.WVector(zero(FT))), - ) - divsinᶠ = divᶠ⁰.(Geometry.WVector.(sin.(centers))) - # DivergenceC2F, SetDivergence # f(z) = cos(z) divᶠ¹ = Operators.DivergenceC2F( @@ -324,8 +298,8 @@ end divcosᶠ = divᶠ¹.(Geometry.WVector.(cos.(centers))) curlᶠ = Operators.CurlC2F( - left = Operators.SetValue(Geometry.Covariant12Vector(zero(FT), zero(FT))), - right = Operators.SetValue(Geometry.Covariant12Vector(zero(FT), zero(FT))), + left = Operators.SetCurl(Geometry.Contravariant12Vector(zero(FT), one(FT))), + right = Operators.SetCurl(Geometry.Contravariant12Vector(zero(FT), -one(FT))), ) curlsinᶠ = curlᶠ.(Geometry.Covariant12Vector.(sin.(centers), zero(FT))) @@ -336,11 +310,7 @@ end # Errors err_grad_sin_c[k] = norm(gradsinᶜ .- Geometry.WVector.(cos.(centers))) err_div_sin_c[k] = norm(divsinᶜ .- cos.(centers)) - err_grad_z_f[k] = norm(∂zᶠ .- Geometry.WVector.(ones(FT, fs))) - err_grad_cos_f1[k] = norm(gradcosᶠ¹ .- Geometry.WVector.(.-sin.(faces))) err_grad_cos_f2[k] = norm(gradcosᶠ² .- Geometry.WVector.(.-sin.(faces))) - err_div_sin_f[k] = - norm(divsinᶠ .- (Geometry.WVector.(cos.(faces))).components.data.:1) err_div_cos_f[k] = norm( divcosᶠ .- (Geometry.WVector.(.-sin.(faces))).components.data.:1, ) @@ -354,15 +324,12 @@ end conv_div_sin_c = convergence_rate(err_div_sin_c, Δh) # GradientC2F conv, with f(z) = z, SetValue conv_grad_z = convergence_rate(err_grad_z_f, Δh) - # GradientC2F conv, with f(z) = cos(z), SetValue - conv_grad_cos_f1 = convergence_rate(err_grad_cos_f1, Δh) # GradientC2F conv, with f(z) = cos(z), SetGradient conv_grad_cos_f2 = convergence_rate(err_grad_cos_f2, Δh) # DivergenceC2F conv, with f(z) = sin(z), SetValue - conv_div_sin_f = convergence_rate(err_div_sin_f, Δh) # DivergenceC2F conv, with f(z) = cos(z), SetDivergence conv_div_cos_f = convergence_rate(err_div_cos_f, Δh) - # CurlC2F with f(z) = sin(z), SetValue + # CurlC2F with f(z) = sin(z), SetCurl conv_curl_sin_f = convergence_rate(err_curl_sin_f, Δh) # GradientF2C conv, with f(z) = sin(z) @@ -383,13 +350,6 @@ end @test norm(err_grad_z_f) ≤ 200 * eps(FT) # Convergence rate for this case is noisy because error very small - # GradientC2F conv, with f(z) = cos(z), SetValue - @test err_grad_cos_f1[3] ≤ err_grad_cos_f1[2] ≤ err_grad_cos_f1[1] ≤ 0.1 - @test conv_grad_cos_f1[1] ≈ 1.5 atol = 0.1 - @test conv_grad_cos_f1[2] ≈ 1.5 atol = 0.1 - @test conv_grad_cos_f1[3] ≈ 1.5 atol = 0.1 - # @test conv_grad_cos_f1[1] ≤ conv_grad_cos_f1[2] ≤ conv_grad_cos_f1[3] - # GradientC2F conv, with f(z) = cos(z), SetGradient @test err_grad_cos_f2[3] ≤ err_grad_cos_f2[2] ≤ err_grad_cos_f2[1] ≤ 0.1 @test conv_grad_cos_f2[1] ≈ 2 atol = 0.1 @@ -397,13 +357,6 @@ end @test conv_grad_cos_f2[3] ≈ 2 atol = 0.1 @test conv_grad_cos_f2[1] ≤ conv_grad_cos_f2[2] ≤ conv_grad_cos_f2[3] - # DivergenceC2F conv, with f(z) = sin(z), SetValue - @test err_div_sin_f[3] ≤ err_div_sin_f[2] ≤ err_div_sin_f[1] ≤ 0.1 - @test conv_div_sin_f[1] ≈ 2 atol = 0.1 - @test conv_div_sin_f[2] ≈ 2 atol = 0.1 - @test conv_div_sin_f[3] ≈ 2 atol = 0.1 - @test conv_div_sin_f[1] ≤ conv_div_sin_f[2] ≤ conv_div_sin_f[3] - # DivergenceC2F conv, with f(z) = cos(z), SetDivergence @test err_div_cos_f[3] ≤ err_div_cos_f[2] ≤ err_div_cos_f[1] ≤ 0.1 @test conv_div_cos_f[1] ≈ 2 atol = 0.1 @@ -411,7 +364,7 @@ end @test conv_div_cos_f[3] ≈ 2 atol = 0.1 @test conv_div_cos_f[1] ≤ conv_div_cos_f[2] ≤ conv_div_cos_f[3] - # CurlC2F with f(z) = sin(z), SetValue + # CurlC2F with f(z) = sin(z), SetCurl @test err_curl_sin_f[3] ≤ err_curl_sin_f[2] ≤ err_curl_sin_f[1] ≤ 0.1 @test conv_curl_sin_f[1] ≈ 2 atol = 0.1 @test conv_curl_sin_f[2] ≈ 2 atol = 0.1 @@ -1066,15 +1019,17 @@ end end end -@testset "Center -> Center Advection" begin +@testset "Center -> Face -> Center Advection" begin function advection(c, f, cs) adv = zeros(eltype(f), cs) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(0.0), - top = Operators.Extrapolate(), + gradc2f = Operators.GradientC2F( + bottom = Operators.SetGradient(Geometry.WVector(FT(1))), + top = Operators.SetGradient(Geometry.WVector(FT(1))), ) - return @. adv = A(c, f) + interpf2c = Operators.InterpolateF2C() + return @. adv = + interpf2c(LinearAlgebra.dot(Geometry.Contravariant3Vector(c), gradc2f(f))) end FT = Float64 diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 33d2cc6f68..6803c8b88f 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -33,19 +33,6 @@ function opt_RightBiasedF2C(face_field) return RB.(identity.(face_field)) end -function opt_AdvectionF2F(face_vel, face_field) - A = Operators.AdvectionF2F() - return A.(face_vel, identity.(face_field)) -end - -function opt_FluxCorrectionF2F_Extrapolate(face_vel, face_field) - FC = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - return FC.(face_vel, identity.(face_field)) -end - function opt_GradientF2C(face_field) ∇ᶜ = Operators.GradientF2C() return Geometry.WVector.(∇ᶜ.(sin.(face_field))) @@ -72,14 +59,6 @@ function opt_InterpolateC2F_SetValue(center_field) return I.(identity.(center_field)) end -function opt_InterpolateC2F_SetGradient(center_field) - I = Operators.InterpolateC2F( - left = Operators.SetGradient(Geometry.WVector(0.0)), - right = Operators.SetGradient(Geometry.WVector(0.0)), - ) - return I.(identity.(center_field)) -end - function opt_InterpolateC2F_Extrapolate(center_field) I = Operators.InterpolateC2F( left = Operators.Extrapolate(), @@ -96,14 +75,6 @@ function opt_WeightedInterpolateC2F_SetValue(weights, center_field) return identity.(WI.(weights, center_field)) end -function opt_WeightedInterpolateC2F_SetGradient(weights, center_field) - WI = Operators.WeightedInterpolateC2F( - left = Operators.SetGradient(Geometry.WVector(0.0)), - right = Operators.SetGradient(Geometry.WVector(0.0)), - ) - return identity.(WI.(weights, center_field)) -end - function opt_WeightedInterpolateC2F_Extrapolate(weights, center_field) WI = Operators.WeightedInterpolateC2F( left = Operators.Extrapolate(), @@ -122,14 +93,6 @@ function opt_RightBiasedC2F(center_field) return RB.(identity.(center_field)) end -function opt_UpwindBiasedProductC2F_SetValue(face_vel, center_field) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(0.0), - right = Operators.SetValue(0.0), - ) - return UB.(face_vel, identity.(center_field)) -end - function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field) UB = Operators.UpwindBiasedProductC2F( left = Operators.Extrapolate(), @@ -138,29 +101,6 @@ function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field) return UB.(face_vel, identity.(center_field)) end -function opt_AdvectionC2C_SetValue(face_vel, center_field) - A = Operators.AdvectionC2C( - left = Operators.SetValue(0.0), - right = Operators.SetValue(0.0), - ) - return A.(face_vel, identity.(center_field)) -end - -function opt_AdvectionC2C_Extrapolate(face_vel, center_field) - A = Operators.AdvectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - return A.(face_vel, identity.(center_field)) -end - -function opt_FluxCorrectionC2C_Extrapolate(face_vel, center_field) - FC = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - return FC.(face_vel, identity.(center_field)) -end function opt_GradientC2F_SetValue(center_field) ∇ᶠ = Operators.GradientC2F( @@ -178,14 +118,6 @@ function opt_GradientC2F_SetGradient(center_field) return Geometry.WVector.(∇ᶠ.(cos.(center_field))) end -function opt_DivergenceC2F_SetValue(center_field) - divᶠ = Operators.DivergenceC2F( - left = Operators.SetValue(Geometry.WVector(0.0)), - right = Operators.SetValue(Geometry.WVector(0.0)), - ) - return divᶠ.(Geometry.WVector.(sin.(center_field))) -end - function opt_DivergenceC2F_SetDivergence(center_field) # DivergenceC2F, SetDivergence divᶠ = Operators.DivergenceC2F( @@ -196,15 +128,6 @@ function opt_DivergenceC2F_SetDivergence(center_field) end -function opt_CurlC2F_SetValue(center_field) - # DivergenceC2F, SetDivergence - curlᶠ = Operators.CurlC2F( - left = Operators.SetValue(Geometry.Covariant1Vector(0.0)), - right = Operators.SetValue(Geometry.Covariant1Vector(0.0)), - ) - return curlᶠ.(Geometry.Covariant1Vector.(cos.(center_field))) -end - # Test that Julia ia able to optimize Stencil operations v1.7+ @static if @isdefined(var"@test_opt") @testset "Scalar Field FiniteDifferenceSpaces optimizations" begin @@ -240,13 +163,6 @@ end @test_opt opt_LeftBiasedF2C(faces) @test_opt opt_RightBiasedF2C(faces) - # @test_opt opt_AdvectionF2F(face_velocities, faces) - - @test_opt opt_FluxCorrectionF2F_Extrapolate( - center_velocities, - faces, - ) - @test_opt opt_GradientF2C(faces) @test_opt opt_DivergenceF2C(faces) @@ -257,17 +173,12 @@ end @test_opt function_filter = filter sum(sin.(centers)) @test_opt opt_InterpolateC2F_SetValue(centers) - @test_opt opt_InterpolateC2F_SetGradient(centers) @test_opt opt_InterpolateC2F_Extrapolate(centers) @test_opt opt_WeightedInterpolateC2F_SetValue( center_values, centers, ) - @test_opt opt_WeightedInterpolateC2F_SetGradient( - center_values, - centers, - ) @test_opt opt_WeightedInterpolateC2F_Extrapolate( center_values, centers, @@ -276,29 +187,15 @@ end @test_opt opt_LeftBiasedC2F(centers) @test_opt opt_RightBiasedC2F(centers) - @test_opt opt_UpwindBiasedProductC2F_SetValue( - face_velocities, - centers, - ) @test_opt opt_UpwindBiasedProductC2F_Extrapolate( face_velocities, centers, ) - @test_opt opt_AdvectionC2C_SetValue(face_velocities, centers) - @test_opt opt_AdvectionC2C_Extrapolate(face_velocities, centers) - - @test_opt opt_FluxCorrectionC2C_Extrapolate( - face_velocities, - centers, - ) - @test_opt opt_GradientC2F_SetValue(centers) @test_opt opt_GradientC2F_SetGradient(centers) - @test_opt opt_DivergenceC2F_SetValue(centers) @test_opt opt_DivergenceC2F_SetDivergence(centers) - @test_opt opt_CurlC2F_SetValue(centers) end end end diff --git a/test/Operators/finitedifference/opt_examples.jl b/test/Operators/finitedifference/opt_examples.jl index 06f98bc08b..f3447b5b26 100644 --- a/test/Operators/finitedifference/opt_examples.jl +++ b/test/Operators/finitedifference/opt_examples.jl @@ -232,8 +232,8 @@ function alloc_test_nested_expressions_4(cfield, ffield) top = Operators.SetValue(0), ) ∇f = Operators.DivergenceC2F(; - bottom = Operators.SetValue(wvec(0)), - top = Operators.SetValue(wvec(0)), + bottom = Operators.SetDivergence(0), + top = Operators.SetDivergence(0), ) LB = Operators.LeftBiasedF2C(; bottom = Operators.SetValue(1)) #! format: off @@ -269,9 +269,9 @@ function alloc_test_nested_expressions_6(cfield, ffield) wvec = Geometry.WVector Ic = Operators.InterpolateF2C() ∇f = Operators.DivergenceC2F(; - bottom = Operators.SetValue(wvec(0)), - top = Operators.SetValue(wvec(0)), - ) + bottom = Operators.SetDivergence(0), + top = Operators.SetDivergence(0), + ) #! format: off @. fz = fx * fy * ∇f(wvec(Ic(fy) * cx)) * ∇f(wvec(Ic(fy) * cx)) * fϕ * fψ # Compile first p = @allocated begin @@ -413,7 +413,10 @@ function alloc_test_nested_expressions_13( zero_bcs = (; bottom = Operators.SetValue(FT(0)), top = Operators.SetValue(FT(0))) I0f = Operators.InterpolateC2F(; zero_bcs...) - ∇f = Operators.DivergenceC2F(; adv_bcs...) + ∇f = Operators.DivergenceC2F(; + bottom = Operators.SetDivergence(0), + top = Operators.SetDivergence(0), + ) # Compile first @inbounds for i in 1:n_tuples @@ -494,14 +497,6 @@ end top = Operators.SetValue(0), ), ) - alloc_test_c2f_interp( - cfield, - ffield, - Operators.InterpolateC2F(; - bottom = Operators.SetGradient(wvec_glob(0)), - top = Operators.SetGradient(wvec_glob(0)), - ), - ) alloc_test_c2f_interp( cfield, ffield, @@ -526,8 +521,8 @@ end ffield, Operators.DivergenceF2C(), Operators.DivergenceC2F(; - bottom = Operators.SetValue(wvec_glob(0)), - top = Operators.SetValue(wvec_glob(0)), + bottom = Operators.SetDivergence(0), + top = Operators.SetDivergence(0), ), ) alloc_test_derivative( @@ -538,8 +533,8 @@ end top = Operators.SetValue(wvec_glob(0)), ), Operators.DivergenceC2F(; - bottom = Operators.SetValue(wvec_glob(0)), - top = Operators.SetValue(wvec_glob(0)), + bottom = Operators.SetDivergence(0), + top = Operators.SetDivergence(0), ), ) alloc_test_derivative( diff --git a/test/Operators/finitedifference/unit_column.jl b/test/Operators/finitedifference/unit_column.jl index deb6bf0fb7..30b996cfcd 100644 --- a/test/Operators/finitedifference/unit_column.jl +++ b/test/Operators/finitedifference/unit_column.jl @@ -48,19 +48,6 @@ device = ClimaComms.device() # Center -> Face operator # first order convergence at boundaries - ∇ᶠ = Operators.GradientC2F( - left = Operators.SetValue(FT(0)), - right = Operators.SetValue(FT(pi)), - ) - ∂z = Geometry.WVector.(∇ᶠ.(centers)) - @test ∂z ≈ Geometry.WVector.(ones(FT, face_space)) rtol = 10 * eps(FT) - - ∇ᶠ = Operators.GradientC2F( - left = Operators.SetValue(FT(1)), - right = Operators.SetValue(FT(-1)), - ) - ∂cos = Geometry.WVector.(∇ᶠ.(cos.(centers))) - @test ∂cos ≈ Geometry.WVector.(.-sin.(faces)) atol = 1e-1 ∇ᶠ = Operators.GradientC2F( left = Operators.SetGradient(Geometry.WVector(FT(0))), @@ -178,15 +165,6 @@ end ∂sin = Geometry.WVector.(∂.(w .* I.(θ))) @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2 - I = Operators.InterpolateC2F( - left = Operators.SetGradient(Geometry.WVector(FT(1))), - right = Operators.SetGradient(Geometry.WVector(FT(-1))), - ) - ∂ = Operators.GradientF2C() - - ∂sin = Geometry.WVector.(∂.(w .* I.(θ))) - @test ∂sin ≈ Geometry.WVector.(cos.(centers)) atol = 1e-2 - # 3) we set boundaries on both: 2nd should take precedence I = Operators.InterpolateC2F( left = Operators.SetValue(FT(NaN)), @@ -298,62 +276,3 @@ end end @test all(cy_ref .== parent(ClimaCore.to_cpu(cy))) end - -# https://github.com/CliMA/ClimaCore.jl/issues/994 -# TODO: make this test more low-level / granular (test `getidx`). -@testset "Spatially varying BC with Grad" begin - FT = Float64 - zmin = FT(1.0) - zmax = FT(2.0) - xlim = FT.((0.0, 10.0)) - ylim = FT.((0.0, 1.0)) - zlim = FT.((zmin, zmax)) - nelements = (1, 1, 5) - npolynomial = 3 - domain_x = Domains.IntervalDomain( - Geometry.XPoint(xlim[1]), - Geometry.XPoint(xlim[2]); - periodic = true, - ) - domain_y = Domains.IntervalDomain( - Geometry.YPoint(ylim[1]), - Geometry.YPoint(ylim[2]); - periodic = true, - ) - plane = Domains.RectangleDomain(domain_x, domain_y) - context = ClimaComms.context() - mesh = Meshes.RectilinearMesh(plane, nelements[1], nelements[2]) - grid_topology = Topologies.Topology2D(context, mesh) - quad = Spaces.Quadratures.GLL{npolynomial + 1}() - horzspace = Spaces.SpectralElementSpace2D(grid_topology, quad) - - - vertdomain = Domains.IntervalDomain( - Geometry.ZPoint(zlim[1]), - Geometry.ZPoint(zlim[2]); - boundary_names = (:bottom, :top), - ) - vertmesh = Meshes.IntervalMesh(vertdomain, nelems = nelements[3]) - vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh) - - hv_center_space = - Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space) - - surface_field = Fields.zeros(horzspace) - - ψ = Fields.ones(hv_center_space) - value = surface_field # same type/instance of underlying space as horizontal space of \psi - - gradc2f_no_bc = Operators.GradientC2F() - divf2c = Operators.DivergenceF2C( - top = Operators.SetValue(Geometry.WVector.(value)), - bottom = Operators.SetValue(Geometry.WVector(FT(0.0))), - ) - @. divf2c(gradc2f_no_bc(ψ)) # runs - - gradc2f = Operators.GradientC2F(; - top = Operators.SetValue(value), - bottom = Operators.SetValue{FT}(0.0), - ) - gradc2f.(ψ) # fails -end diff --git a/test/Operators/hybrid/convergence_2d.jl b/test/Operators/hybrid/convergence_2d.jl index 261ecc9cac..bf015d8452 100644 --- a/test/Operators/hybrid/convergence_2d.jl +++ b/test/Operators/hybrid/convergence_2d.jl @@ -4,11 +4,13 @@ include("utils_2d.jl") function advection(c, f, hv_center_space) adv = zeros(eltype(f), hv_center_space) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(0.0), - top = Operators.Extrapolate(), + gradc2f = Operators.GradientC2F( + bottom = Operators.SetGradient(Geometry.WVector(1.0)), + top = Operators.SetGradient(Geometry.WVector(1.0)), ) - return @. adv = A(c, f) + interpf2c = Operators.InterpolateF2C() + return @. adv = + interpf2c(LinearAlgebra.dot(Geometry.Contravariant3Vector(c), gradc2f(f))) end n_elems_seq = 2 .^ (5, 6, 7, 8) @@ -41,83 +43,6 @@ include("utils_2d.jl") @test conv_adv_c2c[3] ≈ 2 atol = 0.1 end -@testset "1D SE, 1D FD Extruded Domain Discrete Product Rule Operations" begin - - gradc2f = Operators.GradientC2F( - top = Operators.SetValue(0.0), - bottom = Operators.SetValue(0.0), - ) - gradf2c = Operators.GradientF2C() - - n_elems_seq = 2 .^ (5, 6, 7, 8) - err, Δh = zeros(length(n_elems_seq)), zeros(length(n_elems_seq)) - - for (k, n) in enumerate(n_elems_seq) - # Discrete Product Rule Test - # ∂(ab)/∂s = a̅∂b/∂s + b̅∂a/∂s - # a, b are interface variables, and ̅ represents interpolation - # s represents the vertical coordinate, in our case `z` - # For this test, we use a(z) = z and b = sin(z), - hv_center_space, hv_face_space = hvspace_2D(helem = n, velem = n) - ᶠz = Fields.coordinate_field(hv_face_space).z - ᶜz = Fields.coordinate_field(hv_center_space).z - Δh[k] = 1.0 / n - - # advective velocity - # scalar-valued field to be advected - ∂ab_numerical = @. Geometry.WVector(gradf2c(ᶠz * sin(ᶠz))) - ∂ab_analytical = @. Geometry.WVector(ᶜz * cos(ᶜz) + sin(ᶜz)) - - err[k] = norm(∂ab_numerical .- ∂ab_analytical) - end - # Solution convergence rate - grad_pr = convergence_rate(err, Δh) - @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1 - @test grad_pr[1] ≈ 2 atol = 0.1 - @test grad_pr[2] ≈ 2 atol = 0.1 - @test grad_pr[3] ≈ 2 atol = 0.1 -end - -@testset "1D SE, 1D FD Extruded Domain Discrete Product Rule Operations: Stretched" begin - - gradc2f = Operators.GradientC2F( - top = Operators.SetValue(0.0), - bottom = Operators.SetValue(0.0), - ) - gradf2c = Operators.GradientF2C() - - n_elems_seq = 2 .^ (5, 6, 7, 8) - err, Δh = zeros(length(n_elems_seq)), zeros(length(n_elems_seq)) - - for (k, n) in enumerate(n_elems_seq) - # Discrete Product Rule Test - # ∂(ab)/∂s = a̅∂b/∂s + b̅∂a/∂s - # a, b are interface variables, and ̅ represents interpolation - # s represents the vertical coordinate, in our case `z` - # For this test, we use a(z) = z and b = sin(z), - hv_center_space, hv_face_space = hvspace_2D( - helem = n, - velem = n, - stretch = Meshes.ExponentialStretching(2π), - ) - ᶠz = Fields.coordinate_field(hv_face_space).z - ᶜz = Fields.coordinate_field(hv_center_space).z - Δh[k] = 1.0 / n - - # advective velocity - # scalar-valued field to be advected - ∂ab_numerical = @. Geometry.WVector(gradf2c(ᶠz * sin(ᶠz))) - ∂ab_analytical = @. Geometry.WVector(ᶜz * cos(ᶜz) + sin(ᶜz)) - - err[k] = norm(∂ab_numerical .- ∂ab_analytical) - end - # Solution convergence rate - grad_pr = convergence_rate(err, Δh) - @test err[3] ≤ err[2] ≤ err[1] ≤ 0.2 - @test grad_pr[1] ≈ 2 atol = 0.1 - @test grad_pr[2] ≈ 2 atol = 0.1 - @test grad_pr[3] ≈ 2 atol = 0.1 -end @testset "1D SE, 1D FD Extruded Domain horz & vert divergence operator, with Extrapolate BCs" begin diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 574558cde8..6270155b99 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -34,19 +34,6 @@ function opt_RightBiasedF2C(face_field) return RB.(identity.(face_field)) end -function opt_AdvectionF2F(face_vel, face_field) - A = Operators.AdvectionF2F() - return A.(face_vel, identity.(face_field)) -end - -function opt_FluxCorrectionF2F_Extrapolate(face_vel, face_field) - FC = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - return FC.(face_vel, identity.(face_field)) -end - function opt_GradientF2C(face_field) ∇ᶜ = Operators.GradientF2C() return Geometry.WVector.(∇ᶜ.(sin.(face_field))) @@ -73,14 +60,6 @@ function opt_InterpolateC2F_SetValue(center_field) return I.(identity.(center_field)) end -function opt_InterpolateC2F_SetGradient(center_field) - I = Operators.InterpolateC2F( - left = Operators.SetGradient(Geometry.WVector(0.0)), - right = Operators.SetGradient(Geometry.WVector(0.0)), - ) - return I.(identity.(center_field)) -end - function opt_InterpolateC2F_Extrapolate(center_field) I = Operators.InterpolateC2F( left = Operators.Extrapolate(), @@ -97,14 +76,6 @@ function opt_WeightedInterpolateC2F_SetValue(weights, center_field) return identity.(WI.(weights, center_field)) end -function opt_WeightedInterpolateC2F_SetGradient(weights, center_field) - WI = Operators.WeightedInterpolateC2F( - left = Operators.SetGradient(Geometry.WVector(0.0)), - right = Operators.SetGradient(Geometry.WVector(0.0)), - ) - return identity.(WI.(weights, center_field)) -end - function opt_WeightedInterpolateC2F_Extrapolate(weights, center_field) WI = Operators.WeightedInterpolateC2F( left = Operators.Extrapolate(), @@ -123,14 +94,6 @@ function opt_RightBiasedC2F(center_field) return RB.(identity.(center_field)) end -function opt_UpwindBiasedProductC2F_SetValue(face_vel, center_field) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(0.0), - right = Operators.SetValue(0.0), - ) - return UB.(face_vel, identity.(center_field)) -end - function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field) UB = Operators.UpwindBiasedProductC2F( left = Operators.Extrapolate(), @@ -139,30 +102,6 @@ function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field) return UB.(face_vel, identity.(center_field)) end -function opt_AdvectionC2C_SetValue(face_vel, center_field) - A = Operators.AdvectionC2C( - left = Operators.SetValue(0.0), - right = Operators.SetValue(0.0), - ) - return A.(face_vel, identity.(center_field)) -end - -function opt_AdvectionC2C_Extrapolate(face_vel, center_field) - A = Operators.AdvectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - return A.(face_vel, identity.(center_field)) -end - -function opt_FluxCorrectionC2C_Extrapolate(face_vel, center_field) - FC = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - return FC.(face_vel, identity.(center_field)) -end - function opt_GradientC2F_SetValue(center_field) ∇ᶠ = Operators.GradientC2F( left = Operators.SetValue(1.0), @@ -179,14 +118,6 @@ function opt_GradientC2F_SetGradient(center_field) return Geometry.WVector.(∇ᶠ.(cos.(center_field))) end -function opt_DivergenceC2F_SetValue(center_field) - divᶠ = Operators.DivergenceC2F( - left = Operators.SetValue(Geometry.WVector(0.0)), - right = Operators.SetValue(Geometry.WVector(0.0)), - ) - return divᶠ.(Geometry.WVector.(sin.(center_field))) -end - function opt_DivergenceC2F_SetDivergence(center_field) # DivergenceC2F, SetDivergence divᶠ = Operators.DivergenceC2F( @@ -196,15 +127,6 @@ function opt_DivergenceC2F_SetDivergence(center_field) return divᶠ.(Geometry.WVector.(cos.(center_field))) end -function opt_CurlC2F_SetValue(center_field) - # DivergenceC2F, SetDivergence - curlᶠ = Operators.CurlC2F( - left = Operators.SetValue(Geometry.Covariant1Vector(0.0)), - right = Operators.SetValue(Geometry.Covariant1Vector(0.0)), - ) - return curlᶠ.(Geometry.Covariant1Vector.(cos.(center_field))) -end - function hspace1d(FT) hdomain = Domains.IntervalDomain( Geometry.XPoint{FT}(-pi) .. Geometry.XPoint{FT}(pi), @@ -271,13 +193,6 @@ end @test_opt opt_LeftBiasedF2C(faces) @test_opt opt_RightBiasedF2C(faces) - # @test_opt opt_AdvectionF2F(face_velocities, faces) - - @test_opt opt_FluxCorrectionF2F_Extrapolate( - center_velocities, - faces, - ) - @test_opt opt_GradientF2C(faces) @test_opt opt_GradientF2C(faces) @test_opt opt_DivergenceF2C(faces) @@ -289,17 +204,12 @@ end @test_opt function_filter = filter sum(sin.(centers)) @test_opt opt_InterpolateC2F_SetValue(centers) - @test_opt opt_InterpolateC2F_SetGradient(centers) @test_opt opt_InterpolateC2F_Extrapolate(centers) @test_opt opt_WeightedInterpolateC2F_SetValue( center_values, centers, ) - @test_opt opt_WeightedInterpolateC2F_SetGradient( - center_values, - centers, - ) @test_opt opt_WeightedInterpolateC2F_Extrapolate( center_values, centers, @@ -308,29 +218,14 @@ end @test_opt opt_LeftBiasedC2F(centers) @test_opt opt_RightBiasedC2F(centers) - @test_opt opt_UpwindBiasedProductC2F_SetValue( - face_velocities, - centers, - ) @test_opt opt_UpwindBiasedProductC2F_Extrapolate( face_velocities, centers, ) - - @test_opt opt_AdvectionC2C_SetValue(face_velocities, centers) - @test_opt opt_AdvectionC2C_Extrapolate(face_velocities, centers) - - @test_opt opt_FluxCorrectionC2C_Extrapolate( - face_velocities, - centers, - ) - @test_opt opt_GradientC2F_SetValue(centers) @test_opt opt_GradientC2F_SetGradient(centers) - @test_opt opt_DivergenceC2F_SetValue(centers) @test_opt opt_DivergenceC2F_SetDivergence(centers) - @test_opt opt_CurlC2F_SetValue(centers) end end diff --git a/test/Operators/hybrid/simulation_3d.jl b/test/Operators/hybrid/simulation_3d.jl index 1777fa1e4f..bf37adac49 100644 --- a/test/Operators/hybrid/simulation_3d.jl +++ b/test/Operators/hybrid/simulation_3d.jl @@ -15,11 +15,13 @@ device = ClimaComms.device() ) function rhs!(dudt, u, _, t) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(sin(-t)), - top = Operators.Extrapolate(), + gradc2f = Operators.GradientC2F( + bottom = Operators.SetGradient(Geometry.WVector(cos(-t))), + top = Operators.SetGradient(Geometry.WVector(cos(-t))), ) - return @. dudt = -A(V, u) + interpf2c = Operators.InterpolateF2C() + return @. dudt = + -1 * interpf2c(LinearAlgebra.dot(Geometry.Contravariant3Vector(V), gradc2f(u))) end U = sin.(Fields.coordinate_field(hv_center_space).z) @@ -37,7 +39,6 @@ device = ClimaComms.device() @test sol_column_field ≈ ref_column_field rtol = 0.6 end end - @testset "2D SE, 1D FD Extruded Domain ∇ ODE Solve horizontal" begin # Advection Equation @@ -63,7 +64,6 @@ end Δt = 0.01 prob = ODEProblem(rhs!, U, (0.0, 2π)) sol = solve(prob, SSPRK33(), dt = Δt) - @test U ≈ sol.u[end] rtol = 1e-6 end diff --git a/test/Operators/hybrid/unit_2d.jl b/test/Operators/hybrid/unit_2d.jl index f14295126c..20990264af 100644 --- a/test/Operators/hybrid/unit_2d.jl +++ b/test/Operators/hybrid/unit_2d.jl @@ -86,21 +86,11 @@ end # => C ∂_z F = cos(z) hv_center_space, hv_face_space = hvspace_2D() - function advect(f) - advf = zeros(eltype(f), hv_center_space) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(Geometry.UVector(0.0)), # value of f at the boundary (UVector field) - top = Operators.Extrapolate(), - ) - @. advf = A(vC, f) - end - # Vertical advective velocity vC = Geometry.WVector.(ones(Float64, hv_face_space),) # vector-valued field to be advected (one component only, a UVector) f = Geometry.UVector.(sin.(Fields.coordinate_field(hv_center_space).z),) - advf = advect(f) function div!(F) vecdivf = zeros(eltype(F), hv_center_space) @@ -145,10 +135,11 @@ end diffu = zeros(eltype(u), hv_center_space) gradc2f = Operators.GradientC2F( - top = Operators.SetValue(0.0), - bottom = Operators.SetValue(0.0), ) - divf2c = Operators.DivergenceF2C() + divf2c = Operators.DivergenceF2C(; + top = Operators.SetDivergence(0.0), + bottom = Operators.SetDivergence(0.0), + ) @. diffu = divf2c(K * gradc2f(u)) hgrad = Operators.Gradient() @@ -197,11 +188,11 @@ end function vec_diff(K, U) vec_diff = zeros(eltype(U), hv_center_space) - gradc2f = Operators.GradientC2F( - top = Operators.SetValue(Geometry.UVector(0.0)), - bottom = Operators.SetValue(Geometry.UVector(0.0)), + gradc2f = Operators.GradientC2F() + divf2c = Operators.DivergenceF2C(; + top = Operators.SetDivergence(Geometry.UVector(0.0)), + bottom = Operators.SetDivergence(Geometry.UVector(0.0)), ) - divf2c = Operators.DivergenceF2C() @. vec_diff = divf2c(K * gradc2f(U)) hgrad = Operators.Gradient() @@ -262,8 +253,8 @@ end curl = Operators.Curl() curlC2F = Operators.CurlC2F( - bottom = Operators.SetValue(Geometry.Covariant1Vector(0.0)), - top = Operators.SetValue(Geometry.Covariant1Vector(0.0)), + bottom = Operators.SetCurl(Geometry.Contravariant2Vector(2.0)), + top = Operators.SetCurl(Geometry.Contravariant2Vector(2.0)), ) curlu = curlC2F.(u) @@ -297,7 +288,7 @@ end ccoords = Fields.coordinate_field(hv_center_space) fcoords = Fields.coordinate_field(hv_face_space) fcoords_1 = Fields.level(fcoords, ClimaCore.Utilities.half) - curl_bcfield₁ = Geometry.Covariant1Vector.(0.0 .* fcoords_1.z) + curl_bcfield₁ = Geometry.Contravariant2Vector(2.0) curl_bcfield² = Geometry.Contravariant2Vector.(0.0 .* fcoords_1.z) u = @@ -313,8 +304,8 @@ end curl = Operators.Curl() curlC2F = Operators.CurlC2F( - bottom = Operators.SetValue(curl_bcfield₁), - top = Operators.SetValue(Geometry.Covariant1Vector(0.0)), + bottom = Operators.SetCurl(curl_bcfield₁), + top = Operators.SetCurl(Geometry.Contravariant2Vector(2.0)), ) curlu = curlC2F.(u)