Mercurial > repos > public > sbplib_julia
changeset 1026:29ac38ba0744 feature/dissipation_operators
Add StencilOperatorDistinctClosures
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Tue, 22 Mar 2022 08:02:23 +0100 |
parents | e74c41c4b60e |
children | d44cb983cfa4 |
files | src/SbpOperators/SbpOperators.jl src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl |
diffstat | 2 files changed, 59 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/src/SbpOperators/SbpOperators.jl Mon Mar 21 15:12:59 2022 +0100 +++ b/src/SbpOperators/SbpOperators.jl Tue Mar 22 08:02:23 2022 +0100 @@ -23,6 +23,7 @@ include("stencil.jl") include("readoperator.jl") include("volumeops/volume_operator.jl") +include("volumeops/stencil_operator_distinct_closures.jl") include("volumeops/constant_interior_scaling_operator.jl") include("volumeops/derivatives/first_derivative.jl") include("volumeops/derivatives/second_derivative.jl")
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl Tue Mar 22 08:02:23 2022 +0100 @@ -0,0 +1,58 @@ +function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction) + op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure) + return LazyTensors.inflate(op, size(grid), direction) +end + +""" + StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1} + +A one dimensional stencil operator with separate closures for the two boundaries. + +`StencilOperatorDistinctClosures` can be contrasted to `VolumeOperator` in +that it has different closure stencils for the upper and lower boundary. +`VolumeOperator` uses the same closure for both boundaries. Having distinct +closures is useful for representing operators with skewed stencils like upwind +operators. + +See also: [`VolumeOperator`](@ref) +""" +struct StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1,1} + inner_stencil::Stencil{T,K} + lower_closure::NTuple{N,Stencil{T,L}} + upper_closure::NTuple{M,Stencil{T,L}} + size::Tuple{Int} +end + +function StencilOperatorDistinctClosures(grid::EquidistantGrid{1}, inner_stencil, lower_closure, upper_closure) + return StencilOperatorDistinctClosures(inner_stencil, Tuple(lower_closure), Tuple(upper_closure), size(grid)) +end + +lower_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = N +upper_closure_size(::StencilOperatorDistinctClosures{T,K,N,M}) where {T,K,N,M} = M + +LazyTensors.range_size(op::StencilOperatorDistinctClosures) = op.size +LazyTensors.domain_size(op::StencilOperatorDistinctClosures) = op.size + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Lower}) + return @inbounds apply_stencil(op.lower_closure[Int(i)], v, Int(i)) +end + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Interior}) + return apply_stencil(op.inner_stencil, v, Int(i)) +end + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i::Index{Upper}) + stencil_index = Int(i) - (op.size[1]-upper_closure_size(op)) + return @inbounds apply_stencil(op.upper_closure[stencil_index], v, Int(i)) +end + +function LazyTensors.apply(op::StencilOperatorDistinctClosures, v::AbstractVector, i) + if i <= lower_closure_size(op) + LazyTensors.apply(op, v, Index(i, Lower)) + elseif i > op.size[1]-upper_closure_size(op) + LazyTensors.apply(op, v, Index(i, Upper)) + else + LazyTensors.apply(op, v, Index(i, Interior)) + end +end +# TODO: Move this to LazyTensors when we have the region communication down.