From 306d47eb4f07f064c554620ac5773afcbf3f9079 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 29 May 2024 15:15:28 +0100 Subject: [PATCH] Implement the Gauss-Legendre basis for DGSEM --- .../elixir_euler_convergence_pure_fv.jl | 14 +- src/Trixi.jl | 2 +- src/solvers/dgsem/basis_gauss_legendre.jl | 240 ++++++++++++++++++ src/solvers/dgsem/dgsem.jl | 11 +- src/solvers/dgsem_tree/dg_2d.jl | 86 ++++++- 5 files changed, 344 insertions(+), 9 deletions(-) create mode 100644 src/solvers/dgsem/basis_gauss_legendre.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl b/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl index 96184b5ba47..ac4b3709465 100644 --- a/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl +++ b/examples/tree_2d_dgsem/elixir_euler_convergence_pure_fv.jl @@ -11,14 +11,18 @@ initial_condition = initial_condition_convergence_test surface_flux = flux_hllc -basis = LobattoLegendreBasis(3) -volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc) +# basis = LobattoLegendreBasis(3) +# volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc) +# solver = DGSEM(basis, surface_flux, volume_integral) +polydeg = 4 +basis = GaussLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 4) +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (2.0, 2.0) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, + initial_refinement_level = 2, n_cells_max = 10_000) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, @@ -32,7 +36,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 100 * 10 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) @@ -46,7 +50,7 @@ stepsize_callback = StepsizeCallback(cfl = 0.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + # save_solution, stepsize_callback) ############################################################################### diff --git a/src/Trixi.jl b/src/Trixi.jl index d0ee4b6f2a0..392138b0f8d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -230,7 +230,7 @@ export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMe T8codeMesh export DG, - DGSEM, LobattoLegendreBasis, + DGSEM, LobattoLegendreBasis, GaussLegendreBasis, FDSBP, VolumeIntegralWeakForm, VolumeIntegralStrongForm, VolumeIntegralWeakFormProjection, diff --git a/src/solvers/dgsem/basis_gauss_legendre.jl b/src/solvers/dgsem/basis_gauss_legendre.jl new file mode 100644 index 00000000000..b9aafeb9da5 --- /dev/null +++ b/src/solvers/dgsem/basis_gauss_legendre.jl @@ -0,0 +1,240 @@ +# 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}, + BoundaryMatrix <: AbstractMatrix{RealT}, + DerivativeMatrix <: AbstractMatrix{RealT}, + Matrix_MxN <: AbstractMatrix{RealT}, + Matrix_NxM <: AbstractMatrix{RealT}, + Matrix_MxM <: AbstractMatrix{RealT}} <: + AbstractBasisSBP{RealT} + nodes::VectorT + weights::VectorT + inverse_weights::VectorT + + inverse_vandermonde_legendre::InverseVandermondeLegendre + boundary_interpolation::BoundaryMatrix # lhat + + derivative_matrix::DerivativeMatrix # strong form derivative matrix + derivative_split::DerivativeMatrix # strong form derivative matrix minus boundary terms + derivative_split_transpose::DerivativeMatrix # transpose of `derivative_split` + derivative_dhat::DerivativeMatrix # weak form matrix "dhat", + # negative adjoint wrt the SBP dot product + # L2 projection operators + interpolate_N_to_M::Matrix_MxN # interpolates from N nodes to M nodes + project_M_to_N::Matrix_NxM # compute projection via Legendre modes, cut off modes at N, back to N nodes + filter_modal_to_N::Matrix_MxM # compute modal filter via Legendre modes, cut off modes at N, leave it at M nodes + filter_modal_to_cutoff::DerivativeMatrix # compute modal cutoff filter via Legendre modes, cut off modes at polydeg_cutoff +end + +function GaussLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) + nnodes_ = polydeg + 1 + + # compute everything using `Float64` by default + nodes_, weights_ = gauss_nodes_weights(nnodes_) + inverse_weights_ = inv.(weights_) + + _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_) + + boundary_interpolation_ = zeros(nnodes_, 2) + boundary_interpolation_[:, 1] = calc_lhat(-1.0, nodes_, weights_) + boundary_interpolation_[:, 2] = calc_lhat(1.0, 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 + 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_) + + # L2 projection operators + nnodes_projection = polydeg_projection + 1 + nodes_projection, weights_projection = gauss_nodes_weights(nnodes_projection) + interpolate_N_to_M_ = polynomial_interpolation_matrix(nodes_, nodes_projection) + interpolate_N_to_M = Matrix{RealT}(interpolate_N_to_M_) + + project_M_to_N_,filter_modal_to_N_ = calc_projection_matrix(nodes_projection, nodes_) + project_M_to_N = Matrix{RealT}(project_M_to_N_) + filter_modal_to_N = Matrix{RealT}(filter_modal_to_N_) + + nnodes_cutoff = polydeg_cutoff + 1 + nodes_cutoff, weights_cutoff = gauss_nodes_weights(nnodes_cutoff) + _, filter_modal_to_cutoff_ = calc_projection_matrix(nodes_, nodes_cutoff) + filter_modal_to_cutoff = Matrix{RealT}(filter_modal_to_cutoff_) + + return GaussLegendreBasis{RealT, nnodes_, typeof(nodes), + typeof(inverse_vandermonde_legendre), + typeof(boundary_interpolation), + typeof(derivative_matrix), + typeof(interpolate_N_to_M), + typeof(project_M_to_N), + typeof(filter_modal_to_N)}(nodes, weights, + inverse_weights, + inverse_vandermonde_legendre, + boundary_interpolation, + derivative_matrix, + derivative_split, + derivative_split_transpose, + derivative_dhat, + interpolate_N_to_M, + project_M_to_N, + filter_modal_to_N, + filter_modal_to_cutoff) +end + +GaussLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) = GaussLegendreBasis(Float64, polydeg; polydeg_projection, polydeg_cutoff) + +function Base.show(io::IO, basis::GaussLegendreBasis) + @nospecialize basis # reduce precompilation time + + print(io, "GaussLegendreBasis{", real(basis), "}(polydeg=", polydeg(basis), ")") +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)) +end + +function Base.:(==)(b1::GaussLegendreBasis, b2::GaussLegendreBasis) + if typeof(b1) != typeof(b2) + return false + end + + for field in fieldnames(typeof(b1)) + if getfield(b1, field) != getfield(b2, field) + return false + end + end + + return true +end + +@inline Base.real(basis::GaussLegendreBasis{RealT}) where {RealT} = RealT + +@inline function nnodes(basis::GaussLegendreBasis{RealT, NNODES}) where {RealT, NNODES + } + 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 + +# Return the first/last weight of the quadrature associated with `basis`. +# Since the mass matrix for nodal Gauss-Legendre bases is diagonal, +# these weights are the only coefficients necessary for the scaling of +# surface terms/integrals in DGSEM. +left_boundary_weight(basis::GaussLegendreBasis) = first(basis.weights) +right_boundary_weight(basis::GaussLegendreBasis) = last(basis.weights) + +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), ")") +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)) +end + +@inline Base.real(analyzer::GaussLegendreAnalyzer{RealT}) where {RealT} = RealT + +@inline function nnodes(analyzer::GaussLegendreAnalyzer{RealT, NNODES}) where {RealT, + NNODES} + 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/dgsem.jl b/src/solvers/dgsem/dgsem.jl index 27caad4d2dc..d3d7b9f9310 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/src/solvers/dgsem/dgsem.jl @@ -9,6 +9,7 @@ include("interpolation.jl") include("l2projection.jl") include("basis_lobatto_legendre.jl") +include("basis_gauss_legendre.jl") """ DGSEM(; RealT=Float64, polydeg::Integer, @@ -20,7 +21,7 @@ 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} # TODO: Deprecated in v0.3 (no longer documented) function DGSEM(basis::LobattoLegendreBasis, @@ -31,6 +32,14 @@ function DGSEM(basis::LobattoLegendreBasis, return DG{typeof(basis), typeof(mortar), typeof(surface_integral), typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral) end +function DGSEM(basis::GaussLegendreBasis, + surface_flux = flux_central, + volume_integral = VolumeIntegralWeakForm(), + mortar = nothing) + surface_integral = SurfaceIntegralWeakForm(surface_flux) + return DG{typeof(basis), typeof(mortar), typeof(surface_integral), + typeof(volume_integral)}(basis, mortar, surface_integral, volume_integral) +end # TODO: Deprecated in v0.3 (no longer documented) function DGSEM(basis::LobattoLegendreBasis, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 50b8e7eac35..d087e4dcd9d 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -720,7 +720,7 @@ end end function prolong2interfaces!(cache, u, - mesh::TreeMesh{2}, equations, surface_integral, dg::DG) + mesh::TreeMesh{2}, equations, surface_integral, dg::DG{<:LobattoLegendreBasis}) @unpack interfaces = cache @unpack orientations, neighbor_ids = interfaces interfaces_u = interfaces.u @@ -747,6 +747,43 @@ function prolong2interfaces!(cache, u, return nothing end +function prolong2interfaces!(cache, u, + mesh::TreeMesh{2}, equations, surface_integral, dg::DG{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack orientations, neighbor_ids = interfaces + @unpack boundary_interpolation, weights = 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] = 0 + interfaces_u[2, v, j, interface] = 0 + for ii in eachnode(dg) + interfaces_u[1, v, j, interface] += u[v, ii, j, left_element] * weights[ii] * boundary_interpolation[ii, 2] + interfaces_u[2, v, j, interface] += u[v, ii, j, right_element] * weights[ii] * 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] = 0 + interfaces_u[2, v, i, interface] = 0 + for jj in eachnode(dg) + interfaces_u[1, v, i, interface] += u[v, i, jj, left_element] * weights[jj] * boundary_interpolation[jj, 2] + interfaces_u[2, v, i, interface] += u[v, i, jj, right_element] * weights[jj] * boundary_interpolation[jj, 1] + end + end + end + end + + return nothing +end + function calc_interface_flux!(surface_flux_values, mesh::TreeMesh{2}, nonconservative_terms::False, equations, @@ -1271,7 +1308,7 @@ function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DG, cache) + dg::DG{<:LobattoLegendreBasis}, cache) @unpack boundary_interpolation = dg.basis @unpack surface_flux_values = cache.elements @@ -1310,6 +1347,51 @@ 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 = dg.basis + @unpack surface_flux_values = cache.elements + + # Note that all fluxes have been computed with outward-pointing normal vectors. + # Access the factors only once before beginning the loop to increase performance. + # 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[ii, 1]) + + # surface at +x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_values[v, l, 2, element] * + boundary_interpolation[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[jj, 1]) + + # surface at +y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_values[v, l, 4, element] * + boundary_interpolation[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