From 4e59637274065108140f8c4a8807dafefd89ede2 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Wed, 12 Jul 2023 14:24:57 -0700 Subject: [PATCH 1/2] non-product upwinding --- docs/src/operators.md | 16 +- examples/column/bb_fct_advection.jl | 8 +- examples/column/fct_advection.jl | 4 +- examples/column/zalesak_fct_advection.jl | 8 +- examples/hybrid/box/limiters_advection.jl | 4 +- .../hybrid/plane/bubble_2d_invariant_rhoe.jl | 4 +- examples/hybrid/plane/topo_schar_nh.jl | 4 +- examples/hybrid/sphere/deformation_flow.jl | 12 +- examples/hybrid/sphere/hadley_circulation.jl | 8 +- .../hybrid/staggered_nonhydrostatic_model.jl | 4 +- src/Operators/Operators.jl | 4 +- src/Operators/finitedifference/deprecated.jl | 534 ++++++++ .../finitedifference.jl | 1096 +++-------------- src/Operators/finitedifference/upwinding.jl | 582 +++++++++ src/Operators/operator2stencil.jl | 8 +- test/Operators/finitedifference/column.jl | 28 +- .../column_benchmark_kernels.jl | 12 +- .../column_benchmark_utils.jl | 36 +- .../implicit_stencils_utils.jl | 8 +- test/Operators/finitedifference/opt.jl | 12 +- .../finitedifference/opt_examples.jl | 10 +- .../finitedifference/opt_implicit_stencils.jl | 56 +- test/Operators/hybrid/opt.jl | 12 +- 23 files changed, 1388 insertions(+), 1082 deletions(-) create mode 100644 src/Operators/finitedifference/deprecated.jl rename src/Operators/{ => finitedifference}/finitedifference.jl (76%) create mode 100644 src/Operators/finitedifference/upwinding.jl diff --git a/docs/src/operators.md b/docs/src/operators.md index 217a46e460..73808e860b 100644 --- a/docs/src/operators.md +++ b/docs/src/operators.md @@ -65,14 +65,14 @@ InterpolateC2F InterpolateF2C WeightedInterpolateC2F WeightedInterpolateF2C -UpwindBiasedProductC2F -Upwind3rdOrderBiasedProductC2F +LeftBiased1stOrderC2F +RightBiased1stOrderC2F +LeftBiased3rdOrderC2F +RightBiased3rdOrderC2F +Upwind1stOrderC2F +Upwind3rdOrderC2F FCTBorisBook FCTZalesak -LeftBiasedC2F -RightBiasedC2F -LeftBiasedF2C -RightBiasedF2C ``` ### Derivative operators @@ -91,8 +91,6 @@ CurlC2F ```@docs SetBoundaryOperator -FirstOrderOneSided -ThirdOrderOneSided ``` ## Finite difference boundary conditions @@ -103,6 +101,8 @@ SetValue SetGradient SetDivergence Extrapolate +OneSided1stOrder +OneSided3rdOrder ``` ## Internal APIs diff --git a/examples/column/bb_fct_advection.jl b/examples/column/bb_fct_advection.jl index 746ed3792e..33af59f976 100644 --- a/examples/column/bb_fct_advection.jl +++ b/examples/column/bb_fct_advection.jl @@ -48,12 +48,12 @@ function tendency!(yₜ, y, parameters, t) top = Operators.Extrapolate(), ) upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) FCTBorisBook = Operators.FCTBorisBook( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) @. yₜ.q = -divf2c( diff --git a/examples/column/fct_advection.jl b/examples/column/fct_advection.jl index 8b047367f0..86495ed5a5 100644 --- a/examples/column/fct_advection.jl +++ b/examples/column/fct_advection.jl @@ -45,8 +45,8 @@ function f!(dydt, y, parameters, t) top = Operators.Extrapolate(), ) third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) divf2c = Operators.DivergenceF2C( bottom = Operators.SetValue(Geometry.WVector(FT(0.0))), diff --git a/examples/column/zalesak_fct_advection.jl b/examples/column/zalesak_fct_advection.jl index efc20fbf44..5ec508c8ac 100644 --- a/examples/column/zalesak_fct_advection.jl +++ b/examples/column/zalesak_fct_advection.jl @@ -48,12 +48,12 @@ function tendency!(yₜ, y, parameters, t) top = Operators.Extrapolate(), ) upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) FCTZalesak = Operators.FCTZalesak( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) @. yₜ.q = -divf2c( diff --git a/examples/hybrid/box/limiters_advection.jl b/examples/hybrid/box/limiters_advection.jl index daaa8155a7..046ed58cd2 100644 --- a/examples/hybrid/box/limiters_advection.jl +++ b/examples/hybrid/box/limiters_advection.jl @@ -179,8 +179,8 @@ function vertical_tendency!(yₜ, y, cache, t) bottom = Operators.SetValue(Geometry.Contravariant3Vector(FT(0))), ) upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) ax12 = (Geometry.Covariant12Axis(),) ax3 = (Geometry.Covariant3Axis(),) diff --git a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl index 6b3d1107fe..70a68337c4 100644 --- a/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl +++ b/examples/hybrid/plane/bubble_2d_invariant_rhoe.jl @@ -172,8 +172,8 @@ function rhs_invariant!(dY, Y, _, t) ) # 1.c) vertical upwinding third_order_upwind_c2f = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) # 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 diff --git a/examples/hybrid/plane/topo_schar_nh.jl b/examples/hybrid/plane/topo_schar_nh.jl index d4277ef37a..64181f523e 100644 --- a/examples/hybrid/plane/topo_schar_nh.jl +++ b/examples/hybrid/plane/topo_schar_nh.jl @@ -230,8 +230,8 @@ function rhs_invariant!(dY, Y, _, t) f_upwind_product1 = Operators.UpwindBiasedProductC2F() f_upwind_product3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) dρ .= 0 .* cρ diff --git a/examples/hybrid/sphere/deformation_flow.jl b/examples/hybrid/sphere/deformation_flow.jl index db995ae2f8..3697ba667e 100644 --- a/examples/hybrid/sphere/deformation_flow.jl +++ b/examples/hybrid/sphere/deformation_flow.jl @@ -67,16 +67,16 @@ const upwind1 = Operators.UpwindBiasedProductC2F( top = Operators.Extrapolate(), ) const upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) const FCTZalesak = Operators.FCTZalesak( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) const FCTBorisBook = Operators.FCTBorisBook( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) # Reference pressure and density diff --git a/examples/hybrid/sphere/hadley_circulation.jl b/examples/hybrid/sphere/hadley_circulation.jl index 949b2bf9c5..597ad6c4bb 100644 --- a/examples/hybrid/sphere/hadley_circulation.jl +++ b/examples/hybrid/sphere/hadley_circulation.jl @@ -94,12 +94,12 @@ function tendency!(yₜ, y, parameters, t) top = Operators.Extrapolate(), ) upwind3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) FCTZalesak = Operators.FCTZalesak( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) hdiv = Operators.Divergence() hwdiv = Operators.WeakDivergence() diff --git a/examples/hybrid/staggered_nonhydrostatic_model.jl b/examples/hybrid/staggered_nonhydrostatic_model.jl index 302de7696f..1879104818 100644 --- a/examples/hybrid/staggered_nonhydrostatic_model.jl +++ b/examples/hybrid/staggered_nonhydrostatic_model.jl @@ -52,8 +52,8 @@ const ᶜFC = Operators.FluxCorrectionC2C( ) const ᶠupwind_product1 = Operators.UpwindBiasedProductC2F() const ᶠupwind_product3 = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) const ᶜinterp_stencil = Operators.Operator2Stencil(ᶜinterp) diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index 973d59b842..0066f87b3d 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -20,7 +20,9 @@ using ..RecursiveApply include("common.jl") include("spectralelement.jl") include("numericalflux.jl") -include("finitedifference.jl") +include("finitedifference/finitedifference.jl") +include("finitedifference/upwinding.jl") +include("finitedifference/deprecated.jl") include("stencilcoefs.jl") include("operator2stencil.jl") include("pointwisestencil.jl") diff --git a/src/Operators/finitedifference/deprecated.jl b/src/Operators/finitedifference/deprecated.jl new file mode 100644 index 0000000000..470e0ffa3f --- /dev/null +++ b/src/Operators/finitedifference/deprecated.jl @@ -0,0 +1,534 @@ +# to be deprecated in future + +# these were renamed +const LeftBiasedC2F = LeftBiased1stOrderC2F +const RightBiasedC2F = RightBiased1stOrderC2F +const FirstOrderOneSided = OneSided1stOrder +const ThirdOrderOneSided = OneSided3rdOrder + +# F2C aren't used + +""" + L = LeftBiasedF2C(;boundaries) + L.(x) + +Interpolate a face-value field to a center-valued field from the left. +```math +L(x)[i+\\tfrac{1}{2}] = x[i] +``` + +Only the left boundary condition should be set. Currently supported is: +- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. +```math +L(x)[1] = x_0 +``` +""" +struct LeftBiasedF2C{BCS} <: InterpolationOperator + bcs::BCS +end +LeftBiasedF2C(; kwargs...) = LeftBiasedF2C(NamedTuple(kwargs)) + +return_space(::LeftBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) = + Spaces.CenterFiniteDifferenceSpace(space) +return_space(::LeftBiasedF2C, space::Spaces.FaceExtrudedFiniteDifferenceSpace) = + Spaces.CenterExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::LeftBiasedF2C, arg) = ((-half, -half),) +Base.@propagate_inbounds stencil_interior( + ::LeftBiasedF2C, + loc, + space, + idx, + hidx, + arg, +) = getidx(space, arg, loc, idx - half, hidx) +left_interior_idx( + space::AbstractSpace, + ::LeftBiasedF2C, + ::AbstractBoundaryCondition, + arg, +) = left_idx(space) +right_interior_idx( + space::AbstractSpace, + ::LeftBiasedF2C, + ::AbstractBoundaryCondition, + arg, +) = right_idx(space) + +left_interior_idx(space::AbstractSpace, ::LeftBiasedF2C, ::SetValue, arg) = + left_idx(space) + 1 +Base.@propagate_inbounds function stencil_left_boundary( + ::LeftBiasedF2C, + bc::SetValue, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == left_center_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end + + + + +""" + L = LeftBiased3rdOrderF2C(;boundaries) + L.(x) + +Interpolate a face-value field to a center-valued field from the left, using a 3rd-order reconstruction. +```math +L(x)[i+\\tfrac{1}{2}] = \\left(-2 x[i-1] + 10 x[i] + 4 x[i+1] \\right) / 12 +``` + +Only the left boundary condition should be set. Currently supported is: +- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. +```math +L(x)[1] = x_0 +``` +""" +struct LeftBiased3rdOrderF2C{BCS} <: InterpolationOperator + bcs::BCS +end +LeftBiased3rdOrderF2C(; kwargs...) = LeftBiased3rdOrderF2C(NamedTuple(kwargs)) + +return_space(::LeftBiased3rdOrderF2C, space::Spaces.FaceFiniteDifferenceSpace) = + Spaces.CenterFiniteDifferenceSpace(space) +return_space( + ::LeftBiased3rdOrderF2C, + space::Spaces.FaceExtrudedFiniteDifferenceSpace, +) = Spaces.CenterExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::LeftBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),) +Base.@propagate_inbounds stencil_interior( + ::LeftBiased3rdOrderF2C, + loc, + space, + idx, + hidx, + arg, +) = + ( + -2 * getidx(space, arg, loc, idx - 1 - half, hidx) + + 10 * getidx(space, arg, loc, idx - half, hidx) + + 4 * getidx(space, arg, loc, idx + half, hidx) + ) / 12 + +left_interior_idx( + space::AbstractSpace, + ::LeftBiased3rdOrderF2C, + ::AbstractBoundaryCondition, + arg, +) = left_idx(space) + 1 +right_interior_idx( + space::AbstractSpace, + ::LeftBiased3rdOrderF2C, + ::AbstractBoundaryCondition, + arg, +) = right_idx(space) + +Base.@propagate_inbounds function stencil_left_boundary( + ::LeftBiased3rdOrderF2C, + bc::SetValue, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == left_center_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end + + +""" + R = RightBiasedF2C(;boundaries) + R.(x) + +Interpolate a face-valued field to a center-valued field from the right. +```math +R(x)[i] = x[i+\\tfrac{1}{2}] +``` + +Only the right boundary condition should be set. Currently supported is: +- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. +```math +R(x)[n+\\tfrac{1}{2}] = x_0 +``` +""" +struct RightBiasedF2C{BCS} <: InterpolationOperator + bcs::BCS +end +RightBiasedF2C(; kwargs...) = RightBiasedF2C(NamedTuple(kwargs)) + +return_space(::RightBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) = + Spaces.CenterFiniteDifferenceSpace(space) +return_space( + ::RightBiasedF2C, + space::Spaces.FaceExtrudedFiniteDifferenceSpace, +) = Spaces.CenterExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::RightBiasedF2C, arg) = ((half, half),) +Base.@propagate_inbounds stencil_interior( + ::RightBiasedF2C, + loc, + space, + idx, + hidx, + arg, +) = getidx(space, arg, loc, idx + half, hidx) + +left_interior_idx( + space::AbstractSpace, + ::RightBiasedF2C, + ::AbstractBoundaryCondition, + arg, +) = left_idx(space) +right_interior_idx( + space::AbstractSpace, + ::RightBiasedF2C, + ::AbstractBoundaryCondition, + arg, +) = right_idx(space) + +right_interior_idx(space::AbstractSpace, ::RightBiasedF2C, ::SetValue, arg) = + right_idx(space) - 1 +Base.@propagate_inbounds function stencil_right_boundary( + ::RightBiasedF2C, + bc::SetValue, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == right_center_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end + + +""" + R = RightBiased3rdOrderF2C(;boundaries) + R.(x) + +Interpolate a face-valued field to a center-valued field from the right, using a 3rd-order reconstruction. +```math +R(x)[i] = \\left(4 x[i] + 10 x[i+1] -2 x[i+2] \\right) / 12 +``` + +Only the right boundary condition should be set. Currently supported is: +- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. +```math +R(x)[n+\\tfrac{1}{2}] = x_0 +``` +""" +struct RightBiased3rdOrderF2C{BCS} <: InterpolationOperator + bcs::BCS +end +RightBiased3rdOrderF2C(; kwargs...) = RightBiased3rdOrderF2C(NamedTuple(kwargs)) + +return_space( + ::RightBiased3rdOrderF2C, + space::Spaces.FaceFiniteDifferenceSpace, +) = Spaces.CenterFiniteDifferenceSpace(space) +return_space( + ::RightBiased3rdOrderF2C, + space::Spaces.FaceExtrudedFiniteDifferenceSpace, +) = Spaces.CenterExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::RightBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),) +Base.@propagate_inbounds stencil_interior( + ::RightBiased3rdOrderF2C, + loc, + space, + idx, + hidx, + arg, +) = + ( + 4 * getidx(space, arg, loc, idx - half, hidx) + + 10 * getidx(space, arg, loc, idx + half, hidx) - + 2 * getidx(space, arg, loc, idx + half + 1, hidx) + ) / 12 + +boundary_width(::RightBiased3rdOrderF2C, ::SetValue) = 1 +Base.@propagate_inbounds function stencil_right_boundary( + ::RightBiased3rdOrderF2C, + bc::SetValue, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == right_center_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end + + +""" + U = UpwindBiasedProductC2F(;boundaries) + U.(v, x) + +Compute the product of the face-valued vector field `v` and a center-valued +field `x` at cell faces by upwinding `x` according to the direction of `v`. + +More precisely, it is computed based on the sign of the 3rd contravariant +component, and it returns a `Contravariant3Vector`: +```math +U(\\boldsymbol{v},x)[i] = \\begin{cases} + v^3[i] x[i-\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] > 0 \\\\ + v^3[i] x[i+\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] < 0 + \\end{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 + U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = U(\\boldsymbol{v},x)[1 + \\tfrac{1}{2}] + ``` +""" +struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator + bcs::BCS +end +UpwindBiasedProductC2F(; kwargs...) = UpwindBiasedProductC2F(NamedTuple(kwargs)) + +return_eltype(::UpwindBiasedProductC2F, V, A) = + Geometry.Contravariant3Vector{eltype(eltype(V))} + +return_space( + ::UpwindBiasedProductC2F, + velocity_space::Spaces.FaceFiniteDifferenceSpace, + arg_space::Spaces.CenterFiniteDifferenceSpace, +) = velocity_space +return_space( + ::UpwindBiasedProductC2F, + velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace, + arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = velocity_space + + +stencil_interior_width(::UpwindBiasedProductC2F, velocity, arg) = + ((0, 0), (-half, half)) + +Base.@propagate_inbounds function stencil_interior( + ::UpwindBiasedProductC2F, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + u = stencil_interior( + Upwind1stOrderC2F(), + loc, + space, + idx, + hidx, + Geometry.Contravariant3Vector(v³), + arg, + ) + return Geometry.Contravariant3Vector(v³ * u) +end + +boundary_width(::UpwindBiasedProductC2F, ::AbstractBoundaryCondition) = 1 + +Base.@propagate_inbounds function stencil_left_boundary( + ::UpwindBiasedProductC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + u = stencil_left_boundary( + Upwind1stOrderC2F(), + bc, + loc, + space, + idx, + hidx, + Geometry.Contravariant3Vector(v³), + arg, + ) + return Geometry.Contravariant3Vector(v³ * u) +end + +Base.@propagate_inbounds function stencil_right_boundary( + ::UpwindBiasedProductC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + u = stencil_right_boundary( + Upwind1stOrderC2F(), + bc, + loc, + space, + idx, + hidx, + Geometry.Contravariant3Vector(v³), + arg, + ) + return Geometry.Contravariant3Vector(v³ * u) +end + +""" + U = Upwind3rdOrderBiasedProductC2F(;boundaries) + U.(v, x) + +Compute the product of a face-valued vector field `v` and a center-valued field +`x` at cell faces by upwinding `x`, to third-order of accuracy, according to `v` +```math +U(v,x)[i] = \\begin{cases} + v[i] \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12 \\textrm{, if } v[i] > 0 \\\\ + v[i] \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}] \\right) / 12 \\textrm{, if } v[i] < 0 + \\end{cases} +``` +This stencil is based on [WickerSkamarock2002](@cite), eq. 4(a). + +Supported boundary conditions are: +- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind scheme to compute `x` on the left boundary, + and the first-order upwind scheme to compute `x` on the right boundary. +and the third-order upwind reconstruction to compute `x` on the right boundary. + +!!! note + These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a [`DivergenceF2C`](@ref) operator, with a [`SetValue`](@ref) boundary. +""" +struct Upwind3rdOrderBiasedProductC2F{BCS} <: AdvectionOperator + bcs::BCS +end +Upwind3rdOrderBiasedProductC2F(; kwargs...) = + Upwind3rdOrderBiasedProductC2F(NamedTuple(kwargs)) + +return_eltype(::Upwind3rdOrderBiasedProductC2F, V, A) = + Geometry.Contravariant3Vector{eltype(eltype(V))} + +return_space( + ::Upwind3rdOrderBiasedProductC2F, + velocity_space::Spaces.FaceFiniteDifferenceSpace, + arg_space::Spaces.CenterFiniteDifferenceSpace, +) = velocity_space +return_space( + ::Upwind3rdOrderBiasedProductC2F, + velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace, + arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = velocity_space + + +stencil_interior_width(::Upwind3rdOrderBiasedProductC2F, velocity, arg) = + ((0, 0), (-half - 1, half + 1)) + +Base.@propagate_inbounds function stencil_interior( + ::Upwind3rdOrderBiasedProductC2F, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + u = stencil_interior( + Upwind3rdOrderC2F(), + loc, + space, + idx, + hidx, + Geometry.Contravariant3Vector(v³), + arg, + ) + return Geometry.Contravariant3Vector(v³ * u) +end + +boundary_width(::Upwind3rdOrderBiasedProductC2F, ::AbstractBoundaryCondition) = + 2 + + +Base.@propagate_inbounds function stencil_left_boundary( + ::Upwind3rdOrderBiasedProductC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + u = stencil_left_boundary( + Upwind3rdOrderC2F(), + bc, + loc, + space, + idx, + hidx, + Geometry.Contravariant3Vector(v³), + arg, + ) + return Geometry.Contravariant3Vector(v³ * u) +end + +Base.@propagate_inbounds function stencil_right_boundary( + ::Upwind3rdOrderBiasedProductC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + u = stencil_right_boundary( + Upwind3rdOrderC2F(), + bc, + loc, + space, + idx, + hidx, + Geometry.Contravariant3Vector(v³), + arg, + ) + return Geometry.Contravariant3Vector(v³ * u) +end diff --git a/src/Operators/finitedifference.jl b/src/Operators/finitedifference/finitedifference.jl similarity index 76% rename from src/Operators/finitedifference.jl rename to src/Operators/finitedifference/finitedifference.jl index a15cb9cf7e..4fc428a298 100644 --- a/src/Operators/finitedifference.jl +++ b/src/Operators/finitedifference/finitedifference.jl @@ -136,18 +136,19 @@ Set the value at the boundary to be the same as the closest interior point. struct Extrapolate <: AbstractBoundaryCondition end """ - FirstOrderOneSided() + OneSided1stOrder() Use a first-order up/down-wind scheme to compute the value at the boundary. """ -struct FirstOrderOneSided <: AbstractBoundaryCondition end +struct OneSided1stOrder <: AbstractBoundaryCondition end """ - ThirdOrderOneSided() + OneSided3rdOrder() Use a third-order up/down-wind scheme to compute the value at the boundary. """ -struct ThirdOrderOneSided <: AbstractBoundaryCondition end +struct OneSided3rdOrder <: AbstractBoundaryCondition end + abstract type Location end abstract type Boundary <: Location end @@ -514,1038 +515,225 @@ Base.@propagate_inbounds function stencil_right_boundary( a⁻ end -""" - L = LeftBiasedC2F(;boundaries) - L.(x) -Interpolate a center-value field to a face-valued field from the left. -```math -L(x)[i] = x[i-\\tfrac{1}{2}] -``` -Only the left boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -L(x)[\\tfrac{1}{2}] = x_0 -``` -""" -struct LeftBiasedC2F{BCS} <: InterpolationOperator - bcs::BCS -end -LeftBiasedC2F(; kwargs...) = LeftBiasedC2F(NamedTuple(kwargs)) -return_space(::LeftBiasedC2F, space::Spaces.CenterFiniteDifferenceSpace) = - Spaces.FaceFiniteDifferenceSpace(space) -return_space( - ::LeftBiasedC2F, - space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) -stencil_interior_width(::LeftBiasedC2F, arg) = ((-half, -half),) -Base.@propagate_inbounds stencil_interior( - ::LeftBiasedC2F, - loc, - space, - idx, - hidx, - arg, -) = getidx(space, arg, loc, idx - half, hidx) -left_interior_idx( - space::AbstractSpace, - ::LeftBiasedC2F, - ::AbstractBoundaryCondition, - arg, -) = left_idx(space) + 1 -right_interior_idx( - space::AbstractSpace, - ::LeftBiasedC2F, - ::AbstractBoundaryCondition, - arg, -) = right_idx(space) -Base.@propagate_inbounds function stencil_left_boundary( - ::LeftBiasedC2F, - bc::SetValue, - loc, - space, - idx, - hidx, - arg, -) - @assert idx == left_face_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end +abstract type WeightedInterpolationOperator <: InterpolationOperator end +# TODO: this is not in general correct and the return type +# should be based on the component operator types (rdiv, rmul) but we don't have a good way +# of creating ex. one(field_type) for complex fields for inference +return_eltype(::WeightedInterpolationOperator, weights, arg) = eltype(arg) """ - L = LeftBiasedF2C(;boundaries) - L.(x) + WI = WeightedInterpolateF2C(; boundaries) + WI.(w, x) -Interpolate a face-value field to a center-valued field from the left. +Interpolate a face-valued field `x` to centers, weighted by a face-valued field +`w`, using the stencil ```math -L(x)[i+\\tfrac{1}{2}] = x[i] +WI(w, x)[i] = \\frac{ + w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}]) + }{ + w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] + } ``` -Only the left boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -L(x)[1] = x_0 -``` +No boundary conditions are required (or supported) """ -struct LeftBiasedF2C{BCS} <: InterpolationOperator +struct WeightedInterpolateF2C{BCS} <: WeightedInterpolationOperator bcs::BCS end -LeftBiasedF2C(; kwargs...) = LeftBiasedF2C(NamedTuple(kwargs)) - -return_space(::LeftBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) = - Spaces.CenterFiniteDifferenceSpace(space) -return_space(::LeftBiasedF2C, space::Spaces.FaceExtrudedFiniteDifferenceSpace) = - Spaces.CenterExtrudedFiniteDifferenceSpace(space) +WeightedInterpolateF2C(; kwargs...) = WeightedInterpolateF2C(NamedTuple(kwargs)) -stencil_interior_width(::LeftBiasedF2C, arg) = ((-half, -half),) -Base.@propagate_inbounds stencil_interior( - ::LeftBiasedF2C, - loc, - space, - idx, - hidx, - arg, -) = getidx(space, arg, loc, idx - half, hidx) -left_interior_idx( - space::AbstractSpace, - ::LeftBiasedF2C, - ::AbstractBoundaryCondition, - arg, -) = left_idx(space) -right_interior_idx( - space::AbstractSpace, - ::LeftBiasedF2C, - ::AbstractBoundaryCondition, - arg, -) = right_idx(space) +return_space( + ::WeightedInterpolateF2C, + weight_space::Spaces.FaceFiniteDifferenceSpace, + arg_space::Spaces.FaceFiniteDifferenceSpace, +) = Spaces.CenterFiniteDifferenceSpace(arg_space) +return_space( + ::WeightedInterpolateF2C, + weight_space::Spaces.FaceExtrudedFiniteDifferenceSpace, + arg_space::Spaces.FaceExtrudedFiniteDifferenceSpace, +) = Spaces.CenterExtrudedFiniteDifferenceSpace(arg_space) -left_interior_idx(space::AbstractSpace, ::LeftBiasedF2C, ::SetValue, arg) = - left_idx(space) + 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::LeftBiasedF2C, - bc::SetValue, +stencil_interior_width(::WeightedInterpolateF2C, weight, arg) = + ((-half, half), (-half, half)) +Base.@propagate_inbounds function stencil_interior( + ::WeightedInterpolateF2C, loc, space, idx, hidx, + weight, arg, ) - @assert idx == left_center_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) + w⁺ = getidx(space, weight, loc, idx + half, hidx) + w⁻ = getidx(space, weight, loc, idx - half, hidx) + a⁺ = getidx(space, arg, loc, idx + half, hidx) + a⁻ = getidx(space, arg, loc, idx - half, hidx) + RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻)) end +boundary_width(::WeightedInterpolateF2C, ::AbstractBoundaryCondition) = 0 + """ - L = LeftBiased3rdOrderC2F(;boundaries) - L.(x) + WI = WeightedInterpolateC2F(; boundaries) + WI.(w, x) -Interpolate a center-value field to a face-valued field from the left, using a 3rd-order reconstruction. +Interpolate a center-valued field `x` to faces, weighted by a center-valued field +`w`, using the stencil ```math -L(x)[i] = \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12 +WI(w, x)[i] = \\frac{ + w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}]) +}{ + w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] +} ``` -Only the left boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -L(x)[\\tfrac{1}{2}] = x_0 -``` +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). """ -struct LeftBiased3rdOrderC2F{BCS} <: InterpolationOperator +struct WeightedInterpolateC2F{BCS} <: WeightedInterpolationOperator bcs::BCS end -LeftBiased3rdOrderC2F(; kwargs...) = LeftBiased3rdOrderC2F(NamedTuple(kwargs)) +WeightedInterpolateC2F(; kwargs...) = WeightedInterpolateC2F(NamedTuple(kwargs)) return_space( - ::LeftBiased3rdOrderC2F, - space::Spaces.CenterFiniteDifferenceSpace, -) = Spaces.FaceFiniteDifferenceSpace(space) + ::WeightedInterpolateC2F, + weight_space::Spaces.CenterFiniteDifferenceSpace, + arg_space::Spaces.CenterFiniteDifferenceSpace, +) = Spaces.FaceFiniteDifferenceSpace(arg_space) return_space( - ::LeftBiased3rdOrderC2F, - space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) + ::WeightedInterpolateC2F, + weight_space::Spaces.CenterExtrudedFiniteDifferenceSpace, + arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = Spaces.FaceExtrudedFiniteDifferenceSpace(arg_space) -stencil_interior_width(::LeftBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),) -Base.@propagate_inbounds stencil_interior( - ::LeftBiased3rdOrderC2F, +stencil_interior_width(::WeightedInterpolateC2F, weight, arg) = + ((-half, half), (-half, half)) +Base.@propagate_inbounds function stencil_interior( + ::WeightedInterpolateC2F, loc, space, idx, hidx, + weight, arg, -) = - ( - -2 * getidx(space, arg, loc, idx - 1 - half, hidx) + - 10 * getidx(space, arg, loc, idx - half, hidx) + - 4 * getidx(space, arg, loc, idx + half, hidx) - ) / 12 - -left_interior_idx( - space::AbstractSpace, - ::LeftBiased3rdOrderC2F, - ::AbstractBoundaryCondition, - arg, -) = left_idx(space) + 2 -right_interior_idx( - space::AbstractSpace, - ::LeftBiased3rdOrderC2F, - ::AbstractBoundaryCondition, - arg, -) = right_idx(space) - 1 +) + w⁺ = getidx(space, weight, loc, idx + half, hidx) + w⁻ = getidx(space, weight, loc, idx - half, hidx) + a⁺ = getidx(space, arg, loc, idx + half, hidx) + a⁻ = getidx(space, arg, loc, idx - half, hidx) + RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻)) +end +boundary_width(::WeightedInterpolateC2F, ::AbstractBoundaryCondition) = 1 Base.@propagate_inbounds function stencil_left_boundary( - ::LeftBiased3rdOrderC2F, + ::WeightedInterpolateC2F, bc::SetValue, loc, space, idx, hidx, + weight, arg, ) @assert idx == left_face_boundary_idx(space) getidx(space, bc.val, loc, nothing, hidx) end - -""" - L = LeftBiased3rdOrderF2C(;boundaries) - L.(x) - -Interpolate a face-value field to a center-valued field from the left, using a 3rd-order reconstruction. -```math -L(x)[i+\\tfrac{1}{2}] = \\left(-2 x[i-1] + 10 x[i] + 4 x[i+1] \\right) / 12 -``` - -Only the left boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -L(x)[1] = x_0 -``` -""" -struct LeftBiased3rdOrderF2C{BCS} <: InterpolationOperator - bcs::BCS -end -LeftBiased3rdOrderF2C(; kwargs...) = LeftBiased3rdOrderF2C(NamedTuple(kwargs)) - -return_space(::LeftBiased3rdOrderF2C, space::Spaces.FaceFiniteDifferenceSpace) = - Spaces.CenterFiniteDifferenceSpace(space) -return_space( - ::LeftBiased3rdOrderF2C, - space::Spaces.FaceExtrudedFiniteDifferenceSpace, -) = Spaces.CenterExtrudedFiniteDifferenceSpace(space) - -stencil_interior_width(::LeftBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),) -Base.@propagate_inbounds stencil_interior( - ::LeftBiased3rdOrderF2C, +Base.@propagate_inbounds function stencil_right_boundary( + ::WeightedInterpolateC2F, + bc::SetValue, loc, space, idx, hidx, + weight, arg, -) = - ( - -2 * getidx(space, arg, loc, idx - 1 - half, hidx) + - 10 * getidx(space, arg, loc, idx - half, hidx) + - 4 * getidx(space, arg, loc, idx + half, hidx) - ) / 12 - -left_interior_idx( - space::AbstractSpace, - ::LeftBiased3rdOrderF2C, - ::AbstractBoundaryCondition, - arg, -) = left_idx(space) + 1 -right_interior_idx( - space::AbstractSpace, - ::LeftBiased3rdOrderF2C, - ::AbstractBoundaryCondition, - arg, -) = right_idx(space) +) + @assert idx == right_face_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end Base.@propagate_inbounds function stencil_left_boundary( - ::LeftBiased3rdOrderF2C, - bc::SetValue, + ::WeightedInterpolateC2F, + bc::SetGradient, loc, space, idx, hidx, + weight, arg, ) - @assert idx == left_center_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end - -""" - R = RightBiasedC2F(;boundaries) - R.(x) - -Interpolate a center-valued field to a face-valued field from the right. -```math -R(x)[i] = x[i+\\tfrac{1}{2}] -``` - -Only the right boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -R(x)[n+\\tfrac{1}{2}] = x_0 -``` -""" -struct RightBiasedC2F{BCS} <: InterpolationOperator - bcs::BCS + @assert idx == left_face_boundary_idx(space) + a⁺ = getidx(space, arg, loc, idx + half, hidx) + v₃ = Geometry.covariant3( + getidx(space, bc.val, loc, nothing, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ ⊟ RecursiveApply.rdiv(v₃, 2) end -RightBiasedC2F(; kwargs...) = RightBiasedC2F(NamedTuple(kwargs)) - -return_space(::RightBiasedC2F, space::Spaces.CenterFiniteDifferenceSpace) = - Spaces.FaceFiniteDifferenceSpace(space) -return_space( - ::RightBiasedC2F, - space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) - -stencil_interior_width(::RightBiasedC2F, arg) = ((half, half),) -Base.@propagate_inbounds stencil_interior( - ::RightBiasedC2F, +Base.@propagate_inbounds function stencil_right_boundary( + ::WeightedInterpolateC2F, + bc::SetGradient, loc, space, idx, hidx, + weight, arg, -) = getidx(space, arg, loc, idx + half, hidx) +) + @assert idx == right_face_boundary_idx(space) + a⁻ = getidx(space, arg, loc, idx - half, hidx) + v₃ = Geometry.covariant3( + getidx(space, bc.val, loc, nothing, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁻ ⊞ RecursiveApply.rdiv(v₃, 2) +end -left_interior_idx( - space::AbstractSpace, - ::RightBiasedC2F, - ::AbstractBoundaryCondition, - arg, -) = left_idx(space) -right_interior_idx( - space::AbstractSpace, - ::RightBiasedC2F, - ::AbstractBoundaryCondition, +Base.@propagate_inbounds function stencil_left_boundary( + ::WeightedInterpolateC2F, + bc::Extrapolate, + loc, + space, + idx, + hidx, + weight, arg, -) = right_idx(space) - 1 - +) + @assert idx == left_face_boundary_idx(space) + a⁺ = getidx(space, arg, loc, idx + half, hidx) + a⁺ +end Base.@propagate_inbounds function stencil_right_boundary( - ::RightBiasedC2F, - bc::SetValue, + ::WeightedInterpolateC2F, + bc::Extrapolate, loc, space, idx, hidx, + weight, arg, ) @assert idx == right_face_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) + a⁻ = getidx(space, arg, loc, idx - half, hidx) + a⁻ end -""" - R = RightBiasedF2C(;boundaries) - R.(x) - -Interpolate a face-valued field to a center-valued field from the right. -```math -R(x)[i] = x[i+\\tfrac{1}{2}] -``` -Only the right boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -R(x)[n+\\tfrac{1}{2}] = x_0 -``` -""" -struct RightBiasedF2C{BCS} <: InterpolationOperator - bcs::BCS -end -RightBiasedF2C(; kwargs...) = RightBiasedF2C(NamedTuple(kwargs)) -return_space(::RightBiasedF2C, space::Spaces.FaceFiniteDifferenceSpace) = - Spaces.CenterFiniteDifferenceSpace(space) -return_space( - ::RightBiasedF2C, - space::Spaces.FaceExtrudedFiniteDifferenceSpace, -) = Spaces.CenterExtrudedFiniteDifferenceSpace(space) +abstract type AdvectionOperator <: FiniteDifferenceOperator end +return_eltype(::AdvectionOperator, velocity, arg) = eltype(arg) -stencil_interior_width(::RightBiasedF2C, arg) = ((half, half),) -Base.@propagate_inbounds stencil_interior( - ::RightBiasedF2C, - loc, - space, - idx, - hidx, - arg, -) = getidx(space, arg, loc, idx + half, hidx) - -left_interior_idx( - space::AbstractSpace, - ::RightBiasedF2C, - ::AbstractBoundaryCondition, - arg, -) = left_idx(space) -right_interior_idx( - space::AbstractSpace, - ::RightBiasedF2C, - ::AbstractBoundaryCondition, - arg, -) = right_idx(space) - -right_interior_idx(space::AbstractSpace, ::RightBiasedF2C, ::SetValue, arg) = - right_idx(space) - 1 -Base.@propagate_inbounds function stencil_right_boundary( - ::RightBiasedF2C, - bc::SetValue, - loc, - space, - idx, - hidx, - arg, -) - @assert idx == right_center_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end - - -""" - R = RightBiased3rdOrderC2F(;boundaries) - R.(x) - -Interpolate a center-valued field to a face-valued field from the right, using a 3rd-order reconstruction. -```math -R(x)[i] = \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}] \\right) / 12 -``` - -Only the right boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -R(x)[n+\\tfrac{1}{2}] = x_0 -``` -""" -struct RightBiased3rdOrderC2F{BCS} <: InterpolationOperator - bcs::BCS -end -RightBiased3rdOrderC2F(; kwargs...) = RightBiased3rdOrderC2F(NamedTuple(kwargs)) - -return_space( - ::RightBiased3rdOrderC2F, - space::Spaces.CenterFiniteDifferenceSpace, -) = Spaces.FaceFiniteDifferenceSpace(space) -return_space( - ::RightBiased3rdOrderC2F, - space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) - -stencil_interior_width(::RightBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),) -Base.@propagate_inbounds stencil_interior( - ::RightBiased3rdOrderC2F, - loc, - space, - idx, - hidx, - arg, -) = - ( - 4 * getidx(space, arg, loc, idx - half, hidx) + - 10 * getidx(space, arg, loc, idx + half, hidx) - - 2 * getidx(space, arg, loc, idx + half + 1, hidx) - ) / 12 - -boundary_width(::RightBiased3rdOrderC2F, ::SetValue) = 1 -Base.@propagate_inbounds function stencil_right_boundary( - ::RightBiased3rdOrderC2F, - bc::SetValue, - loc, - space, - idx, - hidx, - arg, -) - @assert idx == right_face_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end - -""" - R = RightBiased3rdOrderF2C(;boundaries) - R.(x) - -Interpolate a face-valued field to a center-valued field from the right, using a 3rd-order reconstruction. -```math -R(x)[i] = \\left(4 x[i] + 10 x[i+1] -2 x[i+2] \\right) / 12 -``` - -Only the right boundary condition should be set. Currently supported is: -- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. -```math -R(x)[n+\\tfrac{1}{2}] = x_0 -``` -""" -struct RightBiased3rdOrderF2C{BCS} <: InterpolationOperator - bcs::BCS -end -RightBiased3rdOrderF2C(; kwargs...) = RightBiased3rdOrderF2C(NamedTuple(kwargs)) - -return_space( - ::RightBiased3rdOrderF2C, - space::Spaces.FaceFiniteDifferenceSpace, -) = Spaces.CenterFiniteDifferenceSpace(space) -return_space( - ::RightBiased3rdOrderF2C, - space::Spaces.FaceExtrudedFiniteDifferenceSpace, -) = Spaces.CenterExtrudedFiniteDifferenceSpace(space) - -stencil_interior_width(::RightBiased3rdOrderF2C, arg) = ((-half - 1, half + 1),) -Base.@propagate_inbounds stencil_interior( - ::RightBiased3rdOrderF2C, - loc, - space, - idx, - hidx, - arg, -) = - ( - 4 * getidx(space, arg, loc, idx - half, hidx) + - 10 * getidx(space, arg, loc, idx + half, hidx) - - 2 * getidx(space, arg, loc, idx + half + 1, hidx) - ) / 12 - -boundary_width(::RightBiased3rdOrderF2C, ::SetValue) = 1 -Base.@propagate_inbounds function stencil_right_boundary( - ::RightBiased3rdOrderF2C, - bc::SetValue, - loc, - space, - idx, - hidx, - arg, -) - @assert idx == right_center_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end - -abstract type WeightedInterpolationOperator <: InterpolationOperator end -# TODO: this is not in general correct and the return type -# should be based on the component operator types (rdiv, rmul) but we don't have a good way -# of creating ex. one(field_type) for complex fields for inference -return_eltype(::WeightedInterpolationOperator, weights, arg) = eltype(arg) - -""" - WI = WeightedInterpolateF2C(; boundaries) - WI.(w, x) - -Interpolate a face-valued field `x` to centers, weighted by a face-valued field -`w`, using the stencil -```math -WI(w, x)[i] = \\frac{ - w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}]) - }{ - w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] - } -``` - -No boundary conditions are required (or supported) -""" -struct WeightedInterpolateF2C{BCS} <: WeightedInterpolationOperator - bcs::BCS -end -WeightedInterpolateF2C(; kwargs...) = WeightedInterpolateF2C(NamedTuple(kwargs)) - -return_space( - ::WeightedInterpolateF2C, - weight_space::Spaces.FaceFiniteDifferenceSpace, - arg_space::Spaces.FaceFiniteDifferenceSpace, -) = Spaces.CenterFiniteDifferenceSpace(arg_space) -return_space( - ::WeightedInterpolateF2C, - weight_space::Spaces.FaceExtrudedFiniteDifferenceSpace, - arg_space::Spaces.FaceExtrudedFiniteDifferenceSpace, -) = Spaces.CenterExtrudedFiniteDifferenceSpace(arg_space) - -stencil_interior_width(::WeightedInterpolateF2C, weight, arg) = - ((-half, half), (-half, half)) -Base.@propagate_inbounds function stencil_interior( - ::WeightedInterpolateF2C, - loc, - space, - idx, - hidx, - weight, - arg, -) - w⁺ = getidx(space, weight, loc, idx + half, hidx) - w⁻ = getidx(space, weight, loc, idx - half, hidx) - a⁺ = getidx(space, arg, loc, idx + half, hidx) - a⁻ = getidx(space, arg, loc, idx - half, hidx) - RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻)) -end - -boundary_width(::WeightedInterpolateF2C, ::AbstractBoundaryCondition) = 0 - -""" - WI = WeightedInterpolateC2F(; boundaries) - WI.(w, x) - -Interpolate a center-valued field `x` to faces, weighted by a center-valued field -`w`, using the stencil -```math -WI(w, x)[i] = \\frac{ - w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}]) -}{ - w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}] -} -``` - -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). -""" -struct WeightedInterpolateC2F{BCS} <: WeightedInterpolationOperator - bcs::BCS -end -WeightedInterpolateC2F(; kwargs...) = WeightedInterpolateC2F(NamedTuple(kwargs)) - -return_space( - ::WeightedInterpolateC2F, - weight_space::Spaces.CenterFiniteDifferenceSpace, - arg_space::Spaces.CenterFiniteDifferenceSpace, -) = Spaces.FaceFiniteDifferenceSpace(arg_space) -return_space( - ::WeightedInterpolateC2F, - weight_space::Spaces.CenterExtrudedFiniteDifferenceSpace, - arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = Spaces.FaceExtrudedFiniteDifferenceSpace(arg_space) - -stencil_interior_width(::WeightedInterpolateC2F, weight, arg) = - ((-half, half), (-half, half)) -Base.@propagate_inbounds function stencil_interior( - ::WeightedInterpolateC2F, - loc, - space, - idx, - hidx, - weight, - arg, -) - w⁺ = getidx(space, weight, loc, idx + half, hidx) - w⁻ = getidx(space, weight, loc, idx - half, hidx) - a⁺ = getidx(space, arg, loc, idx + half, hidx) - a⁻ = getidx(space, arg, loc, idx - half, hidx) - RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (w⁺ ⊞ w⁻)) -end - -boundary_width(::WeightedInterpolateC2F, ::AbstractBoundaryCondition) = 1 -Base.@propagate_inbounds function stencil_left_boundary( - ::WeightedInterpolateC2F, - bc::SetValue, - loc, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == left_face_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end -Base.@propagate_inbounds function stencil_right_boundary( - ::WeightedInterpolateC2F, - bc::SetValue, - loc, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == right_face_boundary_idx(space) - getidx(space, bc.val, loc, nothing, hidx) -end - -Base.@propagate_inbounds function stencil_left_boundary( - ::WeightedInterpolateC2F, - bc::SetGradient, - loc, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == left_face_boundary_idx(space) - a⁺ = getidx(space, arg, loc, idx + half, hidx) - v₃ = Geometry.covariant3( - getidx(space, bc.val, loc, nothing, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁺ ⊟ RecursiveApply.rdiv(v₃, 2) -end -Base.@propagate_inbounds function stencil_right_boundary( - ::WeightedInterpolateC2F, - bc::SetGradient, - loc, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == right_face_boundary_idx(space) - a⁻ = getidx(space, arg, loc, idx - half, hidx) - v₃ = Geometry.covariant3( - getidx(space, bc.val, loc, nothing, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁻ ⊞ RecursiveApply.rdiv(v₃, 2) -end - -Base.@propagate_inbounds function stencil_left_boundary( - ::WeightedInterpolateC2F, - bc::Extrapolate, - loc, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == left_face_boundary_idx(space) - a⁺ = getidx(space, arg, loc, idx + half, hidx) - a⁺ -end -Base.@propagate_inbounds function stencil_right_boundary( - ::WeightedInterpolateC2F, - bc::Extrapolate, - loc, - space, - idx, - hidx, - weight, - arg, -) - @assert idx == right_face_boundary_idx(space) - a⁻ = getidx(space, arg, loc, idx - half, hidx) - a⁻ -end - - -abstract type AdvectionOperator <: FiniteDifferenceOperator end -return_eltype(::AdvectionOperator, velocity, arg) = eltype(arg) - -""" - U = UpwindBiasedProductC2F(;boundaries) - U.(v, x) - -Compute the product of the face-valued vector field `v` and a center-valued -field `x` at cell faces by upwinding `x` according to the direction of `v`. - -More precisely, it is computed based on the sign of the 3rd contravariant -component, and it returns a `Contravariant3Vector`: -```math -U(\\boldsymbol{v},x)[i] = \\begin{cases} - v^3[i] x[i-\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] > 0 \\\\ - v^3[i] x[i+\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] < 0 - \\end{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 - U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = U(\\boldsymbol{v},x)[1 + \\tfrac{1}{2}] - ``` -""" -struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator - bcs::BCS -end -UpwindBiasedProductC2F(; kwargs...) = UpwindBiasedProductC2F(NamedTuple(kwargs)) - -return_eltype(::UpwindBiasedProductC2F, V, A) = - Geometry.Contravariant3Vector{eltype(eltype(V))} - -return_space( - ::UpwindBiasedProductC2F, - velocity_space::Spaces.FaceFiniteDifferenceSpace, - arg_space::Spaces.CenterFiniteDifferenceSpace, -) = velocity_space -return_space( - ::UpwindBiasedProductC2F, - velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace, - arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = velocity_space - -function upwind_biased_product(v, a⁻, a⁺) - RecursiveApply.rdiv( - ((v ⊞ RecursiveApply.rmap(abs, v)) ⊠ a⁻) ⊞ - ((v ⊟ RecursiveApply.rmap(abs, v)) ⊠ a⁺), - 2, - ) -end - -stencil_interior_width(::UpwindBiasedProductC2F, velocity, arg) = - ((0, 0), (-half, half)) - -Base.@propagate_inbounds function stencil_interior( - ::UpwindBiasedProductC2F, - loc, - space, - idx, - hidx, - velocity, - arg, -) - a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg) - a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg) - vᶠ = Geometry.contravariant3( - getidx(space, velocity, loc, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, a⁻, a⁺)) -end - -boundary_width(::UpwindBiasedProductC2F, ::AbstractBoundaryCondition) = 1 - -Base.@propagate_inbounds function stencil_left_boundary( - ::UpwindBiasedProductC2F, - bc::SetValue, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_face_boundary_idx(space) - aᴸᴮ = getidx(space, bc.val, loc, nothing, hidx) - a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg) - vᶠ = Geometry.contravariant3( - getidx(space, velocity, loc, 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, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_face_boundary_idx(space) - a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg) - aᴿᴮ = getidx(space, bc.val, loc, nothing, hidx) - vᶠ = Geometry.contravariant3( - getidx(space, velocity, loc, 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, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == left_face_boundary_idx(space) - stencil_interior(op, loc, space, idx + 1, hidx, velocity, arg) -end - -Base.@propagate_inbounds function stencil_right_boundary( - op::UpwindBiasedProductC2F, - ::Extrapolate, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx == right_face_boundary_idx(space) - stencil_interior(op, loc, space, idx - 1, hidx, velocity, arg) -end - -""" - U = Upwind3rdOrderBiasedProductC2F(;boundaries) - U.(v, x) - -Compute the product of a face-valued vector field `v` and a center-valued field -`x` at cell faces by upwinding `x`, to third-order of accuracy, according to `v` -```math -U(v,x)[i] = \\begin{cases} - v[i] \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12 \\textrm{, if } v[i] > 0 \\\\ - v[i] \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}] \\right) / 12 \\textrm{, if } v[i] < 0 - \\end{cases} -``` -This stencil is based on [WickerSkamarock2002](@cite), eq. 4(a). - -Supported boundary conditions are: -- [`FirstOrderOneSided(x₀)`](@ref): uses the first-order downwind scheme to compute `x` on the left boundary, - and the first-order upwind scheme to compute `x` on the right boundary. -- [`ThirdOrderOneSided(x₀)`](@ref): uses the third-order downwind reconstruction to compute `x` on the left boundary, -and the third-order upwind reconstruction to compute `x` on the right boundary. - -!!! note - These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a [`DivergenceF2C`](@ref) operator, with a [`SetValue`](@ref) boundary. -""" -struct Upwind3rdOrderBiasedProductC2F{BCS} <: AdvectionOperator - bcs::BCS -end -Upwind3rdOrderBiasedProductC2F(; kwargs...) = - Upwind3rdOrderBiasedProductC2F(NamedTuple(kwargs)) - -return_eltype(::Upwind3rdOrderBiasedProductC2F, V, A) = - Geometry.Contravariant3Vector{eltype(eltype(V))} - -return_space( - ::Upwind3rdOrderBiasedProductC2F, - velocity_space::Spaces.FaceFiniteDifferenceSpace, - arg_space::Spaces.CenterFiniteDifferenceSpace, -) = velocity_space -return_space( - ::Upwind3rdOrderBiasedProductC2F, - velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace, - arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, -) = velocity_space - -function upwind_3rdorder_biased_product(v, a⁻, a⁻⁻, a⁺, a⁺⁺) - RecursiveApply.rdiv( - (v ⊠ (7 ⊠ (a⁺ + a⁻) ⊟ (a⁺⁺ + a⁻⁻))) ⊟ - (RecursiveApply.rmap(abs, v) ⊠ (3 ⊠ (a⁺ - a⁻) ⊟ (a⁺⁺ - a⁻⁻))), - 12, - ) -end - -stencil_interior_width(::Upwind3rdOrderBiasedProductC2F, velocity, arg) = - ((0, 0), (-half - 1, half + 1)) - -Base.@propagate_inbounds function stencil_interior( - ::Upwind3rdOrderBiasedProductC2F, - loc, - space, - idx, - hidx, - velocity, - arg, -) - a⁻ = getidx(space, arg, loc, idx - half, hidx) - a⁻⁻ = getidx(space, arg, loc, idx - half - 1, hidx) - a⁺ = getidx(space, arg, loc, idx + half, hidx) - a⁺⁺ = getidx(space, arg, loc, idx + half + 1, hidx) - vᶠ = Geometry.contravariant3( - getidx(space, velocity, loc, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - return Geometry.Contravariant3Vector( - upwind_3rdorder_biased_product(vᶠ, a⁻, a⁻⁻, a⁺, a⁺⁺), - ) -end - -boundary_width(::Upwind3rdOrderBiasedProductC2F, ::AbstractBoundaryCondition) = - 2 - -Base.@propagate_inbounds function stencil_left_boundary( - ::Upwind3rdOrderBiasedProductC2F, - bc::FirstOrderOneSided, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx <= left_face_boundary_idx(space) + 1 - v = Geometry.contravariant3( - getidx(space, velocity, loc, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁻ = stencil_interior(LeftBiasedC2F(), loc, space, idx, hidx, arg) - a⁺ = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) - return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺)) -end - -Base.@propagate_inbounds function stencil_right_boundary( - ::Upwind3rdOrderBiasedProductC2F, - bc::FirstOrderOneSided, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx >= right_face_boundary_idx(space) - 1 - v = Geometry.contravariant3( - getidx(space, velocity, loc, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a⁻ = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) - a⁺ = stencil_interior(RightBiasedC2F(), loc, space, idx, hidx, arg) - return Geometry.Contravariant3Vector(upwind_biased_product(v, a⁻, a⁺)) - -end - -Base.@propagate_inbounds function stencil_left_boundary( - ::Upwind3rdOrderBiasedProductC2F, - bc::ThirdOrderOneSided, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx <= left_face_boundary_idx(space) + 1 - - vᶠ = Geometry.contravariant3( - getidx(space, velocity, loc, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) - - return Geometry.Contravariant3Vector(vᶠ * a) -end - -Base.@propagate_inbounds function stencil_right_boundary( - ::Upwind3rdOrderBiasedProductC2F, - bc::ThirdOrderOneSided, - loc, - space, - idx, - hidx, - velocity, - arg, -) - @assert idx <= right_face_boundary_idx(space) - 1 - - vᶠ = Geometry.contravariant3( - getidx(space, velocity, loc, idx, hidx), - Geometry.LocalGeometry(space, idx, hidx), - ) - a = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) - - return Geometry.Contravariant3Vector(vᶠ * a) -end """ U = FCTBorisBook(;boundaries) @@ -1564,7 +752,7 @@ where ``s[i] = +1`` if `` v[i] \\geq 0`` and ``s[i] = -1`` if `` v[i] \\leq 0` This formulation is based on [BorisBook1973](@cite), as reported in [durran2010](@cite) section 5.4.1. Supported boundary conditions are: -- [`FirstOrderOneSided(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary. +- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary. !!! note Similar to the [`Upwind3rdOrderBiasedProductC2F`](@ref) operator, these boundary conditions do not define the value at the actual boundary faces, @@ -1645,7 +833,7 @@ boundary_width(::FCTBorisBook, ::AbstractBoundaryCondition) = 2 Base.@propagate_inbounds function stencil_left_boundary( ::FCTBorisBook, - bc::FirstOrderOneSided, + bc::OneSided1stOrder, loc, space, idx, @@ -1664,7 +852,7 @@ end Base.@propagate_inbounds function stencil_right_boundary( ::FCTBorisBook, - bc::FirstOrderOneSided, + bc::OneSided1stOrder, loc, space, idx, @@ -1701,7 +889,7 @@ This stencil is based on [zalesak1979fully](@cite), as reported in [durran2010]( the corrected antidiffusive flux. Supported boundary conditions are: -- [`FirstOrderOneSided(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary. +- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind reconstruction to compute `x` on the left boundary, and the first-order upwind reconstruction to compute `x` on the right boundary. !!! note Similar to the [`Upwind3rdOrderBiasedProductC2F`](@ref) operator, these boundary conditions do not define @@ -1838,7 +1026,7 @@ boundary_width(::FCTZalesak, ::AbstractBoundaryCondition) = 2 Base.@propagate_inbounds function stencil_left_boundary( ::FCTZalesak, - bc::FirstOrderOneSided, + bc::OneSided1stOrder, loc, space, idx, @@ -1854,7 +1042,7 @@ end Base.@propagate_inbounds function stencil_right_boundary( ::FCTZalesak, - bc::FirstOrderOneSided, + bc::OneSided1stOrder, loc, space, idx, diff --git a/src/Operators/finitedifference/upwinding.jl b/src/Operators/finitedifference/upwinding.jl new file mode 100644 index 0000000000..64bb79e402 --- /dev/null +++ b/src/Operators/finitedifference/upwinding.jl @@ -0,0 +1,582 @@ +""" + L = LeftBiased1stOrderC2F(;boundaries) + L.(x) + +Interpolate a center-value field to a face-valued field from the left. +```math +L(x)[i] = x[i-\\tfrac{1}{2}] +``` + +Only the left boundary condition should be set. Currently supported are: +- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. +```math +L(x)[\\tfrac{1}{2}] = x_0 +``` +- [`Extrapolate()`](@ref): use the right-biased interior value. +```math +L(x)[\\tfrac{1}{2}] = x[1] +``` +""" +struct LeftBiased1stOrderC2F{BCS} <: InterpolationOperator + bcs::BCS +end +LeftBiased1stOrderC2F(; kwargs...) = LeftBiased1stOrderC2F(NamedTuple(kwargs)) + +return_space( + ::LeftBiased1stOrderC2F, + space::Spaces.CenterFiniteDifferenceSpace, +) = Spaces.FaceFiniteDifferenceSpace(space) +return_space( + ::LeftBiased1stOrderC2F, + space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::LeftBiased1stOrderC2F, arg) = ((-half, -half),) +Base.@propagate_inbounds stencil_interior( + ::LeftBiased1stOrderC2F, + loc, + space, + idx, + hidx, + arg, +) = getidx(space, arg, loc, idx - half, hidx) + +left_interior_idx( + space::AbstractSpace, + ::LeftBiased1stOrderC2F, + ::AbstractBoundaryCondition, + arg, +) = left_idx(space) + 1 +right_interior_idx( + space::AbstractSpace, + ::LeftBiased1stOrderC2F, + ::AbstractBoundaryCondition, + arg, +) = right_idx(space) + +Base.@propagate_inbounds function stencil_left_boundary( + ::LeftBiased1stOrderC2F, + bc::SetValue, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == left_face_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end +Base.@propagate_inbounds function stencil_left_boundary( + ::LeftBiased1stOrderC2F, + bc::Extrapolate, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == left_face_boundary_idx(space) + getidx(space, arg, loc, idx + half, hidx) +end + + + +""" + L = LeftBiased3rdOrderC2F(;boundaries) + L.(x) + +Interpolate a center-value field to a face-valued field from the left, using a +3rd-order reconstruction. +```math +L(x)[i] = \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12 +``` + +Only the left boundary condition should be set. Currently supported are: + +- [`OneSided1stOrder()`](@ref): use the one-sided 1st-order reconstruction. +```math +L(x)[1+\\tfrac{1}{2}] = x[1] +``` +- [`OneSided3rdOrder()`](@ref): use the [`RightBiased3rdOrderC2F`](@ref) + reconstruction for ``L(x)[1+\\tfrac{1}{2}]``. + +!!! note + The boundary conditions only specify the first interior point. Actual boundary values + remain undefined and need to be handled separately. +""" +struct LeftBiased3rdOrderC2F{BCS} <: InterpolationOperator + bcs::BCS +end +LeftBiased3rdOrderC2F(; kwargs...) = LeftBiased3rdOrderC2F(NamedTuple(kwargs)) + +return_space( + ::LeftBiased3rdOrderC2F, + space::Spaces.CenterFiniteDifferenceSpace, +) = Spaces.FaceFiniteDifferenceSpace(space) +return_space( + ::LeftBiased3rdOrderC2F, + space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::LeftBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),) +Base.@propagate_inbounds function stencil_interior( + ::LeftBiased3rdOrderC2F, + loc, + space, + idx, + hidx, + arg, +) + FT = Spaces.undertype(space) + + FT(4 / 12) ⊠ getidx(space, arg, loc, idx + half, hidx) ⊞ + FT(10 / 12) ⊠ getidx(space, arg, loc, idx - half, hidx) ⊟ + FT(2 / 12) ⊠ getidx(space, arg, loc, idx - half - 1, hidx) +end + +left_interior_idx( + space::AbstractSpace, + ::LeftBiased3rdOrderC2F, + ::AbstractBoundaryCondition, + arg, +) = left_idx(space) + 2 +right_interior_idx( + space::AbstractSpace, + ::LeftBiased3rdOrderC2F, + ::AbstractBoundaryCondition, + arg, +) = right_idx(space) - 1 + +Base.@propagate_inbounds function stencil_left_boundary( + ::LeftBiased3rdOrderC2F, + bc::OneSided1stOrder, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == left_face_boundary_idx(space) + 1 + stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg) +end +Base.@propagate_inbounds function stencil_left_boundary( + ::LeftBiased3rdOrderC2F, + bc::OneSided3rdOrder, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == left_face_boundary_idx(space) + 1 + stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) +end + + +""" + R = RightBiased1stOrderC2F(;boundaries) + R.(x) + +Interpolate a center-valued field to a face-valued field from the right. +```math +R(x)[i] = x[i+\\tfrac{1}{2}] +``` + +Only the right boundary condition should be set. Currently supported: +- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary. +```math +R(x)[n+\\tfrac{1}{2}] = x_0 +``` +- [`Extrapolate()`](@ref): use the left-biased interior value. +```math +R(x)[n+\\tfrac{1}{2}] = x[n] +``` +""" +struct RightBiased1stOrderC2F{BCS} <: InterpolationOperator + bcs::BCS +end +RightBiased1stOrderC2F(; kwargs...) = RightBiased1stOrderC2F(NamedTuple(kwargs)) + +return_space( + ::RightBiased1stOrderC2F, + space::Spaces.CenterFiniteDifferenceSpace, +) = Spaces.FaceFiniteDifferenceSpace(space) +return_space( + ::RightBiased1stOrderC2F, + space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::RightBiased1stOrderC2F, arg) = ((half, half),) +Base.@propagate_inbounds stencil_interior( + ::RightBiased1stOrderC2F, + loc, + space, + idx, + hidx, + arg, +) = getidx(space, arg, loc, idx + half, hidx) + +left_interior_idx( + space::AbstractSpace, + ::RightBiased1stOrderC2F, + ::AbstractBoundaryCondition, + arg, +) = left_idx(space) +right_interior_idx( + space::AbstractSpace, + ::RightBiased1stOrderC2F, + ::AbstractBoundaryCondition, + arg, +) = right_idx(space) - 1 + +Base.@propagate_inbounds function stencil_right_boundary( + ::RightBiased1stOrderC2F, + bc::SetValue, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == right_face_boundary_idx(space) + getidx(space, bc.val, loc, nothing, hidx) +end +Base.@propagate_inbounds function stencil_right_boundary( + ::RightBiased1stOrderC2F, + bc::Extrapolate, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == right_face_boundary_idx(space) + getidx(space, arg, loc, idx - half, hidx) +end + +""" + R = RightBiased3rdOrderC2F(;boundaries) + R.(x) + +Interpolate a center-valued field to a face-valued field from the right, using a 3rd-order reconstruction. +```math +R(x)[i] = \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}] \\right) / 12 +``` + +Only the left boundary condition should be set. Currently supported are: + +- [`OneSided1stOrder()`](@ref): use the one-sided 1st-order reconstruction. +```math +R(x)[n-\\tfrac{1}{2}] = x[1] +``` +- [`OneSided3rdOrder()`](@ref): use the [`RightBiased3rdOrderC2F`](@ref) + reconstruction for ``R(x)[n-\\tfrac{1}{2}]`` + +!!! note + The boundary conditions only specify the first interior point. Actual boundary values + remain undefined and need to be handled separately. +""" +struct RightBiased3rdOrderC2F{BCS} <: InterpolationOperator + bcs::BCS +end +RightBiased3rdOrderC2F(; kwargs...) = RightBiased3rdOrderC2F(NamedTuple(kwargs)) + +return_space( + ::RightBiased3rdOrderC2F, + space::Spaces.CenterFiniteDifferenceSpace, +) = Spaces.FaceFiniteDifferenceSpace(space) +return_space( + ::RightBiased3rdOrderC2F, + space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = Spaces.FaceExtrudedFiniteDifferenceSpace(space) + +stencil_interior_width(::RightBiased3rdOrderC2F, arg) = ((-half - 1, half + 1),) +Base.@propagate_inbounds function stencil_interior( + ::RightBiased3rdOrderC2F, + loc, + space, + idx, + hidx, + arg, +) + FT = Spaces.undertype(space) + + FT(4 / 12) ⊠ getidx(space, arg, loc, idx - half, hidx) ⊞ + FT(10 / 12) ⊠ getidx(space, arg, loc, idx + half, hidx) ⊟ + FT(2 / 12) ⊠ getidx(space, arg, loc, idx + half + 1, hidx) +end + +Base.@propagate_inbounds function stencil_right_boundary( + ::RightBiased3rdOrderC2F, + bc::OneSided1stOrder, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == right_face_boundary_idx(space) - 1 + stencil_interior(RightBiased1stOrderC2F(), loc, space, idx, hidx, arg) +end +Base.@propagate_inbounds function stencil_right_boundary( + ::RightBiased3rdOrderC2F, + bc::OneSided3rdOrder, + loc, + space, + idx, + hidx, + arg, +) + @assert idx == right_face_boundary_idx(space) - 1 + stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) +end + +""" + U = Upwind1stOrderC2F(;boundaries) + U.(v, x) + +Interpolate the center-valued field `x` to cell faces by 1st-order upwinding `x` +according to the direction of `v`. + +More precisely, it is computed based on the sign of the 3rd contravariant +component, and it returns a `Contravariant3Vector`: +```math +U(\\boldsymbol{v},x)[i] = \\begin{cases} + x[i-\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] > 0 \\\\ + x[i+\\tfrac{1}{2}]\\boldsymbol{e}_3 \\textrm{, if } v^3[i] < 0 + \\end{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} + x_0 \\boldsymbol{e}_3 \\textrm{, if } v^3[\\tfrac{1}{2}] > 0 \\\\ + 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 + U(\\boldsymbol{v},x)[\\tfrac{1}{2}] = U(\\boldsymbol{v},x)[1 + \\tfrac{1}{2}] + ``` +""" +struct Upwind1stOrderC2F{BCS} <: InterpolationOperator + bcs::BCS +end +Upwind1stOrderC2F(; kwargs...) = Upwind1stOrderC2F(NamedTuple(kwargs)) + +return_eltype(::Upwind1stOrderC2F, v, arg) = eltype(arg) + +return_space( + ::Upwind1stOrderC2F, + velocity_space::Spaces.FaceFiniteDifferenceSpace, + arg_space::Spaces.CenterFiniteDifferenceSpace, +) = velocity_space +return_space( + ::Upwind1stOrderC2F, + velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace, + arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = velocity_space + +stencil_interior_width(::Upwind1stOrderC2F, velocity, arg) = + ((0, 0), (-half, half)) + + +Base.@propagate_inbounds function stencil_interior( + ::Upwind1stOrderC2F, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ = stencil_interior(RightBiased1stOrderC2F(), loc, space, idx, hidx, arg) + a⁻ = stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg) + # signbit(v³) == true if v³ < 0 + return signbit(v³) ? a⁺ : a⁻ +end + +boundary_width(::Upwind1stOrderC2F, ::AbstractBoundaryCondition) = 1 + +Base.@propagate_inbounds function stencil_left_boundary( + ::Upwind1stOrderC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + @assert idx == left_face_boundary_idx(space) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ = stencil_interior(RightBiased1stOrderC2F(), loc, space, idx, hidx, arg) + a⁻ = stencil_left_boundary( + LeftBiased1stOrderC2F(), + bc, + loc, + space, + idx, + hidx, + arg, + ) + return signbit(v³) ? a⁺ : a⁻ +end + +Base.@propagate_inbounds function stencil_right_boundary( + ::Upwind1stOrderC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + @assert idx == right_face_boundary_idx(space) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ = stencil_right_boundary( + RightBiased1stOrderC2F(), + bc, + loc, + space, + idx, + hidx, + arg, + ) + a⁻ = stencil_interior(LeftBiased1stOrderC2F(), loc, space, idx, hidx, arg) + return signbit(v³) ? a⁺ : a⁻ +end + + +""" + U = Upwind3rdOrderC2F(;boundaries) + U.(v, x) + +Interpolate the center-valued field `x` to cell faces by 3rd-order upwinding `x` +according to the direction of `v`. The stencil is +```math +U(v,x)[i] = \\begin{cases} + \\left(-2 x[i-\\tfrac{3}{2}] + 10 x[i-\\tfrac{1}{2}] + 4 x[i+\\tfrac{1}{2}] \\right) / 12 \\textrm{, if } v[i] > 0 \\\\ + \\left(4 x[i-\\tfrac{1}{2}] + 10 x[i+\\tfrac{1}{2}] -2 x[i+\\tfrac{3}{2}] \\right) / 12 \\textrm{, if } v[i] < 0 + \\end{cases} +``` +This stencil is based on [WickerSkamarock2002](@cite), eq. 4(a). + +Supported boundary conditions are: +- [`OneSided1stOrder(x₀)`](@ref): uses the first-order downwind scheme to compute `x` on the left boundary, + and the first-order upwind scheme to compute `x` on the right boundary. +- [`OneSided3rdOrder(x₀)`](@ref): uses the third-order downwind reconstruction to compute `x` on the left boundary, +and the third-order upwind reconstruction to compute `x` on the right boundary. + +!!! note + These boundary conditions do not define the value at the actual boundary faces, and so this operator cannot be materialized directly: it needs to be composed with another operator that does not make use of this value, e.g. a [`DivergenceF2C`](@ref) operator, with a [`SetValue`](@ref) boundary. +""" +struct Upwind3rdOrderC2F{BCS} <: InterpolationOperator + bcs::BCS +end +Upwind3rdOrderC2F(; kwargs...) = Upwind3rdOrderC2F(NamedTuple(kwargs)) + +return_eltype(::Upwind3rdOrderC2F, V, arg) = eltype(arg) + +return_space( + ::Upwind3rdOrderC2F, + velocity_space::Spaces.FaceFiniteDifferenceSpace, + arg_space::Spaces.CenterFiniteDifferenceSpace, +) = velocity_space +return_space( + ::Upwind3rdOrderC2F, + velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace, + arg_space::Spaces.CenterExtrudedFiniteDifferenceSpace, +) = velocity_space + +stencil_interior_width(::Upwind3rdOrderC2F, velocity, arg) = + ((0, 0), (-half - 1, half + 1)) + +Base.@propagate_inbounds function stencil_interior( + ::Upwind3rdOrderC2F, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + a⁻ = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + return signbit(v³) ? a⁺ : a⁻ +end + +boundary_width(::Upwind3rdOrderC2F, ::AbstractBoundaryCondition) = 2 + +Base.@propagate_inbounds function stencil_left_boundary( + ::Upwind3rdOrderC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ = stencil_interior(RightBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + a⁻ = stencil_left_boundary( + LeftBiased3rdOrderC2F(), + bc, + loc, + space, + idx, + hidx, + arg, + ) + return signbit(v³) ? a⁺ : a⁻ +end +Base.@propagate_inbounds function stencil_right_boundary( + ::Upwind3rdOrderC2F, + bc::AbstractBoundaryCondition, + loc, + space, + idx, + hidx, + velocity, + arg, +) + v³ = Geometry.contravariant3( + getidx(space, velocity, loc, idx, hidx), + Geometry.LocalGeometry(space, idx, hidx), + ) + a⁺ = stencil_right_boundary( + RightBiased3rdOrderC2F(), + bc, + loc, + space, + idx, + hidx, + arg, + ) + a⁻ = stencil_interior(LeftBiased3rdOrderC2F(), loc, space, idx, hidx, arg) + return signbit(v³) ? a⁺ : a⁻ +end diff --git a/src/Operators/operator2stencil.jl b/src/Operators/operator2stencil.jl index 111a7e4f25..65ebd83f7c 100644 --- a/src/Operators/operator2stencil.jl +++ b/src/Operators/operator2stencil.jl @@ -166,7 +166,7 @@ end # TODO: find out why using Base.@propagate_inbounds blows up compilation time function stencil_interior( - ::Operator2Stencil{<:Union{LeftBiasedF2C, LeftBiasedC2F}}, + ::Operator2Stencil{<:Union{LeftBiasedF2C, LeftBiased1stOrderC2F}}, loc, space, idx, @@ -177,7 +177,7 @@ function stencil_interior( return StencilCoefs{-half, -half}((val⁻,)) end stencil_left_boundary( - ::Operator2Stencil{<:Union{LeftBiasedF2C, LeftBiasedC2F}}, + ::Operator2Stencil{<:Union{LeftBiasedF2C, LeftBiased1stOrderC2F}}, bc::SetValue, loc, space, @@ -189,7 +189,7 @@ stencil_left_boundary( # TODO: find out why using Base.@propagate_inbounds blows up compilation time function stencil_interior( - ::Operator2Stencil{<:Union{RightBiasedF2C, RightBiasedC2F}}, + ::Operator2Stencil{<:Union{RightBiasedF2C, RightBiased1stOrderC2F}}, loc, space, idx, @@ -200,7 +200,7 @@ function stencil_interior( return StencilCoefs{half, half}((val⁺,)) end stencil_right_boundary( - ::Operator2Stencil{<:Union{RightBiasedF2C, RightBiasedC2F}}, + ::Operator2Stencil{<:Union{RightBiasedF2C, RightBiased1stOrderC2F}}, bc::SetValue, loc, space, diff --git a/test/Operators/finitedifference/column.jl b/test/Operators/finitedifference/column.jl index 98d0ef6a24..3f23a42e8a 100644 --- a/test/Operators/finitedifference/column.jl +++ b/test/Operators/finitedifference/column.jl @@ -731,7 +731,7 @@ end @test conv_adv_wc[3] ≈ 2 atol = 0.1 end -@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with FirstOrderOneSided + DivergenceF2C SetValue BCs, constant w" begin +@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided1stOrder + DivergenceF2C SetValue BCs, constant w" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -761,8 +761,8 @@ end s = sin.(centers) third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) divf2c = Operators.DivergenceF2C( @@ -792,7 +792,7 @@ end end end -@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with ThirdOrderOneSided + DivergenceF2C SetValue BCs, varying sign w" begin +@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder + DivergenceF2C SetValue BCs, varying sign w" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -821,8 +821,8 @@ end c = sin.(centers) third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) divf2c = Operators.DivergenceF2C( @@ -903,7 +903,7 @@ end @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2] end -@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with FirstOrderOneSided BCs" begin +@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided1stOrder BCs" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -938,8 +938,8 @@ end top = Operators.Extrapolate(), ) third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided(), + bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder(), ) divf2c = Operators.DivergenceF2C( @@ -973,7 +973,7 @@ end end end -@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with ThirdOrderOneSided BCs" begin +@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder BCs" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -1006,8 +1006,8 @@ end top = Operators.Extrapolate(), ) third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( - bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided(), + bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder(), ) divf2c = Operators.DivergenceF2C( @@ -1071,13 +1071,13 @@ end fyp = parent(fy) # C2F biased operators - LBC2F = Operators.LeftBiasedC2F(; bottom = Operators.SetValue(10)) + LBC2F = Operators.LeftBiased1stOrderC2F(; bottom = Operators.SetValue(10)) @. cy = cos(zc) @. fy = LBC2F(cy) fy_ref = [FT(10), [cyp[i] for i in 1:length(cyp)]...] @test all(fy_ref .== fyp) - RBC2F = Operators.RightBiasedC2F(; top = Operators.SetValue(10)) + RBC2F = Operators.RightBiased1stOrderC2F(; top = Operators.SetValue(10)) @. cy = cos(zc) @. fy = RBC2F(cy) fy_ref = [[cyp[i] for i in 1:length(cyp)]..., FT(10)] diff --git a/test/Operators/finitedifference/column_benchmark_kernels.jl b/test/Operators/finitedifference/column_benchmark_kernels.jl index 57c5fdee7d..d320a52a85 100644 --- a/test/Operators/finitedifference/column_benchmark_kernels.jl +++ b/test/Operators/finitedifference/column_benchmark_kernels.jl @@ -31,8 +31,8 @@ function op_InterpolateC2F!(c, f, bcs) @. f.x = interp(c.y) return nothing end -function op_LeftBiasedC2F!(c, f, bcs) - interp = Operators.LeftBiasedC2F(bcs) +function op_LeftBiased1stOrderC2F!(c, f, bcs) + interp = Operators.LeftBiased1stOrderC2F(bcs) @. f.x = interp(c.y) return nothing end @@ -41,8 +41,8 @@ function op_LeftBiasedF2C!(c, f, bcs = ()) @. c.x = interp(f.y) return nothing end -function op_RightBiasedC2F!(c, f, bcs) - interp = Operators.RightBiasedC2F(bcs) +function op_RightBiased1stOrderC2F!(c, f, bcs) + interp = Operators.RightBiased1stOrderC2F(bcs) @. f.x = interp(c.y) return nothing end @@ -162,9 +162,9 @@ GradientF2C GradientC2F InterpolateF2C InterpolateC2F -LeftBiasedC2F +LeftBiased1stOrderC2F LeftBiasedF2C -RightBiasedC2F +RightBiased1stOrderC2F RightBiasedF2C # Additional operators diff --git a/test/Operators/finitedifference/column_benchmark_utils.jl b/test/Operators/finitedifference/column_benchmark_utils.jl index b9467127f3..7e1bd62606 100644 --- a/test/Operators/finitedifference/column_benchmark_utils.jl +++ b/test/Operators/finitedifference/column_benchmark_utils.jl @@ -187,18 +187,18 @@ function set_vec_value_bcs(c) end function set_upwind_biased_bcs(c) - return (;bottom = Operators.FirstOrderOneSided(), - top = Operators.FirstOrderOneSided()) + return (;bottom = Operators.OneSided1stOrder(), + top = Operators.OneSided1stOrder()) end function set_upwind_biased_bcs(c) - return (;bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided()) + return (;bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder()) end function set_upwind_biased_3_bcs(c) - return (;bottom = Operators.ThirdOrderOneSided(), - top = Operators.ThirdOrderOneSided()) + return (;bottom = Operators.OneSided3rdOrder(), + top = Operators.OneSided3rdOrder()) end function set_top_value_bc(c) @@ -289,9 +289,9 @@ 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!)) = ((), ) bcs_tested(c, ::typeof(op_InterpolateC2F!)) = (set_value_bcs(c), extrapolate_bcs(c)) -bcs_tested(c, ::typeof(op_LeftBiasedC2F!)) = (set_bot_value_bc(c),) +bcs_tested(c, ::typeof(op_LeftBiased1stOrderC2F!)) = (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_RightBiased1stOrderC2F!)) = (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)) @@ -381,9 +381,9 @@ function benchmark_operators_base(trials, t_ave, cfield, ffield, h_space) op_broadcast_example0!, op_broadcast_example1!, op_broadcast_example2!, - op_LeftBiasedC2F!, + op_LeftBiased1stOrderC2F!, op_LeftBiasedF2C!, - op_RightBiasedC2F!, + op_RightBiased1stOrderC2F!, op_RightBiasedF2C!, op_CurlC2F!, #### Mixed / adaptive @@ -428,17 +428,17 @@ function test_results(t_ave) @test t_ave[(:no_h_space, op_InterpolateC2F!, :SetValue, :SetValue)] < 404.5473*ns*buffer @test t_ave[(:no_h_space, op_InterpolateC2F!, :Extrapolate, :Extrapolate)] < 400.4656*ns*buffer @test t_ave[(:no_h_space, op_broadcast_example2!, :none)] < 600*ns*buffer - @test t_ave[(:no_h_space, op_LeftBiasedC2F!, :SetValue)] < 365.2909*ns*buffer + @test t_ave[(:no_h_space, op_LeftBiased1stOrderC2F!, :SetValue)] < 365.2909*ns*buffer @test t_ave[(:no_h_space, op_LeftBiasedF2C!, :none)] < 185.358*ns*buffer @test t_ave[(:no_h_space, op_LeftBiasedF2C!, :SetValue)] < 221.175*ns*buffer - @test t_ave[(:no_h_space, op_RightBiasedC2F!, :SetValue)] < 138.649*ns*buffer + @test t_ave[(:no_h_space, op_RightBiased1stOrderC2F!, :SetValue)] < 138.649*ns*buffer @test t_ave[(:no_h_space, op_RightBiasedF2C!, :none)] < 186.417*ns*buffer @test t_ave[(:no_h_space, op_RightBiasedF2C!, :SetValue)] < 189.139*ns*buffer @test t_ave[(:no_h_space, op_CurlC2F!, :SetCurl, :SetCurl)] < 2.884*μs*buffer @test t_ave[(:no_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 2.926*μs*buffer @test t_ave[(:no_h_space, op_UpwindBiasedProductC2F!, :SetValue, :SetValue)] < 697.341*ns*buffer @test t_ave[(:no_h_space, op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate)] < 659.267*ns*buffer - @test t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 4.483*μs*buffer + @test t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 4.483*μs*buffer @test t_ave[(:no_h_space, op_divgrad_CC!, :SetValue, :SetValue, :none)] < 1.607*μs*buffer @test t_ave[(:no_h_space, op_divgrad_FF!, :none, :SetDivergence, :SetDivergence)] < 1.529*μs*buffer @test t_ave[(:no_h_space, op_div_interp_CC!, :SetValue, :SetValue, :none)] < 1.510*μs*buffer @@ -458,17 +458,17 @@ function test_results(t_ave) @test t_ave[(:has_h_space, op_broadcast_example0!, :none)] < 800*ns*buffer @test t_ave[(:has_h_space, op_broadcast_example1!, :none)] < 39*μs*buffer @test t_ave[(:has_h_space, op_broadcast_example2!, :none)] < 35*μs*buffer - @test t_ave[(:has_h_space, op_LeftBiasedC2F!, :SetValue)] < 619.749*ns*buffer + @test t_ave[(:has_h_space, op_LeftBiased1stOrderC2F!, :SetValue)] < 619.749*ns*buffer @test t_ave[(:has_h_space, op_LeftBiasedF2C!, :none)] < 276.520*ns*buffer @test t_ave[(:has_h_space, op_LeftBiasedF2C!, :SetValue)] < 333.901*ns*buffer - @test t_ave[(:has_h_space, op_RightBiasedC2F!, :SetValue)] < 245.966*ns*buffer + @test t_ave[(:has_h_space, op_RightBiased1stOrderC2F!, :SetValue)] < 245.966*ns*buffer @test t_ave[(:has_h_space, op_RightBiasedF2C!, :none)] < 277.616*ns*buffer @test t_ave[(:has_h_space, op_RightBiasedF2C!, :SetValue)] < 280.969*ns*buffer @test t_ave[(:has_h_space, op_CurlC2F!, :SetCurl, :SetCurl)] < 3.078*μs*buffer @test t_ave[(:has_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 3.159*μs*buffer @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :SetValue, :SetValue)] < 5.197*μs*buffer @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate)] < 5.304*μs*buffer - @test t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 14.304*μs*buffer + @test t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 14.304*μs*buffer @test t_ave[(:has_h_space, op_divgrad_CC!, :SetValue, :SetValue, :none)] < 8.593*μs*buffer @test t_ave[(:has_h_space, op_divgrad_FF!, :none, :SetDivergence, :SetDivergence)] < 8.597*μs*buffer @test t_ave[(:has_h_space, op_div_interp_CC!, :SetValue, :SetValue, :none)] < 1.735*μs*buffer @@ -479,7 +479,7 @@ function test_results(t_ave) # Broken tests @test_broken t_ave[(:no_h_space, op_CurlC2F!, :SetCurl, :SetCurl)] < 500 @test_broken t_ave[(:no_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 500 - @test_broken t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 500 + @test_broken t_ave[(:no_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 500 @test_broken t_ave[(:no_h_space, op_divgrad_uₕ!, :none, :SetValue, :Extrapolate)] < 500 @test_broken t_ave[(:no_h_space, op_divgrad_uₕ!, :none, :SetValue, :SetValue)] < 500 # different with/without h_space @test_broken t_ave[(:has_h_space, op_DivergenceF2C!, :none)] < 500 @@ -489,7 +489,7 @@ function test_results(t_ave) @test_broken t_ave[(:has_h_space, op_CurlC2F!, :SetValue, :SetValue)] < 500 @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :SetValue, :SetValue)] < 800 @test t_ave[(:has_h_space, op_UpwindBiasedProductC2F!, :Extrapolate, :Extrapolate)] < 800 - @test_broken t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :ThirdOrderOneSided, :ThirdOrderOneSided, :SetValue, :SetValue)] < 500 + @test_broken t_ave[(:has_h_space, op_divUpwind3rdOrderBiasedProductC2F!, :OneSided3rdOrder, :OneSided3rdOrder, :SetValue, :SetValue)] < 500 @test_broken t_ave[(:has_h_space, op_divgrad_CC!, :SetValue, :SetValue, :none)] < 500 @test_broken t_ave[(:has_h_space, op_divgrad_FF!, :none, :SetDivergence, :SetDivergence)] < 500 @test t_ave[(:has_h_space, op_div_interp_CC!, :SetValue, :SetValue, :none)] < 800 diff --git a/test/Operators/finitedifference/implicit_stencils_utils.jl b/test/Operators/finitedifference/implicit_stencils_utils.jl index 70f7304f36..d5014de149 100644 --- a/test/Operators/finitedifference/implicit_stencils_utils.jl +++ b/test/Operators/finitedifference/implicit_stencils_utils.jl @@ -118,8 +118,8 @@ function get_all_ops(center_space) OP.InterpolateC2F(bottom = extrap, top = extrap), OP.InterpolateC2F(bottom = zero_scalar, top = zero_scalar), OP.InterpolateC2F(bottom = zero_grad, top = zero_grad), - OP.LeftBiasedC2F(bottom = zero_scalar), - OP.RightBiasedC2F(top = zero_scalar), + OP.LeftBiased1stOrderC2F(bottom = zero_scalar), + OP.RightBiased1stOrderC2F(top = zero_scalar), ) ops_F2C_V2V = ( OP.InterpolateF2C(), @@ -143,8 +143,8 @@ function get_all_ops(center_space) ops_C2F_V2V = ( OP.InterpolateC2F(bottom = extrap, top = extrap), OP.InterpolateC2F(bottom = zero_vector, top = zero_vector), - OP.LeftBiasedC2F(bottom = zero_vector), - OP.RightBiasedC2F(top = zero_vector), + OP.LeftBiased1stOrderC2F(bottom = zero_vector), + OP.RightBiased1stOrderC2F(top = zero_vector), CurriedTwoArgOperator( OP.FluxCorrectionF2F(bottom = extrap, top = extrap), a_FV, diff --git a/test/Operators/finitedifference/opt.jl b/test/Operators/finitedifference/opt.jl index 86f5f2f9ae..e84d7ae345 100644 --- a/test/Operators/finitedifference/opt.jl +++ b/test/Operators/finitedifference/opt.jl @@ -110,13 +110,13 @@ function opt_WeightedInterpolateC2F_Extrapolate(weights, center_field) return identity.(WI.(weights, center_field)) end -function opt_LeftBiasedC2F(center_field) - LB = Operators.LeftBiasedC2F(left = Operators.SetValue(0.0)) +function opt_LeftBiased1stOrderC2F(center_field) + LB = Operators.LeftBiased1stOrderC2F(left = Operators.SetValue(0.0)) return LB.(identity.(center_field)) end -function opt_RightBiasedC2F(center_field) - RB = Operators.RightBiasedC2F(right = Operators.SetValue(0.0)) +function opt_RightBiased1stOrderC2F(center_field) + RB = Operators.RightBiased1stOrderC2F(right = Operators.SetValue(0.0)) return RB.(identity.(center_field)) end @@ -270,8 +270,8 @@ end centers, ) - @test_opt opt_LeftBiasedC2F(centers) - @test_opt opt_RightBiasedC2F(centers) + @test_opt opt_LeftBiased1stOrderC2F(centers) + @test_opt opt_RightBiased1stOrderC2F(centers) @test_opt opt_UpwindBiasedProductC2F_SetValue( face_velocities, diff --git a/test/Operators/finitedifference/opt_examples.jl b/test/Operators/finitedifference/opt_examples.jl index 757a11648f..013fa8bed5 100644 --- a/test/Operators/finitedifference/opt_examples.jl +++ b/test/Operators/finitedifference/opt_examples.jl @@ -176,7 +176,7 @@ function alloc_test_nested_expressions_1(cfield, ffield) (; cx, cy, cz, cϕ, cψ) = cfield ∇c = Operators.DivergenceF2C() wvec = Geometry.WVector - LB = Operators.LeftBiasedC2F(; bottom = Operators.SetValue(1)) + LB = Operators.LeftBiased1stOrderC2F(; bottom = Operators.SetValue(1)) @. cz = cx * cy * ∇c(wvec(LB(cy))) * ∇c(wvec(LB(cx))) * cϕ * cψ # Compile first p = @allocated begin @. cz = cx * cy * ∇c(wvec(LB(cy))) * ∇c(wvec(LB(cx))) * cϕ * cψ @@ -189,7 +189,7 @@ function alloc_test_nested_expressions_2(cfield, ffield) (; cx, cy, cz, cϕ, cψ) = cfield ∇c = Operators.DivergenceF2C() wvec = Geometry.WVector - RB = Operators.RightBiasedC2F(; top = Operators.SetValue(1)) + RB = Operators.RightBiased1stOrderC2F(; top = Operators.SetValue(1)) @. cz = cx * cy * ∇c(wvec(RB(cy))) * ∇c(wvec(RB(cx))) * cϕ * cψ # Compile first p = @allocated begin @. cz = cx * cy * ∇c(wvec(RB(cy))) * ∇c(wvec(RB(cx))) * cϕ * cψ @@ -203,7 +203,7 @@ function alloc_test_nested_expressions_3(cfield, ffield) Ic = Operators.InterpolateF2C() ∇c = Operators.DivergenceF2C() wvec = Geometry.WVector - LB = Operators.LeftBiasedC2F(; bottom = Operators.SetValue(1)) + LB = Operators.LeftBiased1stOrderC2F(; bottom = Operators.SetValue(1)) #! format: off @. cz = cx * cy * ∇c(wvec(LB(Ic(fy) * cx))) * ∇c(wvec(LB(Ic(fy) * cx))) * cϕ * cψ # Compile first p = @allocated begin @@ -497,12 +497,12 @@ end alloc_test_c2f_interp( cfield, ffield, - Operators.LeftBiasedC2F(; bottom = Operators.SetValue(0)), + Operators.LeftBiased1stOrderC2F(; bottom = Operators.SetValue(0)), ) alloc_test_c2f_interp( cfield, ffield, - Operators.RightBiasedC2F(; top = Operators.SetValue(0)), + Operators.RightBiased1stOrderC2F(; top = Operators.SetValue(0)), ) alloc_test_derivative( diff --git a/test/Operators/finitedifference/opt_implicit_stencils.jl b/test/Operators/finitedifference/opt_implicit_stencils.jl index 40b2cf0522..68c91b852a 100644 --- a/test/Operators/finitedifference/opt_implicit_stencils.jl +++ b/test/Operators/finitedifference/opt_implicit_stencils.jl @@ -82,8 +82,8 @@ Operators.Operator2Stencil(op::CurriedTwoArgOperator) = ops_C2F_⍰2⍰..., Operators.InterpolateC2F(bottom = zero_scalar, top = zero_scalar), Operators.InterpolateC2F(bottom = zero_grad, top = zero_grad), - Operators.LeftBiasedC2F(bottom = zero_scalar), - Operators.RightBiasedC2F(top = zero_scalar), + Operators.LeftBiased1stOrderC2F(bottom = zero_scalar), + Operators.RightBiased1stOrderC2F(top = zero_scalar), ) ops_F2C_V2V = ( ops_F2C_⍰2⍰..., @@ -105,8 +105,8 @@ Operators.Operator2Stencil(op::CurriedTwoArgOperator) = ops_C2F_V2V = ( ops_C2F_⍰2⍰..., Operators.InterpolateC2F(bottom = zero_vector, top = zero_vector), - Operators.LeftBiasedC2F(bottom = zero_vector), - Operators.RightBiasedC2F(top = zero_vector), + Operators.LeftBiased1stOrderC2F(bottom = zero_vector), + Operators.RightBiased1stOrderC2F(top = zero_vector), CurriedTwoArgOperator( Operators.FluxCorrectionF2F(bottom = extrap, top = extrap), a_FV, @@ -191,12 +191,12 @@ Operators.Operator2Stencil(op::CurriedTwoArgOperator) = success_tabs = apply_all_stencils(apply_configs) - @test success_tabs[:RightBiasedC2F] == 0 + @test success_tabs[:RightBiased1stOrderC2F] == 0 @test success_tabs[:InterpolateC2F] == 0 @test success_tabs[:CurriedTwoArgOperator] == 0 @test success_tabs[:GradientF2C] == 0 @test success_tabs[:DivergenceF2C] == 0 - @test success_tabs[:LeftBiasedC2F] == 0 + @test success_tabs[:LeftBiased1stOrderC2F] == 0 @test success_tabs[:DivergenceC2F] == 0 @test success_tabs[:GradientC2F] == 0 @test success_tabs[:CurlC2F] == 0 @@ -252,8 +252,8 @@ Operators.Operator2Stencil(op::CurriedTwoArgOperator) = @test success_compose_tabs[(:CurriedTwoArgOperator, :DivergenceF2C)] == 0 @test success_compose_tabs[(:InterpolateC2F, :GradientF2C)] == 0 - @test success_compose_tabs[(:RightBiasedF2C, :LeftBiasedC2F)] == 0 - @test success_compose_tabs[(:DivergenceF2C, :RightBiasedC2F)] == 0 + @test success_compose_tabs[(:RightBiasedF2C, :LeftBiased1stOrderC2F)] == 0 + @test success_compose_tabs[(:DivergenceF2C, :RightBiased1stOrderC2F)] == 0 @test success_compose_tabs[(:DivergenceF2C, :CurriedTwoArgOperator)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :CurriedTwoArgOperator)] == 0 @test success_compose_tabs[(:DivergenceC2F, :LeftBiasedF2C)] == 0 @@ -264,55 +264,55 @@ Operators.Operator2Stencil(op::CurriedTwoArgOperator) = @test success_compose_tabs[(:DivergenceC2F, :GradientF2C)] == 0 @test success_compose_tabs[(:DivergenceC2F, :InterpolateF2C)] == 0 @test success_compose_tabs[(:InterpolateF2C, :InterpolateC2F)] == 0 - @test success_compose_tabs[(:LeftBiasedF2C, :LeftBiasedC2F)] == 0 + @test success_compose_tabs[(:LeftBiasedF2C, :LeftBiased1stOrderC2F)] == 0 @test success_compose_tabs[(:RightBiasedF2C, :CurriedTwoArgOperator)] == 0 - @test success_compose_tabs[(:RightBiasedC2F, :LeftBiasedF2C)] == 0 - @test success_compose_tabs[(:RightBiasedC2F, :CurriedTwoArgOperator)] == 0 + @test success_compose_tabs[(:RightBiased1stOrderC2F, :LeftBiasedF2C)] == 0 + @test success_compose_tabs[(:RightBiased1stOrderC2F, :CurriedTwoArgOperator)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :InterpolateF2C)] == 0 @test success_compose_tabs[(:RightBiasedF2C, :GradientC2F)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :GradientC2F)] == 0 @test success_compose_tabs[(:DivergenceC2F, :RightBiasedF2C)] == 0 - @test success_compose_tabs[(:InterpolateF2C, :LeftBiasedC2F)] == 0 + @test success_compose_tabs[(:InterpolateF2C, :LeftBiased1stOrderC2F)] == 0 @test success_compose_tabs[(:InterpolateC2F, :DivergenceF2C)] == 0 - @test success_compose_tabs[(:CurriedTwoArgOperator, :RightBiasedC2F)] == 0 - @test success_compose_tabs[(:LeftBiasedC2F, :GradientF2C)] == 0 - @test success_compose_tabs[(:LeftBiasedC2F, :InterpolateF2C)] == 0 + @test success_compose_tabs[(:CurriedTwoArgOperator, :RightBiased1stOrderC2F)] == 0 + @test success_compose_tabs[(:LeftBiased1stOrderC2F, :GradientF2C)] == 0 + @test success_compose_tabs[(:LeftBiased1stOrderC2F, :InterpolateF2C)] == 0 @test success_compose_tabs[(:DivergenceF2C, :GradientC2F)] == 0 @test success_compose_tabs[(:DivergenceC2F, :CurriedTwoArgOperator)] == 0 @test success_compose_tabs[(:DivergenceF2C, :DivergenceC2F)] == 0 @test success_compose_tabs[(:InterpolateF2C, :GradientC2F)] == 0 - @test success_compose_tabs[(:LeftBiasedC2F, :RightBiasedF2C)] == 0 + @test success_compose_tabs[(:LeftBiased1stOrderC2F, :RightBiasedF2C)] == 0 @test success_compose_tabs[(:DivergenceF2C, :CurlC2F)] == 0 - @test success_compose_tabs[(:RightBiasedC2F, :DivergenceF2C)] == 0 - @test success_compose_tabs[(:LeftBiasedC2F, :LeftBiasedF2C)] == 0 + @test success_compose_tabs[(:RightBiased1stOrderC2F, :DivergenceF2C)] == 0 + @test success_compose_tabs[(:LeftBiased1stOrderC2F, :LeftBiasedF2C)] == 0 @test success_compose_tabs[(:InterpolateC2F, :LeftBiasedF2C)] == 0 @test success_compose_tabs[(:RightBiasedF2C, :CurlC2F)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :LeftBiasedF2C)] == 0 @test success_compose_tabs[(:LeftBiasedF2C, :GradientC2F)] == 0 - @test success_compose_tabs[(:CurriedTwoArgOperator, :LeftBiasedC2F)] == 0 - @test success_compose_tabs[(:InterpolateF2C, :RightBiasedC2F)] == 0 + @test success_compose_tabs[(:CurriedTwoArgOperator, :LeftBiased1stOrderC2F)] == 0 + @test success_compose_tabs[(:InterpolateF2C, :RightBiased1stOrderC2F)] == 0 @test success_compose_tabs[(:InterpolateF2C, :DivergenceC2F)] == 0 @test success_compose_tabs[(:RightBiasedF2C, :DivergenceC2F)] == 0 - @test success_compose_tabs[(:RightBiasedC2F, :InterpolateF2C)] == 0 + @test success_compose_tabs[(:RightBiased1stOrderC2F, :InterpolateF2C)] == 0 @test success_compose_tabs[(:LeftBiasedF2C, :DivergenceC2F)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :InterpolateC2F)] == 0 @test success_compose_tabs[(:InterpolateC2F, :InterpolateF2C)] == 0 - @test success_compose_tabs[(:LeftBiasedF2C, :RightBiasedC2F)] == 0 + @test success_compose_tabs[(:LeftBiasedF2C, :RightBiased1stOrderC2F)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :DivergenceC2F)] == 0 - @test success_compose_tabs[(:DivergenceF2C, :LeftBiasedC2F)] == 0 - @test success_compose_tabs[(:RightBiasedC2F, :RightBiasedF2C)] == 0 + @test success_compose_tabs[(:DivergenceF2C, :LeftBiased1stOrderC2F)] == 0 + @test success_compose_tabs[(:RightBiased1stOrderC2F, :RightBiasedF2C)] == 0 @test success_compose_tabs[(:InterpolateF2C, :CurlC2F)] == 0 @test success_compose_tabs[(:LeftBiasedF2C, :CurlC2F)] == 0 @test success_compose_tabs[(:InterpolateC2F, :RightBiasedF2C)] == 0 @test success_compose_tabs[(:LeftBiasedF2C, :CurriedTwoArgOperator)] == 0 - @test success_compose_tabs[(:LeftBiasedC2F, :DivergenceF2C)] == 0 + @test success_compose_tabs[(:LeftBiased1stOrderC2F, :DivergenceF2C)] == 0 @test success_compose_tabs[(:DivergenceF2C, :InterpolateC2F)] == 0 @test success_compose_tabs[(:InterpolateF2C, :CurriedTwoArgOperator)] == 0 @test success_compose_tabs[(:InterpolateC2F, :CurriedTwoArgOperator)] == 0 - @test success_compose_tabs[(:LeftBiasedC2F, :CurriedTwoArgOperator)] == 0 - @test success_compose_tabs[(:RightBiasedC2F, :GradientF2C)] == 0 + @test success_compose_tabs[(:LeftBiased1stOrderC2F, :CurriedTwoArgOperator)] == 0 + @test success_compose_tabs[(:RightBiased1stOrderC2F, :GradientF2C)] == 0 @test success_compose_tabs[(:CurriedTwoArgOperator, :GradientF2C)] == 0 - @test success_compose_tabs[(:RightBiasedF2C, :RightBiasedC2F)] == 0 + @test success_compose_tabs[(:RightBiasedF2C, :RightBiased1stOrderC2F)] == 0 @test success_compose_tabs[(:DivergenceC2F, :DivergenceF2C)] == 0 # for k in keys(success_tabs) # for debugging diff --git a/test/Operators/hybrid/opt.jl b/test/Operators/hybrid/opt.jl index 253fe4925e..6b2fd97267 100644 --- a/test/Operators/hybrid/opt.jl +++ b/test/Operators/hybrid/opt.jl @@ -112,13 +112,13 @@ function opt_WeightedInterpolateC2F_Extrapolate(weights, center_field) return identity.(WI.(weights, center_field)) end -function opt_LeftBiasedC2F(center_field) - LB = Operators.LeftBiasedC2F(left = Operators.SetValue(0.0)) +function opt_LeftBiased1stOrderC2F(center_field) + LB = Operators.LeftBiased1stOrderC2F(left = Operators.SetValue(0.0)) return LB.(identity.(center_field)) end -function opt_RightBiasedC2F(center_field) - RB = Operators.RightBiasedC2F(right = Operators.SetValue(0.0)) +function opt_RightBiased1stOrderC2F(center_field) + RB = Operators.RightBiased1stOrderC2F(right = Operators.SetValue(0.0)) return RB.(identity.(center_field)) end @@ -302,8 +302,8 @@ end centers, ) - @test_opt opt_LeftBiasedC2F(centers) - @test_opt opt_RightBiasedC2F(centers) + @test_opt opt_LeftBiased1stOrderC2F(centers) + @test_opt opt_RightBiased1stOrderC2F(centers) @test_opt opt_UpwindBiasedProductC2F_SetValue( face_velocities, From 224eb621ad36de816874d38378ef5c6c27b2317b Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 18 Jul 2023 10:24:12 -0700 Subject: [PATCH 2/2] update some tests --- test/Operators/finitedifference/column.jl | 90 ++++++++----------- .../column_benchmark_kernels.jl | 15 +++- .../column_benchmark_utils.jl | 2 + 3 files changed, 54 insertions(+), 53 deletions(-) diff --git a/test/Operators/finitedifference/column.jl b/test/Operators/finitedifference/column.jl index 3f23a42e8a..bfa7523b97 100644 --- a/test/Operators/finitedifference/column.jl +++ b/test/Operators/finitedifference/column.jl @@ -630,7 +630,7 @@ end @test conv_curl_sin_f[1] ≤ conv_curl_sin_f[2] ≤ conv_curl_sin_f[3] end -@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform) periodic mesh, constant w" begin +@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform) periodic mesh, constant w" begin FT = Float64 n_elems_seq = 2 .^ (5, 6, 7, 8) @@ -651,17 +651,15 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) - # Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind3rdOrderC2F Center -> Face operator # Unitary, constant advective velocity w = Geometry.WVector.(ones(fs)) # c = sin(z), scalar field defined at the centers c = sin.(centers) - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F() - third_order_fluxsinᶠ = third_order_fluxᶠ.(w, c) - + U = Operators.Upwind3rdOrderC2F() divf2c = Operators.DivergenceF2C() - adv_wc = divf2c.(third_order_fluxsinᶠ) + adv_wc = divf2c.(w .* U.(w, c)) Δh[k] = cs.face_local_geometry.J[1] @@ -672,7 +670,7 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z) + # Upwind3rdOrderC2F conv, with f(z) = sin(z) @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 5e-4 @test conv_adv_wc[1] ≈ 3 atol = 0.1 @test conv_adv_wc[2] ≈ 3 atol = 0.1 @@ -680,7 +678,7 @@ end @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2] end -@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform) periodic mesh, varying sign w" begin +@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform) periodic mesh, varying sign w" begin FT = Float64 n_elems_seq = 2 .^ (5, 6, 7, 8) @@ -702,17 +700,15 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) faces = getproperty(Fields.coordinate_field(fs), :z) - # Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind3rdOrderC2F Center -> Face operator # w = cos(z), vertical velocity field defined at the faces w = Geometry.WVector.(cos.(faces)) # c = sin(z), scalar field defined at the centers c = sin.(centers) - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F() - third_order_fluxsinᶠ = third_order_fluxᶠ.(w, c) - + U = Operators.Upwind3rdOrderC2F() divf2c = Operators.DivergenceF2C() - adv_wc = divf2c.(third_order_fluxsinᶠ) + adv_wc = divf2c.(w .* U.(w, c)) Δh[k] = cs.face_local_geometry.J[1] @@ -724,14 +720,14 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z), w(z) = cos(z) + # Upwind3rdOrderC2F conv, with f(z) = sin(z), w(z) = cos(z) @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 4e-3 @test conv_adv_wc[1] ≈ 2 atol = 0.2 @test conv_adv_wc[2] ≈ 2 atol = 0.1 @test conv_adv_wc[3] ≈ 2 atol = 0.1 end -@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided1stOrder + DivergenceF2C SetValue BCs, constant w" begin +@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided1stOrder + DivergenceF2C SetValue BCs, constant w" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -752,7 +748,7 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) - # Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind3rdOrderC2F Center -> Face operator # Unitary, constant advective velocity w = Geometry.WVector.(ones(fs)) # c = sin(z), scalar field defined at the centers @@ -760,7 +756,7 @@ end c = (cos.(centers .- Δz / 2) .- cos.(centers .+ Δz / 2)) ./ Δz s = sin.(centers) - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( + U = Operators.Upwind3rdOrderC2F( bottom = Operators.OneSided1stOrder(), top = Operators.OneSided1stOrder(), ) @@ -774,7 +770,7 @@ end ), ) - adv_wc = divf2c.(third_order_fluxᶠ.(w, c)) + adv_wc = divf2c.(w .* U.(w, c)) Δh[k] = cs.face_local_geometry.J[1] @@ -784,7 +780,7 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z), w(z) = 1 + # Upwind3rdOrderC2F conv, with f(z) = sin(z), w(z) = 1 @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 0.2006 @test conv_adv_wc[1] ≈ 0.5 atol = 0.2 @test conv_adv_wc[2] ≈ 0.5 atol = 0.3 @@ -792,7 +788,7 @@ end end end -@testset "Upwind3rdOrderBiasedProductC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder + DivergenceF2C SetValue BCs, varying sign w" begin +@testset "Upwind3rdOrderC2F + DivergenceF2C on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder + DivergenceF2C SetValue BCs, varying sign w" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -814,13 +810,13 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) faces = getproperty(Fields.coordinate_field(fs), :z) - # Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind3rdOrderC2F Center -> Face operator # w = cos(z), vertical velocity field defined at the faces w = Geometry.WVector.(cos.(faces)) # c = sin(z), scalar field defined at the centers c = sin.(centers) - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( + U = Operators.Upwind3rdOrderC2F( bottom = Operators.OneSided3rdOrder(), top = Operators.OneSided3rdOrder(), ) @@ -829,7 +825,7 @@ end bottom = Operators.SetValue(Geometry.WVector(FT(0.0))), top = Operators.SetValue(Geometry.WVector(FT(0.0))), ) - adv_wc = divf2c.(third_order_fluxᶠ.(w, c)) + adv_wc = divf2c.(w .* U.(w, c)) Δh[k] = cs.face_local_geometry.J[1] # Errors @@ -840,7 +836,7 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z), w(z) = cos(z) + # Upwind3rdOrderC2F conv, with f(z) = sin(z), w(z) = cos(z) @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 2e-1 @test conv_adv_wc[1] ≈ 2 atol = 0.1 @test conv_adv_wc[2] ≈ 2 atol = 0.1 @@ -848,7 +844,7 @@ end end end -@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform) periodic mesh" begin +@testset "Simple FCT: lin combination of Upwind1stOrderC2F + Upwind3rdOrderC2F on (uniform) periodic mesh" begin FT = Float64 n_elems_seq = 2 .^ (5, 6, 7, 8) @@ -870,20 +866,17 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding) - # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind1stOrderC2F & Upwind3rdOrderC2F Center -> Face operator # Unitary, constant advective velocity w = Geometry.WVector.(ones(fs)) # c = sin(z), scalar field defined at the centers c = sin.(centers) - first_order_fluxᶠ = Operators.UpwindBiasedProductC2F() - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F() - first_order_fluxsinᶠ = first_order_fluxᶠ.(w, c) - third_order_fluxsinᶠ = third_order_fluxᶠ.(w, c) + U1 = Operators.Upwind1stOrderC2F() + U3 = Operators.Upwind3rdOrderC2F() divf2c = Operators.DivergenceF2C() - corrected_antidiff_flux = - @. divf2c(C * (third_order_fluxsinᶠ - first_order_fluxsinᶠ)) + corrected_antidiff_flux = @. divf2c(C * w * (U1(w, c) - U3(w, c))) adv_wc = @. divf2c.(first_order_fluxsinᶠ) + corrected_antidiff_flux Δh[k] = cs.face_local_geometry.J[1] @@ -895,7 +888,7 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z) + # Upwind3rdOrderC2F conv, with f(z) = sin(z) @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 5e-4 @test conv_adv_wc[1] ≈ 3 atol = 0.1 @test conv_adv_wc[2] ≈ 3 atol = 0.1 @@ -903,7 +896,7 @@ end @test conv_adv_wc[1] ≤ conv_adv_wc[2] ≤ conv_adv_wc[2] end -@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided1stOrder BCs" begin +@testset "Simple FCT: lin combination of Upwind1stOrderC2F + Upwind3rdOrderC2F on (uniform and stretched) non-periodic mesh, with OneSided1stOrder BCs" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -925,7 +918,7 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding) - # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind1stOrderC2F & Upwind3rdOrderC2F Center -> Face operator # Unitary, constant advective velocity w = Geometry.WVector.(ones(fs)) # c = sin(z), scalar field defined at the centers @@ -933,11 +926,11 @@ end c = (cos.(centers .- Δz / 2) .- cos.(centers .+ Δz / 2)) ./ Δz s = sin.(centers) - first_order_fluxᶠ = Operators.UpwindBiasedProductC2F( + U1 = Operators.Upwind1stOrderC2F( bottom = Operators.Extrapolate(), top = Operators.Extrapolate(), ) - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( + U3 = Operators.Upwind3rdOrderC2F( bottom = Operators.OneSided1stOrder(), top = Operators.OneSided1stOrder(), ) @@ -951,9 +944,7 @@ end ), ) - corrected_antidiff_flux = @. divf2c( - C * (third_order_fluxᶠ(w, c) - first_order_fluxᶠ(w, c)), - ) + corrected_antidiff_flux = @. divf2c(C * w * (U3(w, c) - U1(w, c))) adv_wc = @. divf2c.(first_order_fluxᶠ(w, c)) + corrected_antidiff_flux @@ -965,7 +956,7 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z) + # Upwind3rdOrderC2F conv, with f(z) = sin(z) @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 0.2006 @test conv_adv_wc[1] ≈ 0.5 atol = 0.2 @test conv_adv_wc[2] ≈ 0.5 atol = 0.3 @@ -973,7 +964,7 @@ end end end -@testset "Simple FCT: lin combination of UpwindBiasedProductC2F + Upwind3rdOrderBiasedProductC2F on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder BCs" begin +@testset "Simple FCT: lin combination of Upwind1stOrderC2F + Upwind3rdOrderC2F on (uniform and stretched) non-periodic mesh, with OneSided3rdOrder BCs" begin FT = Float64 n_elems_seq = 2 .^ (4, 6, 8, 10) stretch_fns = (Meshes.Uniform(), Meshes.ExponentialStretching(1.0)) @@ -995,17 +986,17 @@ end centers = getproperty(Fields.coordinate_field(cs), :z) C = FT(1.0) # flux-correction coefficient (falling back to third-order upwinding) - # UpwindBiasedProductC2F & Upwind3rdOrderBiasedProductC2F Center -> Face operator + # Upwind1stOrderC2F & Upwind3rdOrderC2F Center -> Face operator # Unitary, constant advective velocity w = Geometry.WVector.(ones(fs)) # c = sin(z), scalar field defined at the centers c = sin.(centers) - first_order_fluxᶠ = Operators.UpwindBiasedProductC2F( + U1 = Operators.Upwind1stOrderC2F( bottom = Operators.Extrapolate(), top = Operators.Extrapolate(), ) - third_order_fluxᶠ = Operators.Upwind3rdOrderBiasedProductC2F( + U3 = Operators.Upwind3rdOrderC2F( bottom = Operators.OneSided3rdOrder(), top = Operators.OneSided3rdOrder(), ) @@ -1014,11 +1005,8 @@ end bottom = Operators.SetValue(Geometry.WVector(FT(0.0))), top = Operators.SetValue(Geometry.WVector(FT(0.0))), ) - corrected_antidiff_flux = @. divf2c( - C * (third_order_fluxᶠ(w, c) - first_order_fluxᶠ(w, c)), - ) - adv_wc = - @. divf2c.(first_order_fluxᶠ(w, c)) + corrected_antidiff_flux + corrected_antidiff_flux = @. divf2c(C * w * (U3(w, c) - U1(w, c))) + adv_wc = @. divf2c.(U1(w, c)) + corrected_antidiff_flux Δh[k] = cs.face_local_geometry.J[1] # Errors @@ -1028,7 +1016,7 @@ end # Check convergence rate conv_adv_wc = convergence_rate(err_adv_wc, Δh) - # Upwind3rdOrderBiasedProductC2F conv, with f(z) = sin(z) + # Upwind3rdOrderC2F conv, with f(z) = sin(z) @test err_adv_wc[3] ≤ err_adv_wc[2] ≤ err_adv_wc[1] ≤ 5e-1 @test conv_adv_wc[1] ≈ 2.5 atol = 0.1 @test conv_adv_wc[2] ≈ 2.5 atol = 0.1 diff --git a/test/Operators/finitedifference/column_benchmark_kernels.jl b/test/Operators/finitedifference/column_benchmark_kernels.jl index d320a52a85..a5587e2da0 100644 --- a/test/Operators/finitedifference/column_benchmark_kernels.jl +++ b/test/Operators/finitedifference/column_benchmark_kernels.jl @@ -99,13 +99,24 @@ function op_divgrad_uₕ!(c, f, bcs) @. c.uₕ2 = div(f.y * grad(c.uₕ)) return nothing end -function op_divUpwind3rdOrderBiasedProductC2F!(c, f, bcs) +function op_divUpwind3rdOrderBiasedProductC2F_old!(c, f, bcs) upwind = Operators.Upwind3rdOrderBiasedProductC2F(bcs.inner) divf2c = Operators.DivergenceF2C(bcs.outer) @. c.y = divf2c(upwind(f.w, c.x)) return nothing end - +function op_divUpwind1sOrderBiasedProductC2F!(c, f, bcs) + upwind = Operators.Upwind1stOrderC2F(bcs.inner) + divf2c = Operators.DivergenceF2C(bcs.outer) + @. c.y = divf2c(f.w * upwind(f.w, c.x)) + return nothing +end +function op_divUpwind3rdOrderBiasedProductC2F!(c, f, bcs) + upwind = Operators.Upwind3rdOrderC2F(bcs.inner) + divf2c = Operators.DivergenceF2C(bcs.outer) + @. c.y = divf2c(f.w * upwind(f.w, c.x)) + return nothing +end function op_broadcast_example0!(c, f, bcs) Fields.bycolumn(axes(f.ᶠu³)) do colidx CT3 = Geometry.Contravariant3Vector diff --git a/test/Operators/finitedifference/column_benchmark_utils.jl b/test/Operators/finitedifference/column_benchmark_utils.jl index 7e1bd62606..9c9fc14527 100644 --- a/test/Operators/finitedifference/column_benchmark_utils.jl +++ b/test/Operators/finitedifference/column_benchmark_utils.jl @@ -298,6 +298,8 @@ bcs_tested(c, ::typeof(op_UpwindBiasedProductC2F!)) = (set_value_bcs(c), extrapo bcs_tested(c, ::typeof(op_Upwind3rdOrderBiasedProductC2F!)) = (set_upwind_biased_3_bcs(c), extrapolate_bcs(c)) # Composed operators (bcs handled case-by-case) +bcs_tested(c, ::typeof(op_divUpwind1stOrderBiasedProductC2F!)) = + ((; inner = set_upwind_biased_3_bcs(c), outer = set_value_contra3_bcs(c)), ) 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!)) =