Mercurial > repos > public > sbplib_julia

comparison src/SbpOperators/volumeops/derivatives/dissipation.jl @ 1395:bdcdbd4ea9cd feature/boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default. Comment out broken tests for boundary_conditions at sat
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 26 Jul 2023 21:35:50 +0200
parents 08f06bfacd5c
children
comparison
equal deleted inserted replaced
1217:ea2e8254820a 1395:bdcdbd4ea9cd
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