From d03e32101f1c490a72b6c627c3970cc5d4920543 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 27 May 2026 14:08:11 -0700 Subject: [PATCH 01/18] Delete unused test --- .../Operators/finitedifference/unit_column.jl | 72 ------------------- 1 file changed, 72 deletions(-) diff --git a/test/Operators/finitedifference/unit_column.jl b/test/Operators/finitedifference/unit_column.jl index deb6bf0fb7..ebb87eae15 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))), @@ -298,62 +285,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 From befc611aae6fc183d49a94531d503d747341523e Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 9 Jun 2026 13:05:47 -0700 Subject: [PATCH 02/18] rm SetValue bc for gradc2f --- .buildkite/pipeline.yml | 32 ----- docs/Project.toml | 1 + docs/tutorials/introduction.jl | 22 ++-- examples/column/advect_diffusion.jl | 10 +- examples/column/ekman.jl | 6 +- examples/column/heat.jl | 5 +- examples/column/hydrostatic_ekman.jl | 4 +- examples/column/wave.jl | 7 +- src/MatrixFields/operator_matrices.jl | 10 -- src/Operators/finitedifference.jl | 37 +----- test/MatrixFields/operator_matrices.jl | 118 ------------------ .../benchmark_stencils_climacore_kernels.jl | 7 -- .../benchmark_stencils_utils.jl | 5 - .../convergence_advection_diffusion1d.jl | 6 + .../finitedifference/convergence_column.jl | 29 ----- test/Operators/hybrid/convergence_2d.jl | 78 ------------ 16 files changed, 41 insertions(+), 336 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7fdc1b51f8..1ffe1b347c 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -1916,38 +1916,6 @@ steps: artifact_paths: - "examples/column/output/zalesak_fct_advection/*" - - label: ":computer: Column Advection Diffusion Eq" - key: "cpu_column_advect_diff" - retry: *retry_policy - command: - - "julia --color=yes --project=.buildkite examples/column/advect_diffusion.jl" - artifact_paths: - - "examples/column/output/advect_diffusion/*" - - - label: ":computer: Column Ekman Eq" - key: "cpu_column_ekman" - retry: *retry_policy - command: - - "julia --color=yes --project=.buildkite examples/column/ekman.jl" - artifact_paths: - - "examples/column/output/ekman/*" - - - label: ":computer: Column Hydrostatic Ekman Eq" - key: "cpu_column_hydrostatic_ekman" - retry: *retry_policy - command: - - "julia --color=yes --project=.buildkite examples/column/hydrostatic_ekman.jl" - artifact_paths: - - "examples/column/output/hydrostatic_ekman/*" - - - label: ":computer: Column Wave Eq" - key: "cpu_column_wave" - retry: *retry_policy - command: - - "julia --color=yes --project=.buildkite examples/column/wave.jl" - artifact_paths: - - "examples/column/output/wave/*" - - label: ":computer: Column Hydrostatic Balance Eq" key: "cpu_column_hydrostatic" retry: *retry_policy 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/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_diffusion.jl b/examples/column/advect_diffusion.jl index e2bf74f6d6..6bb1ae4def 100644 --- a/examples/column/advect_diffusion.jl +++ b/examples/column/advect_diffusion.jl @@ -77,9 +77,13 @@ function ∑tendencies!(dT, T, z, t) # 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_bottom = Fields.level(T, 1) + bc_vb_T = Operators.SetGradient( + Geometry.Covariant3Vector.( + 2 .* (Fields.level(T, 1) .- FT(gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))) + ), + ) + gradc2f = Operators.GradientC2F(bottom = bc_vb_T, top = bc_gt) divf2c = Operators.DivergenceF2C() return @. dT = divf2c(ν * gradc2f(T)) - A(V, T) diff --git a/examples/column/ekman.jl b/examples/column/ekman.jl index 766ada9cb9..4598c94d45 100644 --- a/examples/column/ekman.jl +++ b/examples/column/ekman.jl @@ -92,14 +92,16 @@ function tendency!(dY, Y, _, t) # 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)) + bcs_top = Operators.SetGradient(@. 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 # 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)) + bcs_top = Operators.SetGradient(@. 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 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..074b4dd63a 100644 --- a/examples/column/hydrostatic_ekman.jl +++ b/examples/column/hydrostatic_ekman.jl @@ -143,7 +143,9 @@ function tendency!(dY, Y, _, t) # 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)[] + bcs_top = Operators.SetGradient(Geometry.Covariant3Vector(1) ⊗ ((uvg - uv_top_val))) ∂c = Operators.DivergenceF2C(bottom = bcs_bottom) ∂f = Operators.GradientC2F(top = bcs_top) duv .= (uv .- Ref(uvg)) .× Ref(Geometry.WVector(f)) 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/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index f84d1ee0b7..9f7ecb4177 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -742,16 +742,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, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 5896625d6b..f24aa20128 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -2829,11 +2829,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 +2838,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 +2862,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( diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index cb6deadee9..c3162fb685 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -164,7 +164,6 @@ end 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,)) @@ -212,7 +211,6 @@ end ᶠ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) @@ -222,7 +220,6 @@ end ᶠ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) @@ -334,119 +331,4 @@ end 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..c0ccd7192c 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) diff --git a/test/Operators/finitedifference/benchmark_stencils_utils.jl b/test/Operators/finitedifference/benchmark_stencils_utils.jl index fc91914d59..2ae050903e 100755 --- a/test/Operators/finitedifference/benchmark_stencils_utils.jl +++ b/test/Operators/finitedifference/benchmark_stencils_utils.jl @@ -216,8 +216,6 @@ bcs_tested(c, ::typeof(op_Upwind3rdOrderBiasedProductC2F!)) = (set_upwind_biased # 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!)) = @@ -406,7 +404,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!, @@ -454,7 +451,6 @@ function test_results_column(t_min) [(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], @@ -500,7 +496,6 @@ function test_results_sphere(t_min) [(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], diff --git a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl index 1ca603b8d2..eb33913b33 100644 --- a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl +++ b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl @@ -73,6 +73,12 @@ convergence_rate(err, Δh) = bottom = bc_vb, top = Operators.Extrapolate(), ) + bc_vb = Operators.SetGradient( + Geometry.WVector( + FT(∇gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ)), + ), + ) + # TODO: Delete this test?? gradc2f = Operators.GradientC2F(; bottom = bc_vb, top = bc_gt) divf2c = Operators.DivergenceF2C() return @. dT = divf2c(ν * gradc2f(T)) - A(V, T) diff --git a/test/Operators/finitedifference/convergence_column.jl b/test/Operators/finitedifference/convergence_column.jl index 0d2c4ff65c..0c8607f726 100644 --- a/test/Operators/finitedifference/convergence_column.jl +++ b/test/Operators/finitedifference/convergence_column.jl @@ -249,7 +249,6 @@ 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)) @@ -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( @@ -336,8 +318,6 @@ 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) @@ -354,8 +334,6 @@ 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 @@ -383,13 +361,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 diff --git a/test/Operators/hybrid/convergence_2d.jl b/test/Operators/hybrid/convergence_2d.jl index 261ecc9cac..30cf689cb4 100644 --- a/test/Operators/hybrid/convergence_2d.jl +++ b/test/Operators/hybrid/convergence_2d.jl @@ -41,84 +41,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 # Divergence operator in 2D Cartesian domain with From e85cf868e93c94b5a32a3e290c6636652f724ebe Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 9 Jun 2026 13:17:06 -0700 Subject: [PATCH 03/18] rm interpc2f setgradient bc --- src/MatrixFields/operator_matrices.jl | 10 ----- src/Operators/finitedifference.jl | 40 +------------------ test/Operators/finitedifference/opt.jl | 9 ----- .../Operators/finitedifference/unit_column.jl | 9 ----- test/Operators/hybrid/opt.jl | 9 ----- 5 files changed, 1 insertion(+), 76 deletions(-) diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index 9f7ecb4177..280354485f 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -390,16 +390,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, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index f24aa20128..0800f68cb1 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, diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 33d2cc6f68..f8c1cbc3e3 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -72,14 +72,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(), @@ -257,7 +249,6 @@ 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( diff --git a/test/Operators/finitedifference/unit_column.jl b/test/Operators/finitedifference/unit_column.jl index ebb87eae15..30b996cfcd 100644 --- a/test/Operators/finitedifference/unit_column.jl +++ b/test/Operators/finitedifference/unit_column.jl @@ -165,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)), diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 574558cde8..6cff697635 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -73,14 +73,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(), @@ -289,7 +281,6 @@ 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( From 1ecc1ab9b7f57d3cf448e7e011579a5bd5c0a908 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 2 Jun 2026 15:14:04 -0700 Subject: [PATCH 04/18] rm weightedinterpc2f setgradient --- src/MatrixFields/operator_matrices.jl | 10 ------- src/Operators/finitedifference.jl | 38 +------------------------- test/Operators/finitedifference/opt.jl | 12 -------- test/Operators/hybrid/opt.jl | 12 -------- 4 files changed, 1 insertion(+), 71 deletions(-) diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index 280354485f..17aa517d03 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -458,16 +458,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, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 0800f68cb1..4e3905c0a2 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -1070,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). @@ -1081,7 +1080,7 @@ struct WeightedInterpolateC2F{BCS} <: WeightedInterpolationOperator assert_valid_bcs( "WeightedInterpolateC2F", kwargs, - (SetValue, SetGradient, Extrapolate), + (SetValue, Extrapolate), ) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end @@ -1137,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, diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index f8c1cbc3e3..5f3a2e7f10 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -88,14 +88,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(), @@ -255,10 +247,6 @@ end center_values, centers, ) - @test_opt opt_WeightedInterpolateC2F_SetGradient( - center_values, - centers, - ) @test_opt opt_WeightedInterpolateC2F_Extrapolate( center_values, centers, diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 6cff697635..c17e761364 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -89,14 +89,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(), @@ -287,10 +279,6 @@ end center_values, centers, ) - @test_opt opt_WeightedInterpolateC2F_SetGradient( - center_values, - centers, - ) @test_opt opt_WeightedInterpolateC2F_Extrapolate( center_values, centers, From 4ea3ac461161038c51d597ae8c7b19913851b3f3 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 9 Jun 2026 14:10:10 -0700 Subject: [PATCH 05/18] rm UpwindBiasedProductC2F setval bc --- .buildkite/pipeline.yml | 16 ------- docs/src/examples.md | 3 -- examples/column/advect.jl | 59 ++++++++++++++++++++++---- examples/column/step.jl | 56 ++++++++++++++++++++---- src/MatrixFields/operator_matrices.jl | 24 ----------- src/Operators/finitedifference.jl | 49 +-------------------- test/MatrixFields/operator_matrices.jl | 2 - test/Operators/finitedifference/opt.jl | 12 ------ test/Operators/hybrid/opt.jl | 12 ------ 9 files changed, 99 insertions(+), 134 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 1ffe1b347c..f9690cf36e 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -1848,22 +1848,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 - command: - - "julia --color=yes --project=.buildkite examples/column/advect.jl" - artifact_paths: - - "examples/column/output/advect/*" - - label: ":computer: Column FCT Advection Eq" key: "cpu_fct_column_advect" retry: *retry_policy diff --git a/docs/src/examples.md b/docs/src/examples.md index a00b801f1a..c8a091c3ad 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -40,9 +40,6 @@ This is discretized using the following This test case is set up in a 1D column domain ``z \in [0, 1]`` and discretized into a mesh of 10 elements. A homogeneous Dirichlet boundary condition is set at the bottom boundary, `bcs_bottom`, setting the temperature to 0. A Neumann boundary condition is applied to the top boundary, `bcs_top`, setting the temperature gradient to 1. -### Advection equation - -The 1D Column advection example in [`examples/column/advect.jl`](https://github.com/CliMA/ClimaCore.jl/blob/main/examples/column/advect.jl). #### Equations and Discretizations diff --git a/examples/column/advect.jl b/examples/column/advect.jl index 5d8aee9a46..411a3ce3e0 100644 --- a/examples/column/advect.jl +++ b/examples/column/advect.jl @@ -8,7 +8,8 @@ import ClimaCore: DataLayouts, Operators, Geometry, - Spaces + Spaces, + Utilities using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 @@ -47,13 +48,33 @@ function tendency1!(dθ, θ, _, t) left = Operators.Extrapolate(), right = Operators.Extrapolate(), ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a - t)), - right = Operators.SetValue(sin(b - 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() - return @. dθ = -∂(UB(V, θ)) + return @. dθ = -∂(set_bcs(UB(V, θ))) end function tendency2!(dθ, θ, _, t) fcc = Operators.FluxCorrectionC2C( @@ -64,12 +85,32 @@ function tendency2!(dθ, θ, _, t) left = Operators.Extrapolate(), right = Operators.Extrapolate(), ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a - t)), - right = Operators.SetValue(sin(b - 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() - return @. dθ = -∂(UB(V, θ)) + fcc(V, θ) + return @. dθ = -∂(set_bcs(UB(V, θ))) + fcc(V, θ) end # use the advection operator function tendency3!(dθ, θ, _, t) diff --git a/examples/column/step.jl b/examples/column/step.jl index 12e68a743f..47da7445f9 100644 --- a/examples/column/step.jl +++ b/examples/column/step.jl @@ -49,13 +49,33 @@ function tendency1!(dθ, θ, _, t) left = Operators.Extrapolate(), right = Operators.Extrapolate(), ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a.z - t)), - right = Operators.SetValue(sin(b.z - 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.z - 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.z - 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() - return @. dθ = -∂(UB(V, θ)) + return @. dθ = -∂(set_bcs(UB(V, θ))) end function tendency2!(dθ, θ, _, t) fcc = Operators.FluxCorrectionC2C( @@ -66,12 +86,32 @@ function tendency2!(dθ, θ, _, t) left = Operators.Extrapolate(), right = Operators.Extrapolate(), ) - UB = Operators.UpwindBiasedProductC2F( - left = Operators.SetValue(sin(a.z - t)), - right = Operators.SetValue(sin(b.z - 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.z - 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.z - 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() - return @. dθ = -∂(UB(V, θ)) + fcc(V, θ) + return @. dθ = -∂(set_bcs(UB(V, θ))) + fcc(V, θ) end # use the advection operator function tendency3!(dθ, θ, _, t) diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index 17aa517d03..4ebf24ec72 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -487,30 +487,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, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 4e3905c0a2..74728b0b14 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -1185,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 @@ -1206,7 +1197,7 @@ struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator assert_valid_bcs( "UpwindBiasedProductC2F", kwargs, - (SetValue, Extrapolate), + (Extrapolate,), ) new{typeof(NamedTuple(kwargs))}(NamedTuple(kwargs)) end @@ -1246,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, diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index c3162fb685..56ea94f46b 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -208,7 +208,6 @@ end ᶜlbias = LeftBiasedF2C() ᶠrbias = RightBiasedC2F(; set_nested_values.top) ᶜwinterp = WeightedInterpolateF2C() - ᶠupwind = UpwindBiasedProductC2F(; set_scalar_values...) ᶜadvect = AdvectionC2C(; extrapolate...) ᶜflux_correct = FluxCorrectionC2C(; extrapolate...) ᶜdiv = DivergenceF2C() @@ -217,7 +216,6 @@ end ᶜ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) ᶜdiv_matrix = MatrixFields.operator_matrix(ᶜdiv) diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 5f3a2e7f10..792b39e6f0 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -106,14 +106,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(), @@ -255,10 +247,6 @@ 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, diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index c17e761364..445b4a3077 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -107,14 +107,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(), @@ -287,10 +279,6 @@ 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, From 33fb619c64b3ed66f13d5f43e50272d3619a57a1 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 9 Jun 2026 16:54:25 -0700 Subject: [PATCH 06/18] rm divc2f setval bc --- .../plane/density_current_2dinvariant_rhoe.jl | 31 +++++++++--- examples/hybrid/plane/topo_agnesi_nh.jl | 24 ++++++++-- examples/hybrid/plane/topo_schar_nh.jl | 12 ----- src/MatrixFields/operator_matrices.jl | 22 --------- src/Operators/finitedifference.jl | 47 ------------------- test/MatrixFields/operator_matrices.jl | 4 -- .../finitedifference/convergence_column.jl | 19 -------- test/Operators/finitedifference/opt.jl | 9 ---- .../finitedifference/opt_examples.jl | 31 +++++------- test/Operators/hybrid/opt.jl | 9 ---- 10 files changed, 58 insertions(+), 150 deletions(-) diff --git a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl index 20df8bb6a8..d8d0810319 100644 --- a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl +++ b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl @@ -16,7 +16,8 @@ import ClimaCore: Spaces, Quadratures, Fields, - Operators + Operators, + Utilities using ClimaCore.Geometry using Logging: global_logger @@ -177,10 +178,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 +244,28 @@ 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 = @. Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * + (2 * lg_bottom_face.invJ) + top_divergence = + @. Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * lg_top_face.invJ) + 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..2f082c4344 100644 --- a/examples/hybrid/plane/topo_agnesi_nh.jl +++ b/examples/hybrid/plane/topo_agnesi_nh.jl @@ -233,10 +233,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,6 +300,25 @@ function rhs_invariant!(dY, Y, _, t) hκ₂∇²uₕ = @. hwdiv(κ₂ * ᶜ∇ₕuₕ) vκ₂∇²uₕ = @. vdivf2c(κ₂ * ᶠ∇ᵥuₕ) hκ₂∇²w = @. hwdiv(κ₂ * ᶠ∇ₕ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 = @. Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * + (2 * lg_bottom_face.invJ) + top_divergence = + @. Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * lg_top_face.invJ) + set_bcs = Operators.SetBoundaryOperator( + bottom = Operators.SetValue(bottom_divergence), + top = Operators.SetValue(top_divergence), + ) vκ₂∇²w = @. 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/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index 4ebf24ec72..7c82ab1745 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -743,28 +743,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, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 74728b0b14..59fd5a39c3 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -3074,11 +3074,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 @@ -3117,48 +3112,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, diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index 56ea94f46b..3e9ee37fd1 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -126,7 +126,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 +138,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, @@ -169,7 +166,6 @@ end 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,)) diff --git a/test/Operators/finitedifference/convergence_column.jl b/test/Operators/finitedifference/convergence_column.jl index 0c8607f726..bf61114117 100644 --- a/test/Operators/finitedifference/convergence_column.jl +++ b/test/Operators/finitedifference/convergence_column.jl @@ -250,7 +250,6 @@ end err_div_sin_c = zeros(FT, length(n_elems_seq)) err_grad_z_f = 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)) @@ -289,14 +288,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( @@ -319,8 +310,6 @@ end err_grad_sin_c[k] = norm(gradsinᶜ .- Geometry.WVector.(cos.(centers))) err_div_sin_c[k] = norm(divsinᶜ .- cos.(centers)) 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, ) @@ -337,7 +326,6 @@ end # 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 @@ -368,13 +356,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 diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 792b39e6f0..bc0bb6ddc1 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -154,14 +154,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( @@ -263,7 +255,6 @@ end @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 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/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 445b4a3077..95e0c5cd6a 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -155,14 +155,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( @@ -295,7 +287,6 @@ end @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 From 1db159aa20e49798b7b76a4384c760de25c53f78 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 9 Jun 2026 16:55:17 -0700 Subject: [PATCH 07/18] undo pipeline changes --- .buildkite/pipeline.yml | 48 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index f9690cf36e..7fdc1b51f8 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -1848,6 +1848,22 @@ 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 + command: + - "julia --color=yes --project=.buildkite examples/column/advect.jl" + artifact_paths: + - "examples/column/output/advect/*" + - label: ":computer: Column FCT Advection Eq" key: "cpu_fct_column_advect" retry: *retry_policy @@ -1900,6 +1916,38 @@ steps: artifact_paths: - "examples/column/output/zalesak_fct_advection/*" + - label: ":computer: Column Advection Diffusion Eq" + key: "cpu_column_advect_diff" + retry: *retry_policy + command: + - "julia --color=yes --project=.buildkite examples/column/advect_diffusion.jl" + artifact_paths: + - "examples/column/output/advect_diffusion/*" + + - label: ":computer: Column Ekman Eq" + key: "cpu_column_ekman" + retry: *retry_policy + command: + - "julia --color=yes --project=.buildkite examples/column/ekman.jl" + artifact_paths: + - "examples/column/output/ekman/*" + + - label: ":computer: Column Hydrostatic Ekman Eq" + key: "cpu_column_hydrostatic_ekman" + retry: *retry_policy + command: + - "julia --color=yes --project=.buildkite examples/column/hydrostatic_ekman.jl" + artifact_paths: + - "examples/column/output/hydrostatic_ekman/*" + + - label: ":computer: Column Wave Eq" + key: "cpu_column_wave" + retry: *retry_policy + command: + - "julia --color=yes --project=.buildkite examples/column/wave.jl" + artifact_paths: + - "examples/column/output/wave/*" + - label: ":computer: Column Hydrostatic Balance Eq" key: "cpu_column_hydrostatic" retry: *retry_policy From 8c23d5ef226b3bf4dbcd08289ab44e83b3663352 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 10 Jun 2026 15:13:22 -0700 Subject: [PATCH 08/18] rm curlc2f setvalue bc --- .buildkite/pipeline.yml | 5 --- examples/column/step.jl | 3 +- src/Operators/finitedifference.jl | 36 ++----------------- test/MatrixFields/operator_matrices.jl | 4 --- .../benchmark_stencils_utils.jl | 25 ++++--------- .../finitedifference/convergence_column.jl | 8 ++--- test/Operators/finitedifference/opt.jl | 10 ------ test/Operators/hybrid/opt.jl | 10 ------ test/Operators/hybrid/unit_2d.jl | 25 ++++++------- 9 files changed, 27 insertions(+), 99 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7fdc1b51f8..50fca7ec1c 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -252,11 +252,6 @@ steps: retry: *retry_policy command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/geometry.jl" - - label: "Unit: tensors" - key: unit_tensors - retry: *retry_policy - command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/tensors.jl" - - label: "Unit: mul_with_projection" key: unit_mul_with_projection retry: *retry_policy diff --git a/examples/column/step.jl b/examples/column/step.jl index 47da7445f9..fc25442c36 100644 --- a/examples/column/step.jl +++ b/examples/column/step.jl @@ -9,7 +9,8 @@ import ClimaCore: DataLayouts, Operators, Geometry, - Spaces + Spaces, + Utilities using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 59fd5a39c3..8ac5a2fabf 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -3082,7 +3082,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...) @@ -3171,19 +3171,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...) @@ -3215,32 +3209,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 3e9ee37fd1..2f358e188a 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -172,7 +172,6 @@ end 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()) @@ -197,7 +196,6 @@ 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...) @@ -207,7 +205,6 @@ end ᶜadvect = AdvectionC2C(; extrapolate...) ᶜflux_correct = FluxCorrectionC2C(; extrapolate...) ᶜ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) @@ -215,7 +212,6 @@ end ᶜadvect_matrix = MatrixFields.operator_matrix(ᶜadvect) ᶜflux_correct_matrix = MatrixFields.operator_matrix(ᶜflux_correct) ᶜdiv_matrix = MatrixFields.operator_matrix(ᶜdiv) - ᶠcurl_matrix = MatrixFields.operator_matrix(ᶠcurl) @test_throws "does not contain any Fields" @. ᶜlbias_matrix() * ᶠinterp_matrix() diff --git a/test/Operators/finitedifference/benchmark_stencils_utils.jl b/test/Operators/finitedifference/benchmark_stencils_utils.jl index 2ae050903e..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,8 +202,8 @@ 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) @@ -221,7 +214,7 @@ bcs_tested(c, ::typeof(op_divgrad_FF!)) = 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)), @@ -432,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], @@ -447,13 +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_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], ] @@ -475,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], @@ -492,13 +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_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_column.jl b/test/Operators/finitedifference/convergence_column.jl index bf61114117..54bbffec51 100644 --- a/test/Operators/finitedifference/convergence_column.jl +++ b/test/Operators/finitedifference/convergence_column.jl @@ -297,8 +297,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))) @@ -328,7 +328,7 @@ end # DivergenceC2F conv, with f(z) = sin(z), SetValue # 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) @@ -363,7 +363,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 diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index bc0bb6ddc1..9de8fd9e8a 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -164,15 +164,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 @@ -256,7 +247,6 @@ end @test_opt opt_GradientC2F_SetGradient(centers) @test_opt opt_DivergenceC2F_SetDivergence(centers) - @test_opt opt_CurlC2F_SetValue(centers) end end end diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 95e0c5cd6a..117dfad102 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -164,15 +164,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), @@ -288,7 +279,6 @@ end @test_opt opt_GradientC2F_SetGradient(centers) @test_opt opt_DivergenceC2F_SetDivergence(centers) - @test_opt opt_CurlC2F_SetValue(centers) end end diff --git a/test/Operators/hybrid/unit_2d.jl b/test/Operators/hybrid/unit_2d.jl index f14295126c..19d465c362 100644 --- a/test/Operators/hybrid/unit_2d.jl +++ b/test/Operators/hybrid/unit_2d.jl @@ -145,10 +145,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 +198,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 +263,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 +298,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 +314,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) From 0dc61f16e47fab120e364074effb8f4ad9ce0f90 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 16 Jun 2026 11:14:04 -0700 Subject: [PATCH 09/18] rm AdvectionC2C --- docs/src/examples.md | 4 +- docs/src/operators.md | 1 - examples/column/advect.jl | 60 +++---- examples/column/advect_diffusion.jl | 34 ++-- examples/column/ekman.jl | 24 ++- examples/column/hydrostatic_ekman.jl | 20 ++- examples/column/step.jl | 60 +++---- src/MatrixFields/operator_matrices.jl | 78 --------- src/Operators/finitedifference.jl | 156 ------------------ test/MatrixFields/operator_matrices.jl | 39 ++--- .../benchmark_stencils_climacore_kernels.jl | 1 - .../convergence_advection_diffusion1d.jl | 33 ++-- .../finitedifference/convergence_column.jl | 49 ------ test/Operators/finitedifference/opt.jl | 19 --- test/Operators/hybrid/convergence_2d.jl | 40 ----- test/Operators/hybrid/opt.jl | 18 -- test/Operators/hybrid/simulation_3d.jl | 34 ---- test/Operators/hybrid/unit_2d.jl | 10 -- 18 files changed, 122 insertions(+), 558 deletions(-) diff --git a/docs/src/examples.md b/docs/src/examples.md index c8a091c3ad..65b720e756 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -81,9 +81,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..24f2d0b9ad 100644 --- a/docs/src/operators.md +++ b/docs/src/operators.md @@ -75,7 +75,6 @@ AbstractTVDSlopeLimiter GradientF2C GradientC2F AdvectionF2F -AdvectionC2C DivergenceF2C DivergenceC2F CurlC2F diff --git a/examples/column/advect.jl b/examples/column/advect.jl index 411a3ce3e0..2bb55f75d6 100644 --- a/examples/column/advect.jl +++ b/examples/column/advect.jl @@ -40,14 +40,6 @@ V = Geometry.WVector.(ones(FT, fs)) # upwinding function tendency1!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) 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)) @@ -81,10 +73,6 @@ function tendency2!(dθ, θ, _, t) left = Operators.Extrapolate(), right = Operators.Extrapolate(), ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) 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)) @@ -114,20 +102,18 @@ function tendency2!(dθ, θ, _, t) 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, θ) + left_center = Fields.level(θ, 1) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) + ) + left_gradient = + Operators.SetGradient(@. 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) @@ -136,15 +122,19 @@ function tendency4!(dθ, θ, _, t) 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, θ) + left_center = Fields.level(θ, 1) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) + ) + left_gradient = + Operators.SetGradient(@. 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(θ))) + + fcc(V, θ) end # use the advection operator diff --git a/examples/column/advect_diffusion.jl b/examples/column/advect_diffusion.jl index 6bb1ae4def..fdab554620 100644 --- a/examples/column/advect_diffusion.jl +++ b/examples/column/advect_diffusion.jl @@ -52,41 +52,31 @@ 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()) - T_bottom = Fields.level(T, 1) - bc_vb_T = Operators.SetGradient( + top_center_left_biased_grad = Geometry.Covariant3Vector.( - 2 .* (Fields.level(T, 1) .- FT(gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))) - ), - ) - gradc2f = Operators.GradientC2F(bottom = bc_vb_T, top = bc_gt) + 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 @show ∑tendencies!(similar(T), T, nothing, 0.0); diff --git a/examples/column/ekman.jl b/examples/column/ekman.jl index 4598c94d45..c6c9a3244e 100644 --- a/examples/column/ekman.jl +++ b/examples/column/ekman.jl @@ -85,26 +85,36 @@ 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 u_top = Fields.level(u, Fields.nlevels(u)) + u_bottom = Fields.level(u, 1) bcs_top = Operators.SetGradient(@. 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(@. Geometry.Covariant3Vector(2 * u_bottom)) + bcs_top_advection = Operators.SetGradient(@. 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 v_top = Fields.level(v, Fields.nlevels(v)) + v_bottom = Fields.level(v, 1) bcs_top = Operators.SetGradient(@. 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(@. Geometry.Covariant3Vector(2 * v_bottom)) + bcs_top_advection = Operators.SetGradient(@. 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/hydrostatic_ekman.jl b/examples/column/hydrostatic_ekman.jl index 074b4dd63a..42b36adc94 100644 --- a/examples/column/hydrostatic_ekman.jl +++ b/examples/column/hydrostatic_ekman.jl @@ -136,20 +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) 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( @@ -177,12 +183,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 index fc25442c36..f8d847ac5c 100644 --- a/examples/column/step.jl +++ b/examples/column/step.jl @@ -42,14 +42,6 @@ V = Geometry.WVector.(ones(FT, fs)) # upwinding function tendency1!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) 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)) @@ -83,10 +75,6 @@ function tendency2!(dθ, θ, _, t) left = Operators.Extrapolate(), right = Operators.Extrapolate(), ) - fcf = Operators.FluxCorrectionF2F( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) 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)) @@ -116,20 +104,18 @@ function tendency2!(dθ, θ, _, t) 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, θ) + left_center = Fields.level(θ, 1) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) + ) + left_gradient = + Operators.SetGradient(@. 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) @@ -138,15 +124,19 @@ function tendency4!(dθ, θ, _, t) 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, θ) + left_center = Fields.level(θ, 1) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) + ) + left_gradient = + Operators.SetGradient(@. 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(θ))) + + fcc(V, θ) end # use the advection operator diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index 7c82ab1745..0114798905 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -19,7 +19,6 @@ const TwoArgFDOperatorWithCenterInput = Union{ Operators.WeightedInterpolateC2F, Operators.UpwindBiasedProductC2F, Operators.Upwind3rdOrderBiasedProductC2F, - Operators.AdvectionC2C, Operators.FluxCorrectionC2C, } const TwoArgFDOperatorWithFaceInput = Union{ @@ -578,63 +577,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, @@ -827,26 +769,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 8ac5a2fabf..a2af5a9bf8 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -2164,162 +2164,6 @@ Base.@propagate_inbounds function stencil_interior( 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) diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index 2f358e188a..bf2a304333 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -27,7 +27,6 @@ import ClimaCore.Operators: Upwind3rdOrderBiasedProductC2F, FCTBorisBook, FCTZalesak, - AdvectionC2C, AdvectionF2F, FluxCorrectionC2C, FluxCorrectionF2F, @@ -154,8 +153,6 @@ 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)) @@ -202,14 +199,12 @@ end ᶜlbias = LeftBiasedF2C() ᶠrbias = RightBiasedC2F(; set_nested_values.top) ᶜwinterp = WeightedInterpolateF2C() - ᶜadvect = AdvectionC2C(; extrapolate...) ᶜflux_correct = FluxCorrectionC2C(; extrapolate...) ᶜdiv = DivergenceF2C() ᶠinterp_matrix = MatrixFields.operator_matrix(ᶠinterp) ᶜlbias_matrix = MatrixFields.operator_matrix(ᶜlbias) ᶠrbias_matrix = MatrixFields.operator_matrix(ᶠrbias) ᶜwinterp_matrix = MatrixFields.operator_matrix(ᶜwinterp) - ᶜadvect_matrix = MatrixFields.operator_matrix(ᶜadvect) ᶜflux_correct_matrix = MatrixFields.operator_matrix(ᶜflux_correct) ᶜdiv_matrix = MatrixFields.operator_matrix(ᶜdiv) @@ -237,7 +232,6 @@ end test_name = "product of six operator matrices", get_result = @lazy( @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * @@ -245,7 +239,6 @@ end ), set_result = @lazy( @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * @@ -258,7 +251,6 @@ end matrices", get_result = @lazy( @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * @@ -267,7 +259,6 @@ end ), set_result = @lazy( @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜadvect_matrix(ᶠuvw) * ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * @@ -277,10 +268,7 @@ end ref_set_result = @lazy( @. ᶜflux_correct( ᶠuvw, - ᶜadvect( - ᶠuvw, - ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))), - ), + ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))), ) ), ) @@ -290,31 +278,28 @@ end 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. diff --git a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl index c0ccd7192c..22b8c16dfd 100644 --- a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl +++ b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl @@ -188,7 +188,6 @@ 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/convergence_advection_diffusion1d.jl b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl index eb33913b33..510e1418e8 100644 --- a/test/Operators/finitedifference/convergence_advection_diffusion1d.jl +++ b/test/Operators/finitedifference/convergence_advection_diffusion1d.jl @@ -61,27 +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(), - ) - bc_vb = Operators.SetGradient( - Geometry.WVector( - FT(∇gaussian(z₀, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ)), - ), - ) - # TODO: Delete this test?? - 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 54bbffec51..886405b8f7 100644 --- a/test/Operators/finitedifference/convergence_column.jl +++ b/test/Operators/finitedifference/convergence_column.jl @@ -1017,52 +1017,3 @@ end @test conv_adv_wc[3] ≈ 2.5 atol = 0.1 end end - -@testset "Center -> Center Advection" begin - - function advection(c, f, cs) - adv = zeros(eltype(f), cs) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(0.0), - top = Operators.Extrapolate(), - ) - return @. adv = A(c, f) - end - - FT = Float64 - n_elems_seq = 2 .^ (5, 6, 7, 8) - err = zeros(FT, length(n_elems_seq)) - Δh = zeros(FT, length(n_elems_seq)) - device = ClimaComms.device() - - for (k, n) in enumerate(n_elems_seq) - domain = Domains.IntervalDomain( - Geometry.ZPoint{FT}(0.0), - Geometry.ZPoint{FT}(4π); - boundary_names = (:bottom, :top), - ) - mesh = Meshes.IntervalMesh(domain; nelems = n) - - cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) - fs = Spaces.face_space(cs) - - # advective velocity - c = Geometry.WVector.(ones(Float64, fs),) - # scalar-valued field to be advected - f = sin.(Fields.coordinate_field(cs).z) - - # Call the advection operator - adv = advection(c, f, cs) - - ClimaComms.allowscalar(device) do - Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)] - end - err[k] = norm(adv .- cos.(Fields.coordinate_field(cs).z)) - end - # AdvectionC2C convergence rate - conv_adv_c2c = convergence_rate(err, Δh) - @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1 - @test conv_adv_c2c[1] ≈ 2 atol = 0.1 - @test conv_adv_c2c[2] ≈ 2 atol = 0.1 - @test conv_adv_c2c[3] ≈ 2 atol = 0.1 -end diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 9de8fd9e8a..73b4c8fc55 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -114,22 +114,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(), @@ -235,9 +219,6 @@ end 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, diff --git a/test/Operators/hybrid/convergence_2d.jl b/test/Operators/hybrid/convergence_2d.jl index 30cf689cb4..b82430bc76 100644 --- a/test/Operators/hybrid/convergence_2d.jl +++ b/test/Operators/hybrid/convergence_2d.jl @@ -1,45 +1,5 @@ include("utils_2d.jl") -@testset "1D SE, 1D FD Extruded Domain vertical advection operator" begin - - 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(), - ) - return @. adv = A(c, f) - end - - 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) - # Advection Operator - # c ∂_z f - # for this test, we use f(z) = sin(z) and c = 1, a WVector - # => c ∂_z f = cos(z) - hv_center_space, hv_face_space = hvspace_2D(helem = n, velem = n) - - Δh[k] = 1.0 / n - - # advective velocity - c = Geometry.WVector.(ones(Float64, hv_face_space),) - # scalar-valued field to be advected - f = sin.(Fields.coordinate_field(hv_center_space).z) - - # Call the advection operator - adv = advection(c, f, hv_center_space) - - err[k] = norm(adv .- cos.(Fields.coordinate_field(hv_center_space).z)) - end - # AdvectionC2C convergence rate - conv_adv_c2c = convergence_rate(err, Δh) - @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1 - @test conv_adv_c2c[1] ≈ 2 atol = 0.1 - @test conv_adv_c2c[2] ≈ 2 atol = 0.1 - @test conv_adv_c2c[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 117dfad102..c6cd3a7fa2 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -115,22 +115,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(), @@ -267,8 +251,6 @@ end centers, ) - @test_opt opt_AdvectionC2C_SetValue(face_velocities, centers) - @test_opt opt_AdvectionC2C_Extrapolate(face_velocities, centers) @test_opt opt_FluxCorrectionC2C_Extrapolate( face_velocities, diff --git a/test/Operators/hybrid/simulation_3d.jl b/test/Operators/hybrid/simulation_3d.jl index 1777fa1e4f..e5bfe25ed0 100644 --- a/test/Operators/hybrid/simulation_3d.jl +++ b/test/Operators/hybrid/simulation_3d.jl @@ -3,40 +3,6 @@ include("utils_3d.jl") using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 device = ClimaComms.device() -@testset "2D SE, 1D FV Extruded Domain ∇ ODE Solve vertical" begin - - hv_center_space, hv_face_space = - hvspace_3D(; device = ClimaComms.CPUSingleThreaded()) - V = - Geometry.UVWVector.( - zeros(Float64, hv_face_space), - zeros(Float64, hv_face_space), - ones(Float64, hv_face_space), - ) - - function rhs!(dudt, u, _, t) - A = Operators.AdvectionC2C( - bottom = Operators.SetValue(sin(-t)), - top = Operators.Extrapolate(), - ) - return @. dudt = -A(V, u) - end - - U = sin.(Fields.coordinate_field(hv_center_space).z) - dudt = zeros(eltype(U), hv_center_space) - rhs!(dudt, U, nothing, 0.0) - - Δt = 0.01 - prob = ODEProblem(rhs!, U, (0.0, 2π)) - sol = solve(prob, SSPRK33(), dt = Δt) - - htopo = Spaces.topology(hv_center_space) - for h in 1:Topologies.nlocalelems(htopo) - sol_column_field = ClimaCore.column(sol.u[end], 1, 1, h) - ref_column_field = ClimaCore.column(U, 1, 1, h) - @test sol_column_field ≈ ref_column_field rtol = 0.6 - end -end @testset "2D SE, 1D FD Extruded Domain ∇ ODE Solve horizontal" begin diff --git a/test/Operators/hybrid/unit_2d.jl b/test/Operators/hybrid/unit_2d.jl index 19d465c362..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) From a3b9e79d1840dcdc8048b8796e4ce3192e15d4e5 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 16 Jun 2026 13:19:52 -0700 Subject: [PATCH 10/18] Add LinVanleerC2F to non-linear ops const --- src/MatrixFields/operator_matrices.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index 0114798905..f1ff79aed5 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -32,7 +32,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} From f40d40342ec99974346b44ab4dab5bc59778be5c Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 16 Jun 2026 13:32:01 -0700 Subject: [PATCH 11/18] rm AdvectionF2F --- docs/src/operators.md | 1 - examples/column/hydrostatic_ekman.jl | 6 ++- src/MatrixFields/operator_matrices.jl | 12 ----- src/Operators/finitedifference.jl | 49 ------------------- test/MatrixFields/operator_matrices.jl | 2 - .../benchmark_stencils_climacore_kernels.jl | 1 - test/Operators/finitedifference/opt.jl | 7 --- test/Operators/hybrid/opt.jl | 7 --- 8 files changed, 4 insertions(+), 81 deletions(-) diff --git a/docs/src/operators.md b/docs/src/operators.md index 24f2d0b9ad..c4066e7176 100644 --- a/docs/src/operators.md +++ b/docs/src/operators.md @@ -74,7 +74,6 @@ AbstractTVDSlopeLimiter ```@docs GradientF2C GradientC2F -AdvectionF2F DivergenceF2C DivergenceC2F CurlC2F diff --git a/examples/column/hydrostatic_ekman.jl b/examples/column/hydrostatic_ekman.jl index 42b36adc94..5449a7bb34 100644 --- a/examples/column/hydrostatic_ekman.jl +++ b/examples/column/hydrostatic_ekman.jl @@ -164,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))), @@ -172,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 diff --git a/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index f1ff79aed5..c7e86b6b17 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -23,7 +23,6 @@ const TwoArgFDOperatorWithCenterInput = Union{ } const TwoArgFDOperatorWithFaceInput = Union{ Operators.WeightedInterpolateF2C, - Operators.AdvectionF2F, Operators.FluxCorrectionF2F, } const ErroneousFDOperator = Union{ @@ -578,17 +577,6 @@ end op_matrix_row_type(::Operators.AdvectionOperator, ::Type{FT}, _) where {FT} = TridiagonalMatrixRow{FT} -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, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index a2af5a9bf8..b59f275dca 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -2115,55 +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 = FluxCorrectionC2C(;boundaries) diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index bf2a304333..25a9535896 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -27,7 +27,6 @@ import ClimaCore.Operators: Upwind3rdOrderBiasedProductC2F, FCTBorisBook, FCTZalesak, - AdvectionF2F, FluxCorrectionC2C, FluxCorrectionF2F, SetBoundaryOperator, @@ -153,7 +152,6 @@ end (ᶠuvw, ᶜscalar), true, ) - 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,)) diff --git a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl index 22b8c16dfd..24086ec7bd 100644 --- a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl +++ b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl @@ -187,7 +187,6 @@ LeftBiased3rdOrderC2F # unused LeftBiased3rdOrderF2C # unused RightBiased3rdOrderC2F # unused RightBiased3rdOrderF2C # unused -AdvectionF2F # only used in ClimaAtmos src/ FluxCorrectionC2C # unused at the moment FluxCorrectionF2F # only used in ClimaAtmos src/ ``` diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 73b4c8fc55..c91022fbf6 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -33,11 +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(), @@ -183,8 +178,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, diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index c6cd3a7fa2..dc63a523c3 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -34,11 +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(), @@ -214,8 +209,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, From 16779312be5161042db5906876e5299fadb24274 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 16 Jun 2026 16:09:38 -0700 Subject: [PATCH 12/18] merge pipeline --- .buildkite/pipeline.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 50fca7ec1c..7fdc1b51f8 100755 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -252,6 +252,11 @@ steps: retry: *retry_policy command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/geometry.jl" + - label: "Unit: tensors" + key: unit_tensors + retry: *retry_policy + command: "julia --color=yes --check-bounds=yes --project=.buildkite test/Geometry/tensors.jl" + - label: "Unit: mul_with_projection" key: unit_mul_with_projection retry: *retry_policy From 2fa7fd7beebbd93d47db9fa8e79031b6ecd70d9a Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Tue, 16 Jun 2026 17:00:41 -0700 Subject: [PATCH 13/18] Re-add deleted tests --- .../finitedifference/convergence_column.jl | 52 +++++++++++++++++++ test/Operators/hybrid/convergence_2d.jl | 43 +++++++++++++++ test/Operators/hybrid/simulation_3d.jl | 36 ++++++++++++- 3 files changed, 130 insertions(+), 1 deletion(-) diff --git a/test/Operators/finitedifference/convergence_column.jl b/test/Operators/finitedifference/convergence_column.jl index 886405b8f7..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 @@ -1017,3 +1018,54 @@ end @test conv_adv_wc[3] ≈ 2.5 atol = 0.1 end end + +@testset "Center -> Face -> Center Advection" begin + + function advection(c, f, cs) + adv = zeros(eltype(f), cs) + gradc2f = Operators.GradientC2F( + bottom = Operators.SetGradient(Geometry.WVector(FT(1))), + top = Operators.SetGradient(Geometry.WVector(FT(1))), + ) + interpf2c = Operators.InterpolateF2C() + return @. adv = + interpf2c(LinearAlgebra.dot(Geometry.Contravariant3Vector(c), gradc2f(f))) + end + + FT = Float64 + n_elems_seq = 2 .^ (5, 6, 7, 8) + err = zeros(FT, length(n_elems_seq)) + Δh = zeros(FT, length(n_elems_seq)) + device = ClimaComms.device() + + for (k, n) in enumerate(n_elems_seq) + domain = Domains.IntervalDomain( + Geometry.ZPoint{FT}(0.0), + Geometry.ZPoint{FT}(4π); + boundary_names = (:bottom, :top), + ) + mesh = Meshes.IntervalMesh(domain; nelems = n) + + cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) + fs = Spaces.face_space(cs) + + # advective velocity + c = Geometry.WVector.(ones(Float64, fs),) + # scalar-valued field to be advected + f = sin.(Fields.coordinate_field(cs).z) + + # Call the advection operator + adv = advection(c, f, cs) + + ClimaComms.allowscalar(device) do + Δh[k] = Spaces.local_geometry_data(fs).J[vindex(1)] + end + err[k] = norm(adv .- cos.(Fields.coordinate_field(cs).z)) + end + # AdvectionC2C convergence rate + conv_adv_c2c = convergence_rate(err, Δh) + @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1 + @test conv_adv_c2c[1] ≈ 2 atol = 0.1 + @test conv_adv_c2c[2] ≈ 2 atol = 0.1 + @test conv_adv_c2c[3] ≈ 2 atol = 0.1 +end diff --git a/test/Operators/hybrid/convergence_2d.jl b/test/Operators/hybrid/convergence_2d.jl index b82430bc76..f5376f62f5 100644 --- a/test/Operators/hybrid/convergence_2d.jl +++ b/test/Operators/hybrid/convergence_2d.jl @@ -1,5 +1,48 @@ include("utils_2d.jl") +@testset "1D SE, 1D FD Extruded Domain vertical advection operator" begin + + function advection(c, f, hv_center_space) + adv = zeros(eltype(f), hv_center_space) + gradc2f = Operators.GradientC2F( + bottom = Operators.SetGradient(Geometry.WVector(FT(1))), + top = Operators.SetGradient(Geometry.WVector(FT(1))), + ) + interpf2c = Operators.InterpolateF2C() + return @. adv = + interpf2c(LinearAlgebra.dot(Geometry.Contravariant3Vector(c), gradc2f(f))) + end + + 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) + # Advection Operator + # c ∂_z f + # for this test, we use f(z) = sin(z) and c = 1, a WVector + # => c ∂_z f = cos(z) + hv_center_space, hv_face_space = hvspace_2D(helem = n, velem = n) + + Δh[k] = 1.0 / n + + # advective velocity + c = Geometry.WVector.(ones(Float64, hv_face_space),) + # scalar-valued field to be advected + f = sin.(Fields.coordinate_field(hv_center_space).z) + + # Call the advection operator + adv = advection(c, f, hv_center_space) + + err[k] = norm(adv .- cos.(Fields.coordinate_field(hv_center_space).z)) + end + # AdvectionC2C convergence rate + conv_adv_c2c = convergence_rate(err, Δh) + @test err[3] ≤ err[2] ≤ err[1] ≤ 0.1 + @test conv_adv_c2c[1] ≈ 2 atol = 0.1 + @test conv_adv_c2c[2] ≈ 2 atol = 0.1 + @test conv_adv_c2c[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/simulation_3d.jl b/test/Operators/hybrid/simulation_3d.jl index e5bfe25ed0..bf37adac49 100644 --- a/test/Operators/hybrid/simulation_3d.jl +++ b/test/Operators/hybrid/simulation_3d.jl @@ -3,7 +3,42 @@ include("utils_3d.jl") using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 device = ClimaComms.device() +@testset "2D SE, 1D FV Extruded Domain ∇ ODE Solve vertical" begin + + hv_center_space, hv_face_space = + hvspace_3D(; device = ClimaComms.CPUSingleThreaded()) + V = + Geometry.UVWVector.( + zeros(Float64, hv_face_space), + zeros(Float64, hv_face_space), + ones(Float64, hv_face_space), + ) + + function rhs!(dudt, u, _, t) + gradc2f = Operators.GradientC2F( + bottom = Operators.SetGradient(Geometry.WVector(cos(-t))), + top = Operators.SetGradient(Geometry.WVector(cos(-t))), + ) + interpf2c = Operators.InterpolateF2C() + return @. dudt = + -1 * interpf2c(LinearAlgebra.dot(Geometry.Contravariant3Vector(V), gradc2f(u))) + end + + U = sin.(Fields.coordinate_field(hv_center_space).z) + dudt = zeros(eltype(U), hv_center_space) + rhs!(dudt, U, nothing, 0.0) + Δt = 0.01 + prob = ODEProblem(rhs!, U, (0.0, 2π)) + sol = solve(prob, SSPRK33(), dt = Δt) + + htopo = Spaces.topology(hv_center_space) + for h in 1:Topologies.nlocalelems(htopo) + sol_column_field = ClimaCore.column(sol.u[end], 1, 1, h) + ref_column_field = ClimaCore.column(U, 1, 1, h) + @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 @@ -29,7 +64,6 @@ device = ClimaComms.device() Δ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 From e8bfb7590c9136713fc4aa03c55254aab951e4fc Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 17 Jun 2026 15:29:34 -0700 Subject: [PATCH 14/18] rm FluxCorrectionC2C --- examples/column/advect.jl | 450 ++++++++++-------- examples/column/step.jl | 254 ---------- .../hybrid/box/bubble_3d_invariant_rhoe.jl | 16 +- .../box/bubble_3d_invariant_rhotheta.jl | 15 +- examples/hybrid/box/bubble_3d_rhotheta.jl | 12 +- examples/hybrid/flux_correction_utils.jl | 56 +++ .../hybrid/plane/bubble_2d_invariant_rhoe.jl | 12 - .../plane/density_current_2d_flux_form.jl | 19 +- .../hybrid/staggered_nonhydrostatic_model.jl | 14 - src/MatrixFields/operator_matrices.jl | 34 -- src/Operators/finitedifference.jl | 98 ---- test/MatrixFields/operator_matrices.jl | 30 +- .../benchmark_stencils_climacore_kernels.jl | 1 - test/Operators/finitedifference/opt.jl | 12 - test/Operators/hybrid/opt.jl | 15 - 15 files changed, 326 insertions(+), 712 deletions(-) delete mode 100644 examples/column/step.jl create mode 100644 examples/hybrid/flux_correction_utils.jl diff --git a/examples/column/advect.jl b/examples/column/advect.jl index 2bb55f75d6..803a720e43 100644 --- a/examples/column/advect.jl +++ b/examples/column/advect.jl @@ -28,217 +28,267 @@ 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) - -V = Geometry.WVector.(ones(FT, fs)) -θ = sin.(Fields.coordinate_field(cs).z) - -# Solve advection Equation: ∂θ/dt = -∂(vθ) - -# 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() +for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step)) - return @. dθ = -∂(set_bcs(UB(V, θ))) -end -function tendency2!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), + cs = Spaces.CenterFiniteDifferenceSpace(device, mesh) + fs = Spaces.FaceFiniteDifferenceSpace(cs) + + 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 + + # Solve advection Equation: ∂θ/dt = -∂(vθ) + + # 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() + + 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) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- + Fields.level(θ, Fields.nlevels(θ) - 1) + ) + right_gradient_extrapolate = Operators.SetGradient(right_center_left_biased_grad) + left_gradient_extrapolate = Operators.SetGradient( + Geometry.Covariant3Vector.(Fields.level(θ, 2) .- left_center), + ) + 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) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- + Fields.level(θ, Fields.nlevels(θ) - 1) + ) + left_gradient = + Operators.SetGradient(@. 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) + right_center_left_biased_grad = + Geometry.Covariant3Vector.( + Fields.level(θ, Fields.nlevels(θ)) .- + Fields.level(θ, Fields.nlevels(θ) - 1) + ) + left_gradient = + Operators.SetGradient(@. 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() + left_gradient_extrapolate = Operators.SetGradient( + Geometry.Covariant3Vector.(Fields.level(θ, 2) .- left_center), + ) + 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, ) - 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ᴸ)), + sol2 = solve( + prob2, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, ) - 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ᴿᴮ)), + sol3 = solve( + prob3, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, + ) + sol4 = solve( + prob4, + SSPRK33(), + dt = Δt, + saveat = collect(0.0:(10 * Δt):t_end), + progress = true, + progress_message = (dt, u, p, t) -> t, ) - set_bcs = Operators.SetBoundaryOperator(; left = left_bc, right = right_bc) - UB = Operators.UpwindBiasedProductC2F() - ∂ = Operators.DivergenceF2C() - return @. dθ = -∂(set_bcs(UB(V, θ))) + fcc(V, θ) -end -# use the advection operator -function tendency3!(dθ, θ, _, t) - left_center = Fields.level(θ, 1) - right_center_left_biased_grad = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) - ) - left_gradient = - Operators.SetGradient(@. 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) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), + ENV["GKSwstype"] = "nul" + using ClimaCorePlots, Plots + 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"), ) - left_center = Fields.level(θ, 1) - right_center_left_biased_grad = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) - ) - left_gradient = - Operators.SetGradient(@. 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(θ))) + - 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/step.jl b/examples/column/step.jl deleted file mode 100644 index f8d847ac5c..0000000000 --- a/examples/column/step.jl +++ /dev/null @@ -1,254 +0,0 @@ -import ClimaComms -ClimaComms.@import_required_backends -using LinearAlgebra -import ClimaCore: - Fields, - Domains, - Topologies, - Meshes, - DataLayouts, - Operators, - Geometry, - Spaces, - Utilities - -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) - 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.z - 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.z - 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() - - return @. dθ = -∂(set_bcs(UB(V, θ))) -end -function tendency2!(dθ, θ, _, t) - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - 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.z - 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.z - 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() - return @. dθ = -∂(set_bcs(UB(V, θ))) + fcc(V, θ) -end -# use the advection operator -function tendency3!(dθ, θ, _, t) - left_center = Fields.level(θ, 1) - right_center_left_biased_grad = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) - ) - left_gradient = - Operators.SetGradient(@. 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) - - fcc = Operators.FluxCorrectionC2C( - left = Operators.Extrapolate(), - right = Operators.Extrapolate(), - ) - left_center = Fields.level(θ, 1) - right_center_left_biased_grad = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- Fields.level(θ, Fields.nlevels(θ) - 1) - ) - left_gradient = - Operators.SetGradient(@. 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(θ))) + - 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/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..ca18dacfc4 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,10 +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(), @@ -272,9 +269,10 @@ 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ₕ) + add_flux_correction_c2c(dρ, w, ρ) + add_flux_correction_c2c(dρθ, w, ρθ) + add_flux_correction_c2c(dρuₕ, w, ρuₕ) + @. dρw += fcf(wc, ρw) @. dρw += fcf(wc, ρw) end diff --git a/examples/hybrid/flux_correction_utils.jl b/examples/hybrid/flux_correction_utils.jl new file mode 100644 index 0000000000..295ba30675 --- /dev/null +++ b/examples/hybrid/flux_correction_utils.jl @@ -0,0 +1,56 @@ +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 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/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 c7e86b6b17..b94ebe5a0a 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -19,7 +19,6 @@ const TwoArgFDOperatorWithCenterInput = Union{ Operators.WeightedInterpolateC2F, Operators.UpwindBiasedProductC2F, Operators.Upwind3rdOrderBiasedProductC2F, - Operators.FluxCorrectionC2C, } const TwoArgFDOperatorWithFaceInput = Union{ Operators.WeightedInterpolateF2C, @@ -577,39 +576,6 @@ end op_matrix_row_type(::Operators.AdvectionOperator, ::Type{FT}, _) where {FT} = TridiagonalMatrixRow{FT} -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)) diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index b59f275dca..14c81dbd63 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -2116,104 +2116,6 @@ Base.@propagate_inbounds function stencil_right_boundary( 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, θ) diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index 25a9535896..af6bbb9e8d 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -27,7 +27,6 @@ import ClimaCore.Operators: Upwind3rdOrderBiasedProductC2F, FCTBorisBook, FCTZalesak, - FluxCorrectionC2C, FluxCorrectionF2F, SetBoundaryOperator, GradientC2F, @@ -152,7 +151,6 @@ end (ᶠuvw, ᶜscalar), 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) @@ -197,13 +195,11 @@ end ᶜlbias = LeftBiasedF2C() ᶠrbias = RightBiasedC2F(; set_nested_values.top) ᶜwinterp = WeightedInterpolateF2C() - ᶜflux_correct = FluxCorrectionC2C(; extrapolate...) ᶜdiv = DivergenceF2C() ᶠinterp_matrix = MatrixFields.operator_matrix(ᶠinterp) ᶜlbias_matrix = MatrixFields.operator_matrix(ᶜlbias) ᶠrbias_matrix = MatrixFields.operator_matrix(ᶠrbias) ᶜwinterp_matrix = MatrixFields.operator_matrix(ᶜwinterp) - ᶜflux_correct_matrix = MatrixFields.operator_matrix(ᶜflux_correct) ᶜdiv_matrix = MatrixFields.operator_matrix(ᶜdiv) @test_throws "does not contain any Fields" @. ᶜlbias_matrix() * @@ -229,15 +225,13 @@ end test_field_broadcast(; test_name = "product of six operator matrices", get_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() ), set_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() @@ -248,26 +242,21 @@ end test_name = "applying six operators to a nested field using operator \ matrices", get_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() * ᶜnested ), set_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * - ᶜwinterp_matrix(ᶠscalar) * + @. ᶜwinterp_matrix(ᶠscalar) * ᶠrbias_matrix() * ᶜlbias_matrix() * ᶠinterp_matrix() * ᶜnested ), ref_set_result = @lazy( - @. ᶜflux_correct( - ᶠuvw, - ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))), - ) + @. ᶜwinterp(ᶠscalar, ᶠrbias(ᶜlbias(ᶠinterp(ᶜnested)))) ), ) # this test is will fail because of incorrect results, not InvalidIRError @@ -275,7 +264,7 @@ 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) * ( + @. ( ( ᶜwinterp_matrix(ᶠscalar) * ( ᶠrbias_matrix() * @@ -285,7 +274,7 @@ end ) ), set_result = @lazy( - @. ᶜflux_correct_matrix(ᶠuvw) * ( + @. ( ( ᶜwinterp_matrix(ᶠscalar) * ( ᶠrbias_matrix() * @@ -295,10 +284,7 @@ end ) ), ref_set_result = @lazy( - @. ᶜflux_correct( - ᶠ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. diff --git a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl index 24086ec7bd..9399aa2afd 100644 --- a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl +++ b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl @@ -187,7 +187,6 @@ LeftBiased3rdOrderC2F # unused LeftBiased3rdOrderF2C # unused RightBiased3rdOrderC2F # unused RightBiased3rdOrderF2C # unused -FluxCorrectionC2C # unused at the moment FluxCorrectionF2F # only used in ClimaAtmos src/ ``` =# diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index c91022fbf6..72f74fb5b6 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -109,13 +109,6 @@ function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field) return UB.(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( @@ -212,11 +205,6 @@ end centers, ) - @test_opt opt_FluxCorrectionC2C_Extrapolate( - face_velocities, - centers, - ) - @test_opt opt_GradientC2F_SetValue(centers) @test_opt opt_GradientC2F_SetGradient(centers) diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index dc63a523c3..79022fb756 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -110,14 +110,6 @@ function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field) return UB.(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), @@ -243,13 +235,6 @@ end face_velocities, centers, ) - - - @test_opt opt_FluxCorrectionC2C_Extrapolate( - face_velocities, - centers, - ) - @test_opt opt_GradientC2F_SetValue(centers) @test_opt opt_GradientC2F_SetGradient(centers) From 9e6bf5380199df770b7bf31fabc0474b3ae2c80d Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 17 Jun 2026 15:36:43 -0700 Subject: [PATCH 15/18] fixes for rebase --- examples/hybrid/plane/density_current_2dinvariant_rhoe.jl | 4 ++-- examples/hybrid/plane/topo_agnesi_nh.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl index d8d0810319..7bdc25d858 100644 --- a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl +++ b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl @@ -257,9 +257,9 @@ function rhs_invariant!(dY, Y, _, t) ᶜ∇ᵥw_bottom = Fields.level(ᶜ∇ᵥw, 1) ᶜ∇ᵥw_top = Fields.level(ᶜ∇ᵥw, Fields.nlevels(ᶜ∇ᵥw)) bottom_divergence = @. Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * - (2 * lg_bottom_face.invJ) + (2 * inv(lg_bottom_face.J)) top_divergence = - @. Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * lg_top_face.invJ) + @. 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), diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl index 2f082c4344..235f50cc49 100644 --- a/examples/hybrid/plane/topo_agnesi_nh.jl +++ b/examples/hybrid/plane/topo_agnesi_nh.jl @@ -312,9 +312,9 @@ function rhs_invariant!(dY, Y, _, t) ᶜ∇ᵥw_bottom = Fields.level(ᶜ∇ᵥw, 1) ᶜ∇ᵥw_top = Fields.level(ᶜ∇ᵥw, Fields.nlevels(ᶜ∇ᵥw)) bottom_divergence = @. Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * - (2 * lg_bottom_face.invJ) + (2 * inv(lg_bottom_face.J)) top_divergence = - @. Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * lg_top_face.invJ) + @. 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), From 50aedaaafa939fbb5e0d80c1903f56b220e4bc00 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 17 Jun 2026 15:55:49 -0700 Subject: [PATCH 16/18] rm FluxCorrectionF2F --- examples/hybrid/box/bubble_3d_rhotheta.jl | 10 +- examples/hybrid/flux_correction_utils.jl | 44 ++++++++ src/MatrixFields/operator_matrices.jl | 1 - src/Operators/finitedifference.jl | 100 ------------------ test/MatrixFields/operator_matrices.jl | 2 - .../benchmark_stencils_climacore_kernels.jl | 1 - test/Operators/finitedifference/opt.jl | 13 --- test/Operators/hybrid/opt.jl | 13 --- 8 files changed, 48 insertions(+), 136 deletions(-) diff --git a/examples/hybrid/box/bubble_3d_rhotheta.jl b/examples/hybrid/box/bubble_3d_rhotheta.jl index ca18dacfc4..af80f02978 100644 --- a/examples/hybrid/box/bubble_3d_rhotheta.jl +++ b/examples/hybrid/box/bubble_3d_rhotheta.jl @@ -214,10 +214,6 @@ function rhs!(dY, Y, _, t) top = Operators.SetValue(Geometry.WVector(0.0)), ) - fcf = Operators.FluxCorrectionF2F( - bottom = Operators.Extrapolate(), - top = Operators.Extrapolate(), - ) uₕ = @. ρuₕ / ρ w = @. ρw / If(ρ) @@ -272,8 +268,10 @@ function rhs!(dY, Y, _, t) add_flux_correction_c2c(dρ, w, ρ) add_flux_correction_c2c(dρθ, w, ρθ) add_flux_correction_c2c(dρuₕ, w, ρuₕ) - @. dρw += fcf(wc, ρw) - @. dρw += fcf(wc, ρw) + add_flux_correction_f2f(dρw, wc, ρw) + add_flux_correction_f2f(dρw, wc, ρw) + # @. dρw += fcf(wc, ρw) + # @. dρw += fcf(wc, ρw) end ### DIFFUSION diff --git a/examples/hybrid/flux_correction_utils.jl b/examples/hybrid/flux_correction_utils.jl index 295ba30675..23aa0e2a99 100644 --- a/examples/hybrid/flux_correction_utils.jl +++ b/examples/hybrid/flux_correction_utils.jl @@ -13,6 +13,7 @@ import ClimaCore: 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)) @@ -54,3 +55,46 @@ function add_flux_correction_c2c(d_, velocity, quantity) ) * 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/src/MatrixFields/operator_matrices.jl b/src/MatrixFields/operator_matrices.jl index b94ebe5a0a..e445e2cdcf 100644 --- a/src/MatrixFields/operator_matrices.jl +++ b/src/MatrixFields/operator_matrices.jl @@ -22,7 +22,6 @@ const TwoArgFDOperatorWithCenterInput = Union{ } const TwoArgFDOperatorWithFaceInput = Union{ Operators.WeightedInterpolateF2C, - Operators.FluxCorrectionF2F, } const ErroneousFDOperator = Union{ Operators.LeftBiased3rdOrderC2F, diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference.jl index 14c81dbd63..1ba1bd2645 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference.jl @@ -2115,106 +2115,6 @@ Base.@propagate_inbounds function stencil_right_boundary( return Geometry.Contravariant3Vector(zero(eltype(eltype(A_field)))) 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 """ diff --git a/test/MatrixFields/operator_matrices.jl b/test/MatrixFields/operator_matrices.jl index af6bbb9e8d..bd426ce662 100644 --- a/test/MatrixFields/operator_matrices.jl +++ b/test/MatrixFields/operator_matrices.jl @@ -27,7 +27,6 @@ import ClimaCore.Operators: Upwind3rdOrderBiasedProductC2F, FCTBorisBook, FCTZalesak, - FluxCorrectionF2F, SetBoundaryOperator, GradientC2F, GradientF2C, @@ -151,7 +150,6 @@ end (ᶠuvw, ᶜscalar), true, ) - test_op_matrix(FluxCorrectionF2F, Extrapolate, (ᶜuvw, ᶠnested)) test_op_matrix(SetBoundaryOperator, SetValue, (ᶠnested,)) test_op_matrix(GradientC2F, Nothing, (ᶜscalar,), true) test_op_matrix(GradientC2F, SetGradient, (ᶜscalar,)) diff --git a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl index 9399aa2afd..ff74e0676e 100644 --- a/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl +++ b/test/Operators/finitedifference/benchmark_stencils_climacore_kernels.jl @@ -187,6 +187,5 @@ LeftBiased3rdOrderC2F # unused LeftBiased3rdOrderF2C # unused RightBiased3rdOrderC2F # unused RightBiased3rdOrderF2C # unused -FluxCorrectionF2F # only used in ClimaAtmos src/ ``` =# diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 72f74fb5b6..6803c8b88f 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -33,14 +33,6 @@ function opt_RightBiasedF2C(face_field) return RB.(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))) @@ -171,11 +163,6 @@ end @test_opt opt_LeftBiasedF2C(faces) @test_opt opt_RightBiasedF2C(faces) - @test_opt opt_FluxCorrectionF2F_Extrapolate( - center_velocities, - faces, - ) - @test_opt opt_GradientF2C(faces) @test_opt opt_DivergenceF2C(faces) diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 79022fb756..6270155b99 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -34,14 +34,6 @@ function opt_RightBiasedF2C(face_field) return RB.(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))) @@ -201,11 +193,6 @@ end @test_opt opt_LeftBiasedF2C(faces) @test_opt opt_RightBiasedF2C(faces) - @test_opt opt_FluxCorrectionF2F_Extrapolate( - center_velocities, - faces, - ) - @test_opt opt_GradientF2C(faces) @test_opt opt_GradientF2C(faces) @test_opt opt_DivergenceF2C(faces) From d329c3f9f05d11cea4eb40c192afc4faf30a3c82 Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Wed, 17 Jun 2026 16:01:25 -0700 Subject: [PATCH 17/18] Up news --- .buildkite/pipeline.yml | 19 ------------------- NEWS.md | 5 +++++ docs/src/examples.md | 3 +++ examples/column/advect.jl | 4 +++- examples/hybrid/box/bubble_3d_rhotheta.jl | 2 -- examples/hybrid/plane/topo_agnesi_nh.jl | 2 +- test/Operators/hybrid/convergence_2d.jl | 4 ++-- 7 files changed, 14 insertions(+), 25 deletions(-) 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/src/examples.md b/docs/src/examples.md index 65b720e756..ebdba0bdd5 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -40,6 +40,9 @@ This is discretized using the following This test case is set up in a 1D column domain ``z \in [0, 1]`` and discretized into a mesh of 10 elements. A homogeneous Dirichlet boundary condition is set at the bottom boundary, `bcs_bottom`, setting the temperature to 0. A Neumann boundary condition is applied to the top boundary, `bcs_top`, setting the temperature gradient to 1. +### Advection equation + +The 1D Column advection example in [`examples/column/advect.jl`](https://github.com/CliMA/ClimaCore.jl/blob/main/examples/column/advect.jl). #### Equations and Discretizations diff --git a/examples/column/advect.jl b/examples/column/advect.jl index 803a720e43..ac8f9ff3f0 100644 --- a/examples/column/advect.jl +++ b/examples/column/advect.jl @@ -12,6 +12,8 @@ import ClimaCore: Utilities using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 +using ClimaCorePlots +using Plots import Logging import TerminalLoggers @@ -228,7 +230,7 @@ for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step)) ) ENV["GKSwstype"] = "nul" - using ClimaCorePlots, Plots + Plots.GRBackend() sim_type = fn == "sin" ? "advect" : "advect_step_function" diff --git a/examples/hybrid/box/bubble_3d_rhotheta.jl b/examples/hybrid/box/bubble_3d_rhotheta.jl index af80f02978..ffd1fa254d 100644 --- a/examples/hybrid/box/bubble_3d_rhotheta.jl +++ b/examples/hybrid/box/bubble_3d_rhotheta.jl @@ -270,8 +270,6 @@ function rhs!(dY, Y, _, t) add_flux_correction_c2c(dρuₕ, w, ρuₕ) add_flux_correction_f2f(dρw, wc, ρw) add_flux_correction_f2f(dρw, wc, ρw) - # @. dρw += fcf(wc, ρw) - # @. dρw += fcf(wc, ρw) end ### DIFFUSION diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl index 235f50cc49..102516762b 100644 --- a/examples/hybrid/plane/topo_agnesi_nh.jl +++ b/examples/hybrid/plane/topo_agnesi_nh.jl @@ -319,7 +319,7 @@ function rhs_invariant!(dY, Y, _, t) bottom = Operators.SetValue(bottom_divergence), top = Operators.SetValue(top_divergence), ) - vκ₂∇²w = @. vdivc2f(κ₂ * ᶜ∇ᵥw) + vκ₂∇²w = @. set_bcs(vdivc2f(κ₂ * ᶜ∇ᵥw)) hκ₂∇²h_tot = @. hwdiv(cρ * κ₂ * ᶜ∇ₕh_tot) vκ₂∇²h_tot = @. vdivf2c(fρ * κ₂ * ᶠ∇ᵥh_tot) diff --git a/test/Operators/hybrid/convergence_2d.jl b/test/Operators/hybrid/convergence_2d.jl index f5376f62f5..bf015d8452 100644 --- a/test/Operators/hybrid/convergence_2d.jl +++ b/test/Operators/hybrid/convergence_2d.jl @@ -5,8 +5,8 @@ include("utils_2d.jl") function advection(c, f, hv_center_space) adv = zeros(eltype(f), hv_center_space) gradc2f = Operators.GradientC2F( - bottom = Operators.SetGradient(Geometry.WVector(FT(1))), - top = Operators.SetGradient(Geometry.WVector(FT(1))), + bottom = Operators.SetGradient(Geometry.WVector(1.0)), + top = Operators.SetGradient(Geometry.WVector(1.0)), ) interpf2c = Operators.InterpolateF2C() return @. adv = From 1c33445ce87237fabe78a632c94e5af6aaf4188b Mon Sep 17 00:00:00 2001 From: imreddyTeja Date: Thu, 18 Jun 2026 11:23:24 -0700 Subject: [PATCH 18/18] lazify bcs --- examples/column/advect.jl | 42 ++++++++++--------- examples/column/advect_diffusion.jl | 8 ++-- examples/column/ekman.jl | 19 ++++++--- .../plane/density_current_2dinvariant_rhoe.jl | 11 +++-- examples/hybrid/plane/topo_agnesi_nh.jl | 11 +++-- 5 files changed, 55 insertions(+), 36 deletions(-) diff --git a/examples/column/advect.jl b/examples/column/advect.jl index ac8f9ff3f0..aeee2e122f 100644 --- a/examples/column/advect.jl +++ b/examples/column/advect.jl @@ -11,6 +11,7 @@ import ClimaCore: Spaces, Utilities +import LazyBroadcast: lazy using OrdinaryDiffEqSSPRK: ODEProblem, solve, SSPRK33 using ClimaCorePlots using Plots @@ -121,15 +122,14 @@ for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step)) 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 = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- - Fields.level(θ, Fields.nlevels(θ) - 1) - ) + @. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1)) right_gradient_extrapolate = Operators.SetGradient(right_center_left_biased_grad) - left_gradient_extrapolate = Operators.SetGradient( - Geometry.Covariant3Vector.(Fields.level(θ, 2) .- left_center), - ) + θ_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, @@ -144,13 +144,14 @@ for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step)) # 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 = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- - Fields.level(θ, Fields.nlevels(θ) - 1) - ) + @. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1)) left_gradient = - Operators.SetGradient(@. Geometry.Covariant3Vector(2 * (left_center - sin(-t)))) + 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() @@ -160,20 +161,21 @@ for (fn, mesh) in zip(("sin", "step"), (mesh_sin, mesh_step)) # 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 = - Geometry.Covariant3Vector.( - Fields.level(θ, Fields.nlevels(θ)) .- - Fields.level(θ, Fields.nlevels(θ) - 1) - ) + @. lazy(Geometry.Covariant3Vector(θ_top - θ_top_m1)) left_gradient = - Operators.SetGradient(@. Geometry.Covariant3Vector(2 * (left_center - sin(-t)))) + 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() - left_gradient_extrapolate = Operators.SetGradient( - Geometry.Covariant3Vector.(Fields.level(θ, 2) .- left_center), - ) + θ_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, diff --git a/examples/column/advect_diffusion.jl b/examples/column/advect_diffusion.jl index fdab554620..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 @@ -63,10 +64,9 @@ function ∑tendencies!(dT, T, z, t) bc_gt = Operators.SetGradient( Geometry.WVector(FT(∇gaussian(z₁, t; ν = ν, δ = δ, 𝓌 = 𝓌, μ = μ))), ) - top_center_left_biased_grad = - Geometry.Covariant3Vector.( - Fields.level(T, Fields.nlevels(T)) .- Fields.level(T, Fields.nlevels(T) - 1) - ) + 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) diff --git a/examples/column/ekman.jl b/examples/column/ekman.jl index c6c9a3244e..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()) @@ -89,11 +90,14 @@ function tendency!(dY, Y, _, t) bcs_bottom = Operators.SetValue(Geometry.WVector(Cd * u_wind * u_1)) # Eq. 4.16 u_top = Fields.level(u, Fields.nlevels(u)) u_bottom = Fields.level(u, 1) - bcs_top = Operators.SetGradient(@. Geometry.Covariant3Vector(2 * (FT(ug) - u_top))) # Eq. 4.18 + 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) - bcs_bottom_advection = Operators.SetGradient(@. Geometry.Covariant3Vector(2 * u_bottom)) - bcs_top_advection = Operators.SetGradient(@. Geometry.Covariant3Vector(-2 * u_top)) + 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) @@ -105,11 +109,14 @@ function tendency!(dY, Y, _, t) bcs_bottom = Operators.SetValue(Geometry.WVector(Cd * u_wind * v_1)) # Eq. 4.17 v_top = Fields.level(v, Fields.nlevels(v)) v_bottom = Fields.level(v, 1) - bcs_top = Operators.SetGradient(@. Geometry.Covariant3Vector(2 * (FT(vg) - v_top))) # Eq. 4.19 + 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) - bcs_bottom_advection = Operators.SetGradient(@. Geometry.Covariant3Vector(2 * v_bottom)) - bcs_top_advection = Operators.SetGradient(@. Geometry.Covariant3Vector(-2 * v_top)) + 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 = diff --git a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl index 7bdc25d858..7fc4e18c4f 100644 --- a/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl +++ b/examples/hybrid/plane/density_current_2dinvariant_rhoe.jl @@ -19,6 +19,7 @@ import ClimaCore: Operators, Utilities using ClimaCore.Geometry +import LazyBroadcast: lazy using Logging: global_logger using TerminalLoggers: TerminalLogger @@ -256,10 +257,14 @@ function rhs_invariant!(dY, Y, _, t) 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 = @. Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * - (2 * inv(lg_bottom_face.J)) + bottom_divergence = @. lazy( + Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * + (2 * inv(lg_bottom_face.J)), + ) top_divergence = - @. Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * inv(lg_top_face.J)) + @. 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), diff --git a/examples/hybrid/plane/topo_agnesi_nh.jl b/examples/hybrid/plane/topo_agnesi_nh.jl index 102516762b..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 @@ -311,10 +312,14 @@ function rhs_invariant!(dY, Y, _, t) 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 = @. Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * - (2 * inv(lg_bottom_face.J)) + bottom_divergence = @. lazy( + Geometry.Jcontravariant3(ᶜ∇ᵥw_bottom, lg_bottom_center) * + (2 * inv(lg_bottom_face.J)), + ) top_divergence = - @. Geometry.Jcontravariant3(ᶜ∇ᵥw_top, lg_top_center) * (-2 * inv(lg_top_face.J)) + @. 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),