Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/derivatives/dissipation.jl @ 1221:b3b4d29b46c3 refactor/grids

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 Fri, 10 Feb 2023 08:36:56 +0100
parents d60a10ad6579
children 356ec6a72974
comparison
equal deleted inserted replaced
1220:93bba649aea2 1221:b3b4d29b46c3
1 """
2 undivided_skewed04(g::EquidistantGrid, p, direction)
3
4 Create undivided difference operators approximating the `p`th derivative. The
5 operators do not satisfy any SBP-property and are meant to be used for
6 building artificial dissipation terms.
7
8 The operators and how they are used to create accurate artifical dissipation
9 is described in "K. Mattsson, M. Svärd, and J. Nordström, “Stable and Accurate
10 Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp.
11 57–79, Aug. 2004"
12 """
13 function undivided_skewed04(g::EquidistantGrid, p, direction)
14 T = eltype(g)
15 interior_weights = T.(dissipation_interior_weights(p))
16
17 D = stencil_operator_distinct_closures(
18 g,
19 dissipation_interior_stencil(interior_weights),
20 dissipation_lower_closure_stencils(interior_weights),
21 dissipation_upper_closure_stencils(interior_weights),
22 direction,
23 )
24 Dᵀ = stencil_operator_distinct_closures(
25 g,
26 dissipation_transpose_interior_stencil(interior_weights),
27 dissipation_transpose_lower_closure_stencils(interior_weights),
28 dissipation_transpose_upper_closure_stencils(interior_weights),
29 direction,
30 )
31
32 return D, Dᵀ
33 end
34
35 undivided_skewed04(g::EquidistantGrid{1}, p) = undivided_skewed04(g, p, 1)
36
37 function dissipation_interior_weights(p)
38 if p == 0
39 return (1,)
40 end
41
42 return (0, dissipation_interior_weights(p-1)...) .- (dissipation_interior_weights(p-1)..., 0)
43 end
44
45 midpoint(weights) = length(weights)÷2 + 1
46 midpoint_transpose(weights) = length(weights)+1 - midpoint(weights)
47
48 function dissipation_interior_stencil(weights)
49 return Stencil(weights..., center=midpoint(weights))
50 end
51 function dissipation_transpose_interior_stencil(weights)
52 if iseven(length(weights))
53 weights = map(-, weights)
54 end
55
56 return Stencil(weights..., center=midpoint_transpose(weights))
57 end
58
59 dissipation_lower_closure_size(weights) = midpoint(weights) - 1
60 dissipation_upper_closure_size(weights) = length(weights) - midpoint(weights)
61
62 function dissipation_lower_closure_stencils(interior_weights)
63 stencil(i) = Stencil(interior_weights..., center=i)
64 return ntuple(i->stencil(i), dissipation_lower_closure_size(interior_weights))
65 end
66
67 function dissipation_upper_closure_stencils(interior_weights)
68 center(i) = length(interior_weights) - dissipation_upper_closure_size(interior_weights) + i
69 stencil(i) = Stencil(interior_weights..., center=center(i))
70 return ntuple(i->stencil(i), dissipation_upper_closure_size(interior_weights))
71 end
72
73 function dissipation_transpose_lower_closure_stencils(interior_weights)
74 closure = ntuple(i->dissipation_transpose_lower_closure_stencil(interior_weights, i), length(interior_weights))
75
76 N = maximum(s->length(s.weights), closure)
77 return right_pad.(closure, N)
78 end
79
80 function dissipation_transpose_upper_closure_stencils(interior_weights)
81 closure = reverse(ntuple(i->dissipation_transpose_upper_closure_stencil(interior_weights, i), length(interior_weights)))
82
83 N = maximum(s->length(s.weights), closure)
84 return left_pad.(closure, N)
85 end
86
87
88 function dissipation_transpose_lower_closure_stencil(interior_weights, i)
89 w = ntuple(k->interior_weights[i], dissipation_lower_closure_size(interior_weights))
90
91 for k ∈ i:-1:1
92 w = (w..., interior_weights[k])
93 end
94
95 return Stencil(w..., center = i)
96 end
97
98 function dissipation_transpose_upper_closure_stencil(interior_weights, i)
99 j = length(interior_weights)+1-i
100 w = ntuple(k->interior_weights[j], dissipation_upper_closure_size(interior_weights))
101
102 for k ∈ j:1:length(interior_weights)
103 w = (interior_weights[k], w...)
104 end
105
106 return Stencil(w..., center = length(interior_weights)-midpoint(interior_weights)+1)
107 end