CoCalc -- test_performance_specializations

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/test/test_performance_specializations_2d.jl
⁵⁵⁸² views
1
module TestPerformanceSpecializations2D
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 2D" begin
16
#! format: noindent
17

18
@timed_testset "TreeMesh2D, flux_shima_etal_turbo" begin
19
    trixi_include(@__MODULE__,
20
                  joinpath(EXAMPLES_DIR, "tree_2d_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 "TreeMesh2D, flux_ranocha_turbo" begin
59
    trixi_include(@__MODULE__,
60
                  joinpath(EXAMPLES_DIR, "tree_2d_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 "StructuredMesh2D, flux_shima_etal_turbo" begin
98
    trixi_include(@__MODULE__,
99
                  joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_euler_ec.jl"),
100
                  cells_per_dimension = (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 "StructuredMesh2D, flux_ranocha_turbo" begin
138
    trixi_include(@__MODULE__,
139
                  joinpath(EXAMPLES_DIR, "structured_2d_dgsem", "elixir_euler_ec.jl"),
140
                  cells_per_dimension = (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 "P4estMesh2D, combine_conservative_and_nonconservative_fluxes" begin
177
    trixi_include(@__MODULE__,
178
                  joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_mhd_rotor.jl"),
179
                  volume_integral = VolumeIntegralFluxDifferencing((flux_hindenlang_gassner,
180
                                                                    flux_nonconservative_powell)),
181
                  amr_callback = nothing, tspan = (0.0, 0.001))
182
    u_ode = copy(sol.u[end])
183

184
    @muladd @inline function flux_hindenlang_gassner_nonconservative_powell(u_ll, u_rr,
185
                                                                            normal_direction,
186
                                                                            equations)
187
        rho_ll, v1_ll, v2_ll, v3_ll, p_ll, B1_ll, B2_ll, B3_ll, psi_ll = cons2prim(u_ll,
188
                                                                                   equations)
189
        rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr, psi_rr = cons2prim(u_rr,
190
                                                                                   equations)
191
        v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
192
        v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
193
        B_dot_n_ll = B1_ll * normal_direction[1] + B2_ll * normal_direction[2]
194
        B_dot_n_rr = B1_rr * normal_direction[1] + B2_rr * normal_direction[2]
195

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

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

251
        v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
252
        v_dot_B_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr
253

254
        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
255
        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
256
        g_left = SVector(0,
257
                         B1_ll * B_dot_n_rr,
258
                         B2_ll * B_dot_n_rr,
259
                         B3_ll * B_dot_n_rr,
260
                         v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr,
261
                         v1_ll * B_dot_n_rr,
262
                         v2_ll * B_dot_n_rr,
263
                         v3_ll * B_dot_n_rr,
264
                         v_dot_n_ll * psi_rr)
265
        g_right = SVector(0,
266
                          B1_rr * B_dot_n_ll,
267
                          B2_rr * B_dot_n_ll,
268
                          B3_rr * B_dot_n_ll,
269
                          v_dot_B_rr * B_dot_n_ll + v_dot_n_rr * psi_rr * psi_ll,
270
                          v1_rr * B_dot_n_ll,
271
                          v2_rr * B_dot_n_ll,
272
                          v3_rr * B_dot_n_ll,
273
                          v_dot_n_rr * psi_ll)
274

275
        f = SVector(f1, f2, f3, f4, f5, f6, f7, f8, f9)
276
        flux_left = f + g_left * 0.5f0
277
        flux_right = f + g_right * 0.5f0
278
        return flux_left, flux_right
279
    end
280

281
    @inline Trixi.combine_conservative_and_nonconservative_fluxes(::typeof(flux_hindenlang_gassner_nonconservative_powell),
282
    equations::IdealGlmMhdEquations2D) = Trixi.True()
283

284
    @muladd @inline function flux_lax_friedrichs_nonconservative_powell(u_ll, u_rr,
285
                                                                        normal_direction,
286
                                                                        equations)
287
        f = FluxLaxFriedrichs(max_abs_speed_naive)(u_ll, u_rr, normal_direction,
288
                                                   equations)
289

290
        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
291
        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
292

293
        inv_rho_ll = 1 / rho_ll
294
        inv_rho_rr = 1 / rho_rr
295

296
        v1_ll = rho_v1_ll * inv_rho_ll
297
        v2_ll = rho_v2_ll * inv_rho_ll
298
        v3_ll = rho_v3_ll * inv_rho_ll
299
        v1_rr = rho_v1_rr * inv_rho_rr
300
        v2_rr = rho_v2_rr * inv_rho_rr
301
        v3_rr = rho_v3_rr * inv_rho_rr
302
        v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
303
        v_dot_B_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr
304

305
        v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
306
        B_dot_n_rr = B1_rr * normal_direction[1] +
307
                     B2_rr * normal_direction[2]
308

309
        B_dot_n_ll = B1_ll * normal_direction[1] +
310
                     B2_ll * normal_direction[2]
311
        v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
312
        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
313
        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
314
        g_left = SVector(0,
315
                         B1_ll * B_dot_n_rr,
316
                         B2_ll * B_dot_n_rr,
317
                         B3_ll * B_dot_n_rr,
318
                         v_dot_B_ll * B_dot_n_rr + v_dot_n_ll * psi_ll * psi_rr,
319
                         v1_ll * B_dot_n_rr,
320
                         v2_ll * B_dot_n_rr,
321
                         v3_ll * B_dot_n_rr,
322
                         v_dot_n_ll * psi_rr)
323

324
        g_right = SVector(0,
325
                          B1_rr * B_dot_n_ll,
326
                          B2_rr * B_dot_n_ll,
327
                          B3_rr * B_dot_n_ll,
328
                          v_dot_B_rr * B_dot_n_ll + v_dot_n_rr * psi_rr * psi_ll,
329
                          v1_rr * B_dot_n_ll,
330
                          v2_rr * B_dot_n_ll,
331
                          v3_rr * B_dot_n_ll,
332
                          v_dot_n_rr * psi_ll)
333
        flux_left = f + 0.5f0 * g_left
334
        flux_right = f + 0.5f0 * g_right
335
        return flux_left, flux_right
336
    end
337

338
    @inline Trixi.combine_conservative_and_nonconservative_fluxes(::typeof(flux_lax_friedrichs_nonconservative_powell),
339
    equations::IdealGlmMhdEquations2D) = Trixi.True()
340

341
    @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
342
                                                                      normal_direction::AbstractVector,
343
                                                                      x, t,
344
                                                                      surface_flux_function::typeof(flux_lax_friedrichs_nonconservative_powell),
345
                                                                      equations)
346

347
        # get the external value of the solution
348
        u_boundary = boundary_condition.boundary_value_function(x, t, equations)
349

350
        # Calculate boundary flux
351
        flux, _ = surface_flux_function(u_inner, u_boundary, normal_direction,
352
                                        equations)
353
        return flux
354
    end
355

356
    trixi_include(@__MODULE__,
357
                  joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_mhd_rotor.jl"),
358
                  surface_flux = flux_lax_friedrichs_nonconservative_powell,
359
                  volume_integral = VolumeIntegralFluxDifferencing(flux_hindenlang_gassner_nonconservative_powell),
360
                  amr_callback = nothing, tspan = (0.0, 0.001))
361

362
    u_ode_specialized = copy(sol.u[end])
363
    @test u_ode_specialized ≈ u_ode
364
end
365
end
366

367
# Clean up afterwards: delete Trixi.jl output directory
368
@test_nowarn rm(outdir, recursive = true)
369

370
end #module
371

372
Product

Resources

Company
