CoCalc -- subcell_limiters

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/solvers/dgsem_p4est/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::P4estMesh{2}, equations)
10
    _, _, dg, cache = mesh_equations_solver_cache(semi)
11

12
    (; neighbor_ids, node_indices) = cache.interfaces
13
    index_range = eachnode(dg)
14

15
    for interface in eachinterface(dg, cache)
16
        # Get element and side index information on the primary element
17
        primary_element = neighbor_ids[1, interface]
18
        primary_indices = node_indices[1, interface]
19

20
        # Get element and side index information on the secondary element
21
        secondary_element = neighbor_ids[2, interface]
22
        secondary_indices = node_indices[2, interface]
23

24
        # Create the local i,j indexing
25
        i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1],
26
                                                                 index_range)
27
        j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2],
28
                                                                 index_range)
29
        i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1],
30
                                                                     index_range)
31
        j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2],
32
                                                                     index_range)
33

34
        i_primary = i_primary_start
35
        j_primary = j_primary_start
36
        i_secondary = i_secondary_start
37
        j_secondary = j_secondary_start
38

39
        for node in eachnode(dg)
40
            var_primary = u[variable, i_primary, j_primary, primary_element]
41
            var_secondary = u[variable, i_secondary, j_secondary, secondary_element]
42

43
            var_min[i_primary, j_primary, primary_element] = min(var_min[i_primary,
44
                                                                         j_primary,
45
                                                                         primary_element],
46
                                                                 var_secondary)
47
            var_max[i_primary, j_primary, primary_element] = max(var_max[i_primary,
48
                                                                         j_primary,
49
                                                                         primary_element],
50
                                                                 var_secondary)
51

52
            var_min[i_secondary, j_secondary, secondary_element] = min(var_min[i_secondary,
53
                                                                               j_secondary,
54
                                                                               secondary_element],
55
                                                                       var_primary)
56
            var_max[i_secondary, j_secondary, secondary_element] = max(var_max[i_secondary,
57
                                                                               j_secondary,
58
                                                                               secondary_element],
59
                                                                       var_primary)
60

61
            # Increment primary element indices
62
            i_primary += i_primary_step
63
            j_primary += j_primary_step
64
            i_secondary += i_secondary_step
65
            j_secondary += j_secondary_step
66
        end
67
    end
68

69
    return nothing
70
end
71

72
@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
73
                                                boundary_conditions::BoundaryConditionPeriodic,
74
                                                mesh::P4estMesh{2},
75
                                                equations, dg, cache)
76
    return nothing
77
end
78

79
@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
80
                                                boundary_conditions,
81
                                                mesh::P4estMesh{2},
82
                                                equations, dg, cache)
83
    (; boundary_condition_types, boundary_indices) = boundary_conditions
84
    (; contravariant_vectors) = cache.elements
85

86
    (; boundaries) = cache
87
    index_range = eachnode(dg)
88

89
    foreach_enumerate(boundary_condition_types) do (i, boundary_condition)
90
        for boundary in boundary_indices[i]
91
            element = boundaries.neighbor_ids[boundary]
92
            node_indices = boundaries.node_indices[boundary]
93
            direction = indices2direction(node_indices)
94

95
            i_node_start, i_node_step = index_to_start_step_2d(node_indices[1],
96
                                                               index_range)
97
            j_node_start, j_node_step = index_to_start_step_2d(node_indices[2],
98
                                                               index_range)
99

100
            i_node = i_node_start
101
            j_node = j_node_start
102
            for i in eachnode(dg)
103
                normal_direction = get_normal_direction(direction,
104
                                                        contravariant_vectors,
105
                                                        i_node, j_node, element)
106

107
                u_inner = get_node_vars(u, equations, dg, i_node, j_node, element)
108

109
                u_outer = get_boundary_outer_state(u_inner, t, boundary_condition,
110
                                                   normal_direction,
111
                                                   mesh, equations, dg, cache,
112
                                                   i_node, j_node, element)
113
                var_outer = u_outer[variable]
114

115
                var_min[i_node, j_node, element] = min(var_min[i_node, j_node, element],
116
                                                       var_outer)
117
                var_max[i_node, j_node, element] = max(var_max[i_node, j_node, element],
118
                                                       var_outer)
119

120
                i_node += i_node_step
121
                j_node += j_node_step
122
            end
123
        end
124
    end
125

126
    return nothing
127
end
128

129
function calc_bounds_onesided_interface!(var_minmax, minmax, variable, u,
130
                                         semi, mesh::P4estMesh{2})
131
    _, equations, dg, cache = mesh_equations_solver_cache(semi)
132

133
    (; neighbor_ids, node_indices) = cache.interfaces
134
    index_range = eachnode(dg)
135

136
    for interface in eachinterface(dg, cache)
137
        # Get element and side index information on the primary element
138
        primary_element = neighbor_ids[1, interface]
139
        primary_indices = node_indices[1, interface]
140

141
        # Get element and side index information on the secondary element
142
        secondary_element = neighbor_ids[2, interface]
143
        secondary_indices = node_indices[2, interface]
144

145
        # Create the local i,j indexing
146
        i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1],
147
                                                                 index_range)
148
        j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2],
149
                                                                 index_range)
150
        i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1],
151
                                                                     index_range)
152
        j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2],
153
                                                                     index_range)
154

155
        i_primary = i_primary_start
156
        j_primary = j_primary_start
157
        i_secondary = i_secondary_start
158
        j_secondary = j_secondary_start
159

160
        for node in eachnode(dg)
161
            var_primary = variable(get_node_vars(u, equations, dg, i_primary, j_primary,
162
                                                 primary_element), equations)
163
            var_secondary = variable(get_node_vars(u, equations, dg, i_secondary,
164
                                                   j_secondary, secondary_element),
165
                                     equations)
166

167
            var_minmax[i_primary, j_primary, primary_element] = minmax(var_minmax[i_primary,
168
                                                                                  j_primary,
169
                                                                                  primary_element],
170
                                                                       var_secondary)
171
            var_minmax[i_secondary, j_secondary, secondary_element] = minmax(var_minmax[i_secondary,
172
                                                                                        j_secondary,
173
                                                                                        secondary_element],
174
                                                                             var_primary)
175

176
            # Increment primary element indices
177
            i_primary += i_primary_step
178
            j_primary += j_primary_step
179
            i_secondary += i_secondary_step
180
            j_secondary += j_secondary_step
181
        end
182
    end
183

184
    return nothing
185
end
186

187
@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
188
                                                boundary_conditions::BoundaryConditionPeriodic,
189
                                                mesh::P4estMesh{2},
190
                                                equations, dg, cache)
191
    return nothing
192
end
193

194
@inline function calc_bounds_onesided_boundary!(var_minmax, minmax, variable, u, t,
195
                                                boundary_conditions,
196
                                                mesh::P4estMesh{2},
197
                                                equations, dg, cache)
198
    (; boundary_condition_types, boundary_indices) = boundary_conditions
199
    (; contravariant_vectors) = cache.elements
200

201
    (; boundaries) = cache
202
    index_range = eachnode(dg)
203

204
    foreach_enumerate(boundary_condition_types) do (i, boundary_condition)
205
        for boundary in boundary_indices[i]
206
            element = boundaries.neighbor_ids[boundary]
207
            node_indices = boundaries.node_indices[boundary]
208
            direction = indices2direction(node_indices)
209

210
            i_node_start, i_node_step = index_to_start_step_2d(node_indices[1],
211
                                                               index_range)
212
            j_node_start, j_node_step = index_to_start_step_2d(node_indices[2],
213
                                                               index_range)
214

215
            i_node = i_node_start
216
            j_node = j_node_start
217
            for node in eachnode(dg)
218
                normal_direction = get_normal_direction(direction,
219
                                                        contravariant_vectors,
220
                                                        i_node, j_node, element)
221

222
                u_inner = get_node_vars(u, equations, dg, i_node, j_node, element)
223

224
                u_outer = get_boundary_outer_state(u_inner, t, boundary_condition,
225
                                                   normal_direction,
226
                                                   mesh, equations, dg, cache,
227
                                                   i_node, j_node, element)
228
                var_outer = variable(u_outer, equations)
229

230
                var_minmax[i_node, j_node, element] = minmax(var_minmax[i_node, j_node,
231
                                                                        element],
232
                                                             var_outer)
233

234
                i_node += i_node_step
235
                j_node += j_node_step
236
            end
237
        end
238
    end
239

240
    return nothing
241
end
242
end # @muladd
243

244
Product

Resources

Company
