CoCalc -- subcell_limiters

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/solvers/dgsem_structured/subcell_limiters_2d.jl
⁵⁵⁹¹ views
1
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
2
# Since these FMAs can increase the performance of many numerical algorithms,
3
# we need to opt-in explicitly.
4
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
5
@muladd begin
6
#! format: noindent
7

8
function calc_bounds_twosided_interface!(var_min, var_max, variable, u,
9
                                         semi, mesh::StructuredMesh{2}, equations)
10
    _, _, dg, cache = mesh_equations_solver_cache(semi)
11

12
    for element in eachelement(dg, cache)
13
        # Get neighboring element ids
14
        left = cache.elements.left_neighbors[1, element]
15
        lower = cache.elements.left_neighbors[2, element]
16

17
        if left != 0
18
            for j in eachnode(dg)
19
                var_left = u[variable, nnodes(dg), j, left]
20
                var_element = u[variable, 1, j, element]
21

22
                var_min[1, j, element] = min(var_min[1, j, element], var_left)
23
                var_max[1, j, element] = max(var_max[1, j, element], var_left)
24

25
                var_min[nnodes(dg), j, left] = min(var_min[nnodes(dg), j, left],
26
                                                   var_element)
27
                var_max[nnodes(dg), j, left] = max(var_max[nnodes(dg), j, left],
28
                                                   var_element)
29
            end
30
        end
31
        if lower != 0
32
            for i in eachnode(dg)
33
                var_lower = u[variable, i, nnodes(dg), lower]
34
                var_element = u[variable, i, 1, element]
35

36
                var_min[i, 1, element] = min(var_min[i, 1, element], var_lower)
37
                var_max[i, 1, element] = max(var_max[i, 1, element], var_lower)
38

39
                var_min[i, nnodes(dg), lower] = min(var_min[i, nnodes(dg), lower],
40
                                                    var_element)
41
                var_max[i, nnodes(dg), lower] = max(var_max[i, nnodes(dg), lower],
42
                                                    var_element)
43
            end
44
        end
45
    end
46

47
    return nothing
48
end
49

50
@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
51
                                                boundary_conditions,
52
                                                mesh::StructuredMesh{2},
53
                                                equations, dg, cache)
54
    if isperiodic(mesh)
55
        return nothing
56
    end
57
    (; contravariant_vectors) = cache.elements
58

59
    linear_indices = LinearIndices(size(mesh))
60
    if !isperiodic(mesh, 1)
61
        # - xi direction
62
        for cell_y in axes(mesh, 2)
63
            element = linear_indices[begin, cell_y]
64
            for j in eachnode(dg)
65
                Ja1 = get_contravariant_vector(1, contravariant_vectors,
66
                                               1, j, element)
67
                u_inner = get_node_vars(u, equations, dg, 1, j, element)
68
                u_outer = get_boundary_outer_state(u_inner, t,
69
                                                   boundary_conditions[1], Ja1, 1,
70
                                                   mesh, equations, dg, cache,
71
                                                   1, j, element)
72
                var_outer = u_outer[variable]
73

74
                var_min[1, j, element] = min(var_min[1, j, element], var_outer)
75
                var_max[1, j, element] = max(var_max[1, j, element], var_outer)
76
            end
77
        end
78
        # + xi direction
79
        for cell_y in axes(mesh, 2)
80
            element = linear_indices[end, cell_y]
81
            for j in eachnode(dg)
82
                Ja1 = get_contravariant_vector(1, contravariant_vectors,
83
                                               nnodes(dg), j, element)
84
                u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element)
85
                u_outer = get_boundary_outer_state(u_inner, t,
86
                                                   boundary_conditions[2], Ja1, 2,
87
                                                   mesh, equations, dg, cache,
88
                                                   nnodes(dg), j, element)
89
                var_outer = u_outer[variable]
90

91
                var_min[nnodes(dg), j, element] = min(var_min[nnodes(dg), j, element],
92
                                                      var_outer)
93
                var_max[nnodes(dg), j, element] = max(var_max[nnodes(dg), j, element],
94
                                                      var_outer)
95
            end
96
        end
97
    end
98
    if !isperiodic(mesh, 2)
99
        # - eta direction
100
        for cell_x in axes(mesh, 1)
101
            element = linear_indices[cell_x, begin]
102
            for i in eachnode(dg)
103
                Ja2 = get_contravariant_vector(2, contravariant_vectors,
104
                                               i, 1, element)
105
                u_inner = get_node_vars(u, equations, dg, i, 1, element)
106
                u_outer = get_boundary_outer_state(u_inner, t,
107
                                                   boundary_conditions[3], Ja2, 3,
108
                                                   mesh, equations, dg, cache,
109
                                                   i, 1, element)
110
                var_outer = u_outer[variable]
111

112
                var_min[i, 1, element] = min(var_min[i, 1, element], var_outer)
113
                var_max[i, 1, element] = max(var_max[i, 1, element], var_outer)
114
            end
115
        end
116
        # - eta direction
117
        for cell_x in axes(mesh, 1)
118
            element = linear_indices[cell_x, end]
119
            for i in eachnode(dg)
120
                Ja2 = get_contravariant_vector(2, contravariant_vectors,
121
                                               i, nnodes(dg), element)
122
                u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element)
123
                u_outer = get_boundary_outer_state(u_inner, t,
124
                                                   boundary_conditions[4], Ja2, 4,
125
                                                   mesh, equations, dg, cache,
126
                                                   i, nnodes(dg), element)
127
                var_outer = u_outer[variable]
128

129
                var_min[i, nnodes(dg), element] = min(var_min[i, nnodes(dg), element],
130
                                                      var_outer)
131
                var_max[i, nnodes(dg), element] = max(var_max[i, nnodes(dg), element],
132
                                                      var_outer)
133
            end
134
        end
135
    end
136

137
    return nothing
138
end
139

140
function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u,
141
                                         semi, mesh::StructuredMesh{2})
142
    _, equations, dg, cache = mesh_equations_solver_cache(semi)
143

144
    for element in eachelement(dg, cache)
145
        # Get neighboring element ids
146
        left = cache.elements.left_neighbors[1, element]
147
        lower = cache.elements.left_neighbors[2, element]
148

149
        if left != 0
150
            for j in eachnode(dg)
151
                var_left = variable(get_node_vars(u, equations, dg, nnodes(dg), j,
152
                                                  left), equations)
153
                var_element = variable(get_node_vars(u, equations, dg, 1, j, element),
154
                                       equations)
155

156
                var_minmax[1, j, element] = minmax(var_minmax[1, j, element], var_left)
157
                var_minmax[nnodes(dg), j, left] = minmax(var_minmax[nnodes(dg), j,
158
                                                                    left], var_element)
159
            end
160
        end
161
        if lower != 0
162
            for i in eachnode(dg)
163
                var_lower = variable(get_node_vars(u, equations, dg, i, nnodes(dg),
164
                                                   lower), equations)
165
                var_element = variable(get_node_vars(u, equations, dg, i, 1, element),
166
                                       equations)
167

168
                var_minmax[i, 1, element] = minmax(var_minmax[i, 1, element], var_lower)
169
                var_minmax[i, nnodes(dg), lower] = minmax(var_minmax[i, nnodes(dg),
170
                                                                     lower],
171
                                                          var_element)
172
            end
173
        end
174
    end
175

176
    return nothing
177
end
178

179
@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
180
                                                boundary_conditions,
181
                                                mesh::StructuredMesh{2},
182
                                                equations, dg, cache)
183
    if isperiodic(mesh)
184
        return nothing
185
    end
186
    (; contravariant_vectors) = cache.elements
187

188
    linear_indices = LinearIndices(size(mesh))
189
    if !isperiodic(mesh, 1)
190
        # - xi direction
191
        for cell_y in axes(mesh, 2)
192
            element = linear_indices[begin, cell_y]
193
            for j in eachnode(dg)
194
                Ja1 = get_contravariant_vector(1, contravariant_vectors,
195
                                               1, j, element)
196
                u_inner = get_node_vars(u, equations, dg, 1, j, element)
197
                u_outer = get_boundary_outer_state(u_inner, t,
198
                                                   boundary_conditions[1], Ja1, 1,
199
                                                   mesh, equations, dg, cache,
200
                                                   1, j, element)
201
                var_outer = variable(u_outer, equations)
202

203
                var_minmax[1, j, element] = minmax(var_minmax[1, j, element], var_outer)
204
            end
205
        end
206
        # + xi direction
207
        for cell_y in axes(mesh, 2)
208
            element = linear_indices[end, cell_y]
209
            for j in eachnode(dg)
210
                Ja1 = get_contravariant_vector(1, contravariant_vectors,
211
                                               nnodes(dg), j, element)
212
                u_inner = get_node_vars(u, equations, dg, nnodes(dg), j, element)
213
                u_outer = get_boundary_outer_state(u_inner, t,
214
                                                   boundary_conditions[2], Ja1, 2,
215
                                                   mesh, equations, dg, cache,
216
                                                   nnodes(dg), j, element)
217
                var_outer = variable(u_outer, equations)
218

219
                var_minmax[nnodes(dg), j, element] = minmax(var_minmax[nnodes(dg), j,
220
                                                                       element],
221
                                                            var_outer)
222
            end
223
        end
224
    end
225
    if !isperiodic(mesh, 2)
226
        # - eta direction
227
        for cell_x in axes(mesh, 1)
228
            element = linear_indices[cell_x, begin]
229
            for i in eachnode(dg)
230
                Ja2 = get_contravariant_vector(2, contravariant_vectors,
231
                                               i, 1, element)
232
                u_inner = get_node_vars(u, equations, dg, i, 1, element)
233
                u_outer = get_boundary_outer_state(u_inner, t,
234
                                                   boundary_conditions[3], Ja2, 3,
235
                                                   mesh, equations, dg, cache,
236
                                                   i, 1, element)
237
                var_outer = variable(u_outer, equations)
238

239
                var_minmax[i, 1, element] = minmax(var_minmax[i, 1, element], var_outer)
240
            end
241
        end
242
        # + eta direction
243
        for cell_x in axes(mesh, 1)
244
            element = linear_indices[cell_x, end]
245
            for i in eachnode(dg)
246
                Ja2 = get_contravariant_vector(2, contravariant_vectors,
247
                                               i, nnodes(dg), element)
248
                u_inner = get_node_vars(u, equations, dg, i, nnodes(dg), element)
249
                u_outer = get_boundary_outer_state(u_inner, t,
250
                                                   boundary_conditions[4], Ja2, 4,
251
                                                   mesh, equations, dg, cache,
252
                                                   i, nnodes(dg), element)
253
                var_outer = variable(u_outer, equations)
254

255
                var_minmax[i, nnodes(dg), element] = minmax(var_minmax[i, nnodes(dg),
256
                                                                       element],
257
                                                            var_outer)
258
            end
259
        end
260
    end
261

262
    return nothing
263
end
264
end # @muladd
265

266
Product

Resources

Company
