From 8f3c9550bcc1af259eb8a4c4233f1fbe07f094aa Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 15:42:45 +0100 Subject: [PATCH 1/7] current status of entropy projection --- src/auxiliary/precompile.jl | 6 +- .../positivity_zhang_shu_dg2d.jl | 23 ++++-- src/solvers/dgsem/basis_lobatto_legendre.jl | 16 ++-- src/solvers/dgsem_tree/dg_2d.jl | 73 +++++++++++++++---- 4 files changed, 86 insertions(+), 32 deletions(-) diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index 70ebf923b87..860dc95a9d9 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -340,9 +340,9 @@ function _precompile_manual_() for RealT in (Float64,), polydeg in 1:7 nnodes_ = polydeg + 1 basis_type = basis_type_dgsem(RealT, nnodes_) - @assert Base.precompile(Tuple{typeof(Trixi.MortarL2), basis_type}) - @assert Base.precompile(Tuple{Type{Trixi.SolutionAnalyzer}, basis_type}) - @assert Base.precompile(Tuple{Type{Trixi.AdaptorL2}, basis_type}) +# @assert Base.precompile(Tuple{typeof(Trixi.MortarL2), basis_type}) +# @assert Base.precompile(Tuple{Type{Trixi.SolutionAnalyzer}, basis_type}) +# @assert Base.precompile(Tuple{Type{Trixi.AdaptorL2}, basis_type}) end # Constructors: callbacks diff --git a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl index e56b900131a..a9aebea8353 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg2d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg2d.jl @@ -56,7 +56,7 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr eq_relax = equations.equations_relaxation # prepare local storage for projection - @unpack interpolate_N_to_M, project_M_to_N = dg.basis + @unpack interpolate_N_to_M, project_M_to_N, filter_modal_to_N = dg.basis nnodes_,nnodes_projection = size(project_M_to_N) nVars = nvariables(eq_relax) RealT = real(dg) @@ -65,10 +65,15 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr f_N = zeros(RealT, nVars, nnodes_, nnodes_) g_N = zeros(RealT, nVars, nnodes_, nnodes_) u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M_raw = zeros(RealT, nVars, nnodes_projection, nnodes_projection) w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxM = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) + #@threaded for element in eachelement(dg, cache) for element in eachelement(dg, cache) @@ -80,18 +85,20 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr end end # bring elemtn u_N to grid (M+1)x(M+1) - multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N) + multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) # compute nodal values of entropy variables w on the M grid for j in 1:nnodes_projection, i in 1:nnodes_projection u_cons = get_node_vars(u_M, eq_relax, dg, i, j) w_ij = cons2entropy(u_cons,eq_relax) - set_node_vars!(w_M,w_ij,eq_relax,dg,i,j) + set_node_vars!(w_M_raw,w_ij,eq_relax,dg,i,j) end # compute projection of w with M values down to N - multiply_dimensionwise!(w_N,project_M_to_N,w_M) - # use w now to compute new u on grid M - multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + multiply_dimensionwise!(w_M,filter_modal_to_N,w_M_raw,tmp_MxM) + + #multiply_dimensionwise!(w_N,project_M_to_N,w_M) + #multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + # compute nodal values of conservative f,g on the M grid for j in 1:nnodes_projection, i in 1:nnodes_projection @@ -103,8 +110,8 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations::JinXinCompr set_node_vars!(g_M,g_cons,eq_relax,dg,i,j) end # compute projection of f with M values down to N, same for g - multiply_dimensionwise!(f_N,project_M_to_N,f_M) - multiply_dimensionwise!(g_N,project_M_to_N,g_M) + multiply_dimensionwise!(f_N,project_M_to_N,f_M,tmp_NxM) + multiply_dimensionwise!(g_N,project_M_to_N,g_M,tmp_NxM) #@assert nnodes_projection == nnodes(dg) #for j in 1:nnodes_projection, i in 1:nnodes_projection # u_cons = get_node_vars(u_N, eq_relax, dg, i, j) diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index d1c9b909c6e..aadc60f18eb 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -20,7 +20,8 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, BoundaryMatrix <: AbstractMatrix{RealT}, DerivativeMatrix <: AbstractMatrix{RealT}, Matrix_MxN <: AbstractMatrix{RealT}, - Matrix_NxM <: AbstractMatrix{RealT}} <: + Matrix_NxM <: AbstractMatrix{RealT}, + Matrix_MxM <: AbstractMatrix{RealT}} <: AbstractBasisSBP{RealT} nodes::VectorT weights::VectorT @@ -37,6 +38,7 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, # 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 end function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) @@ -78,15 +80,17 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ interpolate_N_to_M_ = polynomial_interpolation_matrix(nodes_, nodes_projection) interpolate_N_to_M = Matrix{RealT}(interpolate_N_to_M_) - project_M_to_N_ = calc_projection_matrix(nodes_projection, nodes_) + 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_) return LobattoLegendreBasis{RealT, nnodes_, typeof(nodes), typeof(inverse_vandermonde_legendre), typeof(boundary_interpolation), typeof(derivative_matrix), typeof(interpolate_N_to_M), - typeof(project_M_to_N)}(nodes, weights, + typeof(project_M_to_N), + typeof(filter_modal_to_N)}(nodes, weights, inverse_weights, inverse_vandermonde_legendre, boundary_interpolation, @@ -95,7 +99,8 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ derivative_split_transpose, derivative_dhat, interpolate_N_to_M, - project_M_to_N) + project_M_to_N, + filter_modal_to_N) end LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection) @@ -811,8 +816,9 @@ function calc_projection_matrix(nodes_in,nodes_out) filter_matrix[j,j] = 1 end interpolate_M_to_N = polynomial_interpolation_matrix(nodes_in, nodes_out) + filter_modal = vandermonde_in * filter_matrix * inverse_vandermonde_in projection_matrix = interpolate_M_to_N * vandermonde_in * filter_matrix * inverse_vandermonde_in - return projection_matrix + return projection_matrix, filter_modal end end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index f3a8b35af70..1ea05fbe97a 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -237,31 +237,72 @@ function calc_volume_integral!(du, u, # prepare local storage for projection @unpack interpolate_N_to_M, project_M_to_N = dg.basis nnodes_,nnodes_projection = size(project_M_to_N) - nVars = size(u, 1) + nVars = nvariables(equations) RealT = real(dg) - u_N = zeros(RealT, nVars, nnodes_, nnodes_) - w_N = zeros(RealT, nVars, nnodes_, nnodes_) - f_N = zeros(RealT, nVars, nnodes_, nnodes_) - g_N = zeros(RealT, nVars, nnodes_, nnodes_) - u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) - cache_projection = (; u_N, w_N, f_N, g_N, u_M, w_M, f_M, g_M) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + f_N = zeros(RealT, nVars, nnodes_, nnodes_) + g_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M_raw = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + f_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + g_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + + tmp_MxM = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) @threaded for element in eachelement(dg, cache) - weak_form_kernel_projection!(du, u, element, mesh, + # get element u_N + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + for v in eachvariable(equations) + u_N[v,i,j] = u_node[v] + end + end + # bring elemtn u_N to grid (M+1)x(M+1) + multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) + + # compute nodal values of entropy variables w on the M grid + for j in 1:nnodes_projection, i in 1:nnodes_projection + u_cons = get_node_vars(u_M, equations, dg, i, j) + w_ij = cons2entropy(u_cons,equations) + set_node_vars!(w_M_raw,w_ij,equations,dg,i,j) + end + # compute projection of w with M values down to N + multiply_dimensionwise!(w_M,filter_modal_to_N,w_M_raw,tmp_MxM) + + #multiply_dimensionwise!(w_N,project_M_to_N,w_M) + #multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) + + + # compute nodal values of conservative f,g on the M grid + for j in 1:nnodes_projection, i in 1:nnodes_projection + w_ij = get_node_vars(w_M, equations, dg, i, j) + u_cons = entropy2cons(w_ij, equations) + f_cons = flux(u_cons,1,equations) + set_node_vars!(f_M,f_cons,equations,dg,i,j) + g_cons = flux(u_cons,2,equations) + set_node_vars!(g_M,g_cons,equations,dg,i,j) + end + # compute projection of f with M values down to N, same for g + multiply_dimensionwise!(f_N,project_M_to_N,f_M,tmp_NxM) + multiply_dimensionwise!(g_N,project_M_to_N,g_M,tmp_NxM) + + + weak_form_kernel_projection!(du, u,f_N, g_N, element, mesh, nonconservative_terms, equations, - dg, cache, cache_projection) + dg, cache) end return nothing end -@inline function weak_form_kernel_projection!(du, u, +@inline function weak_form_kernel_projection!(du, u, f_N, g_N, element, mesh::TreeMesh{2}, nonconservative_terms::False, equations, - dg::DGSEM, cache, cache_projection) + dg::DGSEM, cache) # true * [some floating point value] == [exactly the same floating point value] # This can (hopefully) be optimized away due to constant propagation. @unpack derivative_dhat = dg.basis @@ -270,13 +311,13 @@ end for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, element) - flux1 = flux(u_node, 1, equations) + flux1 = get_node_vars(f_N, equations, dg, i,j) for ii in eachnode(dg) multiply_add_to_node_vars!(du, derivative_dhat[ii, i], flux1, equations, dg, ii, j, element) end - flux2 = flux(u_node, 2, equations) + flux2 = get_node_vars(g_N, equations, dg, i,j) for jj in eachnode(dg) multiply_add_to_node_vars!(du, derivative_dhat[jj, j], flux2, equations, dg, i, jj, element) From 811c9637f2471258cdfa4965670ed877ad00ea20 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Tue, 28 May 2024 16:10:58 +0100 Subject: [PATCH 2/7] first draft of projection dgsem --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 12 ++- ...kelvin_helmholtz_instability_projection.jl | 88 +++++++++++++++++++ src/solvers/dgsem_tree/dg_2d.jl | 2 +- 3 files changed, 98 insertions(+), 4 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index e634a383cdf..c32aae817a7 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -8,9 +8,15 @@ equations = CompressibleEulerEquations2D(1.4) initial_condition = initial_condition_weak_blast_wave -volume_flux = flux_ranocha -solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +#volume_flux = flux_ranocha +#solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +# +surface_flux = flux_ranocha +polydeg = 1 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 4 * polydeg) +volume_integral = VolumeIntegralWeakFormProjection() +solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-2.0, -2.0) coordinates_max = (2.0, 2.0) diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl new file mode 100644 index 00000000000..9232c12c6b0 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl @@ -0,0 +1,88 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) + +A version of the classical Kelvin-Helmholtz instability based on +- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021) + A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations + of the Euler Equations + [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017) +""" +function initial_condition_kelvin_helmholtz_instability(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # typical resolution 128^2, 256^2 + # domain size is [-1,+1]^2 + slope = 15 + amplitude = 0.02 + B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) + rho = 0.5 + 0.75 * B + v1 = 0.5 * (B - 1) + v2 = 0.1 * sin(2 * pi * x[1]) + p = 1.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_kelvin_helmholtz_instability + +surface_flux = flux_lax_friedrichs +polydeg = 3 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +#indicator_sc = IndicatorHennemannGassner(equations, basis, +# alpha_max = 0.000, +# alpha_min = 0.0000, +# alpha_smooth = true, +# variable = density_pressure) +# +#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; +# volume_flux_dg = volume_flux, +# volume_flux_fv = surface_flux) +volume_integral = VolumeIntegralWeakFormProjection() +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.7) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval=1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + #save_solution, + 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); +summary_callback() # print the timer summary diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 1ea05fbe97a..7b8208039b1 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -235,7 +235,7 @@ function calc_volume_integral!(du, u, volume_integral::VolumeIntegralWeakFormProjection, dg::DGSEM, cache) # prepare local storage for projection - @unpack interpolate_N_to_M, project_M_to_N = dg.basis + @unpack interpolate_N_to_M, project_M_to_N, filter_modal_to_N = dg.basis nnodes_,nnodes_projection = size(project_M_to_N) nVars = nvariables(equations) RealT = real(dg) From fdc414ee410d6a6d92685c592a8cc662f9a35c31 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 18:01:07 +0100 Subject: [PATCH 3/7] current state --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 4 +- src/callbacks_step/analysis_dg2d.jl | 9 ++- src/solvers/dgsem_tree/dg_2d.jl | 73 +++++++++++++++++------ 3 files changed, 66 insertions(+), 20 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index c32aae817a7..2b273091b6a 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -46,11 +46,11 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 1.0) +stepsize_callback = StepsizeCallback(cfl = 0.1) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + #save_solution, stepsize_callback) ############################################################################### diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index de6b9a2a4a6..5249ba74ea4 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -239,13 +239,20 @@ function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) + u_original = similar(u) + u_original .= u + calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ - integrate_via_indices(u, mesh, equations, dg, cache, + result = integrate_via_indices(u, mesh, equations, dg, cache, du) do u, i, j, element, equations, dg, du u_node = get_node_vars(u, equations, dg, i, j, element) du_node = get_node_vars(du, equations, dg, i, j, element) dot(cons2entropy(u_node, equations), du_node) end + + u .= u_original + return result end function analyze(::Val{:l2_divb}, du, u, t, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 7b8208039b1..5776a2f4d5f 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -113,6 +113,12 @@ function rhs!(du, u, t, mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} + # Copy u for safekeeping + u_original = similar(u) + u_original .= u + + calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) @@ -175,9 +181,55 @@ function rhs!(du, u, t, calc_sources!(du, u, t, source_terms, equations, dg, cache) end + u .= u_original + return nothing end +function calc_entropy_projection!(u_projected, u, mesh, equations, dg, cache) + # prepare local storage for projection + @unpack interpolate_N_to_M, project_M_to_N = dg.basis + nnodes_,nnodes_projection = size(project_M_to_N) + nVars = nvariables(equations) + RealT = real(dg) + u_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + u_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + w_M = zeros(RealT, nVars, nnodes_projection, nnodes_projection) + + tmp_MxN = zeros(RealT, nVars, nnodes_projection, nnodes_) + tmp_NxM = zeros(RealT, nVars, nnodes_, nnodes_projection) + + @threaded for element in eachelement(dg, cache) + # get element u_N + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + for v in eachvariable(equations) + u_N[v,i,j] = u_node[v] + end + end + # bring elemtn u_N to grid (M+1)x(M+1) + multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) + + # compute nodal values of entropy variables w on the M grid + for j in 1:nnodes_projection, i in 1:nnodes_projection + u_cons = get_node_vars(u_M, equations, dg, i, j) + w_ij = cons2entropy(u_cons,equations) + set_node_vars!(w_M,w_ij,equations,dg,i,j) + end + + # compute projection of w with M values down to N + multiply_dimensionwise!(w_N, project_M_to_N ,w_M, tmp_NxM) + + # compute nodal values of conservative variables from the projected entropy variables + for j in eachnode(dg), i in eachnode(dg) + w_ij = get_node_vars(w_N, equations, dg, i, j) + u_cons = entropy2cons(w_ij, equations) + set_node_vars!(u_projected, u_cons, equations, dg, i, j, element) + end + end +end + function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, @@ -257,26 +309,14 @@ function calc_volume_integral!(du, u, # get element u_N for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, j, element) + w_ij = cons2entropy(u_node, equations) for v in eachvariable(equations) - u_N[v,i,j] = u_node[v] + w_N[v,i,j] = w_ij[v] end end # bring elemtn u_N to grid (M+1)x(M+1) - multiply_dimensionwise!(u_M,interpolate_N_to_M,u_N,tmp_MxN) - - # compute nodal values of entropy variables w on the M grid - for j in 1:nnodes_projection, i in 1:nnodes_projection - u_cons = get_node_vars(u_M, equations, dg, i, j) - w_ij = cons2entropy(u_cons,equations) - set_node_vars!(w_M_raw,w_ij,equations,dg,i,j) - end - # compute projection of w with M values down to N - multiply_dimensionwise!(w_M,filter_modal_to_N,w_M_raw,tmp_MxM) - - #multiply_dimensionwise!(w_N,project_M_to_N,w_M) - #multiply_dimensionwise!(w_M,interpolate_N_to_M,w_N) - - + multiply_dimensionwise!(w_M, interpolate_N_to_M, w_N, tmp_MxN) + # compute nodal values of conservative f,g on the M grid for j in 1:nnodes_projection, i in 1:nnodes_projection w_ij = get_node_vars(w_M, equations, dg, i, j) @@ -290,7 +330,6 @@ function calc_volume_integral!(du, u, multiply_dimensionwise!(f_N,project_M_to_N,f_M,tmp_NxM) multiply_dimensionwise!(g_N,project_M_to_N,g_M,tmp_NxM) - weak_form_kernel_projection!(du, u,f_N, g_N, element, mesh, nonconservative_terms, equations, dg, cache) From 559dfa296a930f1719f0c7a3b8b7d17398213b0b Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 19:30:34 +0100 Subject: [PATCH 4/7] Add modal cutoff filter --- examples/tree_2d_dgsem/elixir_euler_ec.jl | 7 +-- ...kelvin_helmholtz_instability_projection.jl | 11 +++-- src/callbacks_step/analysis_dg2d.jl | 8 +-- src/solvers/dgsem/basis_lobatto_legendre.jl | 13 +++-- src/solvers/dgsem_tree/dg_2d.jl | 49 ++++++++++++++++++- 5 files changed, 72 insertions(+), 16 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index 2b273091b6a..5e98f1ce785 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -13,9 +13,10 @@ initial_condition = initial_condition_weak_blast_wave # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 1 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 4 * polydeg) -volume_integral = VolumeIntegralWeakFormProjection() +polydeg = 6 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 1) +# volume_integral = VolumeIntegralWeakFormProjection() +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-2.0, -2.0) diff --git a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl index 9232c12c6b0..77c9572f717 100644 --- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl +++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_projection.jl @@ -33,8 +33,9 @@ end initial_condition = initial_condition_kelvin_helmholtz_instability surface_flux = flux_lax_friedrichs -polydeg = 3 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +# polydeg = 3 +polydeg = 10 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) #indicator_sc = IndicatorHennemannGassner(equations, basis, # alpha_max = 0.000, # alpha_min = 0.0000, @@ -44,7 +45,8 @@ basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) #volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; # volume_flux_dg = volume_flux, # volume_flux_fv = surface_flux) -volume_integral = VolumeIntegralWeakFormProjection() +# volume_integral = VolumeIntegralWeakFormProjection() +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (-1.0, -1.0) @@ -72,7 +74,8 @@ save_solution = SaveSolutionCallback(interval=1000, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.5) +# stepsize_callback = StepsizeCallback(cfl = 0.5) +stepsize_callback = StepsizeCallback(cfl = 1.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 5249ba74ea4..53f586938b1 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -239,9 +239,9 @@ function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) - u_original = similar(u) - u_original .= u - calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # u_original = similar(u) + # u_original .= u + # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ result = integrate_via_indices(u, mesh, equations, dg, cache, @@ -251,7 +251,7 @@ function analyze(::typeof(entropy_timederivative), du, u, t, dot(cons2entropy(u_node, equations), du_node) end - u .= u_original + # u .= u_original return result end diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index aadc60f18eb..1bd3d9efa18 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -39,9 +39,10 @@ struct LobattoLegendreBasis{RealT <: Real, NNODES, 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 LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) +function LobattoLegendreBasis(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 @@ -84,6 +85,11 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ 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_lobatto_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 LobattoLegendreBasis{RealT, nnodes_, typeof(nodes), typeof(inverse_vandermonde_legendre), typeof(boundary_interpolation), @@ -100,10 +106,11 @@ function LobattoLegendreBasis(RealT, polydeg::Integer; polydeg_projection::Integ derivative_dhat, interpolate_N_to_M, project_M_to_N, - filter_modal_to_N) + filter_modal_to_N, + filter_modal_to_cutoff) end -LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection) +LobattoLegendreBasis(polydeg::Integer; polydeg_projection::Integer = 2 * polydeg, polydeg_cutoff::Integer = div(polydeg + 1, 2) - 1) = LobattoLegendreBasis(Float64, polydeg; polydeg_projection, polydeg_cutoff) function Base.show(io::IO, basis::LobattoLegendreBasis) @nospecialize basis # reduce precompilation time diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 5776a2f4d5f..899fda82747 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -116,8 +116,20 @@ function rhs!(du, u, t, # Copy u for safekeeping u_original = similar(u) u_original .= u + u_filter_cons = similar(u) + u_filter_prim = similar(u) - calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) + # calc_filter!(u, u_original, cons2entropy, entropy2cons, mesh, equations, dg, cache) + # calc_filter!(u, u_original, cons2cons, cons2cons, mesh, equations, dg, cache) + # calc_filter!(u, u_original, cons2prim, prim2cons, mesh, equations, dg, cache) + calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) + calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) + + u .= u_filter_cons + @. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] + @. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) + # @autoinfiltrate # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) @@ -181,7 +193,7 @@ function rhs!(du, u, t, calc_sources!(du, u, t, source_terms, equations, dg, cache) end - u .= u_original + # u .= u_original return nothing end @@ -230,6 +242,39 @@ function calc_entropy_projection!(u_projected, u, mesh, equations, dg, cache) end end +function calc_filter!(u_filtered, u, cons2filter, filter2cons, mesh, equations, dg, cache) + # prepare local storage for projection + @unpack filter_modal_to_cutoff = dg.basis + nnodes_ = nnodes(dg) + nVars = nvariables(equations) + RealT = real(dg) + w_N = zeros(RealT, nVars, nnodes_, nnodes_) + w_N_filtered = zeros(RealT, nVars, nnodes_, nnodes_) + + tmp_NxN = zeros(RealT, nVars, nnodes_, nnodes_) + + @threaded for element in eachelement(dg, cache) + # convert u to entropy variables + for j in eachnode(dg), i in eachnode(dg) + u_cons = get_node_vars(u, equations, dg, i, j, element) + w_ij = cons2filter(u_cons, equations) + for v in eachvariable(equations) + w_N[v,i,j] = w_ij[v] + end + end + + # filter entropy variables + multiply_dimensionwise!(w_N_filtered, filter_modal_to_cutoff, w_N, tmp_NxN) + + # compute nodal values of conservative variables from the projected entropy variables + for j in eachnode(dg), i in eachnode(dg) + w_ij = get_node_vars(w_N_filtered, equations, dg, i, j) + u_cons = filter2cons(w_ij, equations) + set_node_vars!(u_filtered, u_cons, equations, dg, i, j, element) + end + end +end + function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, From 7df151648d1819f9511cf34b44ff94c428258f46 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Tue, 28 May 2024 23:48:41 +0100 Subject: [PATCH 5/7] Beginning of using JinXin + filtering --- examples/tree_2d_dgsem/elixir_jin_xin_euler.jl | 8 +++++--- src/solvers/dgsem_tree/dg_2d.jl | 8 ++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl index 6ee8d11df98..e06318a500e 100644 --- a/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl +++ b/examples/tree_2d_dgsem/elixir_jin_xin_euler.jl @@ -29,10 +29,12 @@ end #initial_condition = initial_condition_constant initial_condition = Trixi.InitialConditionJinXin(initial_condition_kelvin_helmholtz_instability) #initial_condition = Trixi.InitialConditionJinXin(initial_condition_density_wave) -polydeg = 3 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg) +polydeg = 5 +polydeg_cutoff = 3 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 1 * polydeg, polydeg_cutoff = polydeg_cutoff) # solver = DGSEM(basis, Trixi.flux_upwind) -solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection()) +#solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakFormProjection()) +solver = DGSEM(basis, Trixi.flux_upwind, VolumeIntegralWeakForm()) #solver = DGSEM(polydeg = 3, surface_flux = Trixi.flux_upwind) #surface_flux = Trixi.flux_upwind diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 899fda82747..1bb0ce91df0 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -117,18 +117,18 @@ function rhs!(du, u, t, u_original = similar(u) u_original .= u u_filter_cons = similar(u) - u_filter_prim = similar(u) + # u_filter_prim = similar(u) # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) # calc_filter!(u, u_original, cons2entropy, entropy2cons, mesh, equations, dg, cache) # calc_filter!(u, u_original, cons2cons, cons2cons, mesh, equations, dg, cache) # calc_filter!(u, u_original, cons2prim, prim2cons, mesh, equations, dg, cache) calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) - calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) + #calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) u .= u_filter_cons - @. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] - @. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) + #@. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] + #@. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) # @autoinfiltrate # Reset du From aabec9c14fac346165be71311d9b7a12808e94f3 Mon Sep 17 00:00:00 2001 From: Gregor Gassner Date: Wed, 29 May 2024 12:37:19 +0100 Subject: [PATCH 6/7] filter testing stuff --- .../elixir_euler_density_wave.jl | 39 ++++++++++++++++--- examples/tree_2d_dgsem/elixir_euler_ec.jl | 18 +++++++-- src/solvers/dgsem_tree/dg_2d.jl | 18 ++++----- 3 files changed, 57 insertions(+), 18 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_euler_density_wave.jl b/examples/tree_2d_dgsem/elixir_euler_density_wave.jl index 0f9e232fa14..e7736252dcf 100644 --- a/examples/tree_2d_dgsem/elixir_euler_density_wave.jl +++ b/examples/tree_2d_dgsem/elixir_euler_density_wave.jl @@ -1,15 +1,32 @@ using OrdinaryDiffEq using Trixi +using LinearAlgebra ############################################################################### # semidiscretization of the compressible Euler equations gamma = 1.4 equations = CompressibleEulerEquations2D(gamma) -initial_condition = initial_condition_density_wave - -solver = DGSEM(polydeg = 5, surface_flux = flux_central) +function initial_condition_density_wave2(x, t, equations::CompressibleEulerEquations2D) + v1 = 0.1 + v2 = 0.2 + rho = 1 + 0.98 * sinpi(2 * (x[1] + x[2] - t * (v1 + v2))) + rho_v1 = rho * v1 + rho_v2 = rho * v2 + p = 20 + rho_e = p / (equations.gamma - 1) + 1 / 2 * rho * (v1^2 + v2^2) + return SVector(rho, rho_v1, rho_v2, rho_e) +end +initial_condition = initial_condition_density_wave2 + +surface_flux = flux_central +polydeg = 17 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 5) +volume_integral = VolumeIntegralWeakForm() +solver = DGSEM(basis, surface_flux, volume_integral) + +# solver = DGSEM(polydeg = 5, surface_flux = flux_central) coordinates_min = (-1.0, -1.0) coordinates_max = (1.0, 1.0) @@ -17,8 +34,20 @@ mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2, n_cells_max = 30_000) +@info "Create semi..." semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +@info "Compute Jacobian..." +J = jacobian_ad_forward(semi) + +@info "Compute eigenvalues..." +λ = eigvals(J) + +@info "max real part" maximum(real.(λ)) +# @info "Plot spectrum..." +# scatter(real.(λ), imag.(λ), label="central flux") +wololo + ############################################################################### # ODE solvers, callbacks etc. @@ -27,7 +56,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 100 * 20 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) @@ -41,7 +70,7 @@ stepsize_callback = StepsizeCallback(cfl = 1.6) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + # save_solution, stepsize_callback) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_euler_ec.jl b/examples/tree_2d_dgsem/elixir_euler_ec.jl index 5e98f1ce785..ced8a5221a2 100644 --- a/examples/tree_2d_dgsem/elixir_euler_ec.jl +++ b/examples/tree_2d_dgsem/elixir_euler_ec.jl @@ -1,20 +1,30 @@ using OrdinaryDiffEq using Trixi - +using Random: seed! ############################################################################### # semidiscretization of the compressible Euler equations equations = CompressibleEulerEquations2D(1.4) +seed!(1) +function initial_condition_random_field(x, t, equations::CompressibleEulerEquations2D) +amplitude = 1.5 +rho = 2 + amplitude * rand() +v1 = -3.1 + amplitude * rand() +v2 = 1.3 + amplitude * rand() +p = 7.54 + amplitude * rand() +return prim2cons(SVector(rho, v1, v2, p), equations) +end initial_condition = initial_condition_weak_blast_wave +# initial_condition = initial_condition_random_field #volume_flux = flux_ranocha #solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha, # volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) # surface_flux = flux_ranocha -polydeg = 6 -basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 1) +polydeg = 11 +basis = LobattoLegendreBasis(polydeg; polydeg_projection = 2 * polydeg, polydeg_cutoff = 3) # volume_integral = VolumeIntegralWeakFormProjection() volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) @@ -47,7 +57,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.1) +stepsize_callback = StepsizeCallback(cfl = 0.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 1bb0ce91df0..50b8e7eac35 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -114,19 +114,19 @@ function rhs!(du, u, t, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Copy u for safekeeping - u_original = similar(u) - u_original .= u + #u_original = similar(u) + #u_original .= u u_filter_cons = similar(u) - # u_filter_prim = similar(u) + u_filter_prim = similar(u) # calc_entropy_projection!(u, u_original, mesh, equations, dg, cache) - # calc_filter!(u, u_original, cons2entropy, entropy2cons, mesh, equations, dg, cache) - # calc_filter!(u, u_original, cons2cons, cons2cons, mesh, equations, dg, cache) - # calc_filter!(u, u_original, cons2prim, prim2cons, mesh, equations, dg, cache) - calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) + #calc_filter!(u, u, cons2entropy, entropy2cons, mesh, equations, dg, cache) + # calc_filter!(u, u, cons2cons, cons2cons, mesh, equations, dg, cache) + calc_filter!(u, u, cons2prim, prim2cons, mesh, equations, dg, cache) + #calc_filter!(u_filter_cons, u, cons2cons, cons2cons, mesh, equations, dg, cache) #calc_filter!(u_filter_prim, u, cons2prim, prim2cons, mesh, equations, dg, cache) - u .= u_filter_cons + #u .= u_filter_cons #@. u_filter_prim[4, ..] = u_filter_prim[4, ..] - 0.5 * (u_filter_prim[2, ..]^2 + u_filter_prim[3, ..]^2) / u_filter_prim[1, ..] #@. u[4, ..] = 0.5 * (u_filter_cons[2, ..]^2 + u_filter_cons[3, ..]^2) / u_filter_cons[1, ..] + u_filter_prim[4, ..]# / (equations.gamma - 1.0) # @autoinfiltrate @@ -247,7 +247,7 @@ function calc_filter!(u_filtered, u, cons2filter, filter2cons, mesh, equations, @unpack filter_modal_to_cutoff = dg.basis nnodes_ = nnodes(dg) nVars = nvariables(equations) - RealT = real(dg) + RealT = eltype(u) w_N = zeros(RealT, nVars, nnodes_, nnodes_) w_N_filtered = zeros(RealT, nVars, nnodes_, nnodes_) From 306d47eb4f07f064c554620ac5773afcbf3f9079 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Wed, 29 May 2024 15:15:28 +0100 Subject: [PATCH 7/7] 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