CoCalc -- subcell_limiters

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/solvers/dgsem_tree/subcell_limiters_3d.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
###############################################################################
9
# IDP Limiting
10
###############################################################################
11

12
###############################################################################
13
# Calculation of local bounds using low-order FV solution
14

15
@inline function calc_bounds_twosided!(var_min, var_max, variable,
16
                                       u::AbstractArray{<:Any, 5}, t,
17
                                       semi, equations)
18
    mesh, _, dg, cache = mesh_equations_solver_cache(semi)
19

20
    # Calc bounds inside elements
21
    @threaded for element in eachelement(dg, cache)
22
        # Calculate bounds at Gauss-Lobatto nodes
23
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
24
            var = u[variable, i, j, k, element]
25
            var_min[i, j, k, element] = var
26
            var_max[i, j, k, element] = var
27
        end
28

29
        # Apply values in x direction
30
        for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg)
31
            var = u[variable, i - 1, j, k, element]
32
            var_min[i, j, k, element] = min(var_min[i, j, k, element], var)
33
            var_max[i, j, k, element] = max(var_max[i, j, k, element], var)
34

35
            var = u[variable, i, j, k, element]
36
            var_min[i - 1, j, k, element] = min(var_min[i - 1, j, k, element], var)
37
            var_max[i - 1, j, k, element] = max(var_max[i - 1, j, k, element], var)
38
        end
39

40
        # Apply values in y direction
41
        for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg)
42
            var = u[variable, i, j - 1, k, element]
43
            var_min[i, j, k, element] = min(var_min[i, j, k, element], var)
44
            var_max[i, j, k, element] = max(var_max[i, j, k, element], var)
45

46
            var = u[variable, i, j, k, element]
47
            var_min[i, j - 1, k, element] = min(var_min[i, j - 1, k, element], var)
48
            var_max[i, j - 1, k, element] = max(var_max[i, j - 1, k, element], var)
49
        end
50

51
        # Apply values in z direction
52
        for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg)
53
            var = u[variable, i, j, k - 1, element]
54
            var_min[i, j, k, element] = min(var_min[i, j, k, element], var)
55
            var_max[i, j, k, element] = max(var_max[i, j, k, element], var)
56

57
            var = u[variable, i, j, k, element]
58
            var_min[i, j, k - 1, element] = min(var_min[i, j, k - 1, element], var)
59
            var_max[i, j, k - 1, element] = max(var_max[i, j, k - 1, element], var)
60
        end
61
    end
62

63
    # Calc bounds at element interfaces and periodic boundaries
64
    calc_bounds_twosided_interface!(var_min, var_max, variable, u,
65
                                    semi, mesh, equations)
66

67
    # Calc bounds at physical boundaries
68
    (; boundary_conditions) = semi
69
    calc_bounds_twosided_boundary!(var_min, var_max, variable, u, t,
70
                                   boundary_conditions,
71
                                   mesh, equations, dg, cache)
72
    return nothing
73
end
74

75
@inline function calc_bounds_twosided_interface!(var_min, var_max, variable, u,
76
                                                 semi, mesh::TreeMesh3D, equations)
77
    _, _, dg, cache = mesh_equations_solver_cache(semi)
78

79
    for interface in eachinterface(dg, cache)
80
        # Get neighboring element ids
81
        left_element = cache.interfaces.neighbor_ids[1, interface]
82
        right_element = cache.interfaces.neighbor_ids[2, interface]
83

84
        orientation = cache.interfaces.orientations[interface]
85

86
        for j in eachnode(dg), i in eachnode(dg)
87
            # Define node indices for left and right element based on the interface orientation
88
            if orientation == 1
89
                # interface in x-direction
90
                index_left = (nnodes(dg), i, j)
91
                index_right = (1, i, j)
92
            elseif orientation == 2
93
                # interface in y-direction
94
                index_left = (i, nnodes(dg), j)
95
                index_right = (i, 1, j)
96
            else # if orientation == 3
97
                # interface in z-direction
98
                index_left = (i, j, nnodes(dg))
99
                index_right = (i, j, 1)
100
            end
101
            var_left = u[variable, index_left..., left_element]
102
            var_right = u[variable, index_right..., right_element]
103

104
            var_min[index_right..., right_element] = min(var_min[index_right...,
105
                                                                 right_element],
106
                                                         var_left)
107
            var_max[index_right..., right_element] = max(var_max[index_right...,
108
                                                                 right_element],
109
                                                         var_left)
110

111
            var_min[index_left..., left_element] = min(var_min[index_left...,
112
                                                               left_element], var_right)
113
            var_max[index_left..., left_element] = max(var_max[index_left...,
114
                                                               left_element], var_right)
115
        end
116
    end
117

118
    return nothing
119
end
120

121
@inline function calc_bounds_twosided_boundary!(var_min, var_max, variable,
122
                                                u, t, boundary_conditions,
123
                                                mesh::TreeMesh{3}, equations,
124
                                                dg, cache)
125
    for boundary in eachboundary(dg, cache)
126
        element = cache.boundaries.neighbor_ids[boundary]
127
        orientation = cache.boundaries.orientations[boundary]
128
        neighbor_side = cache.boundaries.neighbor_sides[boundary]
129

130
        for j in eachnode(dg), i in eachnode(dg)
131
            # Define node indices and boundary index based on the orientation and neighbor_side
132
            if neighbor_side == 2 # Element is on the right, boundary on the left
133
                if orientation == 1 # boundary in x-direction
134
                    node_index = (1, i, j)
135
                    boundary_index = 1
136
                elseif orientation == 2 # boundary in y-direction
137
                    node_index = (i, 1, j)
138
                    boundary_index = 3
139
                else # orientation == 3 # boundary in z-direction
140
                    node_index = (i, j, 1)
141
                    boundary_index = 5
142
                end
143
            else # Element is on the left, boundary on the right
144
                if orientation == 1 # boundary in x-direction
145
                    node_index = (nnodes(dg), i, j)
146
                    boundary_index = 2
147
                elseif orientation == 2 # boundary in y-direction
148
                    node_index = (i, nnodes(dg), j)
149
                    boundary_index = 4
150
                else # orientation == 3 # boundary in z-direction
151
                    node_index = (i, j, nnodes(dg))
152
                    boundary_index = 6
153
                end
154
            end
155
            u_inner = get_node_vars(u, equations, dg, node_index..., element)
156
            u_outer = get_boundary_outer_state(u_inner, t,
157
                                               boundary_conditions[boundary_index],
158
                                               orientation, boundary_index,
159
                                               mesh, equations, dg, cache,
160
                                               node_index..., element)
161
            var_outer = u_outer[variable]
162

163
            var_min[node_index..., element] = min(var_min[node_index..., element],
164
                                                  var_outer)
165
            var_max[node_index..., element] = max(var_max[node_index..., element],
166
                                                  var_outer)
167
        end
168
    end
169

170
    return nothing
171
end
172

173
@inline function calc_bounds_onesided!(var_minmax, min_or_max, variable,
174
                                       u::AbstractArray{<:Any, 5}, t,
175
                                       semi)
176
    mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
177

178
    # The approach used in `calc_bounds_twosided!` is not used here because it requires more
179
    # evaluations of the variable and is therefore slower.
180

181
    # Calc bounds inside elements
182
    @threaded for element in eachelement(dg, cache)
183
        # Reset bounds
184
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
185
            if min_or_max === max
186
                var_minmax[i, j, k, element] = typemin(eltype(var_minmax))
187
            else
188
                var_minmax[i, j, k, element] = typemax(eltype(var_minmax))
189
            end
190
        end
191

192
        # Calculate bounds at Gauss-Lobatto nodes
193
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
194
            var = variable(get_node_vars(u, equations, dg, i, j, k, element), equations)
195
            var_minmax[i, j, k, element] = min_or_max(var_minmax[i, j, k, element], var)
196

197
            if i > 1
198
                var_minmax[i - 1, j, k, element] = min_or_max(var_minmax[i - 1, j, k,
199
                                                                         element], var)
200
            end
201
            if i < nnodes(dg)
202
                var_minmax[i + 1, j, k, element] = min_or_max(var_minmax[i + 1, j, k,
203
                                                                         element], var)
204
            end
205
            if j > 1
206
                var_minmax[i, j - 1, k, element] = min_or_max(var_minmax[i, j - 1, k,
207
                                                                         element], var)
208
            end
209
            if j < nnodes(dg)
210
                var_minmax[i, j + 1, k, element] = min_or_max(var_minmax[i, j + 1, k,
211
                                                                         element], var)
212
            end
213
            if k > 1
214
                var_minmax[i, j, k - 1, element] = min_or_max(var_minmax[i, j, k - 1,
215
                                                                         element], var)
216
            end
217
            if k < nnodes(dg)
218
                var_minmax[i, j, k + 1, element] = min_or_max(var_minmax[i, j, k + 1,
219
                                                                         element], var)
220
            end
221
        end
222
    end
223

224
    # Calc bounds at element interfaces and periodic boundaries
225
    calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u,
226
                                    semi, mesh)
227

228
    # Calc bounds at physical boundaries
229
    (; boundary_conditions) = semi
230
    calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable, u, t,
231
                                   boundary_conditions,
232
                                   mesh, equations, dg, cache)
233

234
    return nothing
235
end
236

237
@inline function calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u,
238
                                                 semi, mesh::TreeMesh{3})
239
    _, equations, dg, cache = mesh_equations_solver_cache(semi)
240

241
    for interface in eachinterface(dg, cache)
242
        # Get neighboring element ids
243
        left_element = cache.interfaces.neighbor_ids[1, interface]
244
        right_element = cache.interfaces.neighbor_ids[2, interface]
245

246
        orientation = cache.interfaces.orientations[interface]
247

248
        for j in eachnode(dg), i in eachnode(dg)
249
            # Define node indices for left and right element based on the interface orientation
250
            if orientation == 1
251
                # interface in x-direction
252
                index_left = (nnodes(dg), i, j)
253
                index_right = (1, i, j)
254
            elseif orientation == 2
255
                # interface in y-direction
256
                index_left = (i, nnodes(dg), j)
257
                index_right = (i, 1, j)
258
            else # if orientation == 3
259
                # interface in z-direction
260
                index_left = (i, j, nnodes(dg))
261
                index_right = (i, j, 1)
262
            end
263
            var_left = variable(get_node_vars(u, equations, dg, index_left...,
264
                                              left_element),
265
                                equations)
266
            var_right = variable(get_node_vars(u, equations, dg, index_right...,
267
                                               right_element),
268
                                 equations)
269

270
            var_minmax[index_right..., right_element] = min_or_max(var_minmax[index_right...,
271
                                                                              right_element],
272
                                                                   var_left)
273
            var_minmax[index_left..., left_element] = min_or_max(var_minmax[index_left...,
274
                                                                            left_element],
275
                                                                 var_right)
276
        end
277
    end
278

279
    return nothing
280
end
281

282
@inline function calc_bounds_onesided_boundary!(var_minmax, min_or_max, variable,
283
                                                u, t, boundary_conditions,
284
                                                mesh::TreeMesh{3}, equations,
285
                                                dg, cache)
286
    for boundary in eachboundary(dg, cache)
287
        element = cache.boundaries.neighbor_ids[boundary]
288
        orientation = cache.boundaries.orientations[boundary]
289
        neighbor_side = cache.boundaries.neighbor_sides[boundary]
290

291
        for j in eachnode(dg), i in eachnode(dg)
292
            # Define node indices and boundary index based on the orientation and neighbor_side
293
            if neighbor_side == 2 # Element is on the right, boundary on the left
294
                if orientation == 1 # boundary in x-direction
295
                    node_index = (1, i, j)
296
                    boundary_index = 1
297
                elseif orientation == 2 # boundary in y-direction
298
                    node_index = (i, 1, j)
299
                    boundary_index = 3
300
                else # orientation == 3 # boundary in z-direction
301
                    node_index = (i, j, 1)
302
                    boundary_index = 5
303
                end
304
            else # Element is on the left, boundary on the right
305
                if orientation == 1 # boundary in x-direction
306
                    node_index = (nnodes(dg), i, j)
307
                    boundary_index = 2
308
                elseif orientation == 2 # boundary in y-direction
309
                    node_index = (i, nnodes(dg), j)
310
                    boundary_index = 4
311
                else # orientation == 3 # boundary in z-direction
312
                    node_index = (i, j, nnodes(dg))
313
                    boundary_index = 6
314
                end
315
            end
316
            u_inner = get_node_vars(u, equations, dg, node_index..., element)
317
            u_outer = get_boundary_outer_state(u_inner, t,
318
                                               boundary_conditions[boundary_index],
319
                                               orientation, boundary_index,
320
                                               mesh, equations, dg, cache,
321
                                               node_index..., element)
322
            var_outer = variable(u_outer, equations)
323

324
            var_minmax[node_index..., element] = min_or_max(var_minmax[node_index...,
325
                                                                       element],
326
                                                            var_outer)
327
        end
328
    end
329

330
    return nothing
331
end
332

333
###############################################################################
334
# Local minimum and maximum limiting of conservative variables
335

336
@inline function idp_local_twosided!(alpha, limiter, u::AbstractArray{<:Any, 5},
337
                                     t, dt, semi, variable)
338
    mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
339
    (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes
340
    (; inverse_weights) = dg.basis
341

342
    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
343
    variable_string = string(variable)
344
    var_min = variable_bounds[Symbol(variable_string, "_min")]
345
    var_max = variable_bounds[Symbol(variable_string, "_max")]
346
    calc_bounds_twosided!(var_min, var_max, variable, u, t, semi, equations)
347

348
    @threaded for element in eachelement(dg, semi.cache)
349
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
350
            inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian,
351
                                                    mesh, i, j, k, element)
352
            var = u[variable, i, j, k, element]
353
            # Real Zalesak type limiter
354
            #   * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids"
355
            #   * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics"
356
            #   Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is
357
            #         for each interface, not each node
358

359
            Qp = max(0, (var_max[i, j, k, element] - var) / dt)
360
            Qm = min(0, (var_min[i, j, k, element] - var) / dt)
361

362
            # Calculate Pp and Pm
363
            # Note: Boundaries of antidiffusive_flux1/2 are constant 0, so they make no difference here.
364
            val_flux1_local = inverse_weights[i] *
365
                              antidiffusive_flux1_R[variable, i, j, k, element]
366
            val_flux1_local_ip1 = -inverse_weights[i] *
367
                                  antidiffusive_flux1_L[variable, i + 1, j, k, element]
368
            val_flux2_local = inverse_weights[j] *
369
                              antidiffusive_flux2_R[variable, i, j, k, element]
370
            val_flux2_local_jp1 = -inverse_weights[j] *
371
                                  antidiffusive_flux2_L[variable, i, j + 1, k, element]
372
            val_flux3_local = inverse_weights[k] *
373
                              antidiffusive_flux3_R[variable, i, j, k, element]
374
            val_flux3_local_jp1 = -inverse_weights[k] *
375
                                  antidiffusive_flux3_L[variable, i, j, k + 1, element]
376

377
            Pp = max(0, val_flux1_local) + max(0, val_flux1_local_ip1) +
378
                 max(0, val_flux2_local) + max(0, val_flux2_local_jp1) +
379
                 max(0, val_flux3_local) + max(0, val_flux3_local_jp1)
380
            Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) +
381
                 min(0, val_flux2_local) + min(0, val_flux2_local_jp1) +
382
                 min(0, val_flux3_local) + min(0, val_flux3_local_jp1)
383

384
            Pp = inverse_jacobian * Pp
385
            Pm = inverse_jacobian * Pm
386

387
            # Compute blending coefficient avoiding division by zero
388
            # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8))
389
            Qp = abs(Qp) /
390
                 (abs(Pp) + eps(typeof(Qp)) * 100 * abs(var_max[i, j, k, element]))
391
            Qm = abs(Qm) /
392
                 (abs(Pm) + eps(typeof(Qm)) * 100 * abs(var_max[i, j, k, element]))
393

394
            # Calculate alpha at nodes
395
            alpha[i, j, k, element] = max(alpha[i, j, k, element], 1 - min(1, Qp, Qm))
396
        end
397
    end
398

399
    return nothing
400
end
401

402
##############################################################################
403
# Local one-sided limiting of nonlinear variables
404

405
@inline function idp_local_onesided!(alpha, limiter, u::AbstractArray{<:Real, 5},
406
                                     t, dt, semi,
407
                                     variable, min_or_max)
408
    mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
409
    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
410
    var_minmax = variable_bounds[Symbol(string(variable), "_", string(min_or_max))]
411
    calc_bounds_onesided!(var_minmax, min_or_max, variable, u, t, semi)
412

413
    # Perform Newton's bisection method to find new alpha
414
    @threaded for element in eachelement(dg, cache)
415
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
416
            inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian,
417
                                                    mesh, i, j, k, element)
418
            u_local = get_node_vars(u, equations, dg, i, j, k, element)
419
            newton_loops_alpha!(alpha, var_minmax[i, j, k, element],
420
                                u_local, i, j, k, element,
421
                                variable, min_or_max,
422
                                initial_check_local_onesided_newton_idp,
423
                                final_check_local_onesided_newton_idp,
424
                                inverse_jacobian, dt, equations, dg, cache, limiter)
425
        end
426
    end
427

428
    return nothing
429
end
430

431
###############################################################################
432
# Global positivity limiting of conservative variables
433

434
@inline function idp_positivity_conservative!(alpha, limiter,
435
                                              u::AbstractArray{<:Real, 5}, dt, semi,
436
                                              variable)
437
    mesh, _, dg, cache = mesh_equations_solver_cache(semi)
438
    (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes
439
    (; inverse_weights) = dg.basis # Plays role of DG subcell sizes
440
    (; positivity_correction_factor) = limiter
441

442
    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
443
    var_min = variable_bounds[Symbol(string(variable), "_min")]
444

445
    @threaded for element in eachelement(dg, semi.cache)
446
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
447
            inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian,
448
                                                    mesh, i, j, k, element)
449
            var = u[variable, i, j, k, element]
450
            if var < 0
451
                error("Safe low-order method produces negative value for conservative variable $variable. Try a smaller time step.")
452
            end
453

454
            # Compute bound
455
            if limiter.local_twosided &&
456
               (variable in limiter.local_twosided_variables_cons) &&
457
               (var_min[i, j, k, element] >= positivity_correction_factor * var)
458
                # Local limiting is more restrictive that positivity limiting
459
                # => Skip positivity limiting for this node
460
                continue
461
            end
462
            var_min[i, j, k, element] = positivity_correction_factor * var
463

464
            # Real one-sided Zalesak-type limiter
465
            # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids"
466
            # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics"
467
            # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is
468
            #       for each interface, not each node
469
            Qm = min(0, (var_min[i, j, k, element] - var) / dt)
470

471
            # Calculate Pm
472
            # Note: Boundaries of antidiffusive_flux1/2/3 are constant 0, so they make no difference here.
473
            val_flux1_local = inverse_weights[i] *
474
                              antidiffusive_flux1_R[variable, i, j, k, element]
475
            val_flux1_local_ip1 = -inverse_weights[i] *
476
                                  antidiffusive_flux1_L[variable, i + 1, j, k, element]
477
            val_flux2_local = inverse_weights[j] *
478
                              antidiffusive_flux2_R[variable, i, j, k, element]
479
            val_flux2_local_jp1 = -inverse_weights[j] *
480
                                  antidiffusive_flux2_L[variable, i, j + 1, k, element]
481
            val_flux3_local = inverse_weights[k] *
482
                              antidiffusive_flux3_R[variable, i, j, k, element]
483
            val_flux3_local_jp1 = -inverse_weights[k] *
484
                                  antidiffusive_flux3_L[variable, i, j, k + 1, element]
485

486
            Pm = min(0, val_flux1_local) + min(0, val_flux1_local_ip1) +
487
                 min(0, val_flux2_local) + min(0, val_flux2_local_jp1) +
488
                 min(0, val_flux3_local) + min(0, val_flux3_local_jp1)
489
            Pm = inverse_jacobian * Pm
490

491
            # Compute blending coefficient avoiding division by zero
492
            # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8))
493
            Qm = abs(Qm) / (abs(Pm) + eps(typeof(Qm)) * 100)
494

495
            # Calculate alpha
496
            alpha[i, j, k, element] = max(alpha[i, j, k, element], 1 - Qm)
497
        end
498
    end
499

500
    return nothing
501
end
502

503
###############################################################################
504
# Global positivity limiting of nonlinear variables
505

506
@inline function idp_positivity_nonlinear!(alpha, limiter,
507
                                           u::AbstractArray{<:Real, 5}, dt, semi,
508
                                           variable)
509
    mesh, equations, dg, cache = mesh_equations_solver_cache(semi)
510
    (; positivity_correction_factor) = limiter
511

512
    (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
513
    var_min = variable_bounds[Symbol(string(variable), "_min")]
514

515
    @threaded for element in eachelement(dg, semi.cache)
516
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
517
            inverse_jacobian = get_inverse_jacobian(cache.elements.inverse_jacobian,
518
                                                    mesh, i, j, k, element)
519

520
            # Compute bound
521
            u_local = get_node_vars(u, equations, dg, i, j, k, element)
522
            var = variable(u_local, equations)
523
            if var < 0
524
                error("Safe low-order method produces negative value for variable $variable. Try a smaller time step.")
525
            end
526
            var_min[i, j, k, element] = positivity_correction_factor * var
527

528
            # Perform Newton's bisection method to find new alpha
529
            newton_loops_alpha!(alpha, var_min[i, j, k, element],
530
                                u_local, i, j, k, element,
531
                                variable, min,
532
                                initial_check_nonnegative_newton_idp,
533
                                final_check_nonnegative_newton_idp,
534
                                inverse_jacobian, dt, equations, dg, cache, limiter)
535
        end
536
    end
537

538
    return nothing
539
end
540

541
###############################################################################
542
# Newton-bisection method
543

544
@inline function newton_loops_alpha!(alpha, bound, u, i, j, k, element,
545
                                     variable, min_or_max,
546
                                     initial_check, final_check,
547
                                     inverse_jacobian, dt,
548
                                     equations::AbstractEquations{3},
549
                                     dg, cache, limiter)
550
    (; inverse_weights) = dg.basis # Plays role of inverse DG-subcell sizes
551
    (; antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R) = cache.antidiffusive_fluxes
552

553
    (; gamma_constant_newton) = limiter
554

555
    indices = (i, j, k, element)
556

557
    # negative xi direction
558
    antidiffusive_flux = gamma_constant_newton * inverse_jacobian *
559
                         inverse_weights[i] *
560
                         get_node_vars(antidiffusive_flux1_R, equations, dg,
561
                                       i, j, k, element)
562
    newton_loop!(alpha, bound, u, indices, variable, min_or_max,
563
                 initial_check, final_check, equations, dt, limiter, antidiffusive_flux)
564

565
    # positive xi direction
566
    antidiffusive_flux = -gamma_constant_newton * inverse_jacobian *
567
                         inverse_weights[i] *
568
                         get_node_vars(antidiffusive_flux1_L, equations, dg,
569
                                       i + 1, j, k, element)
570
    newton_loop!(alpha, bound, u, indices, variable, min_or_max,
571
                 initial_check, final_check, equations, dt, limiter, antidiffusive_flux)
572

573
    # negative eta direction
574
    antidiffusive_flux = gamma_constant_newton * inverse_jacobian *
575
                         inverse_weights[j] *
576
                         get_node_vars(antidiffusive_flux2_R, equations, dg,
577
                                       i, j, k, element)
578
    newton_loop!(alpha, bound, u, indices, variable, min_or_max,
579
                 initial_check, final_check, equations, dt, limiter, antidiffusive_flux)
580

581
    # positive eta direction
582
    antidiffusive_flux = -gamma_constant_newton * inverse_jacobian *
583
                         inverse_weights[j] *
584
                         get_node_vars(antidiffusive_flux2_L, equations, dg,
585
                                       i, j + 1, k, element)
586
    newton_loop!(alpha, bound, u, indices, variable, min_or_max,
587
                 initial_check, final_check, equations, dt, limiter, antidiffusive_flux)
588

589
    # negative zeta direction
590
    antidiffusive_flux = gamma_constant_newton * inverse_jacobian *
591
                         inverse_weights[k] *
592
                         get_node_vars(antidiffusive_flux3_R, equations, dg,
593
                                       i, j, k, element)
594
    newton_loop!(alpha, bound, u, indices, variable, min_or_max,
595
                 initial_check, final_check, equations, dt, limiter, antidiffusive_flux)
596

597
    # positive zeta direction
598
    antidiffusive_flux = -gamma_constant_newton * inverse_jacobian *
599
                         inverse_weights[k] *
600
                         get_node_vars(antidiffusive_flux3_L, equations, dg,
601
                                       i, j, k + 1, element)
602
    newton_loop!(alpha, bound, u, indices, variable, min_or_max,
603
                 initial_check, final_check, equations, dt, limiter, antidiffusive_flux)
604

605
    return nothing
606
end
607
end # @muladd
608

609
Product

Resources

Company
