From 1d7b2561ff656cb9b8c826ee8b7ee389a3370917 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Wed, 9 Oct 2024 16:35:07 +0200 Subject: [PATCH 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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