diff --git a/NEWS.md b/NEWS.md index 8da620e3ae5..2570d6fad24 100644 --- a/NEWS.md +++ b/NEWS.md @@ -26,6 +26,7 @@ This enables in particular adaptive mesh refinement for that solver-mesh combina - The second-order subcell volume integral is no longer limited to reconstruction in primitive variables. Instead, it is possible to reconstruct in custom variables, if functions `cons2recon` and `recon2cons` are provided to `VolumeIntegralPureLGLFiniteVolumeO2` and `VolumeIntegralShockCapturingRRG`([#2817]). +- Add Legendre-Gauss basis for DGSEM and implement solver (`VolumeIntegralWeakForm` and `SurfaceIntegralWeakForm` only) support for conforming 1D & 2D `TreeMesh`es ([#1965]). ## Changes when updating to v0.15 from v0.14.x diff --git a/examples/tree_1d_dgsem/elixir_advection_gauss_legendre.jl b/examples/tree_1d_dgsem/elixir_advection_gauss_legendre.jl new file mode 100644 index 00000000000..fed5af082b1 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_gauss_legendre.jl @@ -0,0 +1,42 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = 1.0 +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +basis = GaussLegendreBasis(3) +solver = DGSEM(basis, flux_godunov) + +coordinates_min = -1.0 +coordinates_max = 1.0 + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000, periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +ode = semidiscretize(semi, (0.0, 1.0)) + +summary_callback = SummaryCallback() + +analysis_callback = AnalysisCallback(semi, interval = 100) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, analysis_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false); dt = 1.0, + ode_default_options()..., callback = callbacks); diff --git a/examples/tree_2d_dgsem/elixir_euler_convergence_gauss.jl b/examples/tree_2d_dgsem/elixir_euler_convergence_gauss.jl new file mode 100644 index 00000000000..4ea4831c6e9 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_convergence_gauss.jl @@ -0,0 +1,52 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_convergence_test + +surface_flux = flux_lax_friedrichs + +polydeg = 3 +basis = GaussLegendreBasis(polydeg) +solver = DGSEM(polydeg = 3, basis = basis, surface_flux = surface_flux) + +coordinates_min = (0.0, 0.0) +coordinates_max = (2.0, 2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + n_cells_max = 100_000, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); diff --git a/src/Trixi.jl b/src/Trixi.jl index 99425ec6323..7707b8434bc 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -263,7 +263,7 @@ export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMe P4estMeshView, P4estMeshCubedSphere, T8codeMesh export DG, - DGSEM, LobattoLegendreBasis, + DGSEM, LobattoLegendreBasis, GaussLegendreBasis, FDSBP, VolumeIntegralWeakForm, VolumeIntegralStrongForm, VolumeIntegralFluxDifferencing, diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 20e19539c10..7f5cec6e731 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -315,7 +315,7 @@ function _precompile_manual_() @assert Base.precompile(Tuple{typeof(Trixi.polynomial_interpolation_matrix), Vector{RealT}, Vector{RealT}}) @assert Base.precompile(Tuple{typeof(Trixi.barycentric_weights), Vector{RealT}}) - @assert Base.precompile(Tuple{typeof(Trixi.calc_Lhat), RealT, Vector{RealT}, + @assert Base.precompile(Tuple{typeof(Trixi.calc_L), RealT, Vector{RealT}, Vector{RealT}}) @assert Base.precompile(Tuple{typeof(Trixi.lagrange_interpolating_polynomials), RealT, Vector{RealT}, Vector{RealT}}) diff --git a/src/solvers/dgsem/basis_gauss_legendre.jl b/src/solvers/dgsem/basis_gauss_legendre.jl new file mode 100644 index 00000000000..4854bff2837 --- /dev/null +++ b/src/solvers/dgsem/basis_gauss_legendre.jl @@ -0,0 +1,292 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +""" + GaussLegendreBasis([RealT=Float64,] polydeg::Integer) + +Create a nodal Gauss-Legendre basis for polynomials of degree `polydeg`. +""" +struct GaussLegendreBasis{RealT <: Real, NNODES, + VectorT <: AbstractVector{RealT}, + InverseVandermondeLegendre <: AbstractMatrix{RealT}, + DerivativeMatrix <: AbstractMatrix{RealT}, + BoundaryMatrix <: AbstractMatrix{RealT}} <: + AbstractBasisSBP{RealT} + nodes::VectorT + weights::VectorT + inverse_weights::VectorT + + inverse_vandermonde_legendre::InverseVandermondeLegendre + + derivative_matrix::DerivativeMatrix # strong form derivative matrix + # `derivative_split` currently not implemented since + # Flux-Differencing is not supported for Gauss-Legendre DGSEM at the moment. + derivative_hat::DerivativeMatrix # weak form matrix "dhat", negative adjoint wrt the SBP dot product + + # Required for Gauss-Legendre nodes (non-trivial interpolation to the boundaries) + boundary_interpolation::BoundaryMatrix # L + boundary_interpolation_inverse_weights::BoundaryMatrix # M^{-1} * L = Lhat +end + +function GaussLegendreBasis(RealT, polydeg::Integer) + nnodes_ = polydeg + 1 + + nodes_, weights_ = gauss_nodes_weights(nnodes_, RealT) + inverse_weights_ = inv.(weights_) + + _, inverse_vandermonde_legendre = vandermonde_legendre(nodes_, RealT) + + derivative_matrix = polynomial_derivative_matrix(nodes_) + derivative_hat = calc_Dhat(derivative_matrix, weights_) + + # 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_) + + boundary_interpolation = zeros(RealT, nnodes_, 2) + boundary_interpolation[:, 1] = calc_L(-one(RealT), nodes_, weights_) + boundary_interpolation[:, 2] = calc_L(one(RealT), nodes_, weights_) + + boundary_interpolation_inverse_weights = copy(boundary_interpolation) + boundary_interpolation_inverse_weights[:, 1] = calc_Lhat(boundary_interpolation[:, + 1], + weights_) + boundary_interpolation_inverse_weights[:, 2] = calc_Lhat(boundary_interpolation[:, + 2], + weights_) + + # 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 GaussLegendreBasis{RealT, nnodes_, typeof(nodes), + typeof(inverse_vandermonde_legendre), + typeof(derivative_matrix), + typeof(boundary_interpolation)}(nodes, weights, + inverse_weights, + inverse_vandermonde_legendre, + derivative_matrix, + derivative_hat, + boundary_interpolation, + boundary_interpolation_inverse_weights) +end + +GaussLegendreBasis(polydeg::Integer) = GaussLegendreBasis(Float64, polydeg) + +function Base.show(io::IO, basis::GaussLegendreBasis) + @nospecialize basis # reduce precompilation time + + print(io, "GaussLegendreBasis{", real(basis), "}(polydeg=", polydeg(basis), ")") + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", basis::GaussLegendreBasis) + @nospecialize basis # reduce precompilation time + + print(io, "GaussLegendreBasis{", real(basis), "} with polynomials of degree ", + polydeg(basis)) + return nothing +end + +@inline Base.real(basis::GaussLegendreBasis{RealT}) where {RealT} = RealT + +@inline function nnodes(basis::GaussLegendreBasis{RealT, NNODES}) where {RealT, NNODES} + return NNODES +end + +""" + eachnode(basis::GaussLegendreBasis) + +Return an iterator over the indices that specify the location in relevant data structures +for the nodes in `basis`. +In particular, not the nodes themselves are returned. +""" +@inline eachnode(basis::GaussLegendreBasis) = Base.OneTo(nnodes(basis)) + +@inline polydeg(basis::GaussLegendreBasis) = nnodes(basis) - 1 + +@inline get_nodes(basis::GaussLegendreBasis) = basis.nodes + +""" + integrate(f, u, basis::GaussLegendreBasis) + +Map the function `f` to the coefficients `u` and integrate with respect to the +quadrature rule given by `basis`. +""" +function integrate(f, u, basis::GaussLegendreBasis) + @unpack weights = basis + + res = zero(f(first(u))) + for i in eachindex(u, weights) + res += f(u[i]) * weights[i] + end + return res +end + +# TODO: Not yet implemented +function MortarL2(basis::GaussLegendreBasis) + return nothing +end + +""" + gauss_nodes_weights(n_nodes::Integer, RealT = Float64) + +Computes nodes ``x_j`` and weights ``w_j`` for the Gauss-Legendre quadrature. +This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book + +- David A. Kopriva, (2009). + Implementing spectral methods for partial differential equations: + Algorithms for scientists and engineers. + [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5) +""" +function gauss_nodes_weights(n_nodes::Integer, RealT = Float64) + n_iterations = 20 + tolerance = 2 * eps(RealT) # Relative tolerance for Newton iteration + + # Initialize output + nodes = ones(RealT, n_nodes) + weights = zeros(RealT, n_nodes) + + # Get polynomial degree for convenience + N = n_nodes - 1 + if N == 0 + nodes .= 0 + weights .= 2 + return nodes, weights + elseif N == 1 + nodes[1] = -sqrt(one(RealT) / 3) + nodes[end] = -nodes[1] + weights .= 1 + return nodes, weights + else # N > 1 + # Use symmetry property of the roots of the Legendre polynomials + for i in 0:(div(N + 1, 2) - 1) + # Starting guess for Newton method + nodes[i + 1] = -cospi(one(RealT) / (2 * N + 2) * (2 * i + 1)) + + # Newton iteration to find root of Legendre polynomial (= integration node) + for k in 0:n_iterations + poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i + 1]) + dx = -poly / deriv + nodes[i + 1] += dx + if abs(dx) < tolerance * abs(nodes[i + 1]) + break + end + + if k == n_iterations + @warn "`gauss_nodes_weights` Newton iteration did not converge" + end + end + + # Calculate weight + poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i + 1]) + weights[i + 1] = (2 * N + 3) / ((1 - nodes[i + 1]^2) * deriv^2) + + # Set nodes and weights according to symmetry properties + nodes[N + 1 - i] = -nodes[i + 1] + weights[N + 1 - i] = weights[i + 1] + end + + # If odd number of nodes, set center node to origin (= 0.0) and calculate weight + if n_nodes % 2 == 1 + poly, deriv = legendre_polynomial_and_derivative(N + 1, zero(RealT)) + nodes[div(N, 2) + 1] = 0 + weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2 + end + + return nodes, weights + end +end + +# L(x), where L(x) is the Lagrange polynomial vector at point x. +function calc_L(x, nodes, weights) + n_nodes = length(nodes) + wbary = barycentric_weights(nodes) + + return lagrange_interpolating_polynomials(x, nodes, wbary) +end + +# Calculate M^{-1} * L(x), where L(x) is the Lagrange polynomial +# vector at point x. +# Not required for the DGSEM with LGL basis, as boundary evaluations +# collapse to boundary node evaluations. +function calc_Lhat(L, weights) + Lhat = copy(L) + for i in 1:length(weights) + Lhat[i] /= weights[i] + end + + return Lhat +end + +struct GaussLegendreAnalyzer{RealT <: Real, NNODES, + VectorT <: AbstractVector{RealT}, + Vandermonde <: AbstractMatrix{RealT}} <: + SolutionAnalyzer{RealT} + nodes::VectorT + weights::VectorT + vandermonde::Vandermonde +end + +function SolutionAnalyzer(basis::GaussLegendreBasis; + analysis_polydeg = 2 * polydeg(basis)) + RealT = real(basis) + nnodes_ = analysis_polydeg + 1 + + # compute everything using `Float64` by default + nodes_, weights_ = gauss_nodes_weights(nnodes_) + vandermonde_ = polynomial_interpolation_matrix(get_nodes(basis), nodes_) + + # type conversions to get the requested real type and enable possible + # optimizations of runtime performance and latency + nodes = SVector{nnodes_, RealT}(nodes_) + weights = SVector{nnodes_, RealT}(weights_) + + vandermonde = Matrix{RealT}(vandermonde_) + + return GaussLegendreAnalyzer{RealT, nnodes_, typeof(nodes), typeof(vandermonde)}(nodes, + weights, + vandermonde) +end + +function Base.show(io::IO, analyzer::GaussLegendreAnalyzer) + @nospecialize analyzer # reduce precompilation time + + print(io, "GaussLegendreAnalyzer{", real(analyzer), "}(polydeg=", + polydeg(analyzer), ")") + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", analyzer::GaussLegendreAnalyzer) + @nospecialize analyzer # reduce precompilation time + + print(io, "GaussLegendreAnalyzer{", real(analyzer), + "} with polynomials of degree ", polydeg(analyzer)) + return nothing +end + +@inline Base.real(analyzer::GaussLegendreAnalyzer{RealT}) where {RealT} = RealT + +@inline function nnodes(analyzer::GaussLegendreAnalyzer{RealT, NNODES}) where {RealT, + NNODES} + return NNODES +end + +""" + eachnode(analyzer::GaussLegendreAnalyzer) + +Return an iterator over the indices that specify the location in relevant data structures +for the nodes in `analyzer`. +In particular, not the nodes themselves are returned. +""" +@inline eachnode(analyzer::GaussLegendreAnalyzer) = Base.OneTo(nnodes(analyzer)) + +@inline polydeg(analyzer::GaussLegendreAnalyzer) = nnodes(analyzer) - 1 +end # @muladd diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index 0d79804e48d..61aa9820a74 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -513,23 +513,6 @@ function barycentric_weights(nodes) return weights end -# Calculate M^{-1} * L(x), where L(x) is the Lagrange polynomial -# vector at point x. -# Not required for the DGSEM with LGL basis, as boundary evaluations -# collapse to boundary node evaluations. -function calc_Lhat(x, nodes, weights) - n_nodes = length(nodes) - wbary = barycentric_weights(nodes) - - Lhat = lagrange_interpolating_polynomials(x, nodes, wbary) - - for i in 1:n_nodes - Lhat[i] /= weights[i] - end - - return Lhat -end - """ lagrange_interpolating_polynomials(x, nodes, wbary) @@ -593,7 +576,8 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer, RealT = Float64) nodes = zeros(RealT, n_nodes) weights = zeros(RealT, n_nodes) - # Special case for polynomial degree zero (first order finite volume) + # Special case for polynomial degree zero (first order finite volume): + # Fall back to Gauss-Legendre if n_nodes == 1 nodes[1] = 0 weights[1] = 2 @@ -679,76 +663,6 @@ function calc_q_and_l(N::Integer, x::Real) return q, qder, L end -""" - gauss_nodes_weights(n_nodes::Integer, RealT = Float64) - -Computes nodes ``x_j`` and weights ``w_j`` for the Gauss-Legendre quadrature. -This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book - -- David A. Kopriva, (2009). - Implementing spectral methods for partial differential equations: - Algorithms for scientists and engineers. - [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5) -""" -function gauss_nodes_weights(n_nodes::Integer, RealT = Float64) - n_iterations = 20 - tolerance = 2 * eps(RealT) # Relative tolerance for Newton iteration - - # Initialize output - nodes = ones(RealT, n_nodes) - weights = zeros(RealT, n_nodes) - - # Get polynomial degree for convenience - N = n_nodes - 1 - if N == 0 - nodes .= 0 - weights .= 2 - return nodes, weights - elseif N == 1 - nodes[1] = -sqrt(one(RealT) / 3) - nodes[end] = -nodes[1] - weights .= 1 - return nodes, weights - else # N > 1 - # Use symmetry property of the roots of the Legendre polynomials - for i in 0:(div(N + 1, 2) - 1) - # Starting guess for Newton method - nodes[i + 1] = -cospi(one(RealT) / (2 * N + 2) * (2 * i + 1)) - - # Newton iteration to find root of Legendre polynomial (= integration node) - for k in 0:n_iterations - poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i + 1]) - dx = -poly / deriv - nodes[i + 1] += dx - if abs(dx) < tolerance * abs(nodes[i + 1]) - break - end - - if k == n_iterations - @warn "`gauss_nodes_weights` Newton iteration did not converge" - end - end - - # Calculate weight - poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i + 1]) - weights[i + 1] = (2 * N + 3) / ((1 - nodes[i + 1]^2) * deriv^2) - - # Set nodes and weights according to symmetry properties - nodes[N + 1 - i] = -nodes[i + 1] - weights[N + 1 - i] = weights[i + 1] - end - - # If odd number of nodes, set center node to origin (= 0.0) and calculate weight - if n_nodes % 2 == 1 - poly, deriv = legendre_polynomial_and_derivative(N + 1, zero(RealT)) - nodes[div(N, 2) + 1] = 0 - weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2 - end - - return nodes, weights - end -end - """ legendre_polynomial_and_derivative(N::Int, x::Real) diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl index e5b55fbe180..925d2770dfc 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/src/solvers/dgsem/dgsem.jl @@ -9,9 +9,11 @@ include("interpolation.jl") include("l2projection.jl") include("basis_lobatto_legendre.jl") +include("basis_gauss_legendre.jl") """ DGSEM(; RealT=Float64, polydeg::Integer, + basis = LobattoLegendreBasis(RealT, polydeg) surface_flux=flux_central, surface_integral=SurfaceIntegralWeakForm(surface_flux), volume_integral=VolumeIntegralWeakForm()) @@ -19,10 +21,10 @@ include("basis_lobatto_legendre.jl") Create a discontinuous Galerkin spectral element method (DGSEM) using a [`LobattoLegendreBasis`](@ref) with polynomials of degree `polydeg`. """ -const DGSEM = DG{Basis} where {Basis <: LobattoLegendreBasis} +const DGSEM = DG{Basis} where {Basis <: AbstractBasisSBP} # This API is no longer documented, and we recommend avoiding its public use. -function DGSEM(basis::LobattoLegendreBasis, +function DGSEM(basis::AbstractBasisSBP, surface_flux = flux_central, volume_integral = VolumeIntegralWeakForm(), mortar = MortarL2(basis)) @@ -32,7 +34,7 @@ function DGSEM(basis::LobattoLegendreBasis, end # This API is no longer documented, and we recommend avoiding its public use. -function DGSEM(basis::LobattoLegendreBasis, +function DGSEM(basis::AbstractBasisSBP, surface_integral::AbstractSurfaceIntegral, volume_integral = VolumeIntegralWeakForm(), mortar = MortarL2(basis)) @@ -61,10 +63,10 @@ end # `trixi_include`. function DGSEM(; RealT = Float64, polydeg::Integer, + basis = LobattoLegendreBasis(RealT, polydeg), surface_flux = flux_central, surface_integral = SurfaceIntegralWeakForm(surface_flux), volume_integral = VolumeIntegralWeakForm()) - basis = LobattoLegendreBasis(RealT, polydeg) return DGSEM(basis, surface_integral, volume_integral) end diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index e3a6c26a9f0..d4cd609c6bd 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -759,15 +759,14 @@ end function calc_surface_integral!(du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, - surface_integral::SurfaceIntegralWeakForm, + equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements # Note that all fluxes have been computed with outward-pointing normal vectors. # This computes the **negative** surface integral contribution, - # i.e., M^{-1} * boundary_interpolation^T (which is for DGSEM just M^{-1} * B) + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) # and the missing "-" is taken care of by `apply_jacobian!`. # # We also use explicit assignments instead of `+=` to let `@muladd` turn these diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index 5fb332b483e..b70661c56ba 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -924,15 +924,14 @@ end function calc_surface_integral!(du, u, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, - equations, - surface_integral::SurfaceIntegralWeakForm, + equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements # Note that all fluxes have been computed with outward-pointing normal vectors. # This computes the **negative** surface integral contribution, - # i.e., M^{-1} * boundary_interpolation^T (which is for DGSEM just M^{-1} * B) + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) # and the missing "-" is taken care of by `apply_jacobian!`. # # We also use explicit assignments instead of `+=` to let `@muladd` turn these diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 95719c5b630..c472e071ab3 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -396,6 +396,36 @@ function prolong2interfaces!(cache, u_or_flux_viscous, return nothing end +function prolong2interfaces!(cache, u_or_flux_viscous, + mesh::TreeMesh{1}, equations, + dg::DG{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack neighbor_ids = interfaces + @unpack boundary_interpolation = dg.basis + interfaces_u = interfaces.u + + @threaded for interface in eachinterface(dg, cache) + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] + + # interface in x-direction + for v in eachvariable(equations) + interfaces_u[1, v, interface] = zero(eltype(interfaces_u)) + interfaces_u[2, v, interface] = zero(eltype(interfaces_u)) + for ii in eachnode(dg) + interfaces_u[1, v, interface] += (u_or_flux_viscous[v, ii, + left_element] * + boundary_interpolation[ii, 2]) + interfaces_u[2, v, interface] += (u_or_flux_viscous[v, ii, + right_element] * + boundary_interpolation[ii, 1]) + end + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{1}, have_nonconservative_terms::False, equations, @@ -587,12 +617,13 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations, surface_integral, dg::DGSEM, cache) + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM, cache) @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements # This computes the **negative** surface integral contribution, - # i.e., M^{-1} * boundary_interpolation^T (which is for DGSEM just M^{-1} * B) + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) # and the missing "-" is taken care of by `apply_jacobian!`. # # We also use explicit assignments instead of `+=` to let `@muladd` turn these @@ -615,6 +646,37 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1 return nothing end +function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DG{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + # This computes the **negative** surface integral contribution, + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Legendre DGSEM M^{-1} * L) + # and the missing "-" is taken care of by `apply_jacobian!`. + # + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + @threaded for element in eachelement(dg, cache) + for v in eachvariable(equations) + for ii in eachnode(dg) + # surface at -x + du[v, ii, element] = (du[v, ii, element] - + surface_flux_values[v, 1, element] * + boundary_interpolation_inverse_weights[ii, 1]) + + # surface at +x + du[v, ii, element] = (du[v, ii, element] + + surface_flux_values[v, 2, element] * + boundary_interpolation_inverse_weights[ii, 2]) + end + end + end + + return nothing +end + function apply_jacobian!(du, mesh::TreeMesh{1}, equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index fe58f057725..c1a845ab06d 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -523,6 +523,47 @@ function prolong2interfaces!(cache, u, mesh::TreeMesh{2}, equations, dg::DG) return nothing end +function prolong2interfaces!(cache, u, mesh::TreeMesh{2}, equations, + dg::DG{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack orientations, neighbor_ids = interfaces + @unpack boundary_interpolation = dg.basis + interfaces_u = interfaces.u + + @threaded for interface in eachinterface(dg, cache) + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] + + if orientations[interface] == 1 + # interface in x-direction + for j in eachnode(dg), v in eachvariable(equations) + interfaces_u[1, v, j, interface] = zero(eltype(interfaces_u)) + interfaces_u[2, v, j, interface] = zero(eltype(interfaces_u)) + for ii in eachnode(dg) + interfaces_u[1, v, j, interface] += (u[v, ii, j, left_element] * + boundary_interpolation[ii, 2]) + interfaces_u[2, v, j, interface] += (u[v, ii, j, right_element] * + boundary_interpolation[ii, 1]) + end + end + else # if orientations[interface] == 2 + # interface in y-direction + for i in eachnode(dg), v in eachvariable(equations) + interfaces_u[1, v, i, interface] = zero(eltype(interfaces_u)) + interfaces_u[2, v, i, interface] = zero(eltype(interfaces_u)) + for jj in eachnode(dg) + interfaces_u[1, v, i, interface] += (u[v, i, jj, left_element] * + boundary_interpolation[jj, 2]) + interfaces_u[2, v, i, interface] += (u[v, i, jj, right_element] * + boundary_interpolation[jj, 1]) + end + end + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, have_nonconservative_terms::False, equations, @@ -1088,7 +1129,7 @@ function calc_surface_integral!(du, u, @unpack surface_flux_values = cache.elements # This computes the **negative** surface integral contribution, - # i.e., M^{-1} * boundary_interpolation^T (which is for DGSEM just M^{-1} * B) + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) # and the missing "-" is taken care of by `apply_jacobian!`. # # We also use explicit assignments instead of `+=` to let `@muladd` turn these @@ -1123,6 +1164,57 @@ function calc_surface_integral!(du, u, return nothing end +function calc_surface_integral!(du, u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DG{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + # This computes the **negative** surface integral contribution, + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Legendre DGSEM M^{-1} * L) + # and the missing "-" is taken care of by `apply_jacobian!`. + # + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + for v in eachvariable(equations) + for ii in eachnode(dg) + # surface at -x + du[v, ii, l, element] = (du[v, ii, l, element] - + surface_flux_values[v, l, 1, element] * + boundary_interpolation_inverse_weights[ii, + 1]) + + # surface at +x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_values[v, l, 2, element] * + boundary_interpolation_inverse_weights[ii, + 2]) + end + + for jj in eachnode(dg) + # surface at -y + du[v, l, jj, element] = (du[v, l, jj, element] - + surface_flux_values[v, l, 3, element] * + boundary_interpolation_inverse_weights[jj, + 1]) + + # surface at +y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_values[v, l, 4, element] * + boundary_interpolation_inverse_weights[jj, + 2]) + end + end + end + end + + return nothing +end + function apply_jacobian!(du, mesh::TreeMesh{2}, equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 247296b2509..41b46f58f95 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -1334,12 +1334,13 @@ end end function calc_surface_integral!(du, u, mesh::Union{TreeMesh{3}, StructuredMesh{3}}, - equations, surface_integral, dg::DGSEM, cache) + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM, cache) @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements # This computes the **negative** surface integral contribution, - # i.e., M^{-1} * boundary_interpolation^T (which is for DGSEM just M^{-1} * B) + # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) # and the missing "-" is taken care of by `apply_jacobian!`. # # We also use explicit assignments instead of `+=` and `-=` to let `@muladd` diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index f54d90d85a4..76a186b9fa6 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -484,7 +484,7 @@ function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, # Note that all fluxes have been computed with outward-pointing normal vectors. # This computes the **negative** surface integral contribution, - # i.e., M^{-1} * boundary_interpolation^T (which is for DGSEM just M^{-1} * B) + # i.e., M^{-1} * boundary_interpolation^T # and the missing "-" is taken care of by `apply_jacobian!`. # # We also use explicit assignments instead of `+=` and `-=` to let `@muladd` diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index 876b326127f..f1a2a7c91c1 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -19,6 +19,14 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_1d_dgsem") @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_advection_gauss_legendre.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_gauss_legendre.jl"), + l2=[2.515203865524688e-6], linf=[8.660338936650191e-6]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_advection_basic.jl (max_abs_speed)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), surface_flux=FluxLaxFriedrichs(max_abs_speed), diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index c8afa8c141a..92da4ec82d9 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -90,6 +90,25 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_euler_convergence_gauss.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_gauss.jl"), + l2=[ + 0.0001657393364512246, + 0.00018603701552875171, + 0.00018603701552861147, + 0.0006686395458793184 + ], + linf=[ + 0.00044692901513743166, + 0.0005793901371469179, + 0.000579390137147362, + 0.002266770997066736 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_euler_density_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index 59060593d64..707f23c2813 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -316,10 +316,12 @@ end nodes = basis.nodes weights = basis.weights - Lhat_minus1 = Trixi.calc_Lhat(-1.0, nodes, weights) + L_minus1 = Trixi.calc_L(-1.0, nodes, weights) + Lhat_minus1 = Trixi.calc_Lhat(L_minus1, weights) @test basis.inverse_weights[1] == Lhat_minus1[1] - Lhat_plus1 = Trixi.calc_Lhat(1.0, nodes, weights) + L_plus1 = Trixi.calc_L(1.0, nodes, weights) + Lhat_plus1 = Trixi.calc_Lhat(L_plus1, weights) @test basis.inverse_weights[p + 1] == Lhat_plus1[p + 1] end end @@ -361,6 +363,16 @@ end end end +@timed_testset "GaussLegendreBasis" begin + basis = GaussLegendreBasis(3) + @test nnodes(basis) == 4 + @test_nowarn show(stdout, "text/plain", basis) + + solution_analyzer = Trixi.SolutionAnalyzer(basis) + @test nnodes(solution_analyzer) == 7 + @test_nowarn show(stdout, "text/plain", solution_analyzer) +end + @testset "containers" begin # Set up mock container mutable struct MyContainer <: Trixi.AbstractContainer