Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
trixi-framework
GitHub Repository: trixi-framework/Trixi.jl
 Path: blob/main/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl
5591 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
# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element

9
# (**without non-conservative terms**).

10
#

11
# See also `flux_differencing_kernel!`.

12
@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u,

13
                                ::Type{<:Union{StructuredMesh{2}, P4estMesh{2}}},

14
                                have_nonconservative_terms::False, equations,

15
                                volume_flux, dg::DGSEM, element, cache)

16
    (; contravariant_vectors) = cache.elements

17
    (; weights, derivative_split) = dg.basis

18
    (; flux_temp_threaded) = cache

19


20
    flux_temp = flux_temp_threaded[Threads.threadid()]

21


22
    # The FV-form fluxes are calculated in a recursive manner, i.e.:

23
    # fhat_(0,1)   = w_0 * FVol_0,

24
    # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j,   for j=1,...,N-1,

25
    # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i).

26


27
    # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated

28
    # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing`

29
    # and saved in in `flux_temp`.

30


31
    # Split form volume flux in orientation 1: x direction

32
    flux_temp .= zero(eltype(flux_temp))

33


34
    for j in eachnode(dg), i in eachnode(dg)

35
        u_node = get_node_vars(u, equations, dg, i, j, element)

36


37
        # pull the contravariant vectors in each coordinate direction

38
        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction

39


40
        # All diagonal entries of `derivative_split` are zero. Thus, we can skip

41
        # the computation of the diagonal terms. In addition, we use the symmetry

42
        # of the `volume_flux` to save half of the possible two-point flux

43
        # computations.

44


45
        # x direction

46
        for ii in (i + 1):nnodes(dg)

47
            u_node_ii = get_node_vars(u, equations, dg, ii, j, element)

48
            # pull the contravariant vectors and compute the average

49
            Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j,

50
                                                   element)

51
            Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)

52


53
            # compute the contravariant sharp flux in the direction of the averaged contravariant vector

54
            fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations)

55
            multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1,

56
                                       equations, dg, i, j)

57
            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,

58
                                       equations, dg, ii, j)

59
        end

60
    end

61


62
    # FV-form flux `fhat` in x direction

63
    for j in eachnode(dg), i in 1:(nnodes(dg) - 1)

64
        for v in eachvariable(equations)

65
            fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j]

66
            fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j]

67
        end

68
    end

69


70
    # Split form volume flux in orientation 2: y direction

71
    flux_temp .= zero(eltype(flux_temp))

72


73
    for j in eachnode(dg), i in eachnode(dg)

74
        u_node = get_node_vars(u, equations, dg, i, j, element)

75


76
        # pull the contravariant vectors in each coordinate direction

77
        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)

78


79
        # y direction

80
        for jj in (j + 1):nnodes(dg)

81
            u_node_jj = get_node_vars(u, equations, dg, i, jj, element)

82
            # pull the contravariant vectors and compute the average

83
            Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj,

84
                                                   element)

85
            Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)

86
            # compute the contravariant sharp flux in the direction of the averaged contravariant vector

87
            fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations)

88
            multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2,

89
                                       equations, dg, i, j)

90
            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,

91
                                       equations, dg, i, jj)

92
        end

93
    end

94


95
    # FV-form flux `fhat` in y direction

96
    for j in 1:(nnodes(dg) - 1), i in eachnode(dg)

97
        for v in eachvariable(equations)

98
            fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j]

99
            fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1]

100
        end

101
    end

102


103
    return nothing

104
end

105


106
# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element

107
# (**with non-conservative terms in "local * symmetric" form**).

108
#

109
# See also `flux_differencing_kernel!`.

110
#

111
# The calculation of the non-conservative staggered "fluxes" requires non-conservative

112
# terms that can be written as a product of local and a symmetric contributions. See, e.g.,

113
#

114
# - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts

115
#   Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.

116
#

117
@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u,

118
                                ::Type{<:Union{StructuredMesh{2}, P4estMesh{2}}},

119
                                have_nonconservative_terms::True, equations,

120
                                volume_flux::Tuple{F_CONS, F_NONCONS}, dg::DGSEM,

121
                                element,

122
                                cache) where {

123
                                              F_CONS <: Function,

124
                                              F_NONCONS <:

125
                                              FluxNonConservative{NonConservativeSymmetric()}

126
                                              }

127
    (; contravariant_vectors) = cache.elements

128
    (; weights, derivative_split) = dg.basis

129
    (; flux_temp_threaded, flux_nonconservative_temp_threaded) = cache

130
    (; fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded) = cache

131


132
    volume_flux_cons, volume_flux_noncons = volume_flux

133


134
    flux_temp = flux_temp_threaded[Threads.threadid()]

135
    flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()]

136


137
    fhat_temp = fhat_temp_threaded[Threads.threadid()]

138
    fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()]

139
    phi = phi_threaded[Threads.threadid()]

140


141
    # The FV-form fluxes are calculated in a recursive manner, i.e.:

142
    # fhat_(0,1)   = w_0 * FVol_0,

143
    # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j,   for j=1,...,N-1,

144
    # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i).

145


146
    # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated

147
    # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing`

148
    # and saved in in `flux_temp`.

149


150
    # Split form volume flux in orientation 1: x direction

151
    flux_temp .= zero(eltype(flux_temp))

152
    flux_noncons_temp .= zero(eltype(flux_noncons_temp))

153


154
    for j in eachnode(dg), i in eachnode(dg)

155
        u_node = get_node_vars(u, equations, dg, i, j, element)

156


157
        # pull the contravariant vectors in each coordinate direction

158
        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction

159


160
        # All diagonal entries of `derivative_split` are zero. Thus, we can skip

161
        # the computation of the diagonal terms. In addition, we use the symmetry

162
        # of `volume_flux_cons` and `volume_flux_noncons` to save half of the possible two-point flux

163
        # computations.

164
        for ii in (i + 1):nnodes(dg)

165
            u_node_ii = get_node_vars(u, equations, dg, ii, j, element)

166
            # pull the contravariant vectors and compute the average

167
            Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j,

168
                                                   element)

169
            Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)

170


171
            # compute the contravariant sharp flux in the direction of the averaged contravariant vector

172
            fluxtilde1 = volume_flux_cons(u_node, u_node_ii, Ja1_avg, equations)

173
            multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1,

174
                                       equations, dg, i, j)

175
            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,

176
                                       equations, dg, ii, j)

177
            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

178
                # We multiply by 0.5 because that is done in other parts of Trixi

179
                flux1_noncons = volume_flux_noncons(u_node, u_node_ii, Ja1_avg,

180
                                                    equations,

181
                                                    NonConservativeSymmetric(), noncons)

182
                multiply_add_to_node_vars!(flux_noncons_temp,

183
                                           0.5f0 * derivative_split[i, ii],

184
                                           flux1_noncons,

185
                                           equations, dg, noncons, i, j)

186
                multiply_add_to_node_vars!(flux_noncons_temp,

187
                                           0.5f0 * derivative_split[ii, i],

188
                                           flux1_noncons,

189
                                           equations, dg, noncons, ii, j)

190
            end

191
        end

192
    end

193


194
    # FV-form flux `fhat` in x direction

195
    fhat_temp[:, 1, :] .= zero(eltype(fhat1_L))

196
    fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L))

197


198
    # Compute local contribution to non-conservative flux

199
    for j in eachnode(dg), i in eachnode(dg)

200
        u_local = get_node_vars(u, equations, dg, i, j, element)

201
        # pull the local contravariant vector

202
        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element)

203
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

204
            set_node_vars!(phi,

205
                           volume_flux_noncons(u_local, Ja1_node, equations,

206
                                               NonConservativeLocal(), noncons),

207
                           equations, dg, noncons, i, j)

208
        end

209
    end

210


211
    for j in eachnode(dg), i in 1:(nnodes(dg) - 1)

212
        # Conservative part

213
        for v in eachvariable(equations)

214
            value = fhat_temp[v, i, j] + weights[i] * flux_temp[v, i, j]

215
            fhat_temp[v, i + 1, j] = value

216
            fhat1_L[v, i + 1, j] = value

217
            fhat1_R[v, i + 1, j] = value

218
        end

219
        # Nonconservative part

220
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),

221
            v in eachvariable(equations)

222


223
            value = fhat_noncons_temp[v, noncons, i, j] +

224
                    weights[i] * flux_noncons_temp[v, noncons, i, j]

225
            fhat_noncons_temp[v, noncons, i + 1, j] = value

226


227
            fhat1_L[v, i + 1, j] = fhat1_L[v, i + 1, j] + phi[v, noncons, i, j] * value

228
            fhat1_R[v, i + 1, j] = fhat1_R[v, i + 1, j] +

229
                                   phi[v, noncons, i + 1, j] * value

230
        end

231
    end

232


233
    # Split form volume flux in orientation 2: y direction

234
    flux_temp .= zero(eltype(flux_temp))

235
    flux_noncons_temp .= zero(eltype(flux_noncons_temp))

236


237
    for j in eachnode(dg), i in eachnode(dg)

238
        u_node = get_node_vars(u, equations, dg, i, j, element)

239


240
        # pull the contravariant vectors in each coordinate direction

241
        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)

242


243
        for jj in (j + 1):nnodes(dg)

244
            u_node_jj = get_node_vars(u, equations, dg, i, jj, element)

245
            # pull the contravariant vectors and compute the average

246
            Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj,

247
                                                   element)

248
            Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)

249
            # compute the contravariant sharp flux in the direction of the averaged contravariant vector

250
            fluxtilde2 = volume_flux_cons(u_node, u_node_jj, Ja2_avg, equations)

251
            multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2,

252
                                       equations, dg, i, j)

253
            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,

254
                                       equations, dg, i, jj)

255
            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

256
                # We multiply by 0.5 because that is done in other parts of Trixi

257
                flux2_noncons = volume_flux_noncons(u_node, u_node_jj, Ja2_avg,

258
                                                    equations,

259
                                                    NonConservativeSymmetric(), noncons)

260
                multiply_add_to_node_vars!(flux_noncons_temp,

261
                                           0.5f0 * derivative_split[j, jj],

262
                                           flux2_noncons,

263
                                           equations, dg, noncons, i, j)

264
                multiply_add_to_node_vars!(flux_noncons_temp,

265
                                           0.5f0 * derivative_split[jj, j],

266
                                           flux2_noncons,

267
                                           equations, dg, noncons, i, jj)

268
            end

269
        end

270
    end

271


272
    # FV-form flux `fhat` in y direction

273
    fhat_temp[:, :, 1] .= zero(eltype(fhat1_L))

274
    fhat_noncons_temp[:, :, :, 1] .= zero(eltype(fhat1_L))

275


276
    # Compute local contribution to non-conservative flux

277
    for j in eachnode(dg), i in eachnode(dg)

278
        u_local = get_node_vars(u, equations, dg, i, j, element)

279
        # pull the local contravariant vector

280
        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)

281
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

282
            set_node_vars!(phi,

283
                           volume_flux_noncons(u_local, Ja2_node, equations,

284
                                               NonConservativeLocal(), noncons),

285
                           equations, dg, noncons, i, j)

286
        end

287
    end

288


289
    for j in 1:(nnodes(dg) - 1), i in eachnode(dg)

290
        # Conservative part

291
        for v in eachvariable(equations)

292
            value = fhat_temp[v, i, j] + weights[j] * flux_temp[v, i, j]

293
            fhat_temp[v, i, j + 1] = value

294
            fhat2_L[v, i, j + 1] = value

295
            fhat2_R[v, i, j + 1] = value

296
        end

297
        # Nonconservative part

298
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),

299
            v in eachvariable(equations)

300


301
            value = fhat_noncons_temp[v, noncons, i, j] +

302
                    weights[j] * flux_noncons_temp[v, noncons, i, j]

303
            fhat_noncons_temp[v, noncons, i, j + 1] = value

304


305
            fhat2_L[v, i, j + 1] = fhat2_L[v, i, j + 1] + phi[v, noncons, i, j] * value

306
            fhat2_R[v, i, j + 1] = fhat2_R[v, i, j + 1] +

307
                                   phi[v, noncons, i, j + 1] * value

308
        end

309
    end

310


311
    return nothing

312
end

313


314
# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element

315
# (**with non-conservative terms in "local * jump" form**).

316
#

317
# See also `flux_differencing_kernel!`.

318
#

319
# The calculation of the non-conservative staggered "fluxes" requires non-conservative

320
# terms that can be written as a product of local and jump contributions.

321
@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u,

322
                                ::Type{<:Union{StructuredMesh{2}, P4estMesh{2}}},

323
                                have_nonconservative_terms::True, equations,

324
                                volume_flux::Tuple{F_CONS, F_NONCONS}, dg::DGSEM,

325
                                element,

326
                                cache) where {

327
                                              F_CONS <: Function,

328
                                              F_NONCONS <:

329
                                              FluxNonConservative{NonConservativeJump()}

330
                                              }

331
    (; contravariant_vectors) = cache.elements

332
    (; weights, derivative_split) = dg.basis

333
    (; flux_temp_threaded, flux_nonconservative_temp_threaded) = cache

334
    (; fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded) = cache

335


336
    volume_flux_cons, volume_flux_noncons = volume_flux

337


338
    flux_temp = flux_temp_threaded[Threads.threadid()]

339
    flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()]

340


341
    fhat_temp = fhat_temp_threaded[Threads.threadid()]

342
    fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()]

343
    phi = phi_threaded[Threads.threadid()]

344


345
    # The FV-form fluxes are calculated in a recursive manner, i.e.:

346
    # fhat_(0,1)   = w_0 * FVol_0,

347
    # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j,   for j=1,...,N-1,

348
    # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i).

349


350
    # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated

351
    # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing`

352
    # and saved in in `flux_temp`.

353


354
    # Split form volume flux in orientation 1: x direction

355
    flux_temp .= zero(eltype(flux_temp))

356
    flux_noncons_temp .= zero(eltype(flux_noncons_temp))

357


358
    for j in eachnode(dg), i in eachnode(dg)

359
        u_node = get_node_vars(u, equations, dg, i, j, element)

360


361
        # pull the contravariant vectors in each coordinate direction

362
        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction

363


364
        # All diagonal entries of `derivative_split` are zero. Thus, we can skip

365
        # the computation of the diagonal terms. In addition, we use the symmetry

366
        # of `volume_flux_cons` and `volume_flux_noncons` to save half of the possible two-point flux

367
        # computations.

368
        for ii in (i + 1):nnodes(dg)

369
            u_node_ii = get_node_vars(u, equations, dg, ii, j, element)

370
            # pull the contravariant vectors and compute the average

371
            Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j,

372
                                                   element)

373
            Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)

374


375
            # compute the contravariant sharp flux in the direction of the averaged contravariant vector

376
            fluxtilde1 = volume_flux_cons(u_node, u_node_ii, Ja1_avg, equations)

377
            multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1,

378
                                       equations, dg, i, j)

379
            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,

380
                                       equations, dg, ii, j)

381
            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

382
                # We multiply by 0.5 because that is done in other parts of Trixi

383
                flux1_noncons = volume_flux_noncons(u_node, u_node_ii, Ja1_avg,

384
                                                    equations,

385
                                                    NonConservativeJump(), noncons)

386
                multiply_add_to_node_vars!(flux_noncons_temp,

387
                                           0.5f0 * derivative_split[i, ii],

388
                                           flux1_noncons,

389
                                           equations, dg, noncons, i, j)

390
                multiply_add_to_node_vars!(flux_noncons_temp,

391
                                           -0.5f0 * derivative_split[ii, i],

392
                                           flux1_noncons,

393
                                           equations, dg, noncons, ii, j)

394
            end

395
        end

396
    end

397


398
    # FV-form flux `fhat` in x direction

399
    fhat_temp[:, 1, :] .= zero(eltype(fhat1_L))

400
    fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L))

401


402
    # Compute local contribution to non-conservative flux

403
    for j in eachnode(dg), i in eachnode(dg)

404
        u_local = get_node_vars(u, equations, dg, i, j, element)

405
        # pull the local contravariant vector

406
        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element)

407
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

408
            set_node_vars!(phi,

409
                           volume_flux_noncons(u_local, Ja1_node, equations,

410
                                               NonConservativeLocal(), noncons),

411
                           equations, dg, noncons, i, j)

412
        end

413
    end

414


415
    for j in eachnode(dg), i in 1:(nnodes(dg) - 1)

416
        # Conservative part

417
        for v in eachvariable(equations)

418
            value = fhat_temp[v, i, j] + weights[i] * flux_temp[v, i, j]

419
            fhat_temp[v, i + 1, j] = value

420
            fhat1_L[v, i + 1, j] = value

421
            fhat1_R[v, i + 1, j] = value

422
        end

423
        # Nonconservative part

424
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),

425
            v in eachvariable(equations)

426


427
            value = fhat_noncons_temp[v, noncons, i, j] +

428
                    weights[i] * flux_noncons_temp[v, noncons, i, j]

429
            fhat_noncons_temp[v, noncons, i + 1, j] = value

430


431
            fhat1_L[v, i + 1, j] = fhat1_L[v, i + 1, j] + phi[v, noncons, i, j] * value

432
            fhat1_R[v, i + 1, j] = fhat1_R[v, i + 1, j] +

433
                                   phi[v, noncons, i + 1, j] * value

434
        end

435
    end

436


437
    # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes.

438
    for j in eachnode(dg)

439
        u_0 = get_node_vars(u, equations, dg, 1, j, element)

440
        Ja1_node_0 = get_contravariant_vector(1, contravariant_vectors, 1, j, element)

441


442
        for i in 2:(nnodes(dg) - 1)

443
            u_i = get_node_vars(u, equations, dg, i, j, element)

444
            Ja1_node_i = get_contravariant_vector(1, contravariant_vectors, i, j,

445
                                                  element)

446
            Ja1_avg = 0.5f0 * (Ja1_node_0 + Ja1_node_i)

447


448
            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

449
                phi_jump = volume_flux_noncons(u_0, u_i, Ja1_avg, equations,

450
                                               NonConservativeJump(), noncons)

451


452
                for v in eachvariable(equations)

453
                    # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5

454
                    fhat1_R[v, i, j] -= phi[v, noncons, i, j] * phi_jump[v]

455
                    fhat1_L[v, i + 1, j] -= phi[v, noncons, i, j] * phi_jump[v]

456
                end

457
            end

458
        end

459
        u_N = get_node_vars(u, equations, dg, nnodes(dg), j, element)

460
        Ja1_node_N = get_contravariant_vector(1, contravariant_vectors, nnodes(dg), j,

461
                                              element)

462
        Ja1_avg = 0.5f0 * (Ja1_node_0 + Ja1_node_N)

463


464
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

465
            phi_jump = volume_flux_noncons(u_0, u_N, Ja1_avg, equations,

466
                                           NonConservativeJump(), noncons)

467


468
            for v in eachvariable(equations)

469
                # The factor of 2 is missing because Trixi multiplies all the non-cons terms with 0.5

470
                fhat1_R[v, nnodes(dg), j] -= phi[v, noncons, nnodes(dg), j] *

471
                                             phi_jump[v]

472
            end

473
        end

474
    end

475


476
    # Split form volume flux in orientation 2: y direction

477
    flux_temp .= zero(eltype(flux_temp))

478
    flux_noncons_temp .= zero(eltype(flux_noncons_temp))

479


480
    for j in eachnode(dg), i in eachnode(dg)

481
        u_node = get_node_vars(u, equations, dg, i, j, element)

482


483
        # pull the contravariant vectors in each coordinate direction

484
        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)

485


486
        for jj in (j + 1):nnodes(dg)

487
            u_node_jj = get_node_vars(u, equations, dg, i, jj, element)

488
            # pull the contravariant vectors and compute the average

489
            Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj,

490
                                                   element)

491
            Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)

492
            # compute the contravariant sharp flux in the direction of the averaged contravariant vector

493
            fluxtilde2 = volume_flux_cons(u_node, u_node_jj, Ja2_avg, equations)

494
            multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2,

495
                                       equations, dg, i, j)

496
            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,

497
                                       equations, dg, i, jj)

498
            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

499
                # We multiply by 0.5 because that is done in other parts of Trixi

500
                flux2_noncons = volume_flux_noncons(u_node, u_node_jj, Ja2_avg,

501
                                                    equations,

502
                                                    NonConservativeJump(), noncons)

503
                multiply_add_to_node_vars!(flux_noncons_temp,

504
                                           0.5f0 * derivative_split[j, jj],

505
                                           flux2_noncons,

506
                                           equations, dg, noncons, i, j)

507
                multiply_add_to_node_vars!(flux_noncons_temp,

508
                                           -0.5f0 * derivative_split[jj, j],

509
                                           flux2_noncons,

510
                                           equations, dg, noncons, i, jj)

511
            end

512
        end

513
    end

514


515
    # FV-form flux `fhat` in y direction

516
    fhat_temp[:, :, 1] .= zero(eltype(fhat1_L))

517
    fhat_noncons_temp[:, :, :, 1] .= zero(eltype(fhat1_L))

518


519
    # Compute local contribution to non-conservative flux

520
    for j in eachnode(dg), i in eachnode(dg)

521
        u_local = get_node_vars(u, equations, dg, i, j, element)

522
        # pull the local contravariant vector

523
        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)

524
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

525
            set_node_vars!(phi,

526
                           volume_flux_noncons(u_local, Ja2_node, equations,

527
                                               NonConservativeLocal(), noncons),

528
                           equations, dg, noncons, i, j)

529
        end

530
    end

531


532
    for j in 1:(nnodes(dg) - 1), i in eachnode(dg)

533
        # Conservative part

534
        for v in eachvariable(equations)

535
            value = fhat_temp[v, i, j] + weights[j] * flux_temp[v, i, j]

536
            fhat_temp[v, i, j + 1] = value

537
            fhat2_L[v, i, j + 1] = value

538
            fhat2_R[v, i, j + 1] = value

539
        end

540
        # Nonconservative part

541
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),

542
            v in eachvariable(equations)

543


544
            value = fhat_noncons_temp[v, noncons, i, j] +

545
                    weights[j] * flux_noncons_temp[v, noncons, i, j]

546
            fhat_noncons_temp[v, noncons, i, j + 1] = value

547


548
            fhat2_L[v, i, j + 1] = fhat2_L[v, i, j + 1] + phi[v, noncons, i, j] * value

549
            fhat2_R[v, i, j + 1] = fhat2_R[v, i, j + 1] +

550
                                   phi[v, noncons, i, j + 1] * value

551
        end

552
    end

553


554
    # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes.

555
    for i in eachnode(dg)

556
        u_0 = get_node_vars(u, equations, dg, i, 1, element)

557
        Ja2_node_0 = get_contravariant_vector(2, contravariant_vectors, i, 1, element)

558


559
        for j in 2:(nnodes(dg) - 1)

560
            u_j = get_node_vars(u, equations, dg, i, j, element)

561
            Ja2_node_j = get_contravariant_vector(2, contravariant_vectors, i, j,

562
                                                  element)

563
            Ja2_avg = 0.5f0 * (Ja2_node_0 + Ja2_node_j)

564


565
            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

566
                phi_jump = volume_flux_noncons(u_0, u_j, Ja2_avg, equations,

567
                                               NonConservativeJump(), noncons)

568


569
                for v in eachvariable(equations)

570
                    # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5

571
                    fhat2_R[v, i, j] -= phi[v, noncons, i, j] * phi_jump[v]

572
                    fhat2_L[v, i, j + 1] -= phi[v, noncons, i, j] * phi_jump[v]

573
                end

574
            end

575
        end

576
        u_N = get_node_vars(u, equations, dg, i, nnodes(dg), element)

577
        Ja2_node_N = get_contravariant_vector(2, contravariant_vectors, i, nnodes(dg),

578
                                              element)

579
        Ja2_avg = 0.5f0 * (Ja2_node_0 + Ja2_node_N)

580


581
        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)

582
            phi_jump = volume_flux_noncons(u_0, u_N, Ja2_avg, equations,

583
                                           NonConservativeJump(), noncons)

584


585
            for v in eachvariable(equations)

586
                # The factor of 2 is missing cause Trixi multiplies all the non-cons terms with 0.5

587
                fhat2_R[v, i, nnodes(dg)] -= phi[v, noncons, i, nnodes(dg)] *

588
                                             phi_jump[v]

589
            end

590
        end

591
    end

592


593
    return nothing

594
end

595
end # @muladd

596


597