Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 889:069e58fb3829 feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor as multidimensional operator
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 07 Feb 2022 20:38:07 +0100
parents d228d1b26729
children f72cc96a58c6
comparison
equal deleted inserted replaced
888:e6d8fd5e8268 889:069e58fb3829
28 """ 28 """
29 SecondDerivativeVariable{T,N,M,K} <: TensorOperator{T,1} 29 SecondDerivativeVariable{T,N,M,K} <: TensorOperator{T,1}
30 30
31 Implements the one-dimensional second derivative with variable coefficients. 31 Implements the one-dimensional second derivative with variable coefficients.
32 """ 32 """
33 struct SecondDerivativeVariable{T,N,M,K,TArray<:AbstractVector} <: TensorMapping{T,1,1} 33 struct SecondDerivativeVariable{Dir,T,D,N,M,K,TArray<:AbstractArray} <: TensorMapping{T,D,D}
34 inner_stencil::NestedStencil{T,N} 34 inner_stencil::NestedStencil{T,N}
35 closure_stencils::NTuple{M,NestedStencil{T,K}} 35 closure_stencils::NTuple{M,NestedStencil{T,K}}
36 size::NTuple{1,Int} 36 size::NTuple{D,Int}
37 coefficient::TArray 37 coefficient::TArray
38
39 function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T,N}, closure_stencils::NTuple{M,NestedStencil{T,K}}, size::NTuple{D,Int}, coefficient::TArray) where {Dir,T,D,N,M,K,TArray<:AbstractArray}
40 return new{Dir,T,D,N,M,K,TArray}(inner_stencil,closure_stencils,size, coefficient)
41 end
42 end
43
44 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir)
45 return SecondDerivativeVariable{dir, dimension(grid)}(inner_stencil, Tuple(closure_stencils), size(grid), coeff)
38 end 46 end
39 47
40 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils) 48 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils)
41 return SecondDerivativeVariable(inner_stencil, Tuple(closure_stencils), size(grid), coeff) 49 return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1)
42 end 50 end
51
52 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
43 53
44 closure_size(::SecondDerivativeVariable{T,N,M}) where {T,N,M} = M 54 closure_size(::SecondDerivativeVariable{T,N,M}) where {T,N,M} = M
45 55
46 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size 56 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size
47 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size 57 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size
48 58
49 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Lower}) where T 59
50 return @inbounds apply_stencil(op.closure_stencils[Int(i)], op.coefficient, v, Int(i)) 60 function derivative_view(op, a, I)
61 d = derivative_direction(op)
62
63 Iview = Base.setindex(I,:,d)
64 return @view a[Iview...]
65
66 # D = domain_dim(op)
67 # Iₗ, _, Iᵣ = split_tuple(I, Val(d-1), Val(1), Val(D-d))
68 # return @view a[Iₗ..., :, Iᵣ...]
51 end 69 end
52 70
53 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Interior}) where T 71 function apply_lower(op::SecondDerivativeVariable, v, I...)
54 return apply_stencil(op.inner_stencil, op.coefficient, v, Int(i)) 72 ṽ = derivative_view(op, v, I)
73 c̃ = derivative_view(op, op.coefficient, I)
74
75 i = I[derivative_direction(op)]
76 return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i)
55 end 77 end
56 78
57 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i::Index{Upper}) where T 79 function apply_interior(op::SecondDerivativeVariable, v, I...)
58 return @inbounds apply_stencil_backwards(op.closure_stencils[op.size[1]-Int(i)+1], op.coefficient, v, Int(i)) 80 ṽ = derivative_view(op, v, I)
81 c̃ = derivative_view(op, op.coefficient, I)
82
83 i = I[derivative_direction(op)]
84 return apply_stencil(op.inner_stencil, c̃, ṽ, i)
59 end 85 end
60 86
61 function LazyTensors.apply(op::SecondDerivativeVariable{T}, v::AbstractVector{T}, i) where T 87 function apply_upper(op::SecondDerivativeVariable, v, I...)
88 ṽ = derivative_view(op, v, I)
89 c̃ = derivative_view(op, op.coefficient, I)
90
91 i = I[derivative_direction(op)]
92 return @inbounds apply_stencil_backwards(op.closure_stencils[op.size[1]-i+1], c̃, ṽ, i)
93 end
94
95 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractVector, I::Vararg{Index})
96 if I[derivative_direction(op)] isa Index{Lower}
97 return apply_lower(op, v, Int.(I)...)
98 elseif I[derivative_direction(op)] isa Index{Upper}
99 return apply_upper(op, v, Int.(I)...)
100 else
101 return apply_interior(op, v, Int.(I)...)
102 end
103 end
104
105 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractVector, I...)
106 i = I[derivative_direction(op)]
62 r = getregion(i, closure_size(op), op.size[1]) 107 r = getregion(i, closure_size(op), op.size[1])
63 return LazyTensors.apply(op, v, Index(i, r)) 108 return LazyTensors.apply(op, v, Index(i, r))
64 end 109 end
65 110
66 # TODO: Rename SecondDerivativeVariable -> VariableSecondDerivative 111 # TODO: Rename SecondDerivativeVariable -> VariableSecondDerivative