Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/derivatives/dissipation.jl @ 1360:f59228534d3a tooling/benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 20 May 2023 15:15:22 +0200
parents 08f06bfacd5c
children
comparison
equal deleted inserted replaced
1321:42738616422e 1360:f59228534d3a
1 """
2 undivided_skewed04(g::TensorGrid, p, direction)
3 undivided_skewed04(g::EquidistantGrid, p)
4
5 Undivided difference operators approximating the `p`th derivative. The
6 operators do not satisfy any SBP property and are meant to be used for
7 building artificial dissipation terms.
8
9 The operators and how they are used to create accurate artificial dissipation
10 is described in "K. Mattsson, M. Svärd, and J. Nordström, “Stable and Accurate
11 Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp.
12 57–79, Aug. 2004"
13 """
14 function undivided_skewed04 end
15
16 function undivided_skewed04(g::TensorGrid, p, direction)
17 D,Dᵀ = undivided_skewed04(g.grids[direction], p)
18 return (
19 LazyTensors.inflate(D, size(g), direction),
20 LazyTensors.inflate(Dᵀ, size(g), direction),
21 )
22 end
23
24 function undivided_skewed04(g::EquidistantGrid, p)
25 T = eltype(g)
26 interior_weights = T.(dissipation_interior_weights(p))
27
28 D = StencilOperatorDistinctClosures(
29 g,
30 dissipation_interior_stencil(interior_weights),
31 dissipation_lower_closure_stencils(interior_weights),
32 dissipation_upper_closure_stencils(interior_weights),
33 )
34 Dᵀ = StencilOperatorDistinctClosures(
35 g,
36 dissipation_transpose_interior_stencil(interior_weights),
37 dissipation_transpose_lower_closure_stencils(interior_weights),
38 dissipation_transpose_upper_closure_stencils(interior_weights),
39 )
40
41 return D, Dᵀ
42 end
43
44 function dissipation_interior_weights(p)
45 if p == 0
46 return (1,)
47 end
48
49 return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0)
50 end
51
52 midpoint(weights) = length(weights)÷2 + 1
53 midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
54
55 function dissipation_interior_stencil(weights)
56 return Stencil(weights..., center=midpoint(weights))
57 end
58 function dissipation_transpose_interior_stencil(weights)
59 if iseven(length(weights))
60 weights = map(-, weights)
61 end
62
63 return Stencil(weights..., center=midpoint_transpose(weights))
64 end
65
66 dissipation_lower_closure_size(weights) = midpoint(weights) - 1
67 dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
68
69 function dissipation_lower_closure_stencils(interior_weights)
70 stencil(i) = Stencil(interior_weights..., center=i)
71 return ntuple(i->stencil(i), dissipation_lower_closure_size(interior_weights))
72 end
73
74 function dissipation_upper_closure_stencils(interior_weights)
75 center(i) = length(interior_weights) - dissipation_upper_closure_size(interior_weights) + i
76 stencil(i) = Stencil(interior_weights..., center=center(i))
77 return ntuple(i->stencil(i), dissipation_upper_closure_size(interior_weights))
78 end
79
80 function dissipation_transpose_lower_closure_stencils(interior_weights)
81 closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
82
83 N = maximum(s->length(s.weights), closure)
84 return right_pad.(closure, N)
85 end
86
87 function dissipation_transpose_upper_closure_stencils(interior_weights)
88 closure = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights)))
89
90 N = maximum(s->length(s.weights), closure)
91 return left_pad.(closure, N)
92 end
93
94
95 function dissipation_transpose_lower_closure_stencil(interior_weights, i)
96 w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights))
97
98 for k ∈ i:-1:1
99 w = (w..., interior_weights[k])
100 end
101
102 return Stencil(w..., center = i)
103 end
104
105 function dissipation_transpose_upper_closure_stencil(interior_weights, i)
106 j = length(interior_weights)+1-i
107 w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights))
108
109 for k ∈ j:1:length(interior_weights)
110 w = (interior_weights[k], w...)
111 end
112
113 return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1)
114 end