From 1d7b2561ff656cb9b8c826ee8b7ee389a3370917 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 9 Oct 2024 16:35:07 +0200 Subject: [PATCH 01/13] Level Info `TreeMesh` without AMR (#2104) * Level Info TreeMesh without AMR * comment * Update src/callbacks_step/analysis.jl * Update src/callbacks_step/analysis.jl --- src/callbacks_step/analysis.jl | 59 ++++++++++++++++++++++++---------- 1 file changed, 42 insertions(+), 17 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 06110d08d28..b6c80f47867 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -314,8 +314,8 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi) mpi_println(" #elements: " * @sprintf("% 14d", nelementsglobal(mesh, solver, cache))) - # Level information (only show for AMR) - print_amr_information(integrator.opts.callback, mesh, solver, cache) + # Level information (only for AMR and/or non-uniform `TreeMesh`es) + print_level_information(integrator.opts.callback, mesh, solver, cache) mpi_println() # Open file for appending and store time step and time information @@ -489,21 +489,7 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi) return nothing end -# Print level information only if AMR is enabled -function print_amr_information(callbacks, mesh, solver, cache) - - # Return early if there is nothing to print - uses_amr(callbacks) || return nothing - - # Get global minimum and maximum level from the AMRController - min_level = max_level = 0 - for cb in callbacks.discrete_callbacks - if cb.affect! isa AMRCallback - min_level = cb.affect!.controller.base_level - max_level = cb.affect!.controller.max_level - end - end - +function print_level_information(mesh, solver, cache, min_level, max_level) # Get local element count per level elements_per_level = get_elements_per_level(min_level, max_level, mesh, solver, cache) @@ -522,6 +508,45 @@ function print_amr_information(callbacks, mesh, solver, cache) return nothing end +function print_level_information(callbacks, mesh::TreeMesh, solver, cache) + if uses_amr(callbacks) + # Get global minimum and maximum level from the AMRController + min_level = max_level = 0 + for cb in callbacks.discrete_callbacks + if cb.affect! isa AMRCallback + min_level = cb.affect!.controller.base_level + max_level = cb.affect!.controller.max_level + end + end + print_level_information(mesh, solver, cache, min_level, max_level) + # Special check for `TreeMesh`es without AMR. + # These meshes may still be non-uniform due to `refinement_patches`, see + # `refine_box!`, `coarsen_box!`, and `refine_sphere!`. + elseif minimum_level(mesh.tree) != maximum_level(mesh.tree) + min_level = minimum_level(mesh.tree) + max_level = maximum_level(mesh.tree) + print_level_information(mesh, solver, cache, min_level, max_level) + else # Uniform mesh + return nothing + end +end + +function print_level_information(callbacks, mesh, solver, cache) + if uses_amr(callbacks) + # Get global minimum and maximum level from the AMRController + min_level = max_level = 0 + for cb in callbacks.discrete_callbacks + if cb.affect! isa AMRCallback + min_level = cb.affect!.controller.base_level + max_level = cb.affect!.controller.max_level + end + end + print_level_information(mesh, solver, cache, min_level, max_level) + else # Uniform mesh + return nothing + end +end + function get_elements_per_level(min_level, max_level, mesh::P4estMesh, solver, cache) elements_per_level = zeros(P4EST_MAXLEVEL + 1) From 1494aa2144995b194985251512ddfe116b78e45e Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Thu, 10 Oct 2024 10:39:29 +0200 Subject: [PATCH 02/13] Remove `normal_direction_ll` for nonconservative terms (#2062) * remove_normal_direction_ll * update testset name StructuredMesh3D * fix BoundaryConditionDoNothing * adjust docstrings * add comment to docstring * update test values * add news item --------- Co-authored-by: Daniel Doehring Co-authored-by: Hendrik Ranocha --- NEWS.md | 16 +++- .../structured_2d_dgsem/elixir_mhd_coupled.jl | 8 -- src/basic_types.jl | 4 +- src/equations/ideal_glm_mhd_2d.jl | 28 +++---- src/equations/ideal_glm_mhd_3d.jl | 27 ++---- src/equations/shallow_water_2d.jl | 30 +++---- src/solvers/dgmulti/dg.jl | 20 ++--- src/solvers/dgmulti/flux_differencing.jl | 10 +-- src/solvers/dgsem_p4est/dg_2d.jl | 18 ++-- src/solvers/dgsem_p4est/dg_3d.jl | 13 +-- src/solvers/dgsem_structured/dg_2d.jl | 30 +++---- src/solvers/dgsem_structured/dg_3d.jl | 40 ++++----- src/solvers/dgsem_unstructured/dg_2d.jl | 17 +--- test/test_dgmulti_2d.jl | 28 +++---- test/test_p4est_2d.jl | 15 ++-- test/test_p4est_3d.jl | 20 ++--- test/test_structured_2d.jl | 62 +++++++------- test/test_structured_3d.jl | 83 ++++++++++--------- test/test_t8code_2d.jl | 13 +-- test/test_threaded.jl | 20 ++--- test/test_type.jl | 25 ++---- test/test_unstructured_2d.jl | 41 ++++----- 22 files changed, 250 insertions(+), 318 deletions(-) diff --git a/NEWS.md b/NEWS.md index 3d64f7b78f7..cba0d3f86d8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,13 +6,27 @@ for human readability. ## Changes when updating to v0.9 from v0.8.x +#### Added + +- Boundary conditions are now supported on nonconservative terms ([#2062]). + #### Changed - We removed the first argument `semi` corresponding to a `Semidiscretization` from the `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]). - The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. ([#2069]) + The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. + ([#2069]) - In functions `rhs!`, `rhs_parabolic!` we removed the unused argument `initial_condition`. ([#2037]) Users should not be affected by this. +- Nonconservative terms depend only on `normal_direction_average` instead of both + `normal_direction_average` and `normal_direction_ll`, such that the function signature is now + identical with conservative fluxes. This required a change of the `normal_direction` in + `flux_nonconservative_powell` ([#2062]). + +#### Deprecated + +#### Removed + ## Changes in the v0.8 lifecycle diff --git a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl index d3aa4ecf582..de274248b45 100644 --- a/examples/structured_2d_dgsem/elixir_mhd_coupled.jl +++ b/examples/structured_2d_dgsem/elixir_mhd_coupled.jl @@ -31,14 +31,6 @@ equations = IdealGlmMhdEquations2D(gamma) cells_per_dimension = (32, 64) -# Extend the definition of the non-conservative Powell flux functions. -import Trixi.flux_nonconservative_powell -function flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll::AbstractVector, - equations::IdealGlmMhdEquations2D) - flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll, normal_direction_ll, - equations) -end volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) solver = DGSEM(polydeg = 3, surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell), diff --git a/src/basic_types.jl b/src/basic_types.jl index ee479a62039..8a2430c1a85 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -71,14 +71,14 @@ struct BoundaryConditionDoNothing end orientation_or_normal_direction, direction::Integer, x, t, surface_flux, equations) - return flux(u_inner, orientation_or_normal_direction, equations) + return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations) end # This version can be called by hyperbolic solvers on unstructured, curved meshes @inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector, x, t, surface_flux, equations) - return flux(u_inner, outward_direction, equations) + return surface_flux(u_inner, u_inner, outward_direction, equations) end # This version can be called by parabolic solvers diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 7a6d19facd1..a43c6b0e619 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -196,18 +196,15 @@ end flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations2D) flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::IdealGlmMhdEquations2D) Non-symmetric two-point flux discretizing the nonconservative (source) term of Powell and the Galilean nonconservative term associated with the GLM multiplier of the [`IdealGlmMhdEquations2D`](@ref). -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. +On curvilinear meshes, the implementation differs from the reference since we use the averaged +normal direction for the metrics dealiasing. ## References - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs, @@ -254,8 +251,7 @@ terms. end @inline function flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::IdealGlmMhdEquations2D) rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr @@ -265,14 +261,9 @@ end v3_ll = rho_v3_ll / rho_ll v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll - # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector - # at the same node location) while `B_dot_n_rr` uses the averaged normal - # direction. The reason for this is that `v_dot_n_ll` depends only on the left - # state and multiplies some gradient while `B_dot_n_rr` is used to compute - # the divergence of B. - v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] - B_dot_n_rr = B1_rr * normal_direction_average[1] + - B2_rr * normal_direction_average[2] + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + B_dot_n_rr = B1_rr * normal_direction[1] + + B2_rr * normal_direction[2] # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) @@ -300,8 +291,9 @@ of the [`IdealGlmMhdEquations2D`](@ref). This implementation uses a non-conservative term that can be written as the product of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm -et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different -results on non-conforming meshes(!). +et al. [`flux_nonconservative_powell`](@ref) for conforming meshes but it yields different +results on non-conforming meshes(!). On curvilinear meshes this formulation applies the +local normal direction compared to the averaged one used in [`flux_nonconservative_powell`](@ref). The two other flux functions with the same name return either the local or symmetric portion of the non-conservative flux based on the type of the diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index e922a2e6fd6..9bb05cb3f89 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -224,18 +224,15 @@ end flux_nonconservative_powell(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations3D) flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::IdealGlmMhdEquations3D) Non-symmetric two-point flux discretizing the nonconservative (source) term of Powell and the Galilean nonconservative term associated with the GLM multiplier of the [`IdealGlmMhdEquations3D`](@ref). -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. +On curvilinear meshes, the implementation differs from the reference since we use the averaged +normal direction for the metrics dealiasing. ## References - Marvin Bohm, Andrew R.Winters, Gregor J. Gassner, Dominik Derigs, @@ -292,8 +289,7 @@ terms. end @inline function flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::IdealGlmMhdEquations3D) rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr @@ -303,16 +299,11 @@ end v3_ll = rho_v3_ll / rho_ll v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll - # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector - # at the same node location) while `B_dot_n_rr` uses the averaged normal - # direction. The reason for this is that `v_dot_n_ll` depends only on the left - # state and multiplies some gradient while `B_dot_n_rr` is used to compute - # the divergence of B. - v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] + - v3_ll * normal_direction_ll[3] - B_dot_n_rr = B1_rr * normal_direction_average[1] + - B2_rr * normal_direction_average[2] + - B3_rr * normal_direction_average[3] + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + B_dot_n_rr = B1_rr * normal_direction[1] + + B2_rr * normal_direction[2] + + B3_rr * normal_direction[3] # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3}) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 4ecaf3b6e14..30eca095aa4 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -267,8 +267,7 @@ end flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::ShallowWaterEquations2D) Non-symmetric two-point volume flux discretizing the nonconservative (source) term @@ -303,8 +302,7 @@ Further details are available in the papers: end @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::ShallowWaterEquations2D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) @@ -312,8 +310,8 @@ end # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) return SVector(0, - normal_direction_average[1] * equations.gravity * h_ll * b_jump, - normal_direction_average[2] * equations.gravity * h_ll * b_jump, + normal_direction[1] * equations.gravity * h_ll * b_jump, + normal_direction[2] * equations.gravity * h_ll * b_jump, 0) end @@ -321,8 +319,7 @@ end flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) flux_nonconservative_fjordholm_etal(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::ShallowWaterEquations2D) Non-symmetric two-point surface flux discretizing the nonconservative (source) term of @@ -365,8 +362,7 @@ and for curvilinear 2D case in the paper: end @inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::ShallowWaterEquations2D) # Pull the necessary left and right state information h_ll, _, _, b_ll = u_ll @@ -376,8 +372,8 @@ end b_jump = b_rr - b_ll # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0) - f2 = normal_direction_average[1] * equations.gravity * h_average * b_jump - f3 = normal_direction_average[2] * equations.gravity * h_average * b_jump + f2 = normal_direction[1] * equations.gravity * h_average * b_jump + f3 = normal_direction[2] * equations.gravity * h_average * b_jump # First and last equations do not have a nonconservative flux f1 = f4 = 0 @@ -424,8 +420,7 @@ end flux_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) flux_nonconservative_audusse_etal(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::ShallowWaterEquations2D) Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. @@ -467,8 +462,7 @@ Further details for the hydrostatic reconstruction and its motivation can be fou end @inline function flux_nonconservative_audusse_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, + normal_direction::AbstractVector, equations::ShallowWaterEquations2D) # Pull the water height and bottom topography on the left h_ll, _, _, b_ll = u_ll @@ -479,8 +473,8 @@ end # Copy the reconstructed water height for easier to read code h_ll_star = u_ll_star[1] - f2 = normal_direction_average[1] * equations.gravity * (h_ll^2 - h_ll_star^2) - f3 = normal_direction_average[2] * equations.gravity * (h_ll^2 - h_ll_star^2) + f2 = normal_direction[1] * equations.gravity * (h_ll^2 - h_ll_star^2) + f3 = normal_direction[2] * equations.gravity * (h_ll^2 - h_ll_star^2) # First and last equations do not have a nonconservative flux f1 = f4 = 0 diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index be7dd3241d3..2d588c5c79d 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -407,11 +407,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm, # Two notes on the use of `flux_nonconservative`: # 1. In contrast to other mesh types, only one nonconservative part needs to be # computed since we loop over the elements, not the unique interfaces. - # 2. In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. However, - # both are the same at watertight interfaces, so we pass `normal` twice. - nonconservative_part = flux_nonconservative(uM, uP, normal, normal, - equations) + nonconservative_part = flux_nonconservative(uM, uP, normal, equations) # The factor 0.5 is necessary for the nonconservative fluxes based on the # interpretation of global SBP operators. flux_face_values[idM] = (conservative_part + 0.5 * nonconservative_part) * @@ -556,14 +552,12 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, surface_flux, equations) # Compute pointwise nonconservative numerical flux at the boundary. - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, there is only one `face_normal` at boundaries, which we pass in twice. - # Note: This does not set any type of boundary condition for the nonconservative term - noncons_flux_at_face_node = nonconservative_flux(u_face_values[i, f], - u_face_values[i, f], - face_normal, face_normal, - equations) + noncons_flux_at_face_node = boundary_condition(u_face_values[i, f], + face_normal, + face_coordinates, + t, + nonconservative_flux, + equations) flux_face_values[i, f] = (cons_flux_at_face_node + 0.5 * noncons_flux_at_face_node) * Jf[i, f] diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl index a2edf4ceca8..88f06607019 100644 --- a/src/solvers/dgmulti/flux_differencing.jl +++ b/src/solvers/dgmulti/flux_differencing.jl @@ -76,10 +76,7 @@ end du_i = du[i] for j in col_ids u_j = u[j] - # The `normal_direction::AbstractVector` has to be passed in twice. - # This is because on curved meshes, nonconservative fluxes are - # evaluated using both the normal and its average at interfaces. - f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations) + f_ij = volume_flux(u_i, u_j, normal_direction, equations) du_i = du_i + 2 * A[i, j] * f_ij end du[i] = du_i @@ -176,11 +173,8 @@ end for id in nzrange(A_base, i) A_ij = vals[id] j = rows[id] - # The `normal_direction::AbstractVector` has to be passed in twice. - # This is because on curved meshes, nonconservative fluxes are - # evaluated using both the normal and its average at interfaces. u_j = u[j] - f_ij = volume_flux(u_i, u_j, normal_direction, normal_direction, equations) + f_ij = volume_flux(u_i, u_j, normal_direction, equations) du_i = du_i + 2 * A_ij * f_ij end du[i] = du_i diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 17b9af04467..3c868289181 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -223,14 +223,8 @@ end flux_ = surface_flux(u_ll, u_rr, normal_direction, equations) # Compute both nonconservative fluxes - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, both are the same at watertight interfaces, so we pass the - # `normal_direction` twice. - noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) - noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations) + noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations) # Store the flux with nonconservative terms on the primary and secondary elements for v in eachvariable(equations) @@ -369,9 +363,8 @@ end flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations) # Compute pointwise nonconservative numerical flux at the boundary. - # Note: This does not set any type of boundary condition for the nonconservative term - noncons_ = nonconservative_flux(u_inner, u_inner, normal_direction, - normal_direction, equations) + noncons_ = boundary_condition(u_inner, normal_direction, x, t, nonconservative_flux, + equations) # Copy flux to element storage in the correct orientation for v in eachvariable(equations) @@ -554,8 +547,7 @@ end # The nonconservative flux is scaled by a factor of 0.5 based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, - equations) + noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations) flux_plus_noncons = flux + 0.5f0 * noncons diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index ece4840b74b..5aabbf7ac60 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -289,14 +289,8 @@ end flux_ = surface_flux(u_ll, u_rr, normal_direction, equations) # Compute both nonconservative fluxes - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, both are the same at watertight interfaces, so we pass the - # `normal_direction` twice. - noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) - noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations) + noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations) # Store the flux with nonconservative terms on the primary and secondary elements for v in eachvariable(equations) @@ -633,8 +627,7 @@ end # Compute nonconservative flux and add it to the flux scaled by a factor of 0.5 based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, - equations) + noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations) flux_plus_noncons = flux + 0.5f0 * noncons # Copy to buffer diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 39bece53ca9..cb4e75149ea 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -195,7 +195,7 @@ end element) Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # Compute the contravariant nonconservative flux. - fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg, + fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_avg, equations) integral_contribution = integral_contribution + derivative_split[i, ii] * fluxtilde1 @@ -210,7 +210,7 @@ end Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector - fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg, + fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_avg, equations) integral_contribution = integral_contribution + derivative_split[j, jj] * fluxtilde2 @@ -337,11 +337,11 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde1_L = ftilde1 + - 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_ll, u_rr, normal_direction, equations) ftilde1_R = ftilde1 + - 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_rr, u_ll, normal_direction, equations) set_node_vars!(fstar1_L, ftilde1_L, equations, dg, i, j) set_node_vars!(fstar1_R, ftilde1_R, equations, dg, i, j) @@ -377,11 +377,11 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde2_L = ftilde2 + - 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_ll, u_rr, normal_direction, equations) ftilde2_R = ftilde2 + - 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_rr, u_ll, normal_direction, equations) set_node_vars!(fstar2_L, ftilde2_L, equations, dg, i, j) set_node_vars!(fstar2_R, ftilde2_R, equations, dg, i, j) @@ -527,18 +527,12 @@ end flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations) # Compute both nonconservative fluxes - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, both are the same at watertight interfaces, so we pass the - # `normal_direction` twice. # Scale with sign_jacobian to ensure that the normal_direction matches that # from the flux above noncons_left = sign_jacobian * - nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + nonconservative_flux(u_ll, u_rr, normal_direction, equations) noncons_right = sign_jacobian * - nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + nonconservative_flux(u_rr, u_ll, normal_direction, equations) for v in eachvariable(equations) # Note the factor 0.5 necessary for the nonconservative fluxes based on diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index 616a6b15131..aeae63183b0 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -228,7 +228,7 @@ end Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector - fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg, + fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_avg, equations) integral_contribution = integral_contribution + derivative_split[i, ii] * fluxtilde1 @@ -243,7 +243,7 @@ end Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector - fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg, + fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_avg, equations) integral_contribution = integral_contribution + derivative_split[j, jj] * fluxtilde2 @@ -258,7 +258,7 @@ end Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector - fluxtilde3 = nonconservative_flux(u_node, u_node_kk, Ja3_node, Ja3_avg, + fluxtilde3 = nonconservative_flux(u_node, u_node_kk, Ja3_avg, equations) integral_contribution = integral_contribution + derivative_split[k, kk] * fluxtilde3 @@ -411,11 +411,11 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde_L = ftilde + - 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_ll, u_rr, normal_direction, equations) ftilde_R = ftilde + - 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_rr, u_ll, normal_direction, equations) set_node_vars!(fstar1_L, ftilde_L, equations, dg, i, j, k) set_node_vars!(fstar1_R, ftilde_R, equations, dg, i, j, k) @@ -449,11 +449,11 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde_L = ftilde + - 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_ll, u_rr, normal_direction, equations) ftilde_R = ftilde + - 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_rr, u_ll, normal_direction, equations) set_node_vars!(fstar2_L, ftilde_L, equations, dg, i, j, k) set_node_vars!(fstar2_R, ftilde_R, equations, dg, i, j, k) @@ -487,11 +487,11 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde_L = ftilde + - 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_ll, u_rr, normal_direction, equations) ftilde_R = ftilde + - 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + 0.5f0 * + nonconservative_flux(u_rr, u_ll, normal_direction, equations) set_node_vars!(fstar3_L, ftilde_L, equations, dg, i, j, k) set_node_vars!(fstar3_R, ftilde_R, equations, dg, i, j, k) @@ -663,18 +663,12 @@ end flux = sign_jacobian * surface_flux(u_ll, u_rr, normal_direction, equations) # Compute both nonconservative fluxes - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, both are the same at watertight interfaces, so we pass the - # `normal_direction` twice. # Scale with sign_jacobian to ensure that the normal_direction matches that # from the flux above noncons_left = sign_jacobian * - nonconservative_flux(u_ll, u_rr, normal_direction, - normal_direction, equations) + nonconservative_flux(u_ll, u_rr, normal_direction, equations) noncons_right = sign_jacobian * - nonconservative_flux(u_rr, u_ll, normal_direction, - normal_direction, equations) + nonconservative_flux(u_rr, u_ll, normal_direction, equations) for v in eachvariable(equations) # Note the factor 0.5 necessary for the nonconservative fluxes based on diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 1f8c1ffa16a..48d4fe153c6 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -247,14 +247,10 @@ function calc_interface_flux!(surface_flux_values, flux = surface_flux(u_ll, u_rr, outward_direction, equations) # Compute both nonconservative fluxes - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, both are the same at watertight interfaces, so we pass the - # `outward_direction` twice. noncons_primary = nonconservative_flux(u_ll, u_rr, outward_direction, - outward_direction, equations) + equations) noncons_secondary = nonconservative_flux(u_rr, u_ll, outward_direction, - outward_direction, equations) + equations) # Copy flux to primary and secondary element storage # Note the sign change for the components in the secondary element! @@ -449,13 +445,8 @@ end flux = boundary_condition(u_inner, outward_direction, x, t, surface_flux, equations) # Compute pointwise nonconservative numerical flux at the boundary. - # In general, nonconservative fluxes can depend on both the contravariant - # vectors (normal direction) at the current node and the averaged ones. - # However, both are the same at watertight interfaces, so we pass the - # `outward_direction` twice. - # Note: This does not set any type of boundary condition for the nonconservative term - noncons_flux = nonconservative_flux(u_inner, u_inner, outward_direction, - outward_direction, equations) + noncons_flux = boundary_condition(u_inner, outward_direction, x, t, + nonconservative_flux, equations) for v in eachvariable(equations) # Note the factor 0.5 necessary for the nonconservative fluxes based on diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 932fb9bc958..2e41591d52c 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -788,26 +788,26 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_wall.jl"), cells_per_dimension=4, l2=[ - 0.0036019536614619687, - 0.001734097206958611, - 0.008375221008997178, + 0.0036019562526881602, + 0.0017340971255535853, + 0.00837522167692243, 0.0, - 0.028596796602124414, - 0.0018573693138866614, - 0.0020807798141551166, + 0.028596802654003512, + 0.0018573697892233679, + 0.0020807798940528956, 0.0, - 5.301188920230166e-5 + 5.301259762428258e-5 ], linf=[ - 0.01692601228199253, - 0.009369662298436778, - 0.04145169295835428, + 0.016925983823703028, + 0.009369659529710701, + 0.04145170727840005, 0.0, - 0.11569908670112738, - 0.00984964453299233, - 0.01141708032148614, + 0.1156990108418654, + 0.009849648257876749, + 0.011417088537145403, 0.0, - 0.0002992631411931389 + 0.0002992621756946904 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index cba79b58b70..46dfa542f34 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -617,17 +617,18 @@ end @trixi_testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.4552094211937344, 0.8918052934760807, 0.8324715234680768, + l2=[0.45520942119628094, 0.8918052934746296, 0.8324715234698744, 0.0, - 0.9801268321975978, 0.10475722739111007, - 0.15551326369033164, + 0.9801268321975156, 0.1047572273917542, + 0.15551326369065485, 0.0, - 2.0602990858239632e-5], - linf=[10.19421969147307, 18.254409357804683, 10.031954811332596, + 2.0604823850230503e-5], + linf=[10.194219681102396, 18.254409373691253, + 10.031954811617876, 0.0, - 19.646870938371492, 1.3938679692894465, 1.8725058401937984, + 19.646870952133547, 1.39386796733875, 1.8725058402095736, 0.0, - 0.0016201762010257296], + 0.001619787048808356], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 869b7554bf2..55f7ab76023 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -571,16 +571,16 @@ end @trixi_testset "elixir_mhd_shockcapturing_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"), - l2=[0.006297229188299052, 0.0064363477630573936, - 0.007109134822960387, 0.0065295379843073945, - 0.02061487028361094, 0.005561406556868266, - 0.007570747563219415, 0.005571060186624124, - 3.910359570546058e-6], - linf=[0.20904050617411984, 0.18630026905465372, - 0.23476537952044518, 0.19430178061639747, - 0.6858488631108304, 0.15169972134884624, - 0.22431157069631724, 0.16823638724229162, - 0.0005352202836463904], + l2=[0.006297229188267704, 0.006436347763092648, + 0.0071091348227321095, 0.00652953798427642, + 0.0206148702828057, 0.005561406556411695, + 0.007570747563696005, 0.005571060186513173, + 3.888176398720913e-6], + linf=[0.20904050630623572, 0.1863002690612441, + 0.2347653795205547, 0.19430178062881898, + 0.6858488630270272, 0.15169972127018583, + 0.22431157058134898, 0.16823638722404644, + 0.0005208971463830214], tspan=(0.0, 0.04), coverage_override=(maxiters = 6, initial_refinement_level = 1, base_level = 1, max_level = 2)) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index ba2565e0dd9..f78f127f434 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -911,16 +911,16 @@ end @trixi_testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2=[0.04937480811868297, 0.06117033019988596, - 0.060998028674664716, 0.03155145889799417, - 0.2319175391388658, 0.02476283192966346, - 0.024483244374818587, 0.035439957899127385, - 0.0016022148194667542], - linf=[0.24749024430983746, 0.2990608279625713, - 0.3966937932860247, 0.22265033744519683, - 0.9757376320946505, 0.12123736788315098, - 0.12837436699267113, 0.17793825293524734, - 0.03460761690059514], + l2=[0.04937478399958968, 0.0611701500558669, + 0.06099805934392425, 0.031551737882277144, + 0.23191853685798858, 0.02476297013104899, + 0.024482975007695532, 0.035440179203707095, + 0.0016002328034991635], + linf=[0.24744671083295033, 0.2990591185187605, + 0.3968520446251412, 0.2226544553988576, + 0.9752669317263143, 0.12117894533967843, + 0.12845218263379432, 0.17795590713819576, + 0.0348517136607105], tspan=(0.0, 0.3)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -934,16 +934,16 @@ end @trixi_testset "elixir_mhd_alfven_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2=[0.02890769490562535, 0.0062599448721613205, - 0.005650300017676721, 0.007334415940022972, - 0.00490446035599909, 0.007202284100220619, - 0.007003258686714405, 0.006734267830082687, - 0.004253003868791559], - linf=[0.17517380432288565, 0.06197353710696667, - 0.038494840938641646, 0.05293345499813148, - 0.03817506476831778, 0.042847170999492534, - 0.03761563456810613, 0.048184237474911844, - 0.04114666955364693], + l2=[0.028905589451357638, 0.006259570019325034, + 0.005649791156739933, 0.0073272570974805004, + 0.004890348793116962, 0.00720944138561451, + 0.0069984328989438115, 0.006729800315219757, + 0.004318314151888631], + linf=[0.17528323378978317, 0.06161030852803388, + 0.0388335541348234, 0.052906440559080926, + 0.0380036034027319, 0.04291841215471082, + 0.03702743958268562, 0.04815794489066357, + 0.0433064571343779], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -1007,16 +1007,18 @@ end @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"), - l2=[0.0364192725149364, 0.0426667193422069, 0.04261673001449095, - 0.025884071405646924, - 0.16181626564020496, 0.017346518770783536, - 0.017291573200291104, 0.026856206495339655, - 0.0007443858043598808], - linf=[0.25144373906033013, 0.32881947152723745, - 0.3053266801502693, 0.20989755319972866, - 0.9927517314507455, 0.1105172121361323, 0.1257708104676617, - 0.1628334844841588, - 0.02624301627479052]) + l2=[0.03641928087745194, 0.04266672246194787, + 0.042616743034675685, + 0.025884076832341982, + 0.16181640309885276, 0.017346521291731105, + 0.017291600359415987, 0.026856207871456043, + 0.0007448774124272682], + linf=[0.25144155032118376, 0.3288086335996786, + 0.30532573631664345, 0.20990150465080706, + 0.9929091025128138, 0.11053858971264774, + 0.12578085409726314, + 0.16283334251103732, + 0.026146463886273865]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_structured_3d.jl b/test/test_structured_3d.jl index ac932b9535e..f4864e45d58 100644 --- a/test/test_structured_3d.jl +++ b/test/test_structured_3d.jl @@ -11,7 +11,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "structured_3d_dgsem") outdir = "out" isdir(outdir) && rm(outdir, recursive = true) -@testset "Structured mesh" begin +@testset "StructuredMesh3D" begin #! format: noindent @trixi_testset "elixir_advection_basic.jl" begin @@ -240,16 +240,17 @@ end @trixi_testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2=[0.009082353036644902, 0.007128360240528109, - 0.006970330025996491, 0.006898850266874514, - 0.03302008823756457, 0.003203389099143526, - 0.003077498677885352, 0.0030740006760477624, - 4.192129696970217e-5], - linf=[0.2883946030582689, 0.25956437344015054, - 0.2614364943543665, 0.24617277938134657, - 1.1370443512475847, 0.1278041831463388, 0.13347391885068594, - 0.1457563463643099, - 0.0021174246048172563], + l2=[0.009082353008355219, 0.007128360330314966, + 0.0069703300260751545, 0.006898850266164216, + 0.033020091335659474, 0.003203389281512797, + 0.0030774985678369746, 0.00307400076520122, + 4.192572922118587e-5], + linf=[0.28839460197220435, 0.25956437090703427, + 0.26143649456148177, 0.24617277684934058, + 1.1370439348603143, 0.12780410700666367, + 0.13347392283166903, + 0.145756208548534, + 0.0021181795153149053], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -263,16 +264,16 @@ end @trixi_testset "elixir_mhd_alfven_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2=[0.003015476175153681, 0.00145499403283373, - 0.0009125744757935803, 0.0017703080480578979, - 0.0013046447673965966, 0.0014564863387645508, - 0.0013332311430907598, 0.001647832598455728, - 0.0013647609788548722], - linf=[0.027510637768610846, 0.02797062834945721, - 0.01274249949295704, 0.038940694415543736, - 0.02200825678588325, 0.03167600959583505, - 0.021420957993862344, 0.03386589835999665, - 0.01888303191983353], + l2=[0.003015390232128414, 0.0014538563096541798, + 0.000912478356719486, 0.0017715065044433436, + 0.0013017575272262197, 0.0014545437537522726, + 0.0013322897333898482, 0.0016493009787844212, + 0.0013747547738038235], + linf=[0.027577067632765795, 0.027912829563483885, + 0.01282206030593043, 0.03911437990598213, + 0.021962225923304324, 0.03169774571258743, + 0.021591564663781426, 0.034028148178115364, + 0.020084593242858988], # Use same polydeg as everything else to prevent long compile times in CI coverage_override=(polydeg = 3,)) # Ensure that we do not have excessive memory allocations @@ -287,16 +288,16 @@ end @trixi_testset "elixir_mhd_alfven_wave.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2=[0.003047854479955232, 0.0014572199588782184, - 0.0009093737183251411, 0.0017937548694553895, - 0.0013010437110755424, 0.0014545607744895874, - 0.001328514015121245, 0.001671342529206066, - 0.0013653963058149186], - linf=[0.027719103797310463, 0.027570111789910784, - 0.012561901006903103, 0.03903568568480584, - 0.021311996934554767, 0.03154849824135775, - 0.020996033645485412, 0.03403185137382961, - 0.019488952445771597], + l2=[0.0030477691235949685, 0.00145609137038748, + 0.0009092809766088607, 0.0017949926915475929, + 0.0012981612165627713, 0.0014525841626158234, + 0.0013275465154956557, 0.0016728767532610933, + 0.0013751925705271012], + linf=[0.02778552932540901, 0.027511633996169835, + 0.012637649797178449, 0.03920805095546112, + 0.02126543791857216, 0.031563506812970266, + 0.02116105422516923, 0.03419432640106229, + 0.020324891223351533], surface_flux=(flux_lax_friedrichs, flux_nonconservative_powell), # Use same polydeg as everything else to prevent long compile times in CI coverage_override=(polydeg = 3,)) @@ -312,16 +313,16 @@ end @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"), - l2=[0.009352631220872144, 0.008058649103542618, - 0.008027041293333663, 0.008071417851552725, - 0.034909149665869485, 0.00393019428600812, - 0.0039219074393817, 0.003906321245184237, - 4.197255300781248e-5], - linf=[0.30749098250807516, 0.2679008863509767, - 0.271243087484388, 0.26545396569129537, - 0.9620950892188596, 0.18163281157498123, - 0.15995708312378454, 0.17918221526906408, - 0.015138346608166353], + l2=[0.009352631216098996, 0.008058649096024162, + 0.00802704129788766, 0.008071417834885589, + 0.03490914976431044, 0.003930194255268652, + 0.003921907459117296, 0.003906321239858786, + 4.1971260184918575e-5], + linf=[0.307491045404509, 0.26790087991041506, + 0.2712430701672931, 0.2654540237991884, + 0.9620943261873176, 0.181632512204141, + 0.15995711137712265, 0.1791807940466812, + 0.015138421396338456], tspan=(0.0, 0.25), # Use same polydeg as everything else to prevent long compile times in CI coverage_override=(polydeg = 3,)) diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index 644995778df..ed9edbee3df 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -305,16 +305,17 @@ end # This test is identical to the one in `test_p4est_2d.jl` besides minor # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.44207324634847545, 0.8804644301177857, 0.8262542320669426, + l2=[0.4420732463420727, 0.8804644301158163, 0.8262542320734158, 0.0, - 0.9615023124189027, 0.10386709616755131, 0.1540308191628843, + 0.9615023124248694, 0.10386709616933161, + 0.15403081916109138, 0.0, - 2.8350276854372125e-5], - linf=[10.04548675437385, 17.998696852394836, 9.575802136190026, + 2.835066224683485e-5], + linf=[10.045486750338348, 17.998696851793447, 9.57580213608948, 0.0, - 19.431290746184473, 1.3821685018474321, 1.8186235976551453, + 19.431290734386764, 1.3821685025605288, 1.8186235976086789, 0.0, - 0.002309422702635547], + 0.0023118793481168537], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_threaded.jl b/test/test_threaded.jl index 760c5ca0d73..a3d52c1923f 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -227,16 +227,16 @@ end @trixi_testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_mhd_ec.jl"), - l2=[0.04937480811868297, 0.06117033019988596, - 0.060998028674664716, 0.03155145889799417, - 0.2319175391388658, 0.02476283192966346, - 0.024483244374818587, 0.035439957899127385, - 0.0016022148194667542], - linf=[0.24749024430983746, 0.2990608279625713, - 0.3966937932860247, 0.22265033744519683, - 0.9757376320946505, 0.12123736788315098, - 0.12837436699267113, 0.17793825293524734, - 0.03460761690059514], + l2=[0.04937478399958968, 0.0611701500558669, + 0.06099805934392425, 0.031551737882277144, + 0.23191853685798858, 0.02476297013104899, + 0.024482975007695532, 0.035440179203707095, + 0.0016002328034991635], + linf=[0.24744671083295033, 0.2990591185187605, + 0.3968520446251412, 0.2226544553988576, + 0.9752669317263143, 0.12117894533967843, + 0.12845218263379432, 0.17795590713819576, + 0.0348517136607105], tspan=(0.0, 0.3)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_type.jl b/test/test_type.jl index d507cded958..62569efe0ec 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1153,8 +1153,7 @@ isdir(outdir) && rm(outdir, recursive = true) zero(RealT)) orientations = [1, 2] directions = [1, 2, 3, 4] - normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), - zero(RealT)) + normal_direction = SVector(one(RealT), zero(RealT)) nonconservative_type_local = Trixi.NonConservativeLocal() nonconservative_type_symmetric = Trixi.NonConservativeSymmetric() nonconservative_terms = [1, 2] @@ -1167,8 +1166,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, + normal_direction, equations)) == RealT @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, equations)) == RealT @@ -1273,9 +1271,7 @@ isdir(outdir) && rm(outdir, recursive = true) zero(RealT)) orientations = [1, 2, 3] directions = [1, 2, 3, 4, 5, 6] - normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), - zero(RealT), - zero(RealT)) + normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == @@ -1285,8 +1281,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, + normal_direction, equations)) == RealT @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, equations)) == RealT @@ -2221,8 +2216,7 @@ isdir(outdir) && rm(outdir, recursive = true) one(RealT)) orientations = [1, 2] directions = [1, 2, 3, 4] - normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), - zero(RealT)) + normal_direction = SVector(one(RealT), zero(RealT)) surface_flux_function = flux_lax_friedrichs dissipation = DissipationLocalLaxFriedrichs() @@ -2242,17 +2236,14 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, + normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_fjordholm_etal(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, + normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, + normal_direction, equations)) == RealT @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction, equations)) == RealT diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 43b18b48e72..c433338a238 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -254,16 +254,17 @@ end @trixi_testset "elixir_mhd_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"), - l2=[0.06418293357851637, 0.12085176618704108, - 0.12085099342419513, 0.07743005602933221, - 0.1622218916638482, 0.04044434425257972, - 0.04044440614962498, 0.05735896706356321, - 0.0020992340041681734], - linf=[0.1417000509328017, 0.3210578460652491, 0.335041095545175, - 0.22500796423572675, - 0.44230628074326406, 0.16743171716317784, - 0.16745989278866702, 0.17700588224362557, - 0.02692320090677309], + l2=[0.06418288595515664, 0.12085170757294698, + 0.12085093463857763, 0.077430018507123, + 0.16221988122574071, 0.040444455755985195, + 0.04044451621612787, 0.05735903066057611, + 0.002095549716217215], + linf=[0.14169585310190325, 0.32104342885987625, + 0.33503526151419405, + 0.22499513309636543, + 0.44231595436029814, 0.16750863202541477, + 0.1675356630213226, 0.1770099359044508, + 0.026783792841168948], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -277,16 +278,16 @@ end @trixi_testset "elixir_mhd_alfven_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave.jl"), - l2=[5.377518922553881e-5, 0.09999999206243514, - 0.09999999206243441, 0.1414213538550799, - 8.770450430886394e-6, 0.0999999926130084, - 0.0999999926130088, 0.14142135396487032, - 1.1553833987291942e-5], - linf=[0.00039334982566352483, 0.14144904937275282, - 0.14144904937277897, 0.20003315928443416, - 6.826863293230012e-5, 0.14146512909995967, - 0.14146512909994702, 0.20006706837452526, - 0.00013645610312810813], + l2=[5.376431895349634e-5, 0.09999999205016862, + 0.09999999205016788, 0.14142135386740418, + 8.767116801867206e-6, 0.09999999259645777, + 0.09999999259645763, 0.14142135397626523, + 1.1559626795684309e-5], + linf=[0.00039380173293024345, 0.14144879547840894, + 0.14144879547843608, 0.2000330663752416, + 7.021503828519293e-5, 0.14146450834000124, + 0.1414645083399998, 0.20006708807562765, + 0.0001375806459241173], tspan=(0.0, 0.5)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 89feae03d9aa2c7004510b895a71bca38dd6580f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 10 Oct 2024 11:04:16 +0200 Subject: [PATCH 03/13] set version to v0.9.0 closes #1997 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 612b65ec6f0..da4be0d79db 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.8.11-DEV" +version = "0.9.0" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From 0c2c36765b137d04ec8634452dc9c3fa96be069f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 10 Oct 2024 11:04:36 +0200 Subject: [PATCH 04/13] set development version to v0.9.1-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index da4be0d79db..395b11be5ab 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.9.0" +version = "0.9.1-DEV" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From a7895796a139c93ce3c30ae7018df23d26cd724e Mon Sep 17 00:00:00 2001 From: Arpit Babbar Date: Thu, 10 Oct 2024 15:43:40 +0200 Subject: [PATCH 05/13] Add velocity functions for different equations (#2086) * Add velocity functions for different equations * Formatter * Fix CI * Format * Export velocity, add docs * Export velocity, add docstring, add unit tests, make MHD velocity like 3D for all dims * Refractor bug fixgit push * Changes as per review * General velocity normal function --------- Co-authored-by: Hendrik Ranocha --- src/equations/compressible_euler_1d.jl | 14 ++- src/equations/compressible_euler_2d.jl | 13 ++ src/equations/compressible_euler_3d.jl | 14 +++ .../compressible_euler_multicomponent_1d.jl | 6 + .../compressible_euler_multicomponent_2d.jl | 14 +++ src/equations/equations.jl | 31 +++++ src/equations/ideal_glm_mhd_1d.jl | 18 ++- src/equations/ideal_glm_mhd_2d.jl | 14 +++ src/equations/ideal_glm_mhd_3d.jl | 14 +++ src/equations/polytropic_euler_2d.jl | 13 ++ src/equations/shallow_water_2d.jl | 12 +- test/test_type.jl | 23 ++++ test/test_unit.jl | 117 ++++++++++++++++++ 13 files changed, 294 insertions(+), 9 deletions(-) diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index a71750ff98c..542a5b7039f 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -412,10 +412,10 @@ See also return SVector(f1, f2, f3) end -# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that -# the normal component is incorporated into the numerical flux. -# -# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a +# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that +# the normal component is incorporated into the numerical flux. +# +# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a # similar implementation. @inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerEquations1D) @@ -943,6 +943,12 @@ end return rho end +@inline function velocity(u, equations::CompressibleEulerEquations1D) + rho = u[1] + v1 = u[2] / rho + return v1 +end + @inline function pressure(u, equations::CompressibleEulerEquations1D) rho, rho_v1, rho_e = u p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index e964d9c2b21..37a4f362101 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1957,6 +1957,19 @@ end return rho end +@inline function velocity(u, equations::CompressibleEulerEquations2D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + return SVector(v1, v2) +end + +@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations2D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + @inline function pressure(u, equations::CompressibleEulerEquations2D) rho, rho_v1, rho_v2, rho_e = u p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 4f4d4553a8f..56a24238fcc 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1692,6 +1692,20 @@ end return rho end +@inline function velocity(u, equations::CompressibleEulerEquations3D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + +@inline function velocity(u, orientation::Int, equations::CompressibleEulerEquations3D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + @inline function pressure(u, equations::CompressibleEulerEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e = u p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index da434579f76..c25a4631185 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -618,4 +618,10 @@ end return SVector{ncomponents(equations), real(equations)}(u[i + 2] * v for i in eachcomponent(equations)) end + +@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations1D) + rho = density(u, equations) + v1 = u[1] / rho + return v1 +end end # @muladd diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 2424ad0301c..2278fe58eb3 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -825,4 +825,18 @@ end return SVector{ncomponents(equations), real(equations)}(u[i + 3] * v for i in eachcomponent(equations)) end + +@inline function velocity(u, equations::CompressibleEulerMulticomponentEquations2D) + rho = density(u, equations) + v1 = u[1] / rho + v2 = u[2] / rho + return SVector(v1, v2) +end + +@inline function velocity(u, orientation::Int, + equations::CompressibleEulerMulticomponentEquations2D) + rho = density(u, equations) + v = u[orientation] / rho + return v +end end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 7875954c5f0..0304ea16cfe 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -307,6 +307,37 @@ The inverse conversion is performed by [`cons2prim`](@ref). """ function prim2cons end +""" + velocity(u, equations) + +Return the velocity vector corresponding to the equations, e.g., fluid velocity for +Euler's equations. The velocity in certain orientation or normal direction (scalar) can be computed +with `velocity(u, orientation, equations)` or `velocity(u, normal_direction, equations)` +respectively. The `velocity(u, normal_direction, equations)` function calls +`velocity(u, equations)` to compute the velocity vector and then normal vector, thus allowing +a general function to be written for the AbstractEquations type. However, the +`velocity(u, orientation, equations)` is written for each equation separately to ensure +only the velocity in the desired direction (orientation) is computed. +`u` is a vector of the conserved variables at a single node, i.e., a vector +of the correct length `nvariables(equations)`. +""" +function velocity end + +@inline function velocity(u, normal_direction::AbstractVector, + equations::AbstractEquations{2}) + vel = velocity(u, equations) + v = vel[1] * normal_direction[1] + vel[2] * normal_direction[2] + return v +end + +@inline function velocity(u, normal_direction::AbstractVector, + equations::AbstractEquations{3}) + vel = velocity(u, equations) + v = vel[1] * normal_direction[1] + vel[2] * normal_direction[2] + + vel[3] * normal_direction[3] + return v +end + """ entropy(u, equations) diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 3253e4410f8..685ad6f8f3f 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -326,7 +326,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, Sdiff = SsR - SsL # Compute HLL values for vStar, BStar - # These correspond to eq. (28) and (30) from the referenced paper + # These correspond to eq. (28) and (30) from the referenced paper # and the classic HLL intermediate state given by (2) rho_HLL = (SsR * rho_rr - SsL * rho_ll - (f_rr[1] - f_ll[1])) / Sdiff @@ -394,7 +394,7 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, B1_rr * (B1_rr * v1_rr + B2_rr * v2_rr + B3_rr * v3_rr))) / SdiffStar # (23) - # Classic HLLC update (32) + # Classic HLLC update (32) f1 = f_rr[1] + SsR * (densStar - u_rr[1]) f2 = f_rr[2] + SsR * (mom_1_Star - u_rr[2]) f3 = f_rr[3] + SsR * (mom_2_Star - u_rr[3]) @@ -555,6 +555,20 @@ end return rho end +@inline function velocity(u, equations::IdealGlmMhdEquations1D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + +@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations1D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + @inline function pressure(u, equations::IdealGlmMhdEquations1D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index a43c6b0e619..9383fdbe961 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -1111,6 +1111,20 @@ end return rho end +@inline function velocity(u, equations::IdealGlmMhdEquations2D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + +@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations2D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + @inline function pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 9bb05cb3f89..8d7d0461398 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -1001,6 +1001,20 @@ end return u[1] end +@inline function velocity(u, equations::IdealGlmMhdEquations3D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + v3 = u[4] / rho + return SVector(v1, v2, v3) +end + +@inline function velocity(u, orientation::Int, equations::IdealGlmMhdEquations3D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + @inline function pressure(u, equations::IdealGlmMhdEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl index 571d4723f6f..98f0f7f1dc1 100644 --- a/src/equations/polytropic_euler_2d.jl +++ b/src/equations/polytropic_euler_2d.jl @@ -389,6 +389,19 @@ end return rho end +@inline function velocity(u, equations::PolytropicEulerEquations2D) + rho = u[1] + v1 = u[2] / rho + v2 = u[3] / rho + return SVector(v1, v2) +end + +@inline function velocity(u, orientation::Int, equations::PolytropicEulerEquations2D) + rho = u[1] + v = u[orientation + 1] / rho + return v +end + @inline function pressure(u, equations::PolytropicEulerEquations2D) rho, rho_v1, rho_v2 = u p = equations.kappa * rho^equations.gamma diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 30eca095aa4..7ba36d6ec92 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -282,7 +282,7 @@ Further details are available in the papers: shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) - Patrick Ersing, Andrew R. Winters (2023) - An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on curvilinear meshes [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ @@ -834,6 +834,12 @@ end return SVector(v1, v2) end +@inline function velocity(u, orientation::Int, equations::ShallowWaterEquations2D) + h = u[1] + v = u[orientation + 1] / h + return v +end + # Convert conservative variables to primitive @inline function cons2prim(u, equations::ShallowWaterEquations2D) h, _, _, b = u @@ -897,11 +903,11 @@ end equations::ShallowWaterEquations2D) Calculate Roe-averaged velocity `v_roe` and wavespeed `c_roe = sqrt{g * h_roe}` depending on direction. -See for instance equation (62) in +See for instance equation (62) in - Paul A. Ullrich, Christiane Jablonowski, and Bram van Leer (2010) High-order finite-volume methods for the shallow-water equations on the sphere [DOI: 10.1016/j.jcp.2010.04.044](https://doi.org/10.1016/j.jcp.2010.04.044) -Or [this slides](https://faculty.washington.edu/rjl/classes/am574w2011/slides/am574lecture20nup3.pdf), +Or [this slides](https://faculty.washington.edu/rjl/classes/am574w2011/slides/am574lecture20nup3.pdf), slides 8 and 9. """ @inline function calc_wavespeed_roe(u_ll, u_rr, orientation::Integer, diff --git a/test/test_type.jl b/test/test_type.jl index 62569efe0ec..92274bf4248 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -167,6 +167,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred cons2entropy(u, equations)) == RealT @test eltype(@inferred entropy2cons(u, equations)) == RealT @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred velocity(u, equations)) == RealT @test typeof(@inferred pressure(u, equations)) == RealT @test typeof(@inferred density_pressure(u, equations)) == RealT @test typeof(@inferred entropy(cons, equations)) == RealT @@ -218,6 +219,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == RealT @@ -257,6 +259,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for orientation in orientations + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) == RealT @@ -300,6 +303,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @@ -367,6 +371,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_shima_etal(u_ll, u_rr, normal_direction, equations)) == RealT @@ -392,6 +397,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT for orientation in orientations + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_shima_etal(u_ll, u_rr, orientation, equations)) == RealT @@ -421,6 +427,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT end + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @@ -456,6 +463,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == RealT + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT @@ -499,11 +507,13 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) == RealT + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, equations)) == RealT for orientation in orientations + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, equations)) == RealT @@ -515,6 +525,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @@ -1094,6 +1105,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, equations)) == @@ -1164,6 +1176,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, normal_direction, @@ -1181,6 +1194,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT for orientation in orientations + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, equations)) == RealT @@ -1219,6 +1233,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT end + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @@ -1279,6 +1294,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, normal_direction, @@ -1296,6 +1312,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT for orientation in orientations + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, equations)) == RealT @@ -1314,6 +1331,7 @@ isdir(outdir) && rm(outdir, recursive = true) equations)) == RealT end + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @@ -2081,6 +2099,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == RealT + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction, equations)) == @@ -2096,6 +2115,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT for orientation in orientations + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation, equations)) == @@ -2108,6 +2128,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end + @test eltype(@inferred velocity(u, equations)) == RealT @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT @test eltype(@inferred cons2prim(u, equations)) == RealT @test eltype(@inferred prim2cons(u, equations)) == RealT @@ -2234,6 +2255,7 @@ isdir(outdir) && rm(outdir, recursive = true) surface_flux_function, equations)) == RealT + @test eltype(@inferred velocity(u, normal_direction, equations)) == RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, normal_direction, @@ -2276,6 +2298,7 @@ isdir(outdir) && rm(outdir, recursive = true) RealT end + @test eltype(@inferred velocity(u, orientation, equations)) == RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT @test eltype(@inferred flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation, diff --git a/test/test_unit.jl b/test/test_unit.jl index 5831122ffe2..bccdcf8faaa 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1737,6 +1737,123 @@ end @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref), 1.803e-5, atol = 5e-8) end + +# Velocity functions are present in many equations and are tested here +@testset "Velocity functions for different equations" begin + gamma = 1.4 + rho = pi * pi + pres = sqrt(pi) + v1, v2, v3 = pi, exp(1.0), exp(pi) # use pi, exp to test with non-trivial numbers + v_vector = SVector(v1, v2, v3) + normal_direction_2d = SVector(pi^2, pi^3) + normal_direction_3d = SVector(normal_direction_2d..., pi^4) + v_normal_1d = v1 * normal_direction_2d[1] + v_normal_2d = v1 * normal_direction_2d[1] + v2 * normal_direction_2d[2] + v_normal_3d = v_normal_2d + v3 * normal_direction_3d[3] + + equations_euler_1d = CompressibleEulerEquations1D(gamma) + u = prim2cons(SVector(rho, v1, pres), equations_euler_1d) + @test isapprox(velocity(u, equations_euler_1d), v1) + + equations_euler_2d = CompressibleEulerEquations2D(gamma) + u = prim2cons(SVector(rho, v1, v2, pres), equations_euler_2d) + @test isapprox(velocity(u, equations_euler_2d), SVector(v1, v2)) + @test isapprox(velocity(u, normal_direction_2d, equations_euler_2d), v_normal_2d) + for orientation in 1:2 + @test isapprox(velocity(u, orientation, equations_euler_2d), + v_vector[orientation]) + end + + equations_euler_3d = CompressibleEulerEquations3D(gamma) + u = prim2cons(SVector(rho, v1, v2, v3, pres), equations_euler_3d) + @test isapprox(velocity(u, equations_euler_3d), SVector(v1, v2, v3)) + @test isapprox(velocity(u, normal_direction_3d, equations_euler_3d), v_normal_3d) + for orientation in 1:3 + @test isapprox(velocity(u, orientation, equations_euler_3d), + v_vector[orientation]) + end + + rho1, rho2 = rho, rho * pi # use pi to test with non-trivial numbers + gammas = (gamma, exp(gamma)) + gas_constants = (0.387, 1.678) # Standard numbers + 0.1 + + equations_multi_euler_1d = CompressibleEulerMulticomponentEquations1D(; gammas, + gas_constants) + u = prim2cons(SVector(v1, pres, rho1, rho2), equations_multi_euler_1d) + @test isapprox(velocity(u, equations_multi_euler_1d), v1) + + equations_multi_euler_2d = CompressibleEulerMulticomponentEquations2D(; gammas, + gas_constants) + u = prim2cons(SVector(v1, v2, pres, rho1, rho2), equations_multi_euler_2d) + @test isapprox(velocity(u, equations_multi_euler_2d), SVector(v1, v2)) + @test isapprox(velocity(u, normal_direction_2d, equations_multi_euler_2d), + v_normal_2d) + for orientation in 1:2 + @test isapprox(velocity(u, orientation, equations_multi_euler_2d), + v_vector[orientation]) + end + + kappa = 0.1 * pi # pi for non-trivial test + equations_polytropic = PolytropicEulerEquations2D(gamma, kappa) + u = prim2cons(SVector(rho, v1, v2), equations_polytropic) + @test isapprox(velocity(u, equations_polytropic), SVector(v1, v2)) + equations_polytropic = CompressibleEulerMulticomponentEquations2D(; gammas, + gas_constants) + u = prim2cons(SVector(v1, v2, pres, rho1, rho2), equations_polytropic) + @test isapprox(velocity(u, equations_polytropic), SVector(v1, v2)) + @test isapprox(velocity(u, normal_direction_2d, equations_polytropic), v_normal_2d) + for orientation in 1:2 + @test isapprox(velocity(u, orientation, equations_polytropic), + v_vector[orientation]) + end + + B1, B2, B3 = pi^3, pi^4, pi^5 + equations_ideal_mhd_1d = IdealGlmMhdEquations1D(gamma) + u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3), equations_ideal_mhd_1d) + @test isapprox(velocity(u, equations_ideal_mhd_1d), SVector(v1, v2, v3)) + for orientation in 1:3 + @test isapprox(velocity(u, orientation, equations_ideal_mhd_1d), + v_vector[orientation]) + end + + psi = exp(0.1) + equations_ideal_mhd_2d = IdealGlmMhdEquations2D(gamma) + u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3, psi), + equations_ideal_mhd_2d) + @test isapprox(velocity(u, equations_ideal_mhd_2d), SVector(v1, v2, v3)) + @test isapprox(velocity(u, normal_direction_2d, equations_ideal_mhd_2d), + v_normal_2d) + for orientation in 1:3 + @test isapprox(velocity(u, orientation, equations_ideal_mhd_2d), + v_vector[orientation]) + end + + equations_ideal_mhd_3d = IdealGlmMhdEquations3D(gamma) + u = prim2cons(SVector(rho, v1, v2, v3, pres, B1, B2, B3, psi), + equations_ideal_mhd_3d) + @test isapprox(velocity(u, equations_ideal_mhd_3d), SVector(v1, v2, v3)) + @test isapprox(velocity(u, normal_direction_3d, equations_ideal_mhd_3d), + v_normal_3d) + for orientation in 1:3 + @test isapprox(velocity(u, orientation, equations_ideal_mhd_3d), + v_vector[orientation]) + end + + H, b = exp(pi), exp(pi^2) + gravity_constant, H0 = 9.91, 0.1 # Standard numbers + 0.1 + shallow_water_1d = ShallowWaterEquations1D(; gravity_constant, H0) + u = prim2cons(SVector(H, v1, b), shallow_water_1d) + @test isapprox(velocity(u, shallow_water_1d), v1) + + shallow_water_2d = ShallowWaterEquations2D(; gravity_constant, H0) + u = prim2cons(SVector(H, v1, v2, b), shallow_water_2d) + @test isapprox(velocity(u, shallow_water_2d), SVector(v1, v2)) + @test isapprox(velocity(u, normal_direction_2d, shallow_water_2d), v_normal_2d) + for orientation in 1:2 + @test isapprox(velocity(u, orientation, shallow_water_2d), + v_vector[orientation]) + end +end end end #module From 0b3d74834f2315b2b033745a56ecf6a85d011141 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 10 Oct 2024 16:56:31 +0200 Subject: [PATCH 06/13] set version to v0.9.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 395b11be5ab..ce6361698db 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.9.1-DEV" +version = "0.9.1" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From efdb5d49c823d689d5014475bbaef09784692158 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 10 Oct 2024 16:56:48 +0200 Subject: [PATCH 07/13] set development version to v0.9.2-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ce6361698db..285d29ce639 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.9.1" +version = "0.9.2-DEV" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" From 44ad3b1bf08e49b6fe4865c258aca1da219c2e51 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Fri, 11 Oct 2024 08:35:30 +0200 Subject: [PATCH 08/13] CompatHelper: bump compat for Trixi to 0.9 for package benchmark, (keep existing compat) (#2106) Co-authored-by: CompatHelper Julia --- benchmark/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/Project.toml b/benchmark/Project.toml index 958b1ca23f6..6647e39381c 100644 --- a/benchmark/Project.toml +++ b/benchmark/Project.toml @@ -8,4 +8,4 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" BenchmarkTools = "0.5, 0.7, 1.0" OrdinaryDiffEq = "5.65, 6" PkgBenchmark = "0.2.10" -Trixi = "0.4, 0.5, 0.6, 0.7, 0.8" +Trixi = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9" From a05a68d91fbe1b0ff769a7cbcf98093893ff51f9 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Sun, 13 Oct 2024 18:33:40 -1000 Subject: [PATCH 09/13] Fix small amount of `Float64` computation --- src/equations/compressible_euler_multicomponent_1d.jl | 11 +++++------ src/equations/compressible_euler_multicomponent_2d.jl | 11 +++++------ 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index c25a4631185..f9413e0f1f9 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -204,13 +204,12 @@ function initial_condition_weak_blast_wave(x, t, prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ? 2^(i - 1) * (1 - 2) / + (RealT(1) - + 2^ncomponents(equations)) : + 2^(i - 1) * (1 - 2) * + RealT(1.1691) / (1 - - 2^ncomponents(equations)) * - one(RealT) : - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - convert(RealT, 1.1691) + 2^ncomponents(equations)) for i in eachcomponent(equations)) v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 2278fe58eb3..bac76c9a021 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -222,13 +222,12 @@ function initial_condition_weak_blast_wave(x, t, prim_rho = SVector{ncomponents(equations), real(equations)}(r > 0.5f0 ? 2^(i - 1) * (1 - 2) / + (RealT(1) - + 2^ncomponents(equations)) : + 2^(i - 1) * (1 - 2) * + RealT(1.1691) / (1 - - 2^ncomponents(equations)) * - one(RealT) : - 2^(i - 1) * (1 - 2) / - (1 - - 2^ncomponents(equations)) * - convert(RealT, 1.1691) + 2^ncomponents(equations)) for i in eachcomponent(equations)) v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi From cdd552ff5ee342e8e9343cdc10a75aebb10df73f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 14 Oct 2024 11:55:47 +0200 Subject: [PATCH 10/13] Entropy bounded limiter (#2078) * start working * get EB limiter to work * Refine EB limiter * fmt * 2D/3D Limiters * introduce parameter of tolerated entropy decrease * change to exp_entropy_diff * fmt * export * docstring * 3d example * exam,ples * examples * examples * fmt * adaptive * examples * fmt * tests + docu * fmt + spelling * comment * Apply suggestions from code review Co-authored-by: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> * Apply suggestions from code review * increase => change * fmt * spelling * revisit * revert * back to decrease * fmt * typo * comments * theta * longer sim * removed unused code * for coverage * Apply suggestions from code review * longer runtime for coverage * vals * reproducible test values * shorten * changes * testvals * test vals * Update src/callbacks_stage/entropy_bounded_limiter_2d.jl * Apply suggestions from code review * Update docs/literate/src/files/shock_capturing.jl * more iters * Apply suggestions from code review * comments * mention step_limiter! capability * Update docs/literate/src/files/shock_capturing.jl * Update src/callbacks_stage/entropy_bounded_limiter.jl * Update src/callbacks_stage/entropy_bounded_limiter.jl * Apply suggestions from code review --------- Co-authored-by: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com> --- README.md | 4 +- docs/literate/src/files/shock_capturing.jl | 141 +++++++++++++++++- docs/src/index.md | 4 + .../elixir_mhd_amr_entropy_bounded.jl | 128 ++++++++++++++++ ...elixir_euler_blast_wave_entropy_bounded.jl | 80 ++++++++++ ...uler_colliding_flow_amr_entropy_bounded.jl | 109 ++++++++++++++ src/Trixi.jl | 2 +- src/callbacks_stage/callbacks_stage.jl | 1 + .../entropy_bounded_limiter.jl | 70 +++++++++ .../entropy_bounded_limiter_1d.jl | 72 +++++++++ .../entropy_bounded_limiter_2d.jl | 76 ++++++++++ .../entropy_bounded_limiter_3d.jl | 76 ++++++++++ test/test_p4est_3d.jl | 37 +++++ test/test_tree_1d_euler.jl | 15 ++ test/test_tree_2d_euler.jl | 28 ++++ 15 files changed, 836 insertions(+), 7 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl create mode 100644 examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl create mode 100644 examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter_1d.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter_2d.jl create mode 100644 src/callbacks_stage/entropy_bounded_limiter_3d.jl diff --git a/README.md b/README.md index 30e9da243d7..70c28799f31 100644 --- a/README.md +++ b/README.md @@ -36,9 +36,11 @@ installation and postprocessing procedures. Its features include: * Discontinuous Galerkin methods * Kinetic energy-preserving and entropy-stable methods based on flux differencing * Entropy-stable shock capturing + * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl) +* Advanced limiting strategies * Positivity-preserving limiting * Subcell invariant domain-preserving (IDP) limiting - * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl) + * Entropy-bounded limiting * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/) * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods) * [Strong stability preserving methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)) diff --git a/docs/literate/src/files/shock_capturing.jl b/docs/literate/src/files/shock_capturing.jl index dd6698c2a86..62388416d80 100644 --- a/docs/literate/src/files/shock_capturing.jl +++ b/docs/literate/src/files/shock_capturing.jl @@ -96,19 +96,19 @@ # [Zhang, Shu (2011)](https://doi.org/10.1098/rspa.2011.0153). # It works the following way. For every passed (scalar) variable and for every DG element we calculate -# the minimal value $value_{min}$. If this value falls below the given threshold $\varepsilon$, +# the minimal value $value_\text{min}$. If this value falls below the given threshold $\varepsilon$, # the approximation is slightly adapted such that the minimal value of the relevant variable lies # now above the threshold. # ```math -# \underline{u}^{new} = \theta * \underline{u} + (1-\theta) * u_{mean} +# \underline{u}^\text{new} = \theta * \underline{u} + (1-\theta) * u_\text{mean} # ``` # where $\underline{u}$ are the collected pointwise evaluation coefficients in element $e$ and -# $u_{mean}$ the integral mean of the quantity in $e$. The new coefficients are a convex combination +# $u_\text{mean}$ the integral mean of the quantity in $e$. The new coefficients are a convex combination # of these two values with factor # ```math -# \theta = \frac{value_{mean} - \varepsilon}{value_{mean} - value_{min}}, +# \theta = \frac{value_\text{mean} - \varepsilon}{value_\text{mean} - value_\text{min}}, # ``` -# where $value_{mean}$ is the relevant variable evaluated for the mean value $u_{mean}$. +# where $value_\text{mean}$ is the relevant variable evaluated for the mean value $u_\text{mean}$. # The adapted approximation keeps the exact same mean value, but the relevant variable is now greater # or equal the threshold $\varepsilon$ at every node in every element. @@ -225,6 +225,137 @@ sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition=false using Plots plot(sol) +# # Entropy bounded limiter + +# As argued in the description of the positivity preserving limiter above it might sometimes be +# necessary to apply advanced techniques to ensure a physically meaningful solution. +# Apart from the positivity of pressure and density, the physical entropy of the system should increase +# over the course of a simulation, see e.g. [this](https://doi.org/10.1016/0168-9274(86)90029-2) paper by Tadmor where this property is +# shown for the compressible Euler equations. +# As this is not necessarily the case for the numerical approximation (especially for the high-order, non-diffusive DG discretizations), +# Lv and Ihme devised an a-posteriori limiter in [this paper](https://doi.org/10.1016/j.jcp.2015.04.026) which can be applied after each Runge-Kutta stage. +# This limiter enforces a non-decrease in the physical, thermodynamic entropy $S$ +# by bounding the entropy decrease (entropy increase is always tolerated) $\Delta S$ in each grid cell. +# +# This translates into a requirement that the entropy of the limited approximation $S\Big(\mathcal{L}\big[\boldsymbol u(\Delta t) \big] \Big)$ should be +# greater or equal than the previous iterates' entropy $S\big(\boldsymbol u(0) \big)$, enforced at each quadrature point: +# ```math +# S\Big(\mathcal{L}\big[\boldsymbol u(\Delta t, \boldsymbol{x}_i) \big] \Big) \overset{!}{\geq} S\big(\boldsymbol u(0, \boldsymbol{x}_i) \big), \quad i = 1, \dots, (k+1)^d +# ``` +# where $k$ denotes the polynomial degree of the element-wise solution and $d$ is the spatial dimension. +# For an ideal gas, the thermodynamic entropy $S(\boldsymbol u) = S(p, \rho)$ is given by +# ```math +# S(p, \rho) = \ln \left( \frac{p}{\rho^\gamma} \right) \: . +# ``` +# Thus, the non-decrease in entropy can be reformulated as +# ```math +# p(\boldsymbol{x}_i) - e^{ S\big(\boldsymbol u(0, \boldsymbol{x}_i) \big)} \cdot \rho(\boldsymbol{x}_i)^\gamma \overset{!}{\geq} 0, \quad i = 1, \dots, (k+1)^d \: . +# ``` +# In a practical simulation, we might tolerate a maximum (exponentiated) entropy decrease per element, i.e., +# ```math +# \Delta e^S \coloneqq \min_{i} \left\{ p(\boldsymbol{x}_i) - e^{ S\big(\boldsymbol u(0, \boldsymbol{x}_i) \big)} \cdot \rho(\boldsymbol{x}_i)^\gamma \right\} < c +# ``` +# with hyper-parameter $c$ which is to be specified by the user. +# The default value for the corresponding parameter $c=$ `exp_entropy_decrease_max` is set to $-10^{-13}$, i.e., slightly less than zero to +# avoid spurious limiter actions for cells in which the entropy remains effectively constant. +# Other values can be specified by setting the `exp_entropy_decrease_max` keyword in the constructor of the limiter: +# ```julia +# stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max=-1e-9) +# ``` +# Smaller values (larger in absolute value) for `exp_entropy_decrease_max` relax the entropy increase requirement and are thus less diffusive. +# On the other hand, for larger values (smaller in absolute value) of `exp_entropy_decrease_max` the limiter acts more often and the solution becomes more diffusive. +# +# In particular, we compute again a limiting parameter $\vartheta \in [0, 1]$ which is then used to blend the +# unlimited nodal values $\boldsymbol u$ with the mean value $\boldsymbol u_{\text{mean}}$ of the element: +# ```math +# \mathcal{L} [\boldsymbol u](\vartheta) \coloneqq (1 - \vartheta) \boldsymbol u + \vartheta \cdot \boldsymbol u_{\text{mean}} +# ``` +# For the exact definition of $\vartheta$ the interested user is referred to section 4.4 of the paper by Lv and Ihme. +# Note that therein the limiting parameter is denoted by $\epsilon$, which is not to be confused with the threshold $\varepsilon$ of the Zhang-Shu limiter. + +# As for the positivity preserving limiter, the entropy bounded limiter may be applied after every Runge-Kutta stage. +# Both fixed timestep methods such as [`CarpenterKennedy2N54`](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) and +# adaptive timestep methods such as [`SSPRK43`](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) are supported. +# We would like to remark that of course every `stage_limiter!` can also be used as a `step_limiter!`, i.e., +# acting only after the full time step has been taken. + +# As an example, we consider a variant of the [1D medium blast wave example](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl) +# wherein the shock capturing method discussed above is employed to handle the shock. +# Here, we use a standard DG solver with HLLC surface flux: +using Trixi + +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc) + +# The remaining setup is the same as in the standard example: + +using OrdinaryDiffEq + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations1D(1.4) + +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) + ## Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + ## Set up polar coordinates + inicenter = SVector(0.0) + x_norm = x[1] - inicenter[1] + r = abs(x_norm) + ## The following code is equivalent to + ## phi = atan(0.0, x_norm) + ## cos_phi = cos(phi) + ## in 1D but faster + cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + + ## Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_blast_wave + +coordinates_min = (-2.0,) +coordinates_max = (2.0,) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +# We specify the `stage_limiter!` supplied to the classic SSPRK33 integrator +# with strict entropy increase enforcement +# (recall that the tolerated exponentiated entropy decrease is set to -1e-13). +stage_limiter! = EntropyBoundedLimiter() + +# We run the simulation with the SSPRK33 method and the entropy bounded limiter: +sol = solve(ode, SSPRK33(stage_limiter!); + dt = 1.0, + callback = callbacks); + +using Plots +plot(sol) # ## Package versions diff --git a/docs/src/index.md b/docs/src/index.md index 0e4749dde33..88752cc7f11 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -32,6 +32,10 @@ installation and postprocessing procedures. Its features include: * Positivity-preserving limiting * Subcell invariant domain-preserving (IDP) limiting * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl) +* Advanced limiting strategies + * Positivity-preserving limiting + * Subcell invariant domain-preserving (IDP) limiting + * Entropy-bounded limiting * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/) * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods) * [Strong stability preserving methods](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Explicit-Strong-Stability-Preserving-Runge-Kutta-Methods-for-Hyperbolic-PDEs-(Conservation-Laws)) diff --git a/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl b/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl new file mode 100644 index 00000000000..00ab84bf4c0 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_mhd_amr_entropy_bounded.jl @@ -0,0 +1,128 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations + +equations = IdealGlmMhdEquations3D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D) + +Weak magnetic blast wave setup taken from Section 6.1 of the paper: +- A. M. Rueda-Ramírez, S. Hennemann, F. J. Hindenlang, A. R. Winters, G. J. Gassner (2021) + An entropy stable nodal discontinuous Galerkin method for the resistive MHD + equations. Part II: Subcell finite volume shock capturing + [doi: 10.1016/j.jcp.2021.110580](https://doi.org/10.1016/j.jcp.2021.110580) +""" +function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D) + # Center of the blast wave is selected for the domain [0, 3]^3 + inicenter = (1.5, 1.5, 1.5) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + delta_0 = 0.1 + r_0 = 0.3 + lambda = exp(5.0 / delta_0 * (r - r_0)) + + prim_inner = SVector(1.2, 0.1, 0.0, 0.1, 0.9, 1.0, 1.0, 1.0, 0.0) + prim_outer = SVector(1.2, 0.2, -0.4, 0.2, 0.3, 1.0, 1.0, 1.0, 0.0) + prim_vars = (prim_inner + lambda * prim_outer) / (1.0 + lambda) + + return prim2cons(prim_vars, equations) +end +initial_condition = initial_condition_blast_wave + +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell) +volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) +polydeg = 3 +volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but with slightly less warping. +# The mapping will be interpolated at tree level, and then refined without changing +# the geometry interpolant. +function mapping(xi_, eta_, zeta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + zeta = 1.5 * zeta_ + 1.5 + + y = eta + + 3 / 11 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + x = xi + + 3 / 11 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + z = zeta + + 3 / 11 * (cos(0.5 * pi * (2 * x - 3) / 3) * + cos(pi * (2 * y - 3) / 3) * + cos(0.5 * pi * (2 * zeta - 3) / 3)) + + return SVector(x, y, z) +end + +trees_per_dimension = (2, 2, 2) +mesh = P4estMesh(trees_per_dimension, + polydeg = 3, + mapping = mapping, + initial_refinement_level = 2, + periodicity = true) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, + variable = density_pressure) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 2, + max_level = 4, max_threshold = 0.15) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +cfl = 3.0 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + stepsize_callback, + glm_speed_callback) + +# Stage limiter required for this high CFL +stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max = -5e-3) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, 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/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl new file mode 100644 index 00000000000..fa81defb6ca --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_blast_wave_entropy_bounded.jl @@ -0,0 +1,80 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations1D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) + +A medium blast wave taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0) + x_norm = x[1] - inicenter[1] + r = abs(x_norm) + # The following code is equivalent to + # phi = atan(0.0, x_norm) + # cos_phi = cos(phi) + # in 1D but faster + cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, p), equations) +end +initial_condition = initial_condition_blast_wave + +# Note: We do not need to use the shock-capturing methodology here, +# in contrast to the standard `euler_blast_wave.jl` example. +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc) + +coordinates_min = (-2.0,) +coordinates_max = (2.0,) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 12.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +stage_limiter! = EntropyBoundedLimiter() + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK33(stage_limiter!); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl new file mode 100644 index 00000000000..6e7cac9a175 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr_entropy_bounded.jl @@ -0,0 +1,109 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.001 # almost isothermal when gamma reaches 1 +equations = CompressibleEulerEquations2D(gamma) + +# This is a hand made colliding flow setup without reference. Features Mach=70 inflow from both +# sides, with relative low temperature, such that pressure keeps relatively small +# Computed with gamma close to 1, to simulate isothermal gas +function initial_condition_colliding_flow_astro(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # resolution 128^2 elements (refined close to the interface) and polydeg=3 (total of 512^2 DOF) + # domain size is [-64,+64]^2 + @unpack gamma = equations + # the quantities are chosen such, that they are as close as possible to the astro examples + # keep in mind, that in the astro example, the physical units are weird (parsec, mega years, ...) + rho = 0.0247 + c = 0.2 + p = c^2 / gamma * rho + vel = 13.907432274789372 + slope = 1.0 + v1 = -vel * tanh(slope * x[1]) + # add small initial disturbance to the field, but only close to the interface + if abs(x[1]) < 10 + v1 = v1 * (1 + 0.01 * sin(pi * x[2])) + end + v2 = 0.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_colliding_flow_astro + +boundary_conditions = (x_neg = BoundaryConditionDirichlet(initial_condition_colliding_flow_astro), + x_pos = BoundaryConditionDirichlet(initial_condition_colliding_flow_astro), + y_neg = boundary_condition_periodic, + y_pos = boundary_condition_periodic) + +surface_flux = flux_lax_friedrichs # HLLC needs more shock capturing (alpha_max) +volume_flux = flux_ranocha # works with Chandrashekar flux as well +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) + +# shock capturing necessary for this tough example, however alpha_max = 0.5 is fine +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.0001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-64.0, -64.0) +coordinates_max = (64.0, 64.0) + +# only refinement in a patch. Needs x=-17/+17 to trigger refinement due to coarse base mesh +refinement_patches = ((type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)), + (type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)), + (type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)), + (type = "box", coordinates_min = (-17, -64), + coordinates_max = (17, 64)) + #(type="box", coordinates_min=(-17, -64), coordinates_max=(17, 64)), # very high resolution, takes about 1000s on 2 cores + ) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + refinement_patches = refinement_patches, + periodicity = (false, true), + n_cells_max = 100_000) +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 25.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) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback) + +stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max = -1.3e-4) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + dt = 1e-2, save_everystep = false, adaptive = true, + callback = callbacks); + +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index d9b9245918e..0cedec782df 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -274,7 +274,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt, export ControllerThreeLevel, ControllerThreeLevelCombined, IndicatorLöhner, IndicatorLoehner, IndicatorMax -export PositivityPreservingLimiterZhangShu +export PositivityPreservingLimiterZhangShu, EntropyBoundedLimiter export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured, ode_default_options diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index d5abc1d227d..f4f6e72e844 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -8,4 +8,5 @@ include("positivity_zhang_shu.jl") include("subcell_limiter_idp_correction.jl") include("subcell_bounds_check.jl") +include("entropy_bounded_limiter.jl") end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter.jl b/src/callbacks_stage/entropy_bounded_limiter.jl new file mode 100644 index 00000000000..f8dc9a612d6 --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter.jl @@ -0,0 +1,70 @@ +# 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 + +mutable struct EntropyBoundedLimiter{RealT <: Real} + exp_entropy_decrease_max::RealT # < 0 +end + +""" + EntropyBoundedLimiter{RealT <: Real}(; exp_entropy_decrease_max = 1f-13) + +Entropy-bounded limiter by +- Lv, Ihme (2015) + Entropy-bounded discontinuous Galerkin scheme for Euler equations + [doi: 10.1016/j.jcp.2015.04.026](https://doi.org/10.1016/j.jcp.2015.04.026) + +This is an ideal-gas specific limiter that bounds the (unphysical) decrease of the thermodynamic entropy +per element from one time step (or Runge-Kutta stage) to the next. +The parameter `exp_entropy_decrease_max` is the maximum allowed exponentiated +entropy decrease per element at each element's node. + +In the original version of the paper, this value is set to zero to +ensure that the entropy does not decrease, i.e., guarantee entropy stability in the sense of +- Tadmor (1986) + A minimum entropy principle in the gas dynamics equations + [doi: 10.1016/0168-9274(86)90029-2](https://doi.org/10.1016/0168-9274(86)90029-2) + +This, however, leads in general to very diffusive solutions for timesteps violating +a CFL condition (Lemma 3 in Lv, Ihme (2015)) which is required for entropy stability in the mean values. +Since most practical simulations will employ a significantly larger timestep, one can relax the +strict entropy increase requirement by setting `exp_entropy_decrease_max` to a negative value. +The limiter acts if the exponentiated entropy decrease on an element is larger than +`exp_entropy_decrease_max`. +This means that if the change in exponentiated entropy lies *below* `exp_entropy_decrease_max` (i.e., larger in absolute value) the limiter takes action. +The choice of the tolerated exponentiated entropy decrease is a problem-specific parameter +which balances the trade-off between accuracy and stability. +""" +function EntropyBoundedLimiter(; + exp_entropy_decrease_max::RealT = -1.0f-13) where {RealT <: + Real} + @assert exp_entropy_decrease_max<0 "Supplied `exp_entropy_decrease_max` expected to be negative" + EntropyBoundedLimiter{RealT}(exp_entropy_decrease_max) +end + +function (limiter!::EntropyBoundedLimiter)(u_ode, integrator, + semi::AbstractSemidiscretization, + t) + @trixi_timeit timer() "entropy-bounded limiter" begin + @assert :uprev in fieldnames(typeof(integrator)) "EntropyBoundedLimiter requires `uprev` for computation of previous entropy" + + u = wrap_array(u_ode, semi) + u_prev = wrap_array(integrator.uprev, semi) + limiter_entropy_bounded!(u, u_prev, limiter!.exp_entropy_decrease_max, + mesh_equations_solver_cache(semi)...) + end +end + +# Exponentiated entropy change for the thermodynamic entropy (see `entropy_thermodynamic`) +# of an ideal gas with constant gamma. +@inline function exp_entropy_change(p, rho, gamma, exp_entropy_prev) + return p - rho^gamma * exp_entropy_prev +end + +include("entropy_bounded_limiter_1d.jl") +include("entropy_bounded_limiter_2d.jl") +include("entropy_bounded_limiter_3d.jl") +end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter_1d.jl b/src/callbacks_stage/entropy_bounded_limiter_1d.jl new file mode 100644 index 00000000000..217261d39f4 --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter_1d.jl @@ -0,0 +1,72 @@ +# 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 + +function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max, + mesh::AbstractMesh{1}, equations, dg::DGSEM, cache) + @unpack weights = dg.basis + + @threaded for element in eachelement(dg, cache) + # Minimum exponentiated entropy within the current `element` + # of the previous iterate `u_prev` + exp_s_min = typemax(eltype(u_prev)) + + # Determine minimum value for entropy difference + # Can use zero here since d_exp_s is defined as min{0, min_x exp_entropy_change(x)} + d_exp_s_min = zero(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_node_prev = get_node_vars(u_prev, equations, dg, i, element) + + exp_s = exp(entropy_thermodynamic(u_node_prev, equations)) + exp_s_min = min(exp_s_min, exp_s) + + d_exp_s = exp_entropy_change(pressure(u_node, equations), + density(u_node, equations), + equations.gamma, + exp_s) + d_exp_s_min = min(d_exp_s_min, d_exp_s) + end + + # Detect if limiting is necessary. + # Limiting only if entropy DECREASE below a user defined threshold is detected. + d_exp_s_min < exp_entropy_decrease_max || continue + + # Compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, element)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + # Note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + + entropy_change_mean = exp_entropy_change(pressure(u_mean, equations), + density(u_mean, equations), + equations.gamma, + exp_s_min) + + epsilon = d_exp_s_min / (d_exp_s_min - entropy_change_mean) + + # In the derivation of the limiter it is assumed that + # entropy_change_mean >= 0 which would imply epsilon <= 1 (maximum limiting). + # However, this might not always be the case in a simulation as + # we usually do not enforce the corresponding CFL condition from Lemma 3. + # Thus, we clip epsilon at 1. + if epsilon > 1 + epsilon = 1 + end + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + set_node_vars!(u, (1 - epsilon) * u_node + epsilon * u_mean, + equations, dg, i, element) + end + end + + return nothing +end +end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter_2d.jl b/src/callbacks_stage/entropy_bounded_limiter_2d.jl new file mode 100644 index 00000000000..36d70e4f2e5 --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter_2d.jl @@ -0,0 +1,76 @@ +# 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 + +function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max, + mesh::AbstractMesh{2}, equations, dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + @threaded for element in eachelement(dg, cache) + # Minimum exponentiated entropy within the current `element` + # of the previous iterate `u_prev` + exp_s_min = typemax(eltype(u_prev)) + + # Determine minimum value for entropy difference + # Can use zero here since d_exp_s is defined as min{0, min_x exp_entropy_change(x)} + d_exp_s_min = zero(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_node_prev = get_node_vars(u_prev, equations, dg, i, j, element) + + exp_s = exp(entropy_thermodynamic(u_node_prev, equations)) + exp_s_min = min(exp_s_min, exp_s) + + d_exp_s = exp_entropy_change(pressure(u_node, equations), + density(u_node, equations), + equations.gamma, + exp_s) + d_exp_s_min = min(d_exp_s_min, d_exp_s) + end + + # Detect if limiting is necessary. + # Limiting only if entropy DECREASE below a user defined threshold is detected. + d_exp_s_min < exp_entropy_decrease_max || continue + # Compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element)) + total_volume = zero(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh, + i, j, element))) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_mean += u_node * weights[i] * weights[j] * volume_jacobian + total_volume += weights[i] * weights[j] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + + entropy_change_mean = exp_entropy_change(pressure(u_mean, equations), + density(u_mean, equations), + equations.gamma, + exp_s_min) + + epsilon = d_exp_s_min / (d_exp_s_min - entropy_change_mean) + + # In the derivation of the limiter it is assumed that + # entropy_change_mean >= 0 which would imply epsilon <= 1 (maximum limiting). + # However, this might not always be the case in a simulation as + # we usually do not enforce the corresponding CFL condition from Lemma 3. + # Thus, we clip epsilon at 1. + if epsilon > 1 + epsilon = 1 + end + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + set_node_vars!(u, (1 - epsilon) * u_node + epsilon * u_mean, + equations, dg, i, j, element) + end + end + + return nothing +end +end # @muladd diff --git a/src/callbacks_stage/entropy_bounded_limiter_3d.jl b/src/callbacks_stage/entropy_bounded_limiter_3d.jl new file mode 100644 index 00000000000..6f194f0814d --- /dev/null +++ b/src/callbacks_stage/entropy_bounded_limiter_3d.jl @@ -0,0 +1,76 @@ +# 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 + +function limiter_entropy_bounded!(u, u_prev, exp_entropy_decrease_max, + mesh::AbstractMesh{3}, equations, dg::DGSEM, cache) + @unpack weights = dg.basis + @unpack inverse_jacobian = cache.elements + + @threaded for element in eachelement(dg, cache) + # Minimum exponentiated entropy within the current `element` + # of the previous iterate `u_prev` + exp_s_min = typemax(eltype(u_prev)) + + # Determine minimum value for entropy difference + # Can use zero here since d_exp_s is defined as min{0, min_x exp_entropy_change(x)} + d_exp_s_min = zero(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + u_node_prev = get_node_vars(u_prev, equations, dg, i, j, k, element) + + exp_s = exp(entropy_thermodynamic(u_node_prev, equations)) + exp_s_min = min(exp_s_min, exp_s) + + d_exp_s = exp_entropy_change(pressure(u_node, equations), + density(u_node, equations), + equations.gamma, + exp_s) + d_exp_s_min = min(d_exp_s_min, d_exp_s) + end + + # Detect if limiting is necessary. Avoid division by ("near") zero + d_exp_s_min < exp_entropy_decrease_max || continue + + # Compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element)) + total_volume = zero(eltype(u)) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh, + i, j, k, element))) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian + total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian + end + # normalize with the total volume + u_mean = u_mean / total_volume + + entropy_change_mean = exp_entropy_change(pressure(u_mean, equations), + density(u_mean, equations), + equations.gamma, + exp_s_min) + + epsilon = d_exp_s_min / (d_exp_s_min - entropy_change_mean) + + # In the derivation of the limiter it is assumed that + # entropy_change_mean >= 0 which would imply epsilon <= 1 (maximum limiting). + # However, this might not always be the case in a simulation as + # we usually do not enforce the corresponding CFL condition from Lemma 3. + # Thus, we clip epsilon at 1. + if epsilon > 1 + epsilon = 1 + end + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + set_node_vars!(u, (1 - epsilon) * u_node + epsilon * u_mean, + equations, dg, i, j, k, element) + end + end + + return nothing +end +end # @muladd diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 55f7ab76023..5cec16300f8 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -594,6 +594,43 @@ end end end +@trixi_testset "elixir_mhd_amr_entropy_bounded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_amr_entropy_bounded.jl"), + l2=[ + 0.005430176785094096, + 0.006185803468926062, + 0.012158513265762224, + 0.006185144232789619, + 0.03509140423905665, + 0.004968215426326584, + 0.006553519141867704, + 0.005008885124643863, + 5.165777182726578e-6 + ], + linf=[ + 0.1864317840224794, + 0.2041246899193812, + 0.36992946717578445, + 0.2327158690965257, + 1.0368624176126007, + 0.1846308291826353, + 0.2062255411778191, + 0.18955666546331185, + 0.0005208969502913304 + ], + tspan=(0.0, 0.04), + coverage_override=(maxiters = 6, initial_refinement_level = 1, + base_level = 1, max_level = 2)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_linearizedeuler_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_convergence.jl"), diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 74908552521..1386c5f5c4d 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -409,6 +409,21 @@ end end end +@trixi_testset "elixir_euler_blast_wave_entropy_bounded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave_entropy_bounded.jl"), + l2=[0.9689207881108007, 0.1617708899929322, 1.3847895715669456], + linf=[2.95591859210077, 0.3135723412205586, 2.3871554358655365]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "test_quasi_1D_entropy" begin a = 0.9 u_1D = SVector(1.1, 0.2, 2.1) diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 6685c289a9b..479a07b97a2 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -732,6 +732,34 @@ end end end +@trixi_testset "elixir_euler_colliding_flow_amr_entropy_bounded.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_colliding_flow_amr_entropy_bounded.jl"), + l2=[ + 0.04120588220419942, + 0.09868046588789257, + 7.446553779796626e-7, + 5.5746513268066105 + ], + linf=[ + 0.3478655090378702, + 0.864011305195807, + 5.419432288048388e-5, + 47.284459667934684 + ], + tspan=(0.0, 1.0), + dt=2.5e-2, adaptive=false, + coverage_override=(maxiters = 10^3,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_astro_jet_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_astro_jet_amr.jl"), l2=[ From ddf24ae02c0e12e619a1f3f2cdaa4caeea735833 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 15 Oct 2024 08:55:57 +0200 Subject: [PATCH 11/13] 1D --- test/test_parabolic_1d.jl | 48 +++++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 062e6363a2f..db3b24397a0 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -67,12 +67,12 @@ end l2=[ 0.0001133835907077494, 6.226282245610444e-5, - 0.0002820171699999139 + 0.0002820171699999139, ], linf=[ 0.0006255102377159538, 0.00036195501456059986, - 0.0016147729485886941 + 0.0016147729485886941, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -94,12 +94,12 @@ end l2=[ 0.00011310615871043463, 6.216495207074201e-5, - 0.00028195843110817814 + 0.00028195843110817814, ], linf=[ 0.0006240837363233886, 0.0003616694320713876, - 0.0016147339542413874 + 0.0016147339542413874, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -116,13 +116,13 @@ end "elixir_navierstokes_convergence_walls.jl"), l2=[ 0.00047023310868269237, - 0.00032181736027057234, - 0.0014966266486095025 + 0.00032181735314911885, + 0.0014966266120305365, ], linf=[ - 0.002996375101363302, - 0.0028639041695096433, - 0.012691132694550689 + 0.0029963751716199916, + 0.0028639040576996177, + 0.012691132369322844, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -144,12 +144,12 @@ end l2=[ 0.0004608500483647771, 0.00032431091222851285, - 0.0015159733360626845 + 0.0015159733360626845, ], linf=[ 0.002754803146635787, - 0.0028567713744625124, - 0.012941793784197131 + 0.0028567714697580906, + 0.012941794048176192, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -168,14 +168,14 @@ end mu = mu(), Prandtl = prandtl_number()), l2=[ - 2.5278845598681636e-5, - 2.5540145802666872e-5, - 0.0001211867535580826 + 2.5279971242617625e-5, + 2.5543758853591306e-5, + 0.00012119767035310137, ], linf=[ - 0.0001466387202588848, - 0.00019422419092429135, - 0.0009556449835592673 + 0.00014663941066395125, + 0.0001942196764966836, + 0.0009556610409156008, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -195,14 +195,14 @@ end Prandtl = prandtl_number(), gradient_variables = GradientVariablesEntropy()), l2=[ - 2.4593521887223632e-5, - 2.3928212900127102e-5, - 0.00011252332663824173 + 2.4594499734628152e-5, + 2.393055300564619e-5, + 0.00011253402930582066, ], linf=[ - 0.00011850494672183132, - 0.00018987676556476442, - 0.0009597461727750556 + 0.00011850830177762006, + 0.00018988219246484224, + 0.0009598012069105266, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 2000103889162a5608700fc931c6f7f8dc047329 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 15 Oct 2024 08:56:41 +0200 Subject: [PATCH 12/13] fmt --- test/test_parabolic_1d.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index db3b24397a0..0069a81b8ec 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -67,12 +67,12 @@ end l2=[ 0.0001133835907077494, 6.226282245610444e-5, - 0.0002820171699999139, + 0.0002820171699999139 ], linf=[ 0.0006255102377159538, 0.00036195501456059986, - 0.0016147729485886941, + 0.0016147729485886941 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -94,12 +94,12 @@ end l2=[ 0.00011310615871043463, 6.216495207074201e-5, - 0.00028195843110817814, + 0.00028195843110817814 ], linf=[ 0.0006240837363233886, 0.0003616694320713876, - 0.0016147339542413874, + 0.0016147339542413874 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -117,12 +117,12 @@ end l2=[ 0.00047023310868269237, 0.00032181735314911885, - 0.0014966266120305365, + 0.0014966266120305365 ], linf=[ 0.0029963751716199916, 0.0028639040576996177, - 0.012691132369322844, + 0.012691132369322844 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -144,12 +144,12 @@ end l2=[ 0.0004608500483647771, 0.00032431091222851285, - 0.0015159733360626845, + 0.0015159733360626845 ], linf=[ 0.002754803146635787, 0.0028567714697580906, - 0.012941794048176192, + 0.012941794048176192 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -170,12 +170,12 @@ end l2=[ 2.5279971242617625e-5, 2.5543758853591306e-5, - 0.00012119767035310137, + 0.00012119767035310137 ], linf=[ 0.00014663941066395125, 0.0001942196764966836, - 0.0009556610409156008, + 0.0009556610409156008 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -197,12 +197,12 @@ end l2=[ 2.4594499734628152e-5, 2.393055300564619e-5, - 0.00011253402930582066, + 0.00011253402930582066 ], linf=[ 0.00011850830177762006, 0.00018988219246484224, - 0.0009598012069105266, + 0.0009598012069105266 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 689296f62a66e592c8823725674dbe8223a9950c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 15 Oct 2024 10:48:19 +0200 Subject: [PATCH 13/13] test vals --- test/test_parabolic_1d.jl | 4 ++-- test/test_parabolic_2d.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 0069a81b8ec..3d59dba4c4f 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -148,8 +148,8 @@ end ], linf=[ 0.002754803146635787, - 0.0028567714697580906, - 0.012941794048176192 + 0.002856771361794534, + 0.012941793749087438 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index ceefb65e99b..10afa608665 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -579,8 +579,8 @@ end @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_nonperiodic_amr.jl"), tspan=(0.0, 0.01), - l2=[0.007934195641974433], - linf=[0.11030265194954081]) + l2=[0.007933686773471083], + linf=[0.11029277682717009]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let