Mercurial > repos > public > sbplib_julia

annotate src/SbpOperators/volumeops/derivatives/second_derivative_variable.jl @ 901:5bbc3ea3821b feature/variable_derivatives

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