Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 2080:0f949681d3d3 refactor/sbp_operators/direction_check tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
author Vidar Stiernström <vidar.stiernstrom@gmail.com>
date Fri, 20 Feb 2026 12:01:05 +0100
parents e9dfc1998d31
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1370
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
1 """
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
2 second_derivative_variable(g, coeff ..., [direction])
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
3
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
4 The variable second derivative operator as a `LazyTensor` on the given grid.
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
5 `coeff` is a grid function of the variable coefficient.
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
6
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
7 Approximates the d/dξ c d/dξ on `g` along the coordinate dimension specified
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
8 by `direction`.
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
9 """
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
10 function second_derivative_variable end
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
11
2080
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
12 function second_derivative_variable(g::TensorGrid, coeff, stencil_set, direction::Int)
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
13 if direction ∉ Interval(0, ndims(g))
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
14 throw(DomainError(direction, "Direction must be inside [0, $(ndims(g))]."))
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
15 end
1370
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
16 inner_stencil = parse_nested_stencil(eltype(coeff), stencil_set["D2variable"]["inner_stencil"])
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
17 closure_stencils = parse_nested_stencil.(eltype(coeff), stencil_set["D2variable"]["closure_stencils"])
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
18
2080
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
19 return second_derivative_variable(g, coeff, inner_stencil, closure_stencils, direction)
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
20 end
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
21
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
22 function second_derivative_variable(g::EquidistantGrid, coeff, stencil_set, direction)
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
23 return second_derivative_variable(TensorGrid(g), coeff, stencil_set, direction)
1370
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
24 end
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
25
1372
3f7f826489a1 Remove type annotation for coeff
Jonatan Werpers <jonatan@werpers.com>
parents: 1371
diff changeset
26 function second_derivative_variable(g::EquidistantGrid, coeff, stencil_set)
2080
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
27 return second_derivative_variable(g::EquidistantGrid, coeff, stencil_set, 1)
1370
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
28 end
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
29
2080
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
30 function second_derivative_variable(g::TensorGrid, coeff, inner_stencil::NestedStencil, closure_stencils, direction)
1371
d0e48c2e6aad Remove stencil input for 1d grid and reorder methods of second_derivative_variable
Jonatan Werpers <jonatan@werpers.com>
parents: 1370
diff changeset
31 check_coefficient(g, coeff)
d0e48c2e6aad Remove stencil input for 1d grid and reorder methods of second_derivative_variable
Jonatan Werpers <jonatan@werpers.com>
parents: 1370
diff changeset
32
2080
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
33 Δxᵢ = spacing(g.grids[direction])
1371
d0e48c2e6aad Remove stencil input for 1d grid and reorder methods of second_derivative_variable
Jonatan Werpers <jonatan@werpers.com>
parents: 1370
diff changeset
34 scaled_inner_stencil = scale(inner_stencil, 1/Δxᵢ^2)
d0e48c2e6aad Remove stencil input for 1d grid and reorder methods of second_derivative_variable
Jonatan Werpers <jonatan@werpers.com>
parents: 1370
diff changeset
35 scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2)
2080
0f949681d3d3 Check that direction of first/second derivative operators is within the dimension of the grid. Add 1D functions for first/second derivative operators that take a direction.
Vidar Stiernström <vidar.stiernstrom@gmail.com>
parents: 1381
diff changeset
36 return SecondDerivativeVariable(coeff, scaled_inner_stencil, scaled_closure_stencils, direction)
1371
d0e48c2e6aad Remove stencil input for 1d grid and reorder methods of second_derivative_variable
Jonatan Werpers <jonatan@werpers.com>
parents: 1370
diff changeset
37 end
1370
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
38
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
39 function check_coefficient(g, coeff)
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
40 if ndims(g) != ndims(coeff)
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
41 throw(ArgumentError("The coefficient has dimension $(ndims(coeff)) while the grid is dimension $(ndims(g))"))
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
42 end
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
43
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
44 if size(g) != size(coeff)
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
45 throw(DimensionMismatch("the size $(size(coeff)) of the coefficient does not match the size $(size(g)) of the grid"))
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
46 end
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
47 end
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
48
4ef8fb75d144 Start splitting out a second_derivative_variable function
Jonatan Werpers <jonatan@werpers.com>
parents: 1368
diff changeset
49
916
fc1f03fc8d65 Remove comments
Jonatan Werpers <jonatan@werpers.com>
parents: 913
diff changeset
50 """
1049
3bb94ce74697 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 940
diff changeset
51 SecondDerivativeVariable{Dir,T,D,...} <: LazyTensor{T,D,D}
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52
916
fc1f03fc8d65 Remove comments
Jonatan Werpers <jonatan@werpers.com>
parents: 913
diff changeset
53 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
54 """
1049
3bb94ce74697 Merge default
Jonatan Werpers <jonatan@werpers.com>
parents: 940
diff changeset
55 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
56 inner_stencil::IStencil
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
57 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
58 coefficient::TArray
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
59
1374
4f3dd84891f6 Simplify constructor for SecondDerivativeVariable
Jonatan Werpers <jonatan@werpers.com>
parents: 1372
diff changeset
60 function SecondDerivativeVariable(coefficient::AbstractArray, inner_stencil::NestedStencil{T}, closure_stencils::NTuple{M,NestedStencil{T}}, dir) where {T,M}
4f3dd84891f6 Simplify constructor for SecondDerivativeVariable
Jonatan Werpers <jonatan@werpers.com>
parents: 1372
diff changeset
61 D = ndims(coefficient)
906
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
62 IStencil = typeof(inner_stencil)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
63 CStencil = eltype(closure_stencils)
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
64 TArray = typeof(coefficient)
1374
4f3dd84891f6 Simplify constructor for SecondDerivativeVariable
Jonatan Werpers <jonatan@werpers.com>
parents: 1372
diff changeset
65 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
66 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
67 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
68
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
69 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
70
891
f72cc96a58c6 Add some size tests
Jonatan Werpers <jonatan@werpers.com>
parents: 889
diff changeset
71 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72
1367
71e89507dd9a Remove size field from SecondDerivativeVariable since it duplicates info from the coefficient field
Jonatan Werpers <jonatan@werpers.com>
parents: 1365
diff changeset
73 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
74 LazyTensors.domain_size(op::SecondDerivativeVariable) = size(op.coefficient)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
76
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
77 function derivative_view(op, a, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
78 d = derivative_direction(op)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
79
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
80 Iview = Base.setindex(I,:,d)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
81 return @view a[Iview...]
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)]
1367
71e89507dd9a Remove size field from SecondDerivativeVariable since it duplicates info from the coefficient field
Jonatan Werpers <jonatan@werpers.com>
parents: 1365
diff changeset
105 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
106 stencil = op.closure_stencils[sz-i+1]
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
107 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
110 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
111 if I[derivative_direction(op)] isa Index{Lower}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
112 return apply_lower(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
113 elseif I[derivative_direction(op)] isa Index{Upper}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
114 return apply_upper(op, v, Int.(I)...)
931
720b1358e06d More strict apply
Jonatan Werpers <jonatan@werpers.com>
parents: 930
diff changeset
115 elseif I[derivative_direction(op)] isa Index{Interior}
720b1358e06d More strict apply
Jonatan Werpers <jonatan@werpers.com>
parents: 930
diff changeset
116 return apply_interior(op, v, Int.(I)...)
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
117 else
931
720b1358e06d More strict apply
Jonatan Werpers <jonatan@werpers.com>
parents: 930
diff changeset
118 error("Invalid region")
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
119 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
120 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
121
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
122 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
123 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
124 sz = domain_size(op)[dir]
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
125
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
126 i = I[dir]
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 = map(i->Index(i, Interior), I)
930
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
129 if 0 < i <= closure_size(op)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
130 I = Base.setindex(I, Index(i, Lower), dir)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
131 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
132 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
133 I = Base.setindex(I, Index(i, Interior), dir)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
134 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
135 elseif sz-closure_size(op) < i <= sz
930
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
136 I = Base.setindex(I, Index(i, Upper), dir)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
137 return LazyTensors.apply(op, v, I...)
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
138 else
1381
e9dfc1998d31 Fix some comments
Jonatan Werpers <jonatan@werpers.com>
parents: 1374
diff changeset
139 error("Bounds error") # This should be `throw(BoundsError())` but the type inference is so fragile that it doesn't work. Needs investigation. / Jonatan 2023-06-08
930
ba5f4a0ec879 Fix type stability for 2d operator
Jonatan Werpers <jonatan@werpers.com>
parents: 917
diff changeset
140 end
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
141 end
932
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
142
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
143
1381
e9dfc1998d31 Fix some comments
Jonatan Werpers <jonatan@werpers.com>
parents: 1374
diff changeset
144 # 2D Specific implementations to avoid type instability
e9dfc1998d31 Fix some comments
Jonatan Werpers <jonatan@werpers.com>
parents: 1374
diff changeset
145 # TBD: Can this be solved by fixing the general methods instead?
932
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
146
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
147
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
148 ## x-direction
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
149 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
150 ṽ = @view v[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
151 c̃ = @view op.coefficient[:,j]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
152
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
153 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
154 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
155
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
156 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
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.inner_stencil, 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_upper(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
1367
71e89507dd9a Remove size field from SecondDerivativeVariable since it duplicates info from the coefficient field
Jonatan Werpers <jonatan@werpers.com>
parents: 1365
diff changeset
167 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
168 stencil = op.closure_stencils[sz-i+1]
932
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
169 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
170 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
171
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 ## y-direction
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
174 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
175 ṽ = @view v[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
176 c̃ = @view op.coefficient[i,:]
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
177
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
178 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
179 end
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
180
863287577ad4 Temporarily add specialized methods for 2D
Jonatan Werpers <jonatan@werpers.com>
parents: 931
diff changeset
181 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
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.inner_stencil, 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_upper(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
1367
71e89507dd9a Remove size field from SecondDerivativeVariable since it duplicates info from the coefficient field
Jonatan Werpers <jonatan@werpers.com>
parents: 1365
diff changeset
192 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
193 stencil = op.closure_stencils[sz-j+1]
932
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