From cd3485bdefaa29a6c664e581378af287ecc7ee69 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 16 Jan 2025 14:55:34 +0100 Subject: [PATCH 1/5] `resize!` only when needed `EulerGravity` (#2224) Co-authored-by: Michael Schlottke-Lakemper --- .../semidiscretization_euler_gravity.jl | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index e268497a9b..6e4f05d78a 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -284,11 +284,6 @@ end function update_gravity!(semi::SemidiscretizationEulerGravity, u_ode) @unpack semi_euler, semi_gravity, parameters, gravity_counter, cache = semi - # Can be changed by AMR - resize!(cache.du_ode, length(cache.u_ode)) - resize!(cache.u_tmp1_ode, length(cache.u_ode)) - resize!(cache.u_tmp2_ode, length(cache.u_ode)) - u_euler = wrap_array(u_ode, semi_euler) u_gravity = wrap_array(cache.u_ode, semi_gravity) du_gravity = wrap_array(cache.du_ode, semi_gravity) @@ -568,7 +563,17 @@ end t, iter; kwargs...) passive_args = ((semi.cache.u_ode, mesh_equations_solver_cache(semi.semi_gravity)...),) - amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)..., semi, t, iter; - kwargs..., passive_args = passive_args) + has_changed = amr_callback(u_ode, mesh_equations_solver_cache(semi.semi_euler)..., + semi, t, iter; + kwargs..., passive_args = passive_args) + + if has_changed + new_length = length(semi.cache.u_ode) + resize!(semi.cache.du_ode, new_length) + resize!(semi.cache.u_tmp1_ode, new_length) + resize!(semi.cache.u_tmp2_ode, new_length) + end + + return has_changed end end # @muladd From 9893b736a88be6b95c82a04b118fc03c228759a3 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 16 Jan 2025 16:16:57 +0100 Subject: [PATCH 2/5] Tidy up comments LBM (#2231) --- src/equations/lattice_boltzmann_2d.jl | 8 -------- src/equations/lattice_boltzmann_3d.jl | 8 -------- src/solvers/dgsem_tree/dg_2d.jl | 2 +- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 2 +- 4 files changed, 2 insertions(+), 18 deletions(-) diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl index b2323789dd..45184b7a70 100644 --- a/src/equations/lattice_boltzmann_2d.jl +++ b/src/equations/lattice_boltzmann_2d.jl @@ -235,9 +235,6 @@ No-slip wall boundary condition using the bounce-back approach. return flux end -# Pre-defined source terms should be implemented as -# function source_terms_WHATEVER(u, x, t, equations::LatticeBoltzmannEquations2D) - # Calculate 1D flux in for a single point @inline function flux(u, orientation::Integer, equations::LatticeBoltzmannEquations2D) if orientation == 1 @@ -248,11 +245,6 @@ end return v_alpha .* u end -# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation -# @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::LatticeBoltzmannEquations2D) -# λ_max = -# end - """ flux_godunov(u_ll, u_rr, orientation, equations::LatticeBoltzmannEquations2D) diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl index 2211ac5bb7..017e0f07fc 100644 --- a/src/equations/lattice_boltzmann_3d.jl +++ b/src/equations/lattice_boltzmann_3d.jl @@ -223,9 +223,6 @@ function initial_condition_constant(x, t, equations::LatticeBoltzmannEquations3D return equilibrium_distribution(rho, v1, v2, v3, equations) end -# Pre-defined source terms should be implemented as -# function source_terms_WHATEVER(u, x, t, equations::LatticeBoltzmannEquations3D) - # Calculate 1D flux in for a single point @inline function flux(u, orientation::Integer, equations::LatticeBoltzmannEquations3D) if orientation == 1 # x-direction @@ -238,11 +235,6 @@ end return v_alpha .* u end -# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation -# @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::LatticeBoltzmannEquations3D) -# λ_max = -# end - """ flux_godunov(u_ll, u_rr, orientation, equations::LatticeBoltzmannEquations3D) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 98bd088920..118b1f75bb 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -1118,7 +1118,7 @@ end # @views mul!(surface_flux_values[v, :, direction, large_element], # mortar_l2.reverse_upper, fstar_upper[v, :]) # @views mul!(surface_flux_values[v, :, direction, large_element], - # mortar_l2.reverse_lower, fstar_lower[v, :], true, true) + # mortar_l2.reverse_lower, fstar_lower[v, :], true, true) # end # The code above could be replaced by the following code. However, the relative efficiency # depends on the types of fstar_upper/fstar_lower and dg.l2mortar_reverse_upper. diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 2123eae545..ce14da8d20 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -725,7 +725,7 @@ end # @views mul!(surface_flux_values[v, :, direction, large_element], # mortar_l2.reverse_upper, fstar_upper[v, :]) # @views mul!(surface_flux_values[v, :, direction, large_element], - # mortar_l2.reverse_lower, fstar_lower[v, :], true, true) + # mortar_l2.reverse_lower, fstar_lower[v, :], true, true) # end # The code above could be replaced by the following code. However, the relative efficiency # depends on the types of fstar_upper/fstar_lower and dg.l2mortar_reverse_upper. From fdd22a8005177bc626f07ebb130709028ccb80e0 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Fri, 17 Jan 2025 09:41:10 -1000 Subject: [PATCH 3/5] Update comments for some `DGSEM` APIs (#2238) * Update * Add --- src/solvers/dgsem/dgsem.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl index 27caad4d2d..41a9722ecf 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/src/solvers/dgsem/dgsem.jl @@ -22,7 +22,7 @@ Create a discontinuous Galerkin spectral element method (DGSEM) using a """ const DGSEM = DG{Basis} where {Basis <: LobattoLegendreBasis} -# TODO: Deprecated in v0.3 (no longer documented) +# This API is no longer documented, and we recommend avoiding its public use. function DGSEM(basis::LobattoLegendreBasis, surface_flux = flux_central, volume_integral = VolumeIntegralWeakForm(), @@ -32,7 +32,7 @@ function DGSEM(basis::LobattoLegendreBasis, typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral) end -# TODO: Deprecated in v0.3 (no longer documented) +# This API is no longer documented, and we recommend avoiding its public use. function DGSEM(basis::LobattoLegendreBasis, surface_integral::AbstractSurfaceIntegral, volume_integral = VolumeIntegralWeakForm(), @@ -41,7 +41,7 @@ function DGSEM(basis::LobattoLegendreBasis, typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral) end -# TODO: Deprecated in v0.3 (no longer documented) +# This API is no longer documented, and we recommend avoiding its public use. function DGSEM(RealT, polydeg::Integer, surface_flux = flux_central, volume_integral = VolumeIntegralWeakForm(), @@ -51,6 +51,7 @@ function DGSEM(RealT, polydeg::Integer, return DGSEM(basis, surface_flux, volume_integral, mortar) end +# This API is no longer documented, and we recommend avoiding its public use. function DGSEM(polydeg, surface_flux = flux_central, volume_integral = VolumeIntegralWeakForm()) DGSEM(Float64, polydeg, surface_flux, volume_integral) From 6349413023e2f8fdb01fd98323c804b7058b2336 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Sat, 18 Jan 2025 07:22:35 -1000 Subject: [PATCH 4/5] Fix redundant type conversion in `LobattoLegendreBasis` (#2239) * Fix * Fix * Fix * Update src/solvers/dgsem/basis_lobatto_legendre.jl Co-authored-by: Hendrik Ranocha * Update src/solvers/dgsem/basis_lobatto_legendre.jl Co-authored-by: Hendrik Ranocha * Update src/solvers/dgsem/basis_lobatto_legendre.jl Co-authored-by: Hendrik Ranocha * Format --------- Co-authored-by: Hendrik Ranocha --- src/solvers/dgsem/basis_lobatto_legendre.jl | 32 +++++++++------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index 3ea668c026..6c8bec9c93 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -40,31 +40,27 @@ function LobattoLegendreBasis(RealT, polydeg::Integer) nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT) inverse_weights_ = inv.(weights_) - _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_, RealT) + _, inverse_vandermonde_legendre = vandermonde_legendre(nodes_, RealT) - boundary_interpolation_ = zeros(RealT, nnodes_, 2) - boundary_interpolation_[:, 1] = calc_lhat(-one(RealT), nodes_, weights_) - boundary_interpolation_[:, 2] = calc_lhat(one(RealT), nodes_, weights_) + boundary_interpolation = zeros(RealT, nnodes_, 2) + boundary_interpolation[:, 1] = calc_lhat(-one(RealT), nodes_, weights_) + boundary_interpolation[:, 2] = calc_lhat(one(RealT), nodes_, weights_) - derivative_matrix_ = polynomial_derivative_matrix(nodes_) - derivative_split_ = calc_dsplit(nodes_, weights_) - derivative_split_transpose_ = Matrix(derivative_split_') - derivative_dhat_ = calc_dhat(nodes_, weights_) + derivative_matrix = polynomial_derivative_matrix(nodes_) + derivative_split = calc_dsplit(nodes_, weights_) + derivative_split_transpose = Matrix(derivative_split') + derivative_dhat = calc_dhat(nodes_, weights_) - # type conversions to get the requested real type and enable possible - # optimizations of runtime performance and latency + # Type conversions to enable possible optimizations of runtime performance + # and latency nodes = SVector{nnodes_, RealT}(nodes_) weights = SVector{nnodes_, RealT}(weights_) inverse_weights = SVector{nnodes_, RealT}(inverse_weights_) - inverse_vandermonde_legendre = convert.(RealT, inverse_vandermonde_legendre_) - boundary_interpolation = convert.(RealT, boundary_interpolation_) - - # Usually as fast as `SMatrix` (when using `let` in the volume integral/`@threaded`) - derivative_matrix = Matrix{RealT}(derivative_matrix_) - derivative_split = Matrix{RealT}(derivative_split_) - derivative_split_transpose = Matrix{RealT}(derivative_split_transpose_) - derivative_dhat = Matrix{RealT}(derivative_dhat_) + # We keep the matrices above stored using the standard `Matrix` type + # since this is usually as fast as `SMatrix` + # (when using `let` in the volume integral/`@threaded`) + # and reduces latency return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes), typeof(inverse_vandermonde_legendre), From 614910491810aa9481165f43084e4aad3d7da39e Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Sun, 19 Jan 2025 21:14:22 -1000 Subject: [PATCH 5/5] Fix type conversion and type instability in `MortarL2` (#2242) * Fix * Fix --- src/solvers/dgsem/basis_lobatto_legendre.jl | 20 ++++++++------------ src/solvers/dgsem/l2projection.jl | 20 ++++++++++---------- 2 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index 6c8bec9c93..048831df55 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -159,17 +159,15 @@ function MortarL2(basis::LobattoLegendreBasis) RealT = real(basis) nnodes_ = nnodes(basis) - forward_upper_ = calc_forward_upper(nnodes_, RealT) - forward_lower_ = calc_forward_lower(nnodes_, RealT) - reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT) - reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT) + forward_upper = calc_forward_upper(nnodes_, RealT) + forward_lower = calc_forward_lower(nnodes_, RealT) + reverse_upper = calc_reverse_upper(nnodes_, Val(:gauss), RealT) + reverse_lower = calc_reverse_lower(nnodes_, Val(:gauss), RealT) - # type conversions to get the requested real type and enable possible - # optimizations of runtime performance and latency - - # Usually as fast as `SMatrix` but better for latency - forward_upper = Matrix{RealT}(forward_upper_) - forward_lower = Matrix{RealT}(forward_lower_) + # We keep the matrices above stored using the standard `Matrix` type + # since this is usually as fast as `SMatrix` + # (when using `let` in the volume integral/`@threaded`) + # and reduces latency # TODO: Taal performance # Check the performance of different implementations of `mortar_fluxes_to_elements!` @@ -179,8 +177,6 @@ function MortarL2(basis::LobattoLegendreBasis) # `@tullio` when the matrix sizes are not necessarily static. # reverse_upper = SMatrix{nnodes_, nnodes_, RealT, nnodes_^2}(reverse_upper_) # reverse_lower = SMatrix{nnodes_, nnodes_, RealT, nnodes_^2}(reverse_lower_) - reverse_upper = Matrix{RealT}(reverse_upper_) - reverse_lower = Matrix{RealT}(reverse_lower_) LobattoLegendreMortarL2{RealT, nnodes_, typeof(forward_upper), typeof(reverse_upper)}(forward_upper, forward_lower, diff --git a/src/solvers/dgsem/l2projection.jl b/src/solvers/dgsem/l2projection.jl index 436c1b9203..7f8f820032 100644 --- a/src/solvers/dgsem/l2projection.jl +++ b/src/solvers/dgsem/l2projection.jl @@ -29,9 +29,9 @@ function calc_forward_upper(n_nodes, RealT = Float64) wbary = barycentric_weights(nodes) # Calculate projection matrix (actually: interpolation) - operator = zeros(n_nodes, n_nodes) + operator = zeros(RealT, n_nodes, n_nodes) for j in 1:n_nodes - poly = lagrange_interpolating_polynomials(1 / 2 * (nodes[j] + 1), nodes, wbary) + poly = lagrange_interpolating_polynomials(0.5f0 * (nodes[j] + 1), nodes, wbary) for i in 1:n_nodes operator[j, i] = poly[i] end @@ -49,9 +49,9 @@ function calc_forward_lower(n_nodes, RealT = Float64) wbary = barycentric_weights(nodes) # Calculate projection matrix (actually: interpolation) - operator = zeros(n_nodes, n_nodes) + operator = zeros(RealT, n_nodes, n_nodes) for j in 1:n_nodes - poly = lagrange_interpolating_polynomials(1 / 2 * (nodes[j] - 1), nodes, wbary) + poly = lagrange_interpolating_polynomials(0.5f0 * (nodes[j] - 1), nodes, wbary) for i in 1:n_nodes operator[j, i] = poly[i] end @@ -70,12 +70,12 @@ function calc_reverse_upper(n_nodes, ::Val{:gauss}, RealT = Float64) gauss_wbary = barycentric_weights(gauss_nodes) # Calculate projection matrix (actually: discrete L2 projection with errors) - operator = zeros(n_nodes, n_nodes) + operator = zeros(RealT, n_nodes, n_nodes) for j in 1:n_nodes - poly = lagrange_interpolating_polynomials(1 / 2 * (gauss_nodes[j] + 1), + poly = lagrange_interpolating_polynomials(0.5f0 * (gauss_nodes[j] + 1), gauss_nodes, gauss_wbary) for i in 1:n_nodes - operator[i, j] = 1 / 2 * poly[i] * gauss_weights[j] / gauss_weights[i] + operator[i, j] = 0.5f0 * poly[i] * gauss_weights[j] / gauss_weights[i] end end @@ -97,12 +97,12 @@ function calc_reverse_lower(n_nodes, ::Val{:gauss}, RealT = Float64) gauss_wbary = barycentric_weights(gauss_nodes) # Calculate projection matrix (actually: discrete L2 projection with errors) - operator = zeros(n_nodes, n_nodes) + operator = zeros(RealT, n_nodes, n_nodes) for j in 1:n_nodes - poly = lagrange_interpolating_polynomials(1 / 2 * (gauss_nodes[j] - 1), + poly = lagrange_interpolating_polynomials(0.5f0 * (gauss_nodes[j] - 1), gauss_nodes, gauss_wbary) for i in 1:n_nodes - operator[i, j] = 1 / 2 * poly[i] * gauss_weights[j] / gauss_weights[i] + operator[i, j] = 0.5f0 * poly[i] * gauss_weights[j] / gauss_weights[i] end end