Mercurial > repos > public > sbplib_julia
changeset 1719:d81d0660b80d
Merge refactor/grids/boundary_identifiers_1d
author | Vidar Stiernström <vidar.stiernstrom@gmail.com> |
---|---|
date | Thu, 05 Sep 2024 08:26:12 -0700 |
parents | fe298044e044 (current diff) 351cf866a3d8 (diff) |
children | 7fdc348f91ff 8317252e4535 |
files | |
diffstat | 12 files changed, 103 insertions(+), 85 deletions(-) [+] |
line wrap: on
line diff
--- a/benchmark/benchmarks.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/benchmark/benchmarks.jl Thu Sep 05 08:26:12 2024 -0700 @@ -3,7 +3,6 @@ using Sbplib using Sbplib.Grids using Sbplib.SbpOperators -using Sbplib.RegionIndices using Sbplib.LazyTensors using LinearAlgebra @@ -168,20 +167,20 @@ Dxx = second_derivative(g2, stencil_set, 1) Dyy = second_derivative(g2, stencil_set, 2) -e₁ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,Lower}()) -e₁ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,Upper}()) -e₂ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,Lower}()) -e₂ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,Upper}()) +e₁ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,LowerBoundary}()) +e₁ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{1,UpperBoundary}()) +e₂ₗ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,LowerBoundary}()) +e₂ᵤ = boundary_restriction(g2, stencil_set, CartesianBoundary{2,UpperBoundary}()) -d₁ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{1,Lower}()) -d₁ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{1,Upper}()) -d₂ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{2,Lower}()) -d₂ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{2,Upper}()) +d₁ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{1,LowerBoundary}()) +d₁ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{1,UpperBoundary}()) +d₂ₗ = normal_derivative(g2, stencil_set, CartesianBoundary{2,LowerBoundary}()) +d₂ᵤ = normal_derivative(g2, stencil_set, CartesianBoundary{2,UpperBoundary}()) -H₁ₗ = inner_product(boundary_grid(g2, CartesianBoundary{1,Lower}()), stencil_set) -H₁ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{1,Upper}()), stencil_set) -H₂ₗ = inner_product(boundary_grid(g2, CartesianBoundary{2,Lower}()), stencil_set) -H₂ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{2,Upper}()), stencil_set) +H₁ₗ = inner_product(boundary_grid(g2, CartesianBoundary{1,LowerBoundary}()), stencil_set) +H₁ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{1,UpperBoundary}()), stencil_set) +H₂ₗ = inner_product(boundary_grid(g2, CartesianBoundary{2,LowerBoundary}()), stencil_set) +H₂ᵤ = inner_product(boundary_grid(g2, CartesianBoundary{2,UpperBoundary}()), stencil_set) SUITE["boundary_terms"]["pre_composition"] = @benchmarkable $u2 .= $(H⁻¹∘e₁ₗ'∘H₁ₗ∘d₁ₗ)*$v2 SUITE["boundary_terms"]["composition"] = @benchmarkable $u2 .= ($H⁻¹∘$e₁ₗ'∘$H₁ₗ∘$d₁ₗ)*$v2
--- a/src/Grids/Grids.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/src/Grids/Grids.jl Thu Sep 05 08:26:12 2024 -0700 @@ -1,6 +1,5 @@ module Grids -using Sbplib.RegionIndices using Sbplib.LazyTensors using StaticArrays @@ -23,6 +22,8 @@ export BoundaryIdentifier export TensorGridBoundary export CartesianBoundary +export LowerBoundary +export UpperBoundary export TensorGrid export ZeroDimGrid @@ -32,9 +33,6 @@ export spacing export equidistant_grid - -abstract type BoundaryIdentifier end - include("grid.jl") include("tensor_grid.jl") include("equidistant_grid.jl")
--- a/src/Grids/equidistant_grid.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/src/Grids/equidistant_grid.jl Thu Sep 05 08:26:12 2024 -0700 @@ -49,12 +49,30 @@ min_spacing(g::EquidistantGrid) = spacing(g) +""" + LowerBoundary <: BoundaryIdentifier -boundary_identifiers(::EquidistantGrid) = (Lower(), Upper()) -boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin]) -boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end]) -boundary_indices(g::EquidistantGrid, id::Lower) = (1,) -boundary_indices(g::EquidistantGrid, id::Upper) = (length(g),) +Boundary identifier for the the lower (left) boundary of a one-dimensional grid. + +See also: [`BoundaryIdentifier`](@ref) +""" +struct LowerBoundary <: BoundaryIdentifier end + +""" + UpperBoundary <: BoundaryIdentifier + +Boundary identifier for the the upper (right) boundary of a one-dimensional grid. + +See also: [`BoundaryIdentifier`](@ref) +""" +struct UpperBoundary <: BoundaryIdentifier end + + +boundary_identifiers(::EquidistantGrid) = (LowerBoundary(), UpperBoundary()) +boundary_grid(g::EquidistantGrid, id::LowerBoundary) = ZeroDimGrid(g[begin]) +boundary_grid(g::EquidistantGrid, id::UpperBoundary) = ZeroDimGrid(g[end]) +boundary_indices(g::EquidistantGrid, id::LowerBoundary) = (firstindex(g),) +boundary_indices(g::EquidistantGrid, id::UpperBoundary) = (lastindex(g),) """ refine(g::EquidistantGrid, r::Int)
--- a/src/Grids/grid.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/src/Grids/grid.jl Thu Sep 05 08:26:12 2024 -0700 @@ -109,6 +109,13 @@ function coarsen end """ + BoundaryIdentifier + +An identifier for a boundary of a grid. +""" +abstract type BoundaryIdentifier end + +""" boundary_identifiers(g::Grid) Identifiers for all the boundaries of `g`.
--- a/src/SbpOperators/boundaryops/boundary_operator.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Thu Sep 05 08:26:12 2024 -0700 @@ -1,27 +1,27 @@ """ - BoundaryOperator{T,R,N} <: LazyTensor{T,0,1} + BoundaryOperator{T,B,N} <: LazyTensor{T,0,1} Implements the boundary operator `op` for 1D as a `LazyTensor` `op` is the restriction of a grid function to the boundary using some closure -`Stencil{T,N}`. The boundary to restrict to is determined by `R`. `op'` is the +`Stencil{T,N}`. The boundary to restrict to is determined by `B`. `op'` is the prolongation of a zero dimensional array to the whole grid using the same closure stencil. """ -struct BoundaryOperator{T,R<:Region,N} <: LazyTensor{T,0,1} +struct BoundaryOperator{T,B<:BoundaryIdentifier,N} <: LazyTensor{T,0,1} stencil::Stencil{T,N} size::Int end """ - BoundaryOperator(grid::EquidistantGrid, closure_stencil, region) + BoundaryOperator(grid::EquidistantGrid, closure_stencil, boundary) Constructs the BoundaryOperator with stencil `closure_stencil` for a `EquidistantGrid` `grid`, restricting to to the boundary specified by -`region`. +`boundary`. """ -function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, region::Region) where {T,N} - return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1]) +function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, boundary::BoundaryIdentifier) where {T,N} + return BoundaryOperator{T,typeof(boundary),N}(closure_stencil,size(grid)[1]) end """ @@ -29,24 +29,24 @@ The size of the closure stencil. """ -closure_size(::BoundaryOperator{T,R,N}) where {T,R,N} = N +closure_size(::BoundaryOperator{T,B,N}) where {T,B,N} = N LazyTensors.range_size(op::BoundaryOperator) = () LazyTensors.domain_size(op::BoundaryOperator) = (op.size,) -function LazyTensors.apply(op::BoundaryOperator{<:Any,Lower}, v::AbstractVector) +function LazyTensors.apply(op::BoundaryOperator{<:Any,LowerBoundary}, v::AbstractVector) apply_stencil(op.stencil,v,1) end -function LazyTensors.apply(op::BoundaryOperator{<:Any,Upper}, v::AbstractVector) +function LazyTensors.apply(op::BoundaryOperator{<:Any,UpperBoundary}, v::AbstractVector) apply_stencil_backwards(op.stencil,v,op.size) end -function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Lower}, v::AbstractArray{<:Any,0}, i::Index{Lower}) +function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,LowerBoundary}, v::AbstractArray{<:Any,0}, i::Index{Lower}) return op.stencil[Int(i)-1]*v[] end -function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,Upper}, v::AbstractArray{<:Any,0}, i::Index{Upper}) +function LazyTensors.apply_transpose(op::BoundaryOperator{<:Any,UpperBoundary}, v::AbstractArray{<:Any,0}, i::Index{Upper}) return op.stencil[op.size[1] - Int(i)]*v[] end
--- a/src/SbpOperators/volumeops/laplace/laplace.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Thu Sep 05 08:26:12 2024 -0700 @@ -69,7 +69,7 @@ Hᵧ = inner_product(boundary_grid(g, id), set) e = boundary_restriction(g, set, id) d = normal_derivative(g, set, id) - B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning) + B = positivity_decomposition(Δ, g, boundary(bc); H_tuning, R_tuning) penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ return penalty_tensor, e end @@ -94,7 +94,7 @@ end """ - positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) + positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning) Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive definite with respect to the boundary quadrature. Here `d` is the normal @@ -104,16 +104,14 @@ penalty and must be greater than 1. For details we refer to <https://doi.org/10.1016/j.jcp.2020.109294> """ -function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning) +function positivity_decomposition(Δ::Laplace, g::Grid, b::BoundaryIdentifier; H_tuning, R_tuning) @assert(H_tuning ≥ 1.) @assert(R_tuning ≥ 1.) - Nτ_H, τ_R = positivity_limits(Δ,g,bc) + Nτ_H, τ_R = positivity_limits(Δ,g,b) return H_tuning*Nτ_H + R_tuning*τ_R end -# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then -# change bc::BoundaryCondition to id::BoundaryIdentifier -function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition) +function positivity_limits(Δ::Laplace, g::EquidistantGrid, b::BoundaryIdentifier) h = spacing(g) θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1]) θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"]) @@ -123,7 +121,7 @@ return τ_H, τ_R end -function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition) - τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc) +function positivity_limits(Δ::Laplace, g::TensorGrid, b::BoundaryIdentifier) + τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(b)], b) return τ_H*ndims(g), τ_R end
--- a/test/Grids/equidistant_grid_test.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/test/Grids/equidistant_grid_test.jl Thu Sep 05 08:26:12 2024 -0700 @@ -1,6 +1,5 @@ using Sbplib.Grids using Test -using Sbplib.RegionIndices using Sbplib.LazyTensors @@ -63,24 +62,24 @@ @testset "boundary_identifiers" begin g = EquidistantGrid(0:0.1:10) - @test boundary_identifiers(g) == (Lower(), Upper()) + @test boundary_identifiers(g) == (LowerBoundary(), UpperBoundary()) @inferred boundary_identifiers(g) end @testset "boundary_grid" begin g = EquidistantGrid(0:0.1:1) - @test boundary_grid(g, Lower()) == ZeroDimGrid(0.0) - @test boundary_grid(g, Upper()) == ZeroDimGrid(1.0) + @test boundary_grid(g, LowerBoundary()) == ZeroDimGrid(0.0) + @test boundary_grid(g, UpperBoundary()) == ZeroDimGrid(1.0) end @testset "boundary_indices" begin g = EquidistantGrid(0:0.1:1) - @test boundary_indices(g, Lower()) == (1,) - @test boundary_indices(g, Upper()) == (11,) + @test boundary_indices(g, LowerBoundary()) == (1,) + @test boundary_indices(g, UpperBoundary()) == (11,) g = EquidistantGrid(2:0.1:10) - @test boundary_indices(g, Lower()) == (1,) - @test boundary_indices(g, Upper()) == (81,) + @test boundary_indices(g, LowerBoundary()) == (1,) + @test boundary_indices(g, UpperBoundary()) == (81,) end
--- a/test/Grids/tensor_grid_test.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/test/Grids/tensor_grid_test.jl Thu Sep 05 08:26:12 2024 -0700 @@ -1,7 +1,6 @@ using Test using Sbplib.Grids using StaticArrays -using Sbplib.RegionIndices @testset "TensorGrid" begin g₁ = EquidistantGrid(range(0,1,length=11)) @@ -170,13 +169,13 @@ end @testset "boundary_identifiers" begin - @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (Lower,Upper,Lower,Upper)) - @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,Lower}(),TensorGridBoundary{1,Upper}()) + @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (LowerBoundary,UpperBoundary,LowerBoundary,UpperBoundary)) + @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,LowerBoundary}(),TensorGridBoundary{1,UpperBoundary}()) end @testset "boundary_grid" begin - @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂) - @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄) + @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, UpperBoundary}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂) + @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, UpperBoundary}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄) end @testset "boundary_indices" begin @@ -184,14 +183,14 @@ g₂ = EquidistantGrid(range(2,3,length=6)) g₄ = ZeroDimGrid(@SVector[1,2]) - @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Lower}()) == (1,:) - @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == (11,:) - @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Lower}()) == (:,1) - @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Upper}()) == (:,6) - @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Lower}()) == (1,) - @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == (11,) - @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Lower}()) == (1,) - @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Upper}()) == (11,) + @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, LowerBoundary}()) == (1,:) + @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, UpperBoundary}()) == (11,:) + @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, LowerBoundary}()) == (:,1) + @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, UpperBoundary}()) == (:,6) + @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, LowerBoundary}()) == (1,) + @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, UpperBoundary}()) == (11,) + @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, LowerBoundary}()) == (1,) + @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, UpperBoundary}()) == (11,) end end
--- a/test/SbpOperators/boundary_conditions/boundary_condition_test.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl Thu Sep 05 08:26:12 2024 -0700 @@ -15,8 +15,8 @@ g = 3.14 f(x,y,z) = x^2+y^2+z^2 @testset "Constructors" begin - @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower} - @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,Lower}} + @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,LowerBoundary} + @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,LowerBoundary}} end @testset "boundary" begin
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Thu Sep 05 08:26:12 2024 -0700 @@ -13,12 +13,12 @@ g_1D = EquidistantGrid(range(0,1,length=11)) @testset "Constructors" begin - @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T - @test BoundaryOperator(g_1D, closure_stencil, Upper()) isa LazyTensor{T,0,1} where T + @test BoundaryOperator(g_1D, closure_stencil, LowerBoundary()) isa LazyTensor{T,0,1} where T + @test BoundaryOperator(g_1D, closure_stencil, UpperBoundary()) isa LazyTensor{T,0,1} where T end - op_l = BoundaryOperator(g_1D, closure_stencil, Lower()) - op_r = BoundaryOperator(g_1D, closure_stencil, Upper()) + op_l = BoundaryOperator(g_1D, closure_stencil, LowerBoundary()) + op_r = BoundaryOperator(g_1D, closure_stencil, UpperBoundary()) @testset "Sizes" begin @test domain_size(op_l) == (11,)
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl Thu Sep 05 08:26:12 2024 -0700 @@ -14,19 +14,19 @@ @testset "boundary_restriction" begin @testset "1D" begin - e_l = boundary_restriction(g_1D,stencil_set,Lower()) - @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Lower()) - @test e_l isa BoundaryOperator{T,Lower} where T + e_l = boundary_restriction(g_1D,stencil_set,LowerBoundary()) + @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),LowerBoundary()) + @test e_l isa BoundaryOperator{T,LowerBoundary} where T @test e_l isa LazyTensor{T,0,1} where T - e_r = boundary_restriction(g_1D,stencil_set,Upper()) - @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Upper()) - @test e_r isa BoundaryOperator{T,Upper} where T + e_r = boundary_restriction(g_1D,stencil_set,UpperBoundary()) + @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),UpperBoundary()) + @test e_r isa BoundaryOperator{T,UpperBoundary} where T @test e_r isa LazyTensor{T,0,1} where T end @testset "2D" begin - e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}()) + e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,UpperBoundary}()) @test e_w isa InflatedTensor @test e_w isa LazyTensor{T,1,2} where T end
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl Thu Sep 05 08:22:47 2024 -0700 +++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl Thu Sep 05 08:26:12 2024 -0700 @@ -12,19 +12,19 @@ @testset "normal_derivative" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) @testset "1D" begin - d_l = normal_derivative(g_1D, stencil_set, Lower()) - @test d_l == normal_derivative(g_1D, stencil_set, Lower()) - @test d_l isa BoundaryOperator{T,Lower} where T + d_l = normal_derivative(g_1D, stencil_set, LowerBoundary()) + @test d_l == normal_derivative(g_1D, stencil_set, LowerBoundary()) + @test d_l isa BoundaryOperator{T,LowerBoundary} where T @test d_l isa LazyTensor{T,0,1} where T end @testset "2D" begin - d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) - d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}()) + d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}()) + d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,UpperBoundary}()) Ix = IdentityTensor{Float64}((size(g_2D)[1],)) Iy = IdentityTensor{Float64}((size(g_2D)[2],)) - d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower()) - d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper()) - @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) + d_l = normal_derivative(g_2D.grids[1], stencil_set, LowerBoundary()) + d_r = normal_derivative(g_2D.grids[2], stencil_set, UpperBoundary()) + @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,LowerBoundary}()) @test d_w == d_l⊗Iy @test d_n == Ix⊗d_r @test d_w isa LazyTensor{T,1,2} where T