CoCalc -- dg_2d_compressible

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/solvers/dgsem_structured/dg_2d_compressible_euler.jl
⁵⁵⁹⁰ views
1

2
# From here on, this file contains specializations of DG methods on the
3
# curved 3D meshes `StructuredMesh{3}, P4estMesh{3}` to the compressible
4
# Euler equations.
5
#
6
# The specialized methods contain relatively verbose and ugly code in the sense
7
# that we do not use the common "pointwise physics" interface of Trixi.jl.
8
# However, this is currently (November 2021) necessary to get a significant
9
# speed-up by using SIMD instructions via LoopVectorization.jl.
10
#
11
# TODO: SIMD. We could get even more speed-up if we did not need to permute
12
#             array dimensions, i.e., if we changed the basic memory layout.
13
#
14
# We do not wrap this code in `@muladd begin ... end` block. Optimizations like
15
# this are handled automatically by LoopVectorization.jl.
16

17
# We specialize on `PtrArray` since these will be returned by `Trixi.wrap_array`
18
# if LoopVectorization.jl can handle the array types. This ensures that `@turbo`
19
# works efficiently here.
20
@inline function flux_differencing_kernel!(_du::PtrArray, u_cons::PtrArray,
21
                                           element,
22
                                           MeshT::Type{<:Union{StructuredMesh{2},
23
                                                               UnstructuredMesh2D,
24
                                                               P4estMesh{2}}},
25
                                           have_nonconservative_terms::False,
26
                                           equations::CompressibleEulerEquations2D,
27
                                           volume_flux::typeof(flux_shima_etal_turbo),
28
                                           dg::DGSEM, cache, alpha)
29
    @unpack derivative_split = dg.basis
30
    @unpack contravariant_vectors = cache.elements
31

32
    # Create a temporary array that will be used to store the RHS with permuted
33
    # indices `[i, j, v]` to allow using SIMD instructions.
34
    # `StrideArray`s with purely static dimensions do not allocate on the heap.
35
    du = StrideArray{eltype(u_cons)}(undef,
36
                                     (ntuple(_ -> StaticInt(nnodes(dg)), ndims(MeshT))...,
37
                                      StaticInt(nvariables(equations))))
38

39
    # Convert conserved to primitive variables on the given `element`.
40
    u_prim = StrideArray{eltype(u_cons)}(undef,
41
                                         (ntuple(_ -> StaticInt(nnodes(dg)),
42
                                                 ndims(MeshT))...,
43
                                          StaticInt(nvariables(equations))))
44

45
    @turbo for j in eachnode(dg), i in eachnode(dg)
46
        rho = u_cons[1, i, j, element]
47
        rho_v1 = u_cons[2, i, j, element]
48
        rho_v2 = u_cons[3, i, j, element]
49
        rho_e_total = u_cons[4, i, j, element]
50

51
        v1 = rho_v1 / rho
52
        v2 = rho_v2 / rho
53
        p = (equations.gamma - 1) * (rho_e_total - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
54

55
        u_prim[i, j, 1] = rho
56
        u_prim[i, j, 2] = v1
57
        u_prim[i, j, 3] = v2
58
        u_prim[i, j, 4] = p
59
    end
60

61
    # x direction
62
    # At first, we create new temporary arrays with permuted memory layout to
63
    # allow using SIMD instructions along the first dimension (which is contiguous
64
    # in memory).
65
    du_permuted = StrideArray{eltype(u_cons)}(undef,
66
                                              (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
67
                                               StaticInt(nvariables(equations))))
68

69
    u_prim_permuted = StrideArray{eltype(u_cons)}(undef,
70
                                                  (StaticInt(nnodes(dg)),
71
                                                   StaticInt(nnodes(dg)),
72
                                                   StaticInt(nvariables(equations))))
73

74
    @turbo for v in eachvariable(equations),
75
               j in eachnode(dg),
76
               i in eachnode(dg)
77

78
        u_prim_permuted[j, i, v] = u_prim[i, j, v]
79
    end
80
    fill!(du_permuted, zero(eltype(du_permuted)))
81

82
    # We must also permute the contravariant vectors.
83
    contravariant_vectors_x = StrideArray{eltype(contravariant_vectors)}(undef,
84
                                                                         (StaticInt(nnodes(dg)),
85
                                                                          StaticInt(nnodes(dg)),
86
                                                                          StaticInt(ndims(MeshT))))
87

88
    @turbo for j in eachnode(dg), i in eachnode(dg)
89
        contravariant_vectors_x[j, i, 1] = contravariant_vectors[1, 1, i, j, element]
90
        contravariant_vectors_x[j, i, 2] = contravariant_vectors[2, 1, i, j, element]
91
    end
92

93
    # Next, we basically inline the volume flux. To allow SIMD vectorization and
94
    # still use the symmetry of the volume flux and the derivative matrix, we
95
    # loop over the triangular part in an outer loop and use a plain inner loop.
96
    for i in eachnode(dg), ii in (i + 1):nnodes(dg)
97
        @turbo for j in eachnode(dg)
98
            rho_ll = u_prim_permuted[j, i, 1]
99
            v1_ll = u_prim_permuted[j, i, 2]
100
            v2_ll = u_prim_permuted[j, i, 3]
101
            p_ll = u_prim_permuted[j, i, 4]
102

103
            rho_rr = u_prim_permuted[j, ii, 1]
104
            v1_rr = u_prim_permuted[j, ii, 2]
105
            v2_rr = u_prim_permuted[j, ii, 3]
106
            p_rr = u_prim_permuted[j, ii, 4]
107

108
            normal_direction_1 = 0.5 * (contravariant_vectors_x[j, i, 1] +
109
                                  contravariant_vectors_x[j, ii, 1])
110
            normal_direction_2 = 0.5 * (contravariant_vectors_x[j, i, 2] +
111
                                  contravariant_vectors_x[j, ii, 2])
112

113
            v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2
114
            v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2
115

116
            # Compute required mean values
117
            rho_avg = 0.5 * (rho_ll + rho_rr)
118
            v1_avg = 0.5 * (v1_ll + v1_rr)
119
            v2_avg = 0.5 * (v2_ll + v2_rr)
120
            v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)
121
            p_avg = 0.5 * (p_ll + p_rr)
122
            velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)
123

124
            # Calculate fluxes depending on normal_direction
125
            f1 = rho_avg * v_dot_n_avg
126
            f2 = f1 * v1_avg + p_avg * normal_direction_1
127
            f3 = f1 * v2_avg + p_avg * normal_direction_2
128
            f4 = (f1 * velocity_square_avg +
129
                  p_avg * v_dot_n_avg * equations.inv_gamma_minus_one
130
                  + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll))
131

132
            # Add scaled fluxes to RHS
133
            factor_i = alpha * derivative_split[i, ii]
134
            du_permuted[j, i, 1] += factor_i * f1
135
            du_permuted[j, i, 2] += factor_i * f2
136
            du_permuted[j, i, 3] += factor_i * f3
137
            du_permuted[j, i, 4] += factor_i * f4
138

139
            factor_ii = alpha * derivative_split[ii, i]
140
            du_permuted[j, ii, 1] += factor_ii * f1
141
            du_permuted[j, ii, 2] += factor_ii * f2
142
            du_permuted[j, ii, 3] += factor_ii * f3
143
            du_permuted[j, ii, 4] += factor_ii * f4
144
        end
145
    end
146

147
    @turbo for v in eachvariable(equations),
148
               j in eachnode(dg),
149
               i in eachnode(dg)
150

151
        du[i, j, v] = du_permuted[j, i, v]
152
    end
153

154
    # y direction
155
    # We must also permute the contravariant vectors.
156
    contravariant_vectors_y = StrideArray{eltype(contravariant_vectors)}(undef,
157
                                                                         (StaticInt(nnodes(dg)),
158
                                                                          StaticInt(nnodes(dg)),
159
                                                                          StaticInt(ndims(MeshT))))
160

161
    @turbo for j in eachnode(dg), i in eachnode(dg)
162
        contravariant_vectors_y[i, j, 1] = contravariant_vectors[1, 2, i, j, element]
163
        contravariant_vectors_y[i, j, 2] = contravariant_vectors[2, 2, i, j, element]
164
    end
165

166
    # The memory layout is already optimal for SIMD vectorization in this loop.
167
    for j in eachnode(dg), jj in (j + 1):nnodes(dg)
168
        @turbo for i in eachnode(dg)
169
            rho_ll = u_prim[i, j, 1]
170
            v1_ll = u_prim[i, j, 2]
171
            v2_ll = u_prim[i, j, 3]
172
            p_ll = u_prim[i, j, 4]
173

174
            rho_rr = u_prim[i, jj, 1]
175
            v1_rr = u_prim[i, jj, 2]
176
            v2_rr = u_prim[i, jj, 3]
177
            p_rr = u_prim[i, jj, 4]
178

179
            normal_direction_1 = 0.5 * (contravariant_vectors_y[i, j, 1] +
180
                                  contravariant_vectors_y[i, jj, 1])
181
            normal_direction_2 = 0.5 * (contravariant_vectors_y[i, j, 2] +
182
                                  contravariant_vectors_y[i, jj, 2])
183

184
            v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2
185
            v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2
186

187
            # Compute required mean values
188
            rho_avg = 0.5 * (rho_ll + rho_rr)
189
            v1_avg = 0.5 * (v1_ll + v1_rr)
190
            v2_avg = 0.5 * (v2_ll + v2_rr)
191
            v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)
192
            p_avg = 0.5 * (p_ll + p_rr)
193
            velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)
194

195
            # Calculate fluxes depending on normal_direction
196
            f1 = rho_avg * v_dot_n_avg
197
            f2 = f1 * v1_avg + p_avg * normal_direction_1
198
            f3 = f1 * v2_avg + p_avg * normal_direction_2
199
            f4 = (f1 * velocity_square_avg +
200
                  p_avg * v_dot_n_avg * equations.inv_gamma_minus_one
201
                  + 0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll))
202

203
            # Add scaled fluxes to RHS
204
            factor_j = alpha * derivative_split[j, jj]
205
            du[i, j, 1] += factor_j * f1
206
            du[i, j, 2] += factor_j * f2
207
            du[i, j, 3] += factor_j * f3
208
            du[i, j, 4] += factor_j * f4
209

210
            factor_jj = alpha * derivative_split[jj, j]
211
            du[i, jj, 1] += factor_jj * f1
212
            du[i, jj, 2] += factor_jj * f2
213
            du[i, jj, 3] += factor_jj * f3
214
            du[i, jj, 4] += factor_jj * f4
215
        end
216
    end
217

218
    # Finally, we add the temporary RHS computed here to the global RHS in the
219
    # given `element`.
220
    @turbo for v in eachvariable(equations),
221
               j in eachnode(dg),
222
               i in eachnode(dg)
223

224
        _du[v, i, j, element] += du[i, j, v]
225
    end
226
end
227

228
@inline function flux_differencing_kernel!(_du::PtrArray, u_cons::PtrArray,
229
                                           element,
230
                                           MeshT::Type{<:Union{StructuredMesh{2},
231
                                                               UnstructuredMesh2D,
232
                                                               P4estMesh{2}}},
233
                                           have_nonconservative_terms::False,
234
                                           equations::CompressibleEulerEquations2D,
235
                                           volume_flux::typeof(flux_ranocha_turbo),
236
                                           dg::DGSEM, cache, alpha)
237
    @unpack derivative_split = dg.basis
238
    @unpack contravariant_vectors = cache.elements
239

240
    # Create a temporary array that will be used to store the RHS with permuted
241
    # indices `[i, j, v]` to allow using SIMD instructions.
242
    # `StrideArray`s with purely static dimensions do not allocate on the heap.
243
    du = StrideArray{eltype(u_cons)}(undef,
244
                                     (ntuple(_ -> StaticInt(nnodes(dg)), ndims(MeshT))...,
245
                                      StaticInt(nvariables(equations))))
246

247
    # Convert conserved to primitive variables on the given `element`. In addition
248
    # to the usual primitive variables, we also compute logarithms of the density
249
    # and pressure to increase the performance of the required logarithmic mean
250
    # values.
251
    u_prim = StrideArray{eltype(u_cons)}(undef,
252
                                         (ntuple(_ -> StaticInt(nnodes(dg)),
253
                                                 ndims(MeshT))...,
254
                                          StaticInt(nvariables(equations) + 2))) # We also compute "+ 2" logs
255

256
    @turbo for j in eachnode(dg), i in eachnode(dg)
257
        rho = u_cons[1, i, j, element]
258
        rho_v1 = u_cons[2, i, j, element]
259
        rho_v2 = u_cons[3, i, j, element]
260
        rho_e_total = u_cons[4, i, j, element]
261

262
        v1 = rho_v1 / rho
263
        v2 = rho_v2 / rho
264
        p = (equations.gamma - 1) * (rho_e_total - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
265

266
        u_prim[i, j, 1] = rho
267
        u_prim[i, j, 2] = v1
268
        u_prim[i, j, 3] = v2
269
        u_prim[i, j, 4] = p
270
        u_prim[i, j, 5] = log(rho)
271
        u_prim[i, j, 6] = log(p)
272
    end
273

274
    # x direction
275
    # At first, we create new temporary arrays with permuted memory layout to
276
    # allow using SIMD instructions along the first dimension (which is contiguous
277
    # in memory).
278
    du_permuted = StrideArray{eltype(u_cons)}(undef,
279
                                              (StaticInt(nnodes(dg)), StaticInt(nnodes(dg)),
280
                                               StaticInt(nvariables(equations))))
281

282
    u_prim_permuted = StrideArray{eltype(u_cons)}(undef,
283
                                                  (StaticInt(nnodes(dg)),
284
                                                   StaticInt(nnodes(dg)),
285
                                                   StaticInt(nvariables(equations) + 2)))
286

287
    @turbo for v in indices(u_prim, 3), # v in eachvariable(equations) misses +2 logs
288
               j in eachnode(dg),
289
               i in eachnode(dg)
290

291
        u_prim_permuted[j, i, v] = u_prim[i, j, v]
292
    end
293
    fill!(du_permuted, zero(eltype(du_permuted)))
294

295
    # We must also permute the contravariant vectors.
296
    contravariant_vectors_x = StrideArray{eltype(contravariant_vectors)}(undef,
297
                                                                         (StaticInt(nnodes(dg)),
298
                                                                          StaticInt(nnodes(dg)),
299
                                                                          StaticInt(ndims(MeshT))))
300

301
    @turbo for j in eachnode(dg), i in eachnode(dg)
302
        contravariant_vectors_x[j, i, 1] = contravariant_vectors[1, 1, i, j, element]
303
        contravariant_vectors_x[j, i, 2] = contravariant_vectors[2, 1, i, j, element]
304
    end
305

306
    # Next, we basically inline the volume flux. To allow SIMD vectorization and
307
    # still use the symmetry of the volume flux and the derivative matrix, we
308
    # loop over the triangular part in an outer loop and use a plain inner loop.
309
    for i in eachnode(dg), ii in (i + 1):nnodes(dg)
310
        @turbo for j in eachnode(dg)
311
            rho_ll = u_prim_permuted[j, i, 1]
312
            v1_ll = u_prim_permuted[j, i, 2]
313
            v2_ll = u_prim_permuted[j, i, 3]
314
            p_ll = u_prim_permuted[j, i, 4]
315
            log_rho_ll = u_prim_permuted[j, i, 5]
316
            log_p_ll = u_prim_permuted[j, i, 6]
317

318
            rho_rr = u_prim_permuted[j, ii, 1]
319
            v1_rr = u_prim_permuted[j, ii, 2]
320
            v2_rr = u_prim_permuted[j, ii, 3]
321
            p_rr = u_prim_permuted[j, ii, 4]
322
            log_rho_rr = u_prim_permuted[j, ii, 5]
323
            log_p_rr = u_prim_permuted[j, ii, 6]
324

325
            normal_direction_1 = 0.5 * (contravariant_vectors_x[j, i, 1] +
326
                                  contravariant_vectors_x[j, ii, 1])
327
            normal_direction_2 = 0.5 * (contravariant_vectors_x[j, i, 2] +
328
                                  contravariant_vectors_x[j, ii, 2])
329

330
            v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2
331
            v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2
332

333
            # Compute required mean values
334
            # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
335
            # it efficiently. This is equivalent to
336
            #   rho_mean = ln_mean(rho_ll, rho_rr)
337
            x1 = rho_ll
338
            log_x1 = log_rho_ll
339
            y1 = rho_rr
340
            log_y1 = log_rho_rr
341
            x1_plus_y1 = x1 + y1
342
            y1_minus_x1 = y1 - x1
343
            z1 = y1_minus_x1^2 / x1_plus_y1^2
344
            special_path1 = x1_plus_y1 / (2 + z1 * (2 / 3 + z1 * (2 / 5 + 2 / 7 * z1)))
345
            regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
346
            rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)
347

348
            # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
349
            # in exact arithmetic since
350
            #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
351
            #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
352
            # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
353
            x2 = rho_ll * p_rr
354
            log_x2 = log_rho_ll + log_p_rr
355
            y2 = rho_rr * p_ll
356
            log_y2 = log_rho_rr + log_p_ll
357
            x2_plus_y2 = x2 + y2
358
            y2_minus_x2 = y2 - x2
359
            z2 = y2_minus_x2^2 / x2_plus_y2^2
360
            special_path2 = (2 + z2 * (2 / 3 + z2 * (2 / 5 + 2 / 7 * z2))) / x2_plus_y2
361
            regular_path2 = (log_y2 - log_x2) / y2_minus_x2
362
            inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)
363

364
            v1_avg = 0.5 * (v1_ll + v1_rr)
365
            v2_avg = 0.5 * (v2_ll + v2_rr)
366
            p_avg = 0.5 * (p_ll + p_rr)
367
            velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)
368

369
            # Calculate fluxes depending on normal_direction
370
            f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
371
            f2 = f1 * v1_avg + p_avg * normal_direction_1
372
            f3 = f1 * v2_avg + p_avg * normal_direction_2
373
            f4 = (f1 *
374
                  (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one)
375
                  +
376
                  0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll))
377

378
            # Add scaled fluxes to RHS
379
            factor_i = alpha * derivative_split[i, ii]
380
            du_permuted[j, i, 1] += factor_i * f1
381
            du_permuted[j, i, 2] += factor_i * f2
382
            du_permuted[j, i, 3] += factor_i * f3
383
            du_permuted[j, i, 4] += factor_i * f4
384

385
            factor_ii = alpha * derivative_split[ii, i]
386
            du_permuted[j, ii, 1] += factor_ii * f1
387
            du_permuted[j, ii, 2] += factor_ii * f2
388
            du_permuted[j, ii, 3] += factor_ii * f3
389
            du_permuted[j, ii, 4] += factor_ii * f4
390
        end
391
    end
392

393
    @turbo for v in eachvariable(equations),
394
               j in eachnode(dg),
395
               i in eachnode(dg)
396

397
        du[i, j, v] = du_permuted[j, i, v]
398
    end
399

400
    # y direction
401
    # We must also permute the contravariant vectors.
402
    contravariant_vectors_y = StrideArray{eltype(contravariant_vectors)}(undef,
403
                                                                         (StaticInt(nnodes(dg)),
404
                                                                          StaticInt(nnodes(dg)),
405
                                                                          StaticInt(ndims(MeshT))))
406

407
    @turbo for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
408
        contravariant_vectors_y[i, j, 1] = contravariant_vectors[1, 2, i, j, element]
409
        contravariant_vectors_y[i, j, 2] = contravariant_vectors[2, 2, i, j, element]
410
    end
411

412
    # The memory layout is already optimal for SIMD vectorization in this loop.
413
    for j in eachnode(dg), jj in (j + 1):nnodes(dg)
414
        @turbo for i in eachnode(dg)
415
            rho_ll = u_prim[i, j, 1]
416
            v1_ll = u_prim[i, j, 2]
417
            v2_ll = u_prim[i, j, 3]
418
            p_ll = u_prim[i, j, 4]
419
            log_rho_ll = u_prim[i, j, 5]
420
            log_p_ll = u_prim[i, j, 6]
421

422
            rho_rr = u_prim[i, jj, 1]
423
            v1_rr = u_prim[i, jj, 2]
424
            v2_rr = u_prim[i, jj, 3]
425
            p_rr = u_prim[i, jj, 4]
426
            log_rho_rr = u_prim[i, jj, 5]
427
            log_p_rr = u_prim[i, jj, 6]
428

429
            normal_direction_1 = 0.5 * (contravariant_vectors_y[i, j, 1] +
430
                                  contravariant_vectors_y[i, jj, 1])
431
            normal_direction_2 = 0.5 * (contravariant_vectors_y[i, j, 2] +
432
                                  contravariant_vectors_y[i, jj, 2])
433

434
            v_dot_n_ll = v1_ll * normal_direction_1 + v2_ll * normal_direction_2
435
            v_dot_n_rr = v1_rr * normal_direction_1 + v2_rr * normal_direction_2
436

437
            # Compute required mean values
438
            # We inline the logarithmic mean to allow LoopVectorization.jl to optimize
439
            # it efficiently. This is equivalent to
440
            #   rho_mean = ln_mean(rho_ll, rho_rr)
441
            x1 = rho_ll
442
            log_x1 = log_rho_ll
443
            y1 = rho_rr
444
            log_y1 = log_rho_rr
445
            x1_plus_y1 = x1 + y1
446
            y1_minus_x1 = y1 - x1
447
            z1 = y1_minus_x1^2 / x1_plus_y1^2
448
            special_path1 = x1_plus_y1 / (2 + z1 * (2 / 3 + z1 * (2 / 5 + 2 / 7 * z1)))
449
            regular_path1 = y1_minus_x1 / (log_y1 - log_x1)
450
            rho_mean = ifelse(z1 < 1.0e-4, special_path1, regular_path1)
451

452
            # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
453
            # in exact arithmetic since
454
            #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
455
            #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
456
            # inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
457
            x2 = rho_ll * p_rr
458
            log_x2 = log_rho_ll + log_p_rr
459
            y2 = rho_rr * p_ll
460
            log_y2 = log_rho_rr + log_p_ll
461
            x2_plus_y2 = x2 + y2
462
            y2_minus_x2 = y2 - x2
463
            z2 = y2_minus_x2^2 / x2_plus_y2^2
464
            special_path2 = (2 + z2 * (2 / 3 + z2 * (2 / 5 + 2 / 7 * z2))) / x2_plus_y2
465
            regular_path2 = (log_y2 - log_x2) / y2_minus_x2
466
            inv_rho_p_mean = p_ll * p_rr * ifelse(z2 < 1.0e-4, special_path2, regular_path2)
467

468
            v1_avg = 0.5 * (v1_ll + v1_rr)
469
            v2_avg = 0.5 * (v2_ll + v2_rr)
470
            p_avg = 0.5 * (p_ll + p_rr)
471
            velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr)
472

473
            # Calculate fluxes depending on normal_direction
474
            f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
475
            f2 = f1 * v1_avg + p_avg * normal_direction_1
476
            f3 = f1 * v2_avg + p_avg * normal_direction_2
477
            f4 = (f1 *
478
                  (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one)
479
                  +
480
                  0.5 * (p_ll * v_dot_n_rr + p_rr * v_dot_n_ll))
481

482
            # Add scaled fluxes to RHS
483
            factor_j = alpha * derivative_split[j, jj]
484
            du[i, j, 1] += factor_j * f1
485
            du[i, j, 2] += factor_j * f2
486
            du[i, j, 3] += factor_j * f3
487
            du[i, j, 4] += factor_j * f4
488

489
            factor_jj = alpha * derivative_split[jj, j]
490
            du[i, jj, 1] += factor_jj * f1
491
            du[i, jj, 2] += factor_jj * f2
492
            du[i, jj, 3] += factor_jj * f3
493
            du[i, jj, 4] += factor_jj * f4
494
        end
495
    end
496

497
    # Finally, we add the temporary RHS computed here to the global RHS in the
498
    # given `element`.
499
    @turbo for v in eachvariable(equations),
500
               j in eachnode(dg),
501
               i in eachnode(dg)
502

503
        _du[v, i, j, element] += du[i, j, v]
504
    end
505
end
506

507
Product

Resources

Company
