Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/derivatives/dissipation.jl @ 1073:5a3281429a48 feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/variable_derivatives
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 24 Mar 2022 12:35:14 +0100
parents c89c6b63c7f4
children 423a6442efc3
comparison
equal deleted inserted replaced
1068:0b0444adacd3 1073:5a3281429a48
1 function undivided_dissipation(g::EquidistantGrid, p, direction)
2 T = eltype(g)
3 interior_weights = T.(dissipation_interior_weights(p))
4
5 D = stencil_operator_distinct_closures(
6 g,
7 dissipation_interior_stencil(interior_weights),
8 dissipation_lower_closure_stencils(interior_weights),
9 dissipation_upper_closure_stencils(interior_weights),
10 direction,
11 )
12 Dᵀ = stencil_operator_distinct_closures(
13 g,
14 dissipation_transpose_interior_stencil(interior_weights),
15 dissipation_transpose_lower_closure_stencils(interior_weights),
16 dissipation_transpose_upper_closure_stencils(interior_weights),
17 direction,
18 )
19
20 return D, Dᵀ
21 end
22
23 undivided_dissipation(g::EquidistantGrid{1}, p) = undivided_dissipation(g, p, 1)
24
25 function dissipation_interior_weights(p)
26 if p == 0
27 return (1,)
28 end
29
30 return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0)
31 end
32
33 midpoint(weights) = length(weights)÷2 + 1
34 midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
35
36 function dissipation_interior_stencil(weights)
37 return Stencil(weights..., center=midpoint(weights))
38 end
39 function dissipation_transpose_interior_stencil(weights)
40 if iseven(length(weights))
41 weights = map(-, weights)
42 end
43
44 return Stencil(weights..., center=midpoint_transpose(weights))
45 end
46
47 dissipation_lower_closure_size(weights) = midpoint(weights) - 1
48 dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
49
50 dissipation_lower_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=i ), dissipation_lower_closure_size(interior_weights))
51 dissipation_upper_closure_stencils(interior_weights) = ntuple(i->Stencil(interior_weights..., center=length(interior_weights)-dissipation_upper_closure_size(interior_weights)+i), dissipation_upper_closure_size(interior_weights))
52
53 function dissipation_transpose_lower_closure_stencils(interior_weights)
54 closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
55
56 N = maximum(s->length(s.weights), closure)
57 return right_pad.(closure, N)
58 end
59
60 function dissipation_transpose_upper_closure_stencils(interior_weights)
61 closure = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights)))
62
63 N = maximum(s->length(s.weights), closure)
64 return left_pad.(closure, N)
65 end
66
67
68 function dissipation_transpose_lower_closure_stencil(interior_weights, i)
69 w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights))
70
71 for k ∈ i:-1:1
72 w = (w..., interior_weights[k])
73 end
74
75 return Stencil(w..., center = i)
76 end
77
78 function dissipation_transpose_upper_closure_stencil(interior_weights, i)
79 j = length(interior_weights)+1-i
80 w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights))
81
82 for k ∈ j:1:length(interior_weights)
83 w = (interior_weights[k], w...)
84 end
85
86 return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1)
87 end