Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 1367:71e89507dd9a feature/variable_derivatives

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