Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 1355:102ebdaf7c11 feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 08 Feb 2023 21:21:28 +0100
parents 3bb94ce74697
children 4684c7f1c4cb
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
916
fc1f03fc8d65 Remove comments
Jonatan Werpers <jonatan@werpers.com>
parents: 913
diff changeset
3 # REVIEW: Fixa docs
fc1f03fc8d65 Remove comments
Jonatan Werpers <jonatan@werpers.com>
parents: 913
diff changeset
4 """
1049
3bb94ce74697 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 940
diff changeset
5 SecondDerivativeVariable{Dir,T,D,...} <: LazyTensor{T,D,D}
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6
916
fc1f03fc8d65 Remove comments
Jonatan Werpers <jonatan@werpers.com>
parents: 913
diff changeset
7 A second derivative operator in direction `Dir` with a variable coefficient.
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 """
1049
3bb94ce74697 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 940
diff changeset
9 struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: LazyTensor{T,D,D}
906
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
10 inner_stencil::IStencil
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
11 closure_stencils::NTuple{M,CStencil}
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
12 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
13 coefficient::TArray
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
14
906
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
15 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
16 IStencil = typeof(inner_stencil)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
17 CStencil = eltype(closure_stencils)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
18 TArray = typeof(coefficient)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
19 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
20 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
21 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
22
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
23 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, inner_stencil, closure_stencils, dir)
940
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
24 check_coefficient(grid, coeff)
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
25
901
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
26 Δxᵢ = spacing(grid)[dir]
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
27 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
28 scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2)
1355
102ebdaf7c11 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 1049
diff changeset
29 return SecondDerivativeVariable{dir, ndims(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
30 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31
937
025a506ca2fa Fix dispatch bug for 1D grids
Jonatan Werpers <jonatan@werpers.com>
parents: 932
diff changeset
32 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil::NestedStencil, closure_stencils)
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
33 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
34 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35
917
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
36 @doc raw"""
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
37 SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir)
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
38
1049
3bb94ce74697 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 940
diff changeset
39 Create a `LazyTensor` for the second derivative with a variable coefficient
917
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
40 `coeff` on `grid` from the stencils in `stencil_set`. The direction is
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
41 determined by `dir`.
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
42
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
43 `coeff` is a grid function on `grid`.
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
44
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
45 # Example
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
46 With
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
47 ```
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
48 D = SecondDerivativeVariable(g, c, stencil_set, 2)
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
49 ```
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
50 then `D*u` approximates
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
51 ```math
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
52 \frac{\partial}{\partial y} c(x,y) \frac{\partial u}{\partial y},
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
53 ```
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
54 on ``(0,1)⨯(0,1)`` represented by `g`.
eb054537fc63 More docs
Jonatan Werpers <jonatan@werpers.com>
parents: 916
diff changeset
55 """
940
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
56 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir::Int)
901
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
57 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
58 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
59
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
60 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
61 end
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
62
940
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
63 function check_coefficient(grid, coeff)
1355
102ebdaf7c11 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 1049
diff changeset
64 if ndims(grid) != ndims(coeff)
102ebdaf7c11 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 1049
diff changeset
65 throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(ndims(grid))"))
940
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
66 end
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
67
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
68 if size(grid) != size(coeff)
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
69 throw(DimensionMismatch("the size $(size(coeff)) of the coefficient does not match the size $(size(grid)) of the grid"))
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
70 end
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
71 end
b15f39ae1643 Add size checking for the coefficient
Jonatan Werpers <jonatan@werpers.com>
parents: 937
diff changeset
72
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
73 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
74
891
f72cc96a58c6 Add some size tests
Jonatan Werpers <jonatan@werpers.com>
parents: 889
diff changeset
75 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
80
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
81 function derivative_view(op, a, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
82 d = derivative_direction(op)
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 Iview = Base.setindex(I,:,d)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
85 return @view a[Iview...]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
86 end
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 function apply_lower(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
89 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
90 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
91
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
92 i = I[derivative_direction(op)]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
93 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
94 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
96 function apply_interior(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
97 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
98 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
99
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
100 i = I[derivative_direction(op)]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
101 return apply_stencil(op.inner_stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
104 function apply_upper(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
105 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
106 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
107
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
108 i = I[derivative_direction(op)]
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
109 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
110 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
113 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
114 if I[derivative_direction(op)] isa Index{Lower}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
115 return apply_lower(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
116 elseif I[derivative_direction(op)] isa Index{Upper}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
117 return apply_upper(op, v, Int.(I)...)
931
720b1358e06d More strict apply
Jonatan Werpers <jonatan@werpers.com>
parents: 930
diff changeset
118 elseif I[derivative_direction(op)] isa Index{Interior}
720b1358e06d More strict apply
Jonatan Werpers <jonatan@werpers.com>
parents: 930
diff changeset
119 return apply_interior(op, v, Int.(I)...)
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
120 else
931
720b1358e06d More strict apply
Jonatan Werpers <jonatan@werpers.com>
parents: 930
diff changeset
121 error("Invalid region")
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
122 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
123 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
124
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
125 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
126 dir = derivative_direction(op)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
127
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
128 i = I[dir]
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
129
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
130 I = map(i->Index(i, Interior), I)
930
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
131 if 0 < i <= closure_size(op)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
132 I = Base.setindex(I, Index(i, Lower), dir)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
133 return LazyTensors.apply(op, v, I...)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
134 elseif closure_size(op) < i <= op.size[dir]-closure_size(op)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
135 I = Base.setindex(I, Index(i, Interior), dir)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
136 return LazyTensors.apply(op, v, I...)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
137 elseif op.size[dir]-closure_size(op) < i <= op.size[dir]
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
138 I = Base.setindex(I, Index(i, Upper), dir)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
139 return LazyTensors.apply(op, v, I...)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
140 else
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
141 error("Bounds error") # TODO: Make this more standard
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
142 end
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
143 end
932
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
144
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
145
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
146 ## 2D Specific implementations to avoid instability
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
147 ## TODO: Should really be solved by fixing the general methods instead
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
148
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
149
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
150 ## x-direction
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
151 function apply_lower(op::SecondDerivativeVariable{1}, v, i, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
152 ṽ = @view v[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
153 c̃ = @view op.coefficient[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
154
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
155 return @inbounds apply_stencil(op.closure_stencils[i], c̃, ṽ, i)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
156 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
157
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
158 function apply_interior(op::SecondDerivativeVariable{1}, v, i, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
159 ṽ = @view v[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
160 c̃ = @view op.coefficient[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
161
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
162 return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, i)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
163 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
164
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
165 function apply_upper(op::SecondDerivativeVariable{1}, v, i, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
166 ṽ = @view v[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
167 c̃ = @view op.coefficient[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
168
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
169 stencil = op.closure_stencils[op.size[derivative_direction(op)]-i+1]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
170 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
171 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
172
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
173
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
174 ## y-direction
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
175 function apply_lower(op::SecondDerivativeVariable{2}, v, i, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
176 ṽ = @view v[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
177 c̃ = @view op.coefficient[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
178
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
179 return @inbounds apply_stencil(op.closure_stencils[j], c̃, ṽ, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
180 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
181
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
182 function apply_interior(op::SecondDerivativeVariable{2}, v, i, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
183 ṽ = @view v[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
184 c̃ = @view op.coefficient[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
185
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
186 return @inbounds apply_stencil(op.inner_stencil, c̃, ṽ, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
187 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
188
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
189 function apply_upper(op::SecondDerivativeVariable{2}, v, i, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
190 ṽ = @view v[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
191 c̃ = @view op.coefficient[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
192
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
193 stencil = op.closure_stencils[op.size[derivative_direction(op)]-j+1]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
194 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, j)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
195 end