CoCalc -- test_performance_specializations

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/test/test_performance_specializations_3d.jl
⁵⁵⁸² views
1
module TestPerformanceSpecializations3D
2

3
using Test
4
using Trixi
5
using Trixi: @muladd
6

7
include("test_trixi.jl")
8

9
EXAMPLES_DIR = examples_dir()
10

11
# Start with a clean environment: remove Trixi.jl output directory if it exists
12
outdir = "out"
13
isdir(outdir) && rm(outdir, recursive = true)
14

15
@testset "Performance specializations 3D" begin
16
#! format: noindent
17

18
@timed_testset "TreeMesh3D, flux_shima_etal_turbo" begin
19
    trixi_include(@__MODULE__,
20
                  joinpath(EXAMPLES_DIR, "tree_3d_dgsem", "elixir_euler_ec.jl"),
21
                  initial_refinement_level = 0, tspan = (0.0, 0.0), polydeg = 3,
22
                  volume_flux = flux_shima_etal_turbo,
23
                  surface_flux = flux_shima_etal_turbo)
24
    u_ode = copy(sol.u[end])
25
    du_ode = zero(u_ode)
26

27
    # Preserve original memory since it will be `unsafe_wrap`ped and might
28
    # thus otherwise be garbage collected
29
    GC.@preserve u_ode du_ode begin
30
        u = Trixi.wrap_array(u_ode, semi)
31
        du = Trixi.wrap_array(du_ode, semi)
32
        have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)
33

34
        # Call the optimized default version
35
        du .= 0
36
        Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh),
37
                                        have_nonconservative_terms, semi.equations,
38
                                        semi.solver.volume_integral.volume_flux,
39
                                        semi.solver, semi.cache, true)
40
        du_specialized = du[:, :, :, :, 1]
41

42
        # Call the plain version - note the argument type `Function` of
43
        # `semi.solver.volume_integral.volume_flux`
44
        du .= 0
45
        invoke(Trixi.flux_differencing_kernel!,
46
               Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)},
47
                     typeof(have_nonconservative_terms), typeof(semi.equations),
48
                     Function, typeof(semi.solver), typeof(semi.cache), Bool},
49
               du, u, 1, typeof(semi.mesh),
50
               have_nonconservative_terms, semi.equations,
51
               semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
52
        du_baseline = du[:, :, :, :, 1]
53

54
        @test du_specialized ≈ du_baseline
55
    end
56
end
57

58
@timed_testset "TreeMesh3D, flux_ranocha_turbo" begin
59
    trixi_include(@__MODULE__,
60
                  joinpath(EXAMPLES_DIR, "tree_3d_dgsem", "elixir_euler_ec.jl"),
61
                  initial_refinement_level = 0, tspan = (0.0, 0.0), polydeg = 3,
62
                  volume_flux = flux_ranocha_turbo, surface_flux = flux_ranocha_turbo)
63
    u_ode = copy(sol.u[end])
64
    du_ode = zero(u_ode)
65

66
    # Preserve original memory since it will be `unsafe_wrap`ped and might
67
    # thus otherwise be garbage collected
68
    GC.@preserve u_ode du_ode begin
69
        u = Trixi.wrap_array(u_ode, semi)
70
        du = Trixi.wrap_array(du_ode, semi)
71
        have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)
72

73
        # Call the optimized default version
74
        du .= 0
75
        Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh),
76
                                        have_nonconservative_terms, semi.equations,
77
                                        semi.solver.volume_integral.volume_flux,
78
                                        semi.solver, semi.cache, true)
79
        du_specialized = du[:, :, :, :, 1]
80

81
        # Call the plain version - note the argument type `Function` of
82
        # `semi.solver.volume_integral.volume_flux`
83
        du .= 0
84
        invoke(Trixi.flux_differencing_kernel!,
85
               Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)},
86
                     typeof(have_nonconservative_terms), typeof(semi.equations),
87
                     Function, typeof(semi.solver), typeof(semi.cache), Bool},
88
               du, u, 1, typeof(semi.mesh),
89
               have_nonconservative_terms, semi.equations,
90
               semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
91
        du_baseline = du[:, :, :, :, 1]
92

93
        @test du_specialized ≈ du_baseline
94
    end
95
end
96

97
@timed_testset "StructuredMesh3D, flux_shima_etal_turbo" begin
98
    trixi_include(@__MODULE__,
99
                  joinpath(EXAMPLES_DIR, "structured_3d_dgsem", "elixir_euler_ec.jl"),
100
                  cells_per_dimension = (1, 1, 1), tspan = (0.0, 0.0), polydeg = 3,
101
                  volume_flux = flux_shima_etal_turbo,
102
                  surface_flux = flux_shima_etal_turbo)
103
    u_ode = copy(sol.u[end])
104
    du_ode = zero(u_ode)
105

106
    # Preserve original memory since it will be `unsafe_wrap`ped and might
107
    # thus otherwise be garbage collected
108
    GC.@preserve u_ode du_ode begin
109
        u = Trixi.wrap_array(u_ode, semi)
110
        du = Trixi.wrap_array(du_ode, semi)
111
        have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)
112

113
        # Call the optimized default version
114
        du .= 0
115
        Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh),
116
                                        have_nonconservative_terms, semi.equations,
117
                                        semi.solver.volume_integral.volume_flux,
118
                                        semi.solver, semi.cache, true)
119
        du_specialized = du[:, :, :, :, 1]
120

121
        # Call the plain version - note the argument type `Function` of
122
        # `semi.solver.volume_integral.volume_flux`
123
        du .= 0
124
        invoke(Trixi.flux_differencing_kernel!,
125
               Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)},
126
                     typeof(have_nonconservative_terms), typeof(semi.equations),
127
                     Function, typeof(semi.solver), typeof(semi.cache), Bool},
128
               du, u, 1, typeof(semi.mesh),
129
               have_nonconservative_terms, semi.equations,
130
               semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
131
        du_baseline = du[:, :, :, :, 1]
132

133
        @test du_specialized ≈ du_baseline
134
    end
135
end
136

137
@timed_testset "StructuredMesh3D, flux_ranocha_turbo" begin
138
    trixi_include(@__MODULE__,
139
                  joinpath(EXAMPLES_DIR, "structured_3d_dgsem", "elixir_euler_ec.jl"),
140
                  cells_per_dimension = (1, 1, 1), tspan = (0.0, 0.0), polydeg = 3,
141
                  volume_flux = flux_ranocha_turbo, surface_flux = flux_ranocha_turbo)
142
    u_ode = copy(sol.u[end])
143
    du_ode = zero(u_ode)
144

145
    # Preserve original memory since it will be `unsafe_wrap`ped and might
146
    # thus otherwise be garbage collected
147
    GC.@preserve u_ode du_ode begin
148
        u = Trixi.wrap_array(u_ode, semi)
149
        du = Trixi.wrap_array(du_ode, semi)
150
        have_nonconservative_terms = Trixi.have_nonconservative_terms(semi.equations)
151

152
        # Call the optimized default version
153
        du .= 0
154
        Trixi.flux_differencing_kernel!(du, u, 1, typeof(semi.mesh),
155
                                        have_nonconservative_terms, semi.equations,
156
                                        semi.solver.volume_integral.volume_flux,
157
                                        semi.solver, semi.cache, true)
158
        du_specialized = du[:, :, :, :, 1]
159

160
        # Call the plain version - note the argument type `Function` of
161
        # `semi.solver.volume_integral.volume_flux`
162
        du .= 0
163
        invoke(Trixi.flux_differencing_kernel!,
164
               Tuple{typeof(du), typeof(u), Integer, Type{typeof(semi.mesh)},
165
                     typeof(have_nonconservative_terms), typeof(semi.equations),
166
                     Function, typeof(semi.solver), typeof(semi.cache), Bool},
167
               du, u, 1, typeof(semi.mesh),
168
               have_nonconservative_terms, semi.equations,
169
               semi.solver.volume_integral.volume_flux, semi.solver, semi.cache, true)
170
        du_baseline = du[:, :, :, :, 1]
171

172
        @test du_specialized ≈ du_baseline
173
    end
174
end
175

176
@timed_testset "P4estMesh3D, combine_conservative_and_nonconservative_fluxes" begin
177
    trixi_include(@__MODULE__,
178
                  joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
179
                           "elixir_mhd_alfven_wave_nonperiodic.jl"),
180
                  volume_integral = VolumeIntegralFluxDifferencing((flux_hindenlang_gassner,
181
                                                                    flux_nonconservative_powell)),
182
                  tspan = (0.0, 0.1))
183
    u_ode = copy(sol.u[end])
184

185
    @muladd @inline function flux_hindenlang_gassner_nonconservative_powell(u_ll, u_rr,
186
                                                                            normal_direction::AbstractVector,
187
                                                                            equations::IdealGlmMhdEquations3D)
188
        # Unpack left and right states
189
        rho_ll, v1_ll, v2_ll, v3_ll, p_ll, B1_ll, B2_ll, B3_ll, psi_ll = cons2prim(u_ll,
190
                                                                                   equations)
191
        rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr, psi_rr = cons2prim(u_rr,
192
                                                                                   equations)
193
        v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
194
                     v3_ll * normal_direction[3]
195
        v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] +
196
                     v3_rr * normal_direction[3]
197
        B_dot_n_ll = B1_ll * normal_direction[1] + B2_ll * normal_direction[2] +
198
                     B3_ll * normal_direction[3]
199
        B_dot_n_rr = B1_rr * normal_direction[1] + B2_rr * normal_direction[2] +
200
                     B3_rr * normal_direction[3]
201

202
        # Compute the necessary mean values needed for either direction
203
        rho_mean = Trixi.ln_mean(rho_ll, rho_rr)
204
        # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
205
        # in exact arithmetic since
206
        #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
207
        #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
208
        inv_rho_p_mean = p_ll * p_rr * Trixi.inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
209
        v1_avg = 0.5f0 * (v1_ll + v1_rr)
210
        v2_avg = 0.5f0 * (v2_ll + v2_rr)
211
        v3_avg = 0.5f0 * (v3_ll + v3_rr)
212
        p_avg = 0.5f0 * (p_ll + p_rr)
213
        psi_avg = 0.5f0 * (psi_ll + psi_rr)
214
        velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr)
215
        magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr)
216

217
        # Calculate fluxes depending on normal_direction
218
        f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
219
        f2 = (f1 * v1_avg + (p_avg + magnetic_square_avg) * normal_direction[1]
220
              -
221
              0.5f0 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll))
222
        f3 = (f1 * v2_avg + (p_avg + magnetic_square_avg) * normal_direction[2]
223
              -
224
              0.5f0 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll))
225
        f4 = (f1 * v3_avg + (p_avg + magnetic_square_avg) * normal_direction[3]
226
              -
227
              0.5f0 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll))
228
        #f5 below
229
        f6 = (equations.c_h * psi_avg * normal_direction[1]
230
              +
231
              0.5f0 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll +
232
               v_dot_n_rr * B1_rr - v1_rr * B_dot_n_rr))
233
        f7 = (equations.c_h * psi_avg * normal_direction[2]
234
              +
235
              0.5f0 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll +
236
               v_dot_n_rr * B2_rr - v2_rr * B_dot_n_rr))
237
        f8 = (equations.c_h * psi_avg * normal_direction[3]
238
              +
239
              0.5f0 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll +
240
               v_dot_n_rr * B3_rr - v3_rr * B_dot_n_rr))
241
        f9 = equations.c_h * 0.5f0 * (B_dot_n_ll + B_dot_n_rr)
242
        # total energy flux is complicated and involves the previous components
243
        f5 = (f1 *
244
              (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one)
245
              +
246
              0.5f0 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll
247
               + (v_dot_n_ll * B1_ll * B1_rr + v_dot_n_rr * B1_rr * B1_ll)
248
               + (v_dot_n_ll * B2_ll * B2_rr + v_dot_n_rr * B2_rr * B2_ll)
249
               + (v_dot_n_ll * B3_ll * B3_rr + v_dot_n_rr * B3_rr * B3_ll)
250
               -
251
               (v1_ll * B_dot_n_ll * B1_rr + v1_rr * B_dot_n_rr * B1_ll)
252
               -
253
               (v2_ll * B_dot_n_ll * B2_rr + v2_rr * B_dot_n_rr * B2_ll)
254
               -
255
               (v3_ll * B_dot_n_ll * B3_rr + v3_rr * B_dot_n_rr * B3_ll)
256
               +
257
               equations.c_h * (B_dot_n_ll * psi_rr + B_dot_n_rr * psi_ll)))
258

259
        v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
260

261
        v_dot_B_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr
262
        f = SVector(f1, f2, f3, f4, f5, f6, f7, f8, f9)
263
        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
264
        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3})
265
        g_left = SVector(0,
266
                         B1_ll * B_dot_n_rr,
267
                         B2_ll * B_dot_n_rr,
268
                         B3_ll * B_dot_n_rr,
269
                         v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr,
270
                         v1_ll * B_dot_n_rr,
271
                         v2_ll * B_dot_n_rr,
272
                         v3_ll * B_dot_n_rr,
273
                         v_dot_n_ll * psi_rr)
274

275
        g_right = SVector(0,
276
                          B1_rr * B_dot_n_ll,
277
                          B2_rr * B_dot_n_ll,
278
                          B3_rr * B_dot_n_ll,
279
                          v_dot_B_rr * B_dot_n_ll + v_dot_n_rr * psi_rr * psi_ll,
280
                          v1_rr * B_dot_n_ll,
281
                          v2_rr * B_dot_n_ll,
282
                          v3_rr * B_dot_n_ll,
283
                          v_dot_n_rr * psi_ll)
284
        flux_left = f + 0.5f0 * g_left
285
        flux_right = f + 0.5f0 * g_right
286
        return flux_left, flux_right
287
    end
288

289
    @inline Trixi.combine_conservative_and_nonconservative_fluxes(::typeof(flux_hindenlang_gassner_nonconservative_powell),
290
    equations::IdealGlmMhdEquations3D) = Trixi.True()
291

292
    @muladd @inline function flux_hlle_nonconservative_powell(u_ll, u_rr,
293
                                                              normal_direction::AbstractVector,
294
                                                              equations::IdealGlmMhdEquations3D)
295
        f = flux_hlle(u_ll, u_rr, normal_direction, equations)
296

297
        rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
298
        rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
299

300
        v1_ll = rho_v1_ll / rho_ll
301
        v2_ll = rho_v2_ll / rho_ll
302
        v3_ll = rho_v3_ll / rho_ll
303
        v1_rr = rho_v1_rr / rho_rr
304
        v2_rr = rho_v2_rr / rho_rr
305
        v3_rr = rho_v3_rr / rho_rr
306
        v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
307
        v_dot_B_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr
308

309
        v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
310
                     v3_ll * normal_direction[3]
311
        B_dot_n_rr = B1_rr * normal_direction[1] +
312
                     B2_rr * normal_direction[2] +
313
                     B3_rr * normal_direction[3]
314
        v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] +
315
                     v3_rr * normal_direction[3]
316
        B_dot_n_ll = B1_ll * normal_direction[1] +
317
                     B2_ll * normal_direction[2] +
318
                     B3_ll * normal_direction[3]
319

320
        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
321
        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3})
322
        g_left = SVector(0,
323
                         B1_ll * B_dot_n_rr,
324
                         B2_ll * B_dot_n_rr,
325
                         B3_ll * B_dot_n_rr,
326
                         v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr,
327
                         v1_ll * B_dot_n_rr,
328
                         v2_ll * B_dot_n_rr,
329
                         v3_ll * B_dot_n_rr,
330
                         v_dot_n_ll * psi_rr)
331

332
        g_right = SVector(0,
333
                          B1_rr * B_dot_n_ll,
334
                          B2_rr * B_dot_n_ll,
335
                          B3_rr * B_dot_n_ll,
336
                          v_dot_B_rr * B_dot_n_ll + v_dot_n_rr * psi_rr * psi_ll,
337
                          v1_rr * B_dot_n_ll,
338
                          v2_rr * B_dot_n_ll,
339
                          v3_rr * B_dot_n_ll,
340
                          v_dot_n_rr * psi_ll)
341
        flux_left = f + 0.5f0 * g_left
342
        flux_right = f + 0.5f0 * g_right
343

344
        return flux_left, flux_right
345
    end
346

347
    @inline Trixi.combine_conservative_and_nonconservative_fluxes(::typeof(flux_hlle_nonconservative_powell),
348
    equations::IdealGlmMhdEquations3D) = Trixi.True()
349

350
    @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
351
                                                                      normal_direction::AbstractVector,
352
                                                                      x, t,
353
                                                                      surface_flux_function::typeof(flux_hlle_nonconservative_powell),
354
                                                                      equations)
355

356
        # get the external value of the solution
357
        u_boundary = boundary_condition.boundary_value_function(x, t, equations)
358

359
        # Calculate boundary flux
360
        flux, _ = surface_flux_function(u_inner, u_boundary, normal_direction,
361
                                        equations)
362
        return flux
363
    end
364

365
    trixi_include(@__MODULE__,
366
                  joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
367
                           "elixir_mhd_alfven_wave_nonperiodic.jl"),
368
                  surface_flux = flux_hlle_nonconservative_powell,
369
                  volume_integral = VolumeIntegralFluxDifferencing(flux_hindenlang_gassner_nonconservative_powell),
370
                  tspan = (0.0, 0.1))
371

372
    u_ode_specialized = copy(sol.u[end])
373
    @test u_ode_specialized ≈ u_ode
374
end
375
end
376

377
# Clean up afterwards: delete Trixi.jl output directory
378
@test_nowarn rm(outdir, recursive = true)
379

380
end #module
381

382
Product

Resources

Company
