Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 916:fc1f03fc8d65 feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove comments
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 18 Feb 2022 07:19:08 +0100
parents cad3a9f82009
children eb054537fc63
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 """
fc1f03fc8d65 Remove comments
Jonatan Werpers <jonatan@werpers.com>
parents: 913
diff changeset
5 SecondDerivativeVariable{Dir,T,D,...} <: TensorMapping{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 """
906
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
9 struct SecondDerivativeVariable{Dir,T,D,M,IStencil<:NestedStencil{T},CStencil<:NestedStencil{T},TArray<:AbstractArray} <: TensorMapping{T,D,D}
33c7e266e1a9 Add fix for type instability in SecondDerivativeVariable along with tests
Jonatan Werpers <jonatan@werpers.com>
parents: 901
diff changeset
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)
901
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
24 Δxᵢ = spacing(grid)[dir]
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
25 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
26 scaled_closure_stencils = scale.(Tuple(closure_stencils), 1/Δxᵢ^2)
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
27 return SecondDerivativeVariable{dir, dimension(grid)}(scaled_inner_stencil, scaled_closure_stencils, size(grid), coeff)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
882
9098fc936776 Add the coefficient as a part of the struct. Wrap tests in testsets
Jonatan Werpers <jonatan@werpers.com>
parents: 881
diff changeset
30 function SecondDerivativeVariable(grid::EquidistantGrid{1}, coeff::AbstractVector, inner_stencil, closure_stencils)
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
31 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
32 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33
901
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
34 function SecondDerivativeVariable(grid::EquidistantGrid, coeff::AbstractArray, stencil_set, dir)
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
35 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
36 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
37
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
38 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
39 end
5bbc3ea3821b Implement building from stencil set. Add tests for real operators
Jonatan Werpers <jonatan@werpers.com>
parents: 893
diff changeset
40
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
41 derivative_direction(::SecondDerivativeVariable{Dir}) where {Dir} = Dir
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
42
891
f72cc96a58c6 Add some size tests
Jonatan Werpers <jonatan@werpers.com>
parents: 889
diff changeset
43 closure_size(op::SecondDerivativeVariable) = length(op.closure_stencils)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 LazyTensors.range_size(op::SecondDerivativeVariable) = op.size
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 LazyTensors.domain_size(op::SecondDerivativeVariable) = op.size
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
48
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
49 function derivative_view(op, a, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
50 d = derivative_direction(op)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
51
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
52 Iview = Base.setindex(I,:,d)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
53 return @view a[Iview...]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
54 end
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
55
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
56 function apply_lower(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
57 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
58 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
59
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
60 i = I[derivative_direction(op)]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
61 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
62 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
64 function apply_interior(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
65 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
66 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
67
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
68 i = I[derivative_direction(op)]
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
69 return apply_stencil(op.inner_stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71
889
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
72 function apply_upper(op::SecondDerivativeVariable, v, I...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
73 ṽ = derivative_view(op, v, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
74 c̃ = derivative_view(op, op.coefficient, I)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
75
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
76 i = I[derivative_direction(op)]
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
77 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
78 return @inbounds apply_stencil_backwards(stencil, c̃, ṽ, i)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79 end
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
81 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
82 if I[derivative_direction(op)] isa Index{Lower}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
83 return apply_lower(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
84 elseif I[derivative_direction(op)] isa Index{Upper}
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
85 return apply_upper(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
86 else
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
87 return apply_interior(op, v, Int.(I)...)
069e58fb3829 Refactor as multidimensional operator
Jonatan Werpers <jonatan@werpers.com>
parents: 884
diff changeset
88 end
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
893
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
91 function LazyTensors.apply(op::SecondDerivativeVariable, v::AbstractArray, I...)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
92 dir = derivative_direction(op)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
93
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
94 i = I[dir]
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
95 r = getregion(i, closure_size(op), op.size[dir])
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
96
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
97 I = map(i->Index(i, Interior), I)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
98 I = Base.setindex(I, Index(i, r), dir)
422c9f22cf92 Make higher dimensions work
Jonatan Werpers <jonatan@werpers.com>
parents: 891
diff changeset
99 return LazyTensors.apply(op, v, I...)
881
aa4875f9a530 Start implementing the variable second derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 end