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.