annotate src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 906:33c7e266e1a9 feature/variable_derivatives

Add fix for type instability in SecondDerivativeVariable along with tests
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 15 Feb 2022 09:56:07 +0100
parents 5bbc3ea3821b
children cad3a9f82009
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 export SecondDerivativeVariable
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 # """
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 # SecondDerivativeVariable(grid, inner_stencil, closure_stencils, parity, direction)
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 # Creates a volume operator on a `Dim`-dimensional grid acting along the
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 # specified coordinate `direction`. The action of the operator is determined by
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 # the stencils `inner_stencil` and `closure_stencils`. When `Dim=1`, the
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 # corresponding `SecondDerivativeVariable` tensor mapping is returned. When `Dim>1`, the
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 # returned operator is the appropriate outer product of a one-dimensional
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 # operators and `IdentityMapping`s, e.g for `Dim=3` the volume operator in the
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 # y-direction is `I⊗op⊗I`.
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 # """
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 # function volume_operator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity, direction)
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 # #TODO: Check that direction <= Dim?
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 # # Create 1D volume operator in along coordinate direction
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 # op = SecondDerivativeVariable(restrict(grid, direction), inner_stencil, closure_stencils, parity)
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 # # Create 1D IdentityMappings for each coordinate direction
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 # one_d_grids = restrict.(Ref(grid), Tuple(1:dimension(grid)))
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 # Is = IdentityMapping{eltype(grid)}.(size.(one_d_grids))
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 # # Formulate the correct outer product sequence of the identity mappings and
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 # # the volume operator
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 # parts = Base.setindex(Is, op, direction)
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 # return foldl(⊗, parts)
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 # end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 """
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29 SecondDerivativeVariable{T,N,M,K} <: TensorOperator{T,1}
882
9098fc936776 Add the coefficient as a part of the struct. Wrap tests in testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 881
diff changeset
30
9098fc936776 Add the coefficient as a part of the struct. Wrap tests in testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 881
diff changeset
31 Implements the one-dimensional second derivative with variable coefficients.
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 """
906
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
33 struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: TensorMapping{T,D,D}
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
34 inner_stencil::IStencil
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
35 closure_stencils::NTuple{M,CStencil}
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
36 size::NTuple{D,Int}
882
9098fc936776 Add the coefficient as a part of the struct. Wrap tests in testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 881
diff changeset
37 coefficient::TArray
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
38
906
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
39 function SecondDerivativeVariable{Dir, D}(inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, size::NTuple{D,Int}, coefficient::AbstractArray) where {Dir,T,D,M}
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
40 IStencil = typeof(inner_stencil)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
41 CStencil = eltype(closure_stencils)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
42 TArray = typeof(coefficient)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
43 return new{Dir,T,D,M,IStencil,CStencil,TArray}(inner_stencil,closure_stencils,size, coefficient)
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
44 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
45 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
46
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
47 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir)
901
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
48 Δxᵢ = spacing(grid)[dir]
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
49 scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2)
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
50 scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2)
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
51 return SecondDerivativeVariable{dir, dimension(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53
882
9098fc936776 Add the coefficient as a part of the struct. Wrap tests in testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 881
diff changeset
54 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils)
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
55 return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, 1)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57
901
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
58 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir)
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
59 inner_stencil = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"])
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
60 closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"])
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
61
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
62 return SecondDerivativeVariable(grid, coeff, inner_stencil, closure_stencils, dir)
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
63 end
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
64
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
65 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
66
891
f72cc96a58c6 Add some size tests
Jonatan Werpers <jonatan@werpers.com>
parents: 889
diff changeset
67 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
72
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
73 function derivative_view(op, a, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
74 d = derivative_direction(op)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
75
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
76 Iview = Base.setindex(I,:,d)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
77 return @view a[Iview...]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
78
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
79 # D = domain_dim(op)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
80 # Iₗ, _, Iᵣ = split_tuple(I, Val(d-1), Val(1), Val(D-d))
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
81 # return @view a[Iₗ..., :, Iᵣ...]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
82 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
83
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
84 function apply_lower(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
85 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
86 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
87
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
88 i = I[derivative_direction(op)]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
89 return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
92 function apply_interior(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
93 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
94 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
95
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
96 i = I[derivative_direction(op)]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
97 return apply_stencil(op.inner_stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
100 function apply_upper(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
101 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
102 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
103
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
104 i = I[derivative_direction(op)]
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
105 stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1]
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
106 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
109 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I::Vararg{Index})
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
110 if I[derivative_direction(op)] isa Index{Lower}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
111 return apply_lower(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
112 elseif I[derivative_direction(op)] isa Index{Upper}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
113 return apply_upper(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
114 else
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
115 return apply_interior(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
116 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
117 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
118
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
119 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
120 dir = derivative_direction(op)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
121
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
122 i = I[dir]
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
123 r = getregion(i, closure_size(op), op.size[dir])
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
124
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
125 I = map(i->Index(i, Interior), I)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
126 I = Base.setindex(I, Index(i, r), dir)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
127 return LazyTensors.apply(op, v, I...)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
128 end
884
d228d1b26729 Add tests for application
Jonatan Werpers <jonatan@werpers.com>
parents: 882
diff changeset
129
d228d1b26729 Add tests for application
Jonatan Werpers <jonatan@werpers.com>
parents: 882
diff changeset
130 # TODO: Rename SecondDerivativeVariable -> VariableSecondDerivative