CoCalc -- analysis

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/callbacks_step/analysis_dg3d.jl
²⁰⁵⁵ views
1
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
2
# Since these FMAs can increase the performance of many numerical algorithms,
3
# we need to opt-in explicitly.
4
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
5
@muladd begin
6
#! format: noindent
7

8
function create_cache_analysis(analyzer, mesh::TreeMesh{3},
9
                               equations, dg::DG, cache,
10
                               RealT, uEltype)
11

12
    # pre-allocate buffers
13
    # We use `StrideArray`s here since these buffers are used in performance-critical
14
    # places and the additional information passed to the compiler makes them faster
15
    # than native `Array`s.
16
    u_local = StrideArray(undef, uEltype,
17
                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
18
                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
19
    u_tmp1 = StrideArray(undef, uEltype,
20
                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
21
                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
22
    u_tmp2 = StrideArray(undef, uEltype,
23
                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
24
                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
25
    x_local = StrideArray(undef, RealT,
26
                          StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)),
27
                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
28
    x_tmp1 = StrideArray(undef, RealT,
29
                         StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)),
30
                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
31
    x_tmp2 = StrideArray(undef, RealT,
32
                         StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)),
33
                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
34

35
    return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2)
36
end
37

38
# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
39
function create_cache_analysis(analyzer,
40
                               mesh::P4estMesh{3, NDIMS_AMBIENT},
41
                               equations, dg::DG, cache,
42
                               RealT, uEltype) where {NDIMS_AMBIENT}
43

44
    # pre-allocate buffers
45
    # We use `StrideArray`s here since these buffers are used in performance-critical
46
    # places and the additional information passed to the compiler makes them faster
47
    # than native `Array`s.
48
    u_local = StrideArray(undef, uEltype,
49
                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
50
                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
51
    u_tmp1 = StrideArray(undef, uEltype,
52
                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
53
                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
54
    u_tmp2 = StrideArray(undef, uEltype,
55
                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
56
                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
57
    x_local = StrideArray(undef, RealT,
58
                          StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
59
                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
60
    x_tmp1 = StrideArray(undef, RealT,
61
                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
62
                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
63
    x_tmp2 = StrideArray(undef, RealT,
64
                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
65
                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
66
    jacobian_local = StrideArray(undef, RealT,
67
                                 StaticInt(nnodes(analyzer)),
68
                                 StaticInt(nnodes(analyzer)),
69
                                 StaticInt(nnodes(analyzer)))
70
    jacobian_tmp1 = StrideArray(undef, RealT,
71
                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)),
72
                                StaticInt(nnodes(dg)))
73
    jacobian_tmp2 = StrideArray(undef, RealT,
74
                                StaticInt(nnodes(analyzer)),
75
                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
76

77
    return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local,
78
            jacobian_tmp1, jacobian_tmp2)
79
end
80

81
function create_cache_analysis(analyzer,
82
                               mesh::Union{StructuredMesh{3}, T8codeMesh{3}},
83
                               equations, dg::DG, cache,
84
                               RealT, uEltype)
85

86
    # pre-allocate buffers
87
    # We use `StrideArray`s here since these buffers are used in performance-critical
88
    # places and the additional information passed to the compiler makes them faster
89
    # than native `Array`s.
90
    u_local = StrideArray(undef, uEltype,
91
                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
92
                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
93
    u_tmp1 = StrideArray(undef, uEltype,
94
                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
95
                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
96
    u_tmp2 = StrideArray(undef, uEltype,
97
                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
98
                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
99
    x_local = StrideArray(undef, RealT,
100
                          StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)),
101
                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
102
    x_tmp1 = StrideArray(undef, RealT,
103
                         StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)),
104
                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
105
    x_tmp2 = StrideArray(undef, RealT,
106
                         StaticInt(ndims(equations)), StaticInt(nnodes(analyzer)),
107
                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
108
    jacobian_local = StrideArray(undef, RealT,
109
                                 StaticInt(nnodes(analyzer)),
110
                                 StaticInt(nnodes(analyzer)),
111
                                 StaticInt(nnodes(analyzer)))
112
    jacobian_tmp1 = StrideArray(undef, RealT,
113
                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)),
114
                                StaticInt(nnodes(dg)))
115
    jacobian_tmp2 = StrideArray(undef, RealT,
116
                                StaticInt(nnodes(analyzer)),
117
                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
118

119
    return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local,
120
            jacobian_tmp1, jacobian_tmp2)
121
end
122

123
function calc_error_norms(func, u, t, analyzer,
124
                          mesh::TreeMesh{3}, equations, initial_condition,
125
                          dg::DGSEM, cache, cache_analysis)
126
    @unpack vandermonde, weights = analyzer
127
    @unpack node_coordinates = cache.elements
128
    @unpack u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2 = cache_analysis
129

130
    # Set up data structures
131
    l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1, 1), equations))
132
    linf_error = copy(l2_error)
133

134
    # Iterate over all elements for error calculations
135
    for element in eachelement(dg, cache)
136
        # Interpolate solution and node locations to analysis nodes
137
        multiply_dimensionwise!(u_local, vandermonde, view(u, :, :, :, :, element),
138
                                u_tmp1, u_tmp2)
139
        multiply_dimensionwise!(x_local, vandermonde,
140
                                view(node_coordinates, :, :, :, :, element), x_tmp1,
141
                                x_tmp2)
142

143
        # Calculate errors at each analysis node
144
        volume_jacobian_ = volume_jacobian(element, mesh, cache)
145

146
        for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer)
147
            u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j,
148
                                                        k), t, equations)
149
            diff = func(u_exact, equations) -
150
                   func(get_node_vars(u_local, equations, dg, i, j, k), equations)
151
            l2_error += diff .^ 2 *
152
                        (weights[i] * weights[j] * weights[k] * volume_jacobian_)
153
            linf_error = @. max(linf_error, abs(diff))
154
        end
155
    end
156

157
    # For L2 error, divide by total volume
158
    total_volume_ = total_volume(mesh)
159
    l2_error = @. sqrt(l2_error / total_volume_)
160

161
    return l2_error, linf_error
162
end
163

164
function calc_error_norms(func, u, t, analyzer,
165
                          mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
166
                          equations, initial_condition,
167
                          dg::DGSEM, cache, cache_analysis)
168
    @unpack vandermonde, weights = analyzer
169
    @unpack node_coordinates, inverse_jacobian = cache.elements
170
    @unpack u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local, jacobian_tmp1, jacobian_tmp2 = cache_analysis
171

172
    # Set up data structures
173
    l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1, 1), equations))
174
    linf_error = copy(l2_error)
175
    total_volume = zero(real(mesh))
176

177
    # Iterate over all elements for error calculations
178
    for element in eachelement(dg, cache)
179
        # Interpolate solution and node locations to analysis nodes
180
        multiply_dimensionwise!(u_local, vandermonde, view(u, :, :, :, :, element),
181
                                u_tmp1, u_tmp2)
182
        multiply_dimensionwise!(x_local, vandermonde,
183
                                view(node_coordinates, :, :, :, :, element), x_tmp1,
184
                                x_tmp2)
185
        multiply_scalar_dimensionwise!(jacobian_local, vandermonde,
186
                                       inv.(view(inverse_jacobian, :, :, :, element)),
187
                                       jacobian_tmp1, jacobian_tmp2)
188

189
        # Calculate errors at each analysis node
190
        for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer)
191
            u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j,
192
                                                        k), t, equations)
193
            diff = func(u_exact, equations) -
194
                   func(get_node_vars(u_local, equations, dg, i, j, k), equations)
195
            # We take absolute value as we need the Jacobian here for the volume calculation
196
            abs_jacobian_local_ijk = abs(jacobian_local[i, j, k])
197

198
            l2_error += diff .^ 2 *
199
                        (weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk)
200
            linf_error = @. max(linf_error, abs(diff))
201
            total_volume += (weights[i] * weights[j] * weights[k] *
202
                             abs_jacobian_local_ijk)
203
        end
204
    end
205

206
    # For L2 error, divide by total volume
207
    l2_error = @. sqrt(l2_error / total_volume)
208

209
    return l2_error, linf_error
210
end
211

212
function integrate_via_indices(func::Func, u,
213
                               mesh::TreeMesh{3}, equations, dg::DGSEM, cache,
214
                               args...; normalize = true) where {Func}
215
    @unpack weights = dg.basis
216

217
    # Initialize integral with zeros of the right shape
218
    integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))
219

220
    # Use quadrature to numerically integrate over entire domain
221
    @batch reduction=(+, integral) for element in eachelement(dg, cache)
222
        volume_jacobian_ = volume_jacobian(element, mesh, cache)
223
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
224
            integral += volume_jacobian_ * weights[i] * weights[j] * weights[k] *
225
                        func(u, i, j, k, element, equations, dg, args...)
226
        end
227
    end
228

229
    # Normalize with total volume
230
    if normalize
231
        integral = integral / total_volume(mesh)
232
    end
233

234
    return integral
235
end
236

237
function integrate_via_indices(func::Func, u,
238
                               mesh::Union{StructuredMesh{3}, P4estMesh{3},
239
                                           T8codeMesh{3}},
240
                               equations, dg::DGSEM, cache,
241
                               args...; normalize = true) where {Func}
242
    @unpack weights = dg.basis
243

244
    # Initialize integral with zeros of the right shape
245
    integral = zero(func(u, 1, 1, 1, 1, equations, dg, args...))
246
    total_volume = zero(real(mesh))
247

248
    # Use quadrature to numerically integrate over entire domain
249
    @batch reduction=((+, integral), (+, total_volume)) for element in eachelement(dg,
250
                                                                                   cache)
251
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
252
            volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, k, element]))
253
            integral += volume_jacobian * weights[i] * weights[j] * weights[k] *
254
                        func(u, i, j, k, element, equations, dg, args...)
255
            total_volume += volume_jacobian * weights[i] * weights[j] * weights[k]
256
        end
257
    end
258

259
    # Normalize with total volume
260
    if normalize
261
        integral = integral / total_volume
262
    end
263

264
    return integral
265
end
266

267
function integrate(func::Func, u,
268
                   mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3},
269
                               T8codeMesh{3}},
270
                   equations, dg::DG, cache; normalize = true) where {Func}
271
    integrate_via_indices(u, mesh, equations, dg, cache;
272
                          normalize = normalize) do u, i, j, k, element, equations, dg
273
        u_local = get_node_vars(u, equations, dg, i, j, k, element)
274
        return func(u_local, equations)
275
    end
276
end
277

278
function integrate(func::Func, u,
279
                   mesh::Union{TreeMesh{3}, P4estMesh{3}},
280
                   equations, equations_parabolic,
281
                   dg::DGSEM,
282
                   cache, cache_parabolic; normalize = true) where {Func}
283
    gradients_x, gradients_y, gradients_z = cache_parabolic.viscous_container.gradients
284
    integrate_via_indices(u, mesh, equations, dg, cache;
285
                          normalize = normalize) do u, i, j, k, element, equations, dg
286
        u_local = get_node_vars(u, equations, dg, i, j, k, element)
287
        gradients_1_local = get_node_vars(gradients_x, equations_parabolic, dg, i, j, k,
288
                                          element)
289
        gradients_2_local = get_node_vars(gradients_y, equations_parabolic, dg, i, j, k,
290
                                          element)
291
        gradients_3_local = get_node_vars(gradients_z, equations_parabolic, dg, i, j, k,
292
                                          element)
293
        return func(u_local, (gradients_1_local, gradients_2_local, gradients_3_local),
294
                    equations_parabolic)
295
    end
296
end
297

298
function analyze(::typeof(entropy_timederivative), du, u, t,
299
                 mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3},
300
                             T8codeMesh{3}},
301
                 equations, dg::DG, cache)
302
    # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ
303
    integrate_via_indices(u, mesh, equations, dg, cache,
304
                          du) do u, i, j, k, element, equations, dg, du
305
        u_node = get_node_vars(u, equations, dg, i, j, k, element)
306
        du_node = get_node_vars(du, equations, dg, i, j, k, element)
307
        dot(cons2entropy(u_node, equations), du_node)
308
    end
309
end
310

311
function analyze(::Val{:l2_divb}, du, u, t,
312
                 mesh::TreeMesh{3}, equations,
313
                 dg::DGSEM, cache)
314
    integrate_via_indices(u, mesh, equations, dg, cache, cache,
315
                          dg.basis.derivative_matrix) do u, i, j, k, element, equations,
316
                                                         dg, cache, derivative_matrix
317
        divb = zero(eltype(u))
318
        for l in eachnode(dg)
319
            u_ljk = get_node_vars(u, equations, dg, l, j, k, element)
320
            u_ilk = get_node_vars(u, equations, dg, i, l, k, element)
321
            u_ijl = get_node_vars(u, equations, dg, i, j, l, element)
322

323
            B_ljk = magnetic_field(u_ljk, equations)
324
            B_ilk = magnetic_field(u_ilk, equations)
325
            B_ijl = magnetic_field(u_ijl, equations)
326

327
            divb += (derivative_matrix[i, l] * B_ljk[1] +
328
                     derivative_matrix[j, l] * B_ilk[2] +
329
                     derivative_matrix[k, l] * B_ijl[3])
330
        end
331
        divb *= cache.elements.inverse_jacobian[element]
332
        divb^2
333
    end |> sqrt
334
end
335

336
function analyze(::Val{:l2_divb}, du, u, t,
337
                 mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
338
                 equations,
339
                 dg::DGSEM, cache)
340
    @unpack contravariant_vectors = cache.elements
341
    integrate_via_indices(u, mesh, equations, dg, cache, cache,
342
                          dg.basis.derivative_matrix) do u, i, j, k, element, equations,
343
                                                         dg, cache, derivative_matrix
344
        divb = zero(eltype(u))
345
        # Get the contravariant vectors Ja^1, Ja^2, and Ja^3
346
        Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors,
347
                                                    i, j, k, element)
348
        Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors,
349
                                                    i, j, k, element)
350
        Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors,
351
                                                    i, j, k, element)
352
        # Compute the transformed divergence
353
        for l in eachnode(dg)
354
            u_ljk = get_node_vars(u, equations, dg, l, j, k, element)
355
            u_ilk = get_node_vars(u, equations, dg, i, l, k, element)
356
            u_ijl = get_node_vars(u, equations, dg, i, j, l, element)
357

358
            B_ljk = magnetic_field(u_ljk, equations)
359
            B_ilk = magnetic_field(u_ilk, equations)
360
            B_ijl = magnetic_field(u_ijl, equations)
361

362
            divb += (derivative_matrix[i, l] *
363
                     (Ja11 * B_ljk[1] + Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
364
                     derivative_matrix[j, l] *
365
                     (Ja21 * B_ilk[1] + Ja22 * B_ilk[2] + Ja23 * B_ilk[3]) +
366
                     derivative_matrix[k, l] *
367
                     (Ja31 * B_ijl[1] + Ja32 * B_ijl[2] + Ja33 * B_ijl[3]))
368
        end
369
        divb *= cache.elements.inverse_jacobian[i, j, k, element]
370
        divb^2
371
    end |> sqrt
372
end
373

374
function analyze(::Val{:linf_divb}, du, u, t,
375
                 mesh::TreeMesh{3}, equations,
376
                 dg::DGSEM, cache)
377
    @unpack derivative_matrix, weights = dg.basis
378

379
    # integrate over all elements to get the divergence-free condition errors
380
    linf_divb = zero(eltype(u))
381
    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
382
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
383
            divb = zero(eltype(u))
384
            for l in eachnode(dg)
385
                u_ljk = get_node_vars(u, equations, dg, l, j, k, element)
386
                u_ilk = get_node_vars(u, equations, dg, i, l, k, element)
387
                u_ijl = get_node_vars(u, equations, dg, i, j, l, element)
388

389
                B_ljk = magnetic_field(u_ljk, equations)
390
                B_ilk = magnetic_field(u_ilk, equations)
391
                B_ijl = magnetic_field(u_ijl, equations)
392

393
                divb += (derivative_matrix[i, l] * B_ljk[1] +
394
                         derivative_matrix[j, l] * B_ilk[2] +
395
                         derivative_matrix[k, l] * B_ijl[3])
396
            end
397
            divb *= cache.elements.inverse_jacobian[element]
398
            linf_divb = max(linf_divb, abs(divb))
399
        end
400
    end
401

402
    if mpi_isparallel()
403
        # Base.max instead of max needed, see comment in src/auxiliary/math.jl
404
        linf_divb = MPI.Allreduce!(Ref(linf_divb), Base.max, mpi_comm())[]
405
    end
406

407
    return linf_divb
408
end
409

410
function analyze(::Val{:linf_divb}, du, u, t,
411
                 mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
412
                 equations,
413
                 dg::DGSEM, cache)
414
    @unpack derivative_matrix, weights = dg.basis
415
    @unpack contravariant_vectors = cache.elements
416

417
    # integrate over all elements to get the divergence-free condition errors
418
    linf_divb = zero(eltype(u))
419
    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
420
        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
421
            divb = zero(eltype(u))
422
            # Get the contravariant vectors Ja^1, Ja^2, and Ja^3
423
            Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors,
424
                                                        i, j, k, element)
425
            Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors,
426
                                                        i, j, k, element)
427
            Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors,
428
                                                        i, j, k, element)
429
            # Compute the transformed divergence
430
            for l in eachnode(dg)
431
                u_ljk = get_node_vars(u, equations, dg, l, j, k, element)
432
                u_ilk = get_node_vars(u, equations, dg, i, l, k, element)
433
                u_ijl = get_node_vars(u, equations, dg, i, j, l, element)
434

435
                B_ljk = magnetic_field(u_ljk, equations)
436
                B_ilk = magnetic_field(u_ilk, equations)
437
                B_ijl = magnetic_field(u_ijl, equations)
438

439
                divb += (derivative_matrix[i, l] * (Ja11 * B_ljk[1] +
440
                          Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
441
                         derivative_matrix[j, l] * (Ja21 * B_ilk[1] +
442
                          Ja22 * B_ilk[2] + Ja23 * B_ilk[3]) +
443
                         derivative_matrix[k, l] * (Ja31 * B_ijl[1] +
444
                          Ja32 * B_ijl[2] + Ja33 * B_ijl[3]))
445
            end
446
            divb *= cache.elements.inverse_jacobian[i, j, k, element]
447
            linf_divb = max(linf_divb, abs(divb))
448
        end
449
    end
450

451
    if mpi_isparallel()
452
        # Base.max instead of max needed, see comment in src/auxiliary/math.jl
453
        linf_divb = MPI.Allreduce!(Ref(linf_divb), Base.max, mpi_comm())[]
454
    end
455

456
    return linf_divb
457
end
458
end # @muladd
459

460
Product

Resources

Company
