From 034b1587807b5098e64043b4711604b481ba2186 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 13 Dec 2023 09:15:37 +0100 Subject: [PATCH 01/17] bump lower compat bound of TimerOutputs.jl to v0.5.7 (#1772) * bump lower compat bound of TimerOutputs.jl to v0.5.7 * format --- Project.toml | 2 +- .../semidiscretization_euler_acoustics.jl | 10 +++++----- .../semidiscretization_euler_gravity.jl | 10 +++++----- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index ace71a50d8b..65074c5a1ca 100644 --- a/Project.toml +++ b/Project.toml @@ -85,7 +85,7 @@ StrideArrays = "0.1.18" StructArrays = "0.6" SummationByPartsOperators = "0.5.41" T8code = "0.4.3" -TimerOutputs = "0.5" +TimerOutputs = "0.5.7" Triangulate = "2.0" TriplotBase = "0.1" TriplotRecipes = "0.1" diff --git a/src/semidiscretization/semidiscretization_euler_acoustics.jl b/src/semidiscretization/semidiscretization_euler_acoustics.jl index e49fe81177a..173523ff892 100644 --- a/src/semidiscretization/semidiscretization_euler_acoustics.jl +++ b/src/semidiscretization/semidiscretization_euler_acoustics.jl @@ -51,11 +51,11 @@ function SemidiscretizationEulerAcoustics(semi_acoustics::SemiAcoustics, semi_euler::SemiEuler; source_region = x -> true, weights = x -> 1.0) where - {Mesh, - SemiAcoustics <: - SemidiscretizationHyperbolic{Mesh, <:AbstractAcousticPerturbationEquations}, - SemiEuler <: - SemidiscretizationHyperbolic{Mesh, <:AbstractCompressibleEulerEquations}} + {Mesh, + SemiAcoustics <: + SemidiscretizationHyperbolic{Mesh, <:AbstractAcousticPerturbationEquations}, + SemiEuler <: + SemidiscretizationHyperbolic{Mesh, <:AbstractCompressibleEulerEquations}} cache = create_cache(SemidiscretizationEulerAcoustics, source_region, weights, mesh_equations_solver_cache(semi_acoustics)...) diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index a9a60a4ff04..4201344df80 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -117,11 +117,11 @@ Construct a semidiscretization of the compressible Euler equations with self-gra function SemidiscretizationEulerGravity(semi_euler::SemiEuler, semi_gravity::SemiGravity, parameters) where - {Mesh, - SemiEuler <: - SemidiscretizationHyperbolic{Mesh, <:AbstractCompressibleEulerEquations}, - SemiGravity <: - SemidiscretizationHyperbolic{Mesh, <:AbstractHyperbolicDiffusionEquations}} + {Mesh, + SemiEuler <: + SemidiscretizationHyperbolic{Mesh, <:AbstractCompressibleEulerEquations}, + SemiGravity <: + SemidiscretizationHyperbolic{Mesh, <:AbstractHyperbolicDiffusionEquations}} u_ode = compute_coefficients(zero(real(semi_gravity)), semi_gravity) du_ode = similar(u_ode) u_tmp1_ode = similar(u_ode) From 6acc09eaa815ad8cbcc0adaf5d160b3299daf071 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 13 Dec 2023 14:05:14 +0100 Subject: [PATCH 02/17] set version to v0.6.4 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 65074c5a1ca..8be83c7924f 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.6.4-pre" +version = "0.6.4" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 766c7981fdfcbfcbb034888bb627a8f4668259eb Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 13 Dec 2023 14:05:26 +0100 Subject: [PATCH 03/17] set development version to v0.6.5-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8be83c7924f..539dafc3034 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.6.4" +version = "0.6.5-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 14796ea607ad49f49e2258dac87ff30a6530595d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 16:40:21 +0100 Subject: [PATCH 04/17] Central SBP finite difference solver for `UnstructuredMesh2D` (#1773) * containers and kernels for FDSBP solver on UnstructuredMesh2D * add elixirs and corresponding tests * apply formatter to new and edited files * add advection equation test to up coverage * Apply suggestions from code review Co-authored-by: Hendrik Ranocha * update variable name to N --------- Co-authored-by: Hendrik Ranocha --- .../elixir_advection_basic.jl | 69 ++++++ .../elixir_euler_free_stream.jl | 77 ++++++ .../elixir_euler_source_terms.jl | 65 ++++++ src/solvers/dg.jl | 8 +- .../dgsem_unstructured/containers_2d.jl | 15 +- src/solvers/fdsbp_tree/fdsbp_2d.jl | 2 +- .../fdsbp_unstructured/containers_2d.jl | 124 ++++++++++ src/solvers/fdsbp_unstructured/fdsbp.jl | 14 ++ src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 219 ++++++++++++++++++ test/test_unstructured_2d.jl | 61 +++++ 10 files changed, 643 insertions(+), 11 deletions(-) create mode 100644 examples/unstructured_2d_fdsbp/elixir_advection_basic.jl create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl create mode 100644 src/solvers/fdsbp_unstructured/containers_2d.jl create mode 100644 src/solvers/fdsbp_unstructured/fdsbp.jl create mode 100644 src/solvers/fdsbp_unstructured/fdsbp_2d.jl diff --git a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl new file mode 100644 index 00000000000..c181203e7a4 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl @@ -0,0 +1,69 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +############################################################################### +# Get the FDSBP approximation operator + +D_SBP = derivative_operator(SummationByPartsOperators.MattssonAlmquistVanDerWeide2018Accurate(), + derivative_order = 1, accuracy_order = 4, + xmin = -1.0, xmax = 1.0, N = 15) +solver = FDSBP(D_SBP, + surface_integral = SurfaceIntegralStrongForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralStrongForm()) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +############################################################################### +# create the semidiscretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0, 1.0)) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 100, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 1.6) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl new file mode 100644 index 00000000000..7ada50c0c65 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl @@ -0,0 +1,77 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +# Free-stream initial condition +initial_condition = initial_condition_constant + +# Boundary conditions for free-stream testing +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:Body => boundary_condition_free_stream, + :Button1 => boundary_condition_free_stream, + :Button2 => boundary_condition_free_stream, + :Eye1 => boundary_condition_free_stream, + :Eye2 => boundary_condition_free_stream, + :Smile => boundary_condition_free_stream, + :Bowtie => boundary_condition_free_stream) + +############################################################################### +# Get the FDSBP approximation space + +D_SBP = derivative_operator(SummationByPartsOperators.MattssonAlmquistVanDerWeide2018Accurate(), + derivative_order = 1, accuracy_order = 4, + xmin = -1.0, xmax = 1.0, N = 12) +solver = FDSBP(D_SBP, + surface_integral = SurfaceIntegralStrongForm(flux_hll), + volume_integral = VolumeIntegralStrongForm()) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, + alive_callback, save_solution) + +############################################################################### +# run the simulation + +# set small tolerances for the free-stream preservation test +sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12, + save_everystep = false, callback = callbacks) +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl new file mode 100644 index 00000000000..edcd221bf59 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl @@ -0,0 +1,65 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the FDSBP approximation operator + +D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(), + derivative_order = 1, accuracy_order = 4, + xmin = -1.0, xmax = 1.0, N = 10) +solver = FDSBP(D_SBP, + surface_integral = SurfaceIntegralStrongForm(flux_lax_friedrichs), + volume_integral = VolumeIntegralStrongForm()) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, + alive_callback, save_solution) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(), abstol = 1.0e-9, reltol = 1.0e-9, + save_everystep = false, callback = callbacks) +summary_callback() # print the timer summary diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 9e5ebc7f9b5..9b61df62cc3 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -41,8 +41,8 @@ standard textbooks. Applications [doi: 10.1007/978-0-387-72067-8](https://doi.org/10.1007/978-0-387-72067-8) -`VolumeIntegralWeakForm()` is only implemented for conserved terms as -non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +`VolumeIntegralWeakForm()` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see [`VolumeIntegralFluxDifferencing`](@ref). This treatment is required to achieve, e.g., entropy-stability or well-balancedness. """ @@ -415,7 +415,8 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG) summary_line(io, "surface integral", dg.surface_integral |> typeof |> nameof) show(increment_indent(io), mime, dg.surface_integral) summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof) - if !(dg.volume_integral isa VolumeIntegralWeakForm) + if !(dg.volume_integral isa VolumeIntegralWeakForm) && + !(dg.volume_integral isa VolumeIntegralStrongForm) show(increment_indent(io), mime, dg.volume_integral) end summary_footer(io) @@ -598,6 +599,7 @@ include("dgsem/dgsem.jl") # and boundary conditions weakly. Thus, these methods can re-use a lot of # functionality implemented for DGSEM. include("fdsbp_tree/fdsbp.jl") +include("fdsbp_unstructured/fdsbp.jl") function allocate_coefficients(mesh::AbstractMesh, equations, dg::DG, cache) # We must allocate a `Vector` in order to be able to `resize!` it (AMR). diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl index 13eeaeabffb..f51dd09801b 100644 --- a/src/solvers/dgsem_unstructured/containers_2d.jl +++ b/src/solvers/dgsem_unstructured/containers_2d.jl @@ -45,7 +45,7 @@ end eachelement(elements::UnstructuredElementContainer2D) Return an iterator over the indices that specify the location in relevant data structures -for the elements in `elements`. +for the elements in `elements`. In particular, not the elements themselves are returned. """ @inline function eachelement(elements::UnstructuredElementContainer2D) @@ -84,24 +84,25 @@ function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis) # loop through elements and call the correct constructor based on whether the element is curved for element in eachelement(elements) if mesh.element_is_curved[element] - init_element!(elements, element, basis.nodes, + init_element!(elements, element, basis, view(mesh.surface_curves, :, element)) else # straight sided element for i in 1:4, j in 1:2 # pull the (x,y) values of these corners out of the global corners array four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]] end - init_element!(elements, element, basis.nodes, four_corners) + init_element!(elements, element, basis, four_corners) end end end # initialize all the values in the container of a general element (either straight sided or curved) -function init_element!(elements, element, nodes, corners_or_surface_curves) - calc_node_coordinates!(elements.node_coordinates, element, nodes, +function init_element!(elements, element, basis::LobattoLegendreBasis, + corners_or_surface_curves) + calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), corners_or_surface_curves) - calc_metric_terms!(elements.jacobian_matrix, element, nodes, + calc_metric_terms!(elements.jacobian_matrix, element, get_nodes(basis), corners_or_surface_curves) calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix) @@ -109,7 +110,7 @@ function init_element!(elements, element, nodes, corners_or_surface_curves) calc_contravariant_vectors!(elements.contravariant_vectors, element, elements.jacobian_matrix) - calc_normal_directions!(elements.normal_directions, element, nodes, + calc_normal_directions!(elements.normal_directions, element, get_nodes(basis), corners_or_surface_curves) return elements diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl index beff605629a..09d18cecd75 100644 --- a/src/solvers/fdsbp_tree/fdsbp_2d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_2d.jl @@ -9,7 +9,7 @@ #! format: noindent # 2D caches -function create_cache(mesh::TreeMesh{2}, equations, +function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations, volume_integral::VolumeIntegralStrongForm, dg, uEltype) prototype = Array{SVector{nvariables(equations), uEltype}, ndims(mesh)}(undef, ntuple(_ -> nnodes(dg), diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl new file mode 100644 index 00000000000..3857c2d8a20 --- /dev/null +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -0,0 +1,124 @@ +# !!! warning "Experimental implementation (curvilinear FDSBP)" +# This is an experimental feature and may change in future releases. + +# 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 + +# initialize all the values in the container of a general FD block (either straight sided or curved) +# OBS! Requires the SBP derivative matrix in order to compute metric terms that are free-stream preserving +function init_element!(elements, element, basis::AbstractDerivativeOperator, + corners_or_surface_curves) + calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis), + corners_or_surface_curves) + + calc_metric_terms!(elements.jacobian_matrix, element, basis, + elements.node_coordinates) + + calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix) + + calc_contravariant_vectors!(elements.contravariant_vectors, element, + elements.jacobian_matrix) + + calc_normal_directions!(elements.normal_directions, element, + elements.jacobian_matrix) + + return elements +end + +# construct the metric terms for a FDSBP element "block". Directly use the derivative matrix +# applied to the node coordinates. +# TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices +function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator, + node_coordinates) + + # storage format: + # jacobian_matrix[1,1,:,:,:] <- X_xi + # jacobian_matrix[1,2,:,:,:] <- X_eta + # jacobian_matrix[2,1,:,:,:] <- Y_xi + # jacobian_matrix[2,2,:,:,:] <- Y_eta + + # Compute the xi derivatives by applying D on the left + # This is basically the same as + # jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element] + # but uses only matrix-vector products instead of a matrix-matrix product. + for j in eachnode(D_SBP) + mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP, + view(node_coordinates, 1, :, j, element)) + end + # jacobian_matrix[2, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[2, :, :, element] + for j in eachnode(D_SBP) + mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP, + view(node_coordinates, 2, :, j, element)) + end + + # Compute the eta derivatives by applying transpose of D on the right + # jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)' + for i in eachnode(D_SBP) + mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP, + view(node_coordinates, 1, i, :, element)) + end + # jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(D_SBP)' + for i in eachnode(D_SBP) + mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP, + view(node_coordinates, 2, i, :, element)) + end + + return jacobian_matrix +end + +# construct the normal direction vectors (but not actually normalized) for a curved sided FDSBP element "block" +# normalization occurs on the fly during the surface flux computation +# OBS! This assumes that the boundary points are included. +function calc_normal_directions!(normal_directions, element, jacobian_matrix) + + # normal directions on the boundary for the left (local side 4) and right (local side 2) + N = size(jacobian_matrix, 4) + for j in 1:N + # +x side or side 2 in the local indexing + X_xi = jacobian_matrix[1, 1, N, j, element] + X_eta = jacobian_matrix[1, 2, N, j, element] + Y_xi = jacobian_matrix[2, 1, N, j, element] + Y_eta = jacobian_matrix[2, 2, N, j, element] + Jtemp = X_xi * Y_eta - X_eta * Y_xi + normal_directions[1, j, 2, element] = sign(Jtemp) * (Y_eta) + normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta) + + # -x side or side 4 in the local indexing + X_xi = jacobian_matrix[1, 1, 1, j, element] + X_eta = jacobian_matrix[1, 2, 1, j, element] + Y_xi = jacobian_matrix[2, 1, 1, j, element] + Y_eta = jacobian_matrix[2, 2, 1, j, element] + Jtemp = X_xi * Y_eta - X_eta * Y_xi + normal_directions[1, j, 4, element] = -sign(Jtemp) * (Y_eta) + normal_directions[2, j, 4, element] = -sign(Jtemp) * (-X_eta) + end + + # normal directions on the boundary for the top (local side 3) and bottom (local side 1) + N = size(jacobian_matrix, 3) + for i in 1:N + # -y side or side 1 in the local indexing + X_xi = jacobian_matrix[1, 1, i, 1, element] + X_eta = jacobian_matrix[1, 2, i, 1, element] + Y_xi = jacobian_matrix[2, 1, i, 1, element] + Y_eta = jacobian_matrix[2, 2, i, 1, element] + Jtemp = X_xi * Y_eta - X_eta * Y_xi + normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi) + normal_directions[2, i, 1, element] = -sign(Jtemp) * (X_xi) + + # +y side or side 3 in the local indexing + X_xi = jacobian_matrix[1, 1, i, N, element] + X_eta = jacobian_matrix[1, 2, i, N, element] + Y_xi = jacobian_matrix[2, 1, i, N, element] + Y_eta = jacobian_matrix[2, 2, i, N, element] + Jtemp = X_xi * Y_eta - X_eta * Y_xi + normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi) + normal_directions[2, i, 3, element] = sign(Jtemp) * (X_xi) + end + + return normal_directions +end +end # @muladd diff --git a/src/solvers/fdsbp_unstructured/fdsbp.jl b/src/solvers/fdsbp_unstructured/fdsbp.jl new file mode 100644 index 00000000000..dee9776abb7 --- /dev/null +++ b/src/solvers/fdsbp_unstructured/fdsbp.jl @@ -0,0 +1,14 @@ +# !!! warning "Experimental implementation (curvilinear FDSBP)" +# This is an experimental feature and may change in future releases. + +# 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 + +# dimension specific curvilinear implementations and data structures +include("containers_2d.jl") +include("fdsbp_2d.jl") +end # @muladd diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl new file mode 100644 index 00000000000..b459f4c42cc --- /dev/null +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -0,0 +1,219 @@ +# !!! warning "Experimental implementation (curvilinear FDSBP)" +# This is an experimental feature and may change in future releases. + +# 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 + +# 2D unstructured cache +function create_cache(mesh::UnstructuredMesh2D, equations, dg::FDSBP, RealT, uEltype) + elements = init_elements(mesh, equations, dg.basis, RealT, uEltype) + + interfaces = init_interfaces(mesh, elements) + + boundaries = init_boundaries(mesh, elements) + + cache = (; elements, interfaces, boundaries) + + # Add specialized parts of the cache required to for efficient flux computations + cache = (; cache..., + create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) + + return cache +end + +# TODO: FD; Upwind versions of surface / volume integral + +# 2D volume integral contributions for `VolumeIntegralStrongForm` +# OBS! This is the standard (not de-aliased) form of the volume integral. +# So it is not provably stable for variable coefficients due to the the metric terms. +@inline function calc_volume_integral!(du, u, + mesh::UnstructuredMesh2D, + nonconservative_terms::False, equations, + volume_integral::VolumeIntegralStrongForm, + dg::FDSBP, cache) + D = dg.basis # SBP derivative operator + @unpack f_threaded = cache + @unpack contravariant_vectors = cache.elements + + # SBP operators from SummationByPartsOperators.jl implement the basic interface + # of matrix-vector multiplication. Thus, we pass an "array of structures", + # packing all variables per node in an `SVector`. + if nvariables(equations) == 1 + # `reinterpret(reshape, ...)` removes the leading dimension only if more + # than one variable is used. + u_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, + du), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + else + u_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u) + du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, + du) + end + + # Use the tensor product structure to compute the discrete derivatives of + # the contravariant fluxes line-by-line and add them to `du` for each element. + @threaded for element in eachelement(dg, cache) + f_element = f_threaded[Threads.threadid()] + u_element = view(u_vectors, :, :, element) + + # x direction + for j in eachnode(dg) + for i in eachnode(dg) + Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + f_element[i, j] = flux(u_element[i, j], Ja1, equations) + end + mul!(view(du_vectors, :, j, element), D, view(f_element, :, j), + one(eltype(du)), one(eltype(du))) + end + + # y direction + for i in eachnode(dg) + for j in eachnode(dg) + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + f_element[i, j] = flux(u_element[i, j], Ja2, equations) + end + mul!(view(du_vectors, i, :, element), D, view(f_element, i, :), + one(eltype(du)), one(eltype(du))) + end + end + + return nothing +end + +# Note! The local side numbering for the unstructured quadrilateral element implementation differs +# from the structured TreeMesh or StructuredMesh local side numbering: +# +# TreeMesh/StructuredMesh sides versus UnstructuredMesh sides +# 4 3 +# ----------------- ----------------- +# | | | | +# | ^ eta | | ^ eta | +# 1 | | | 2 4 | | | 2 +# | | | | | | +# | ---> xi | | ---> xi | +# ----------------- ----------------- +# 3 1 +# Therefore, we require a different surface integral routine here despite their similar structure. +# Also, the normal directions are already outward pointing for `UnstructuredMesh2D` so all the +# surface contributions are added. +function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, + equations, surface_integral::SurfaceIntegralStrongForm, + dg::DG, cache) + inv_weight_left = inv(left_boundary_weight(dg.basis)) + inv_weight_right = inv(right_boundary_weight(dg.basis)) + @unpack normal_directions, surface_flux_values = cache.elements + + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + # surface at -x + u_node = get_node_vars(u, equations, dg, 1, l, element) + # compute internal flux in normal direction on side 4 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) + multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), + equations, dg, 1, l, element) + + # surface at +x + u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) + # compute internal flux in normal direction on side 2 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) + multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), + equations, dg, nnodes(dg), l, element) + + # surface at -y + u_node = get_node_vars(u, equations, dg, l, 1, element) + # compute internal flux in normal direction on side 1 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) + multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), + equations, dg, l, 1, element) + + # surface at +y + u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) + # compute internal flux in normal direction on side 3 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) + multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), + equations, dg, l, nnodes(dg), element) + end + end + + return nothing +end + +# AnalysisCallback +function integrate_via_indices(func::Func, u, + mesh::UnstructuredMesh2D, equations, + dg::FDSBP, cache, args...; normalize = true) where {Func} + # TODO: FD. This is rather inefficient right now and allocates... + weights = diag(SummationByPartsOperators.mass_matrix(dg.basis)) + + # Initialize integral with zeros of the right shape + integral = zero(func(u, 1, 1, 1, equations, dg, args...)) + total_volume = zero(real(mesh)) + + # Use quadrature to numerically integrate over entire domain + for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element])) + integral += volume_jacobian * weights[i] * weights[j] * + func(u, i, j, element, equations, dg, args...) + total_volume += volume_jacobian * weights[i] * weights[j] + end + end + + # Normalize with total volume + if normalize + integral = integral / total_volume + end + + return integral +end + +function calc_error_norms(func, u, t, analyzer, + mesh::UnstructuredMesh2D, equations, initial_condition, + dg::FDSBP, cache, cache_analysis) + # TODO: FD. This is rather inefficient right now and allocates... + weights = diag(SummationByPartsOperators.mass_matrix(dg.basis)) + @unpack node_coordinates, inverse_jacobian = cache.elements + + # Set up data structures + l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations)) + linf_error = copy(l2_error) + total_volume = zero(real(mesh)) + + # Iterate over all elements for error calculations + for element in eachelement(dg, cache) + for j in eachnode(analyzer), i in eachnode(analyzer) + volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element])) + u_exact = initial_condition(get_node_coords(node_coordinates, equations, dg, + i, j, element), t, equations) + diff = func(u_exact, equations) - + func(get_node_vars(u, equations, dg, i, j, element), equations) + l2_error += diff .^ 2 * (weights[i] * weights[j] * volume_jacobian) + linf_error = @. max(linf_error, abs(diff)) + total_volume += weights[i] * weights[j] * volume_jacobian + end + end + + # For L2 error, divide by total volume + l2_error = @. sqrt(l2_error / total_volume) + + return l2_error, linf_error +end +end # @muladd diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 5341d86a7d1..139b423ead1 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -664,6 +664,67 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +# TODO: FD; for now put the unstructured tests for the 2D FDSBP here. +@trixi_testset "FDSBP (central): elixir_advection_basic.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_advection_basic.jl"), + l2=[0.0001105211407319266], + linf=[0.0004199363734466166]) + # 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 "FDSBP (central): elixir_euler_source_terms.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_source_terms.jl"), + l2=[8.155544666380138e-5, + 0.0001477863788446318, + 0.00014778637884460072, + 0.00045584189984542687], + linf=[0.0002670775876922882, + 0.0005683064706873964, + 0.0005683064706762941, + 0.0017770812025146299], + tspan=(0.0, 0.05)) + # 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 "FDSBP (central): elixir_euler_free_stream.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_free_stream.jl"), + l2=[5.4329175009362306e-14, + 1.0066867437607972e-13, + 6.889210012578449e-14, + 1.568290814572709e-13], + linf=[5.963762816918461e-10, + 5.08869890669672e-11, + 1.1581377523661729e-10, + 4.61017890529547e-11], + tspan=(0.0, 0.1), + atol=1.0e-11) + # 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 end # Clean up afterwards: delete Trixi.jl output directory From c7c4cf7c827177fa3e32942aa58b316468ac0295 Mon Sep 17 00:00:00 2001 From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com> Date: Thu, 14 Dec 2023 03:54:19 -0600 Subject: [PATCH 05/17] Compressible Euler Quasi-1D (#1757) * implementation of quasi 1d compressible Euler Equation * added example elixir for quasi 1d compressible Euler equations * added example elixir with a discontinuous initial condition * including and exported CompressibleEulerEquationsQuasi1D * formatting * added entropy conservative test * fixing spelling * formatting * Update examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * updated test_tree_1d_euler.jl and formatting * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * adding consistency check for flux * Update compressible_euler_quasi_1d.jl * formatting * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update compressible_euler_quasi_1d and test_tree_1d_euler * formatting * Update examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Daniel Doehring * update boundary condition slip wall * update compressible_euler_quasi_1d.jl * Update src/equations/compressible_euler_quasi_1d.jl * Update src/equations/compressible_euler_quasi_1d.jl * remove boundary_condition_slip_wall --------- Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> Co-authored-by: Daniel Doehring Co-authored-by: Daniel Doehring Co-authored-by: Hendrik Ranocha --- .../elixir_euler_quasi_1d_discontinuous.jl | 85 +++++ .../tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 73 ++++ .../elixir_euler_quasi_1d_source_terms.jl | 60 ++++ src/Trixi.jl | 1 + src/equations/compressible_euler_quasi_1d.jl | 328 ++++++++++++++++++ src/equations/equations.jl | 1 + test/test_tree_1d_euler.jl | 73 ++++ test/test_unit.jl | 13 + 8 files changed, 634 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl create mode 100644 examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl create mode 100644 examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl create mode 100644 src/equations/compressible_euler_quasi_1d.jl diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl new file mode 100644 index 00000000000..cc4535be028 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl @@ -0,0 +1,85 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the quasi 1d compressible Euler equations +# See Chan et al. https://doi.org/10.48550/arXiv.2307.12089 for details + +equations = CompressibleEulerEquationsQuasi1D(1.4) + +""" + initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D) + +A discontinuous initial condition taken from +- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) + High order entropy stable schemes for the quasi-one-dimensional + shallow water and compressible Euler equations + [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089) +""" +function initial_condition_discontinuity(x, t, + equations::CompressibleEulerEquationsQuasi1D) + rho = (x[1] < 0) ? 3.4718 : 2.0 + v1 = (x[1] < 0) ? -2.5923 : -3.0 + p = (x[1] < 0) ? 5.7118 : 2.639 + a = (x[1] < 0) ? 1.0 : 1.5 + + return prim2cons(SVector(rho, v1, p, a), equations) +end + +initial_condition = initial_condition_discontinuity + +surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal) +volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal) + +basis = LobattoLegendreBasis(3) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + 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 = (-1.0,) +coordinates_max = (1.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, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl new file mode 100644 index 00000000000..ae1b2b24b62 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl @@ -0,0 +1,73 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the quasi 1d compressible Euler equations with a discontinuous nozzle width function. +# See Chan et al. https://doi.org/10.48550/arXiv.2307.12089 for details + +equations = CompressibleEulerEquationsQuasi1D(1.4) + +# Setup a truly discontinuous density function and nozzle width for +# this academic testcase of entropy conservation. The errors from the analysis +# callback are not important but the entropy error for this test case +# `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=6`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D) + v1 = 0.1 + rho = 2.0 + 0.1 * x[1] + p = 3.0 + a = 2.0 + x[1] + + return prim2cons(SVector(rho, v1, p, a), equations) +end + +initial_condition = initial_condition_ec + +surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal) +volume_flux = surface_flux +solver = DGSEM(polydeg = 4, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0,) +coordinates_max = (1.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, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 + +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl new file mode 100644 index 00000000000..91bb1ba6e8c --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl @@ -0,0 +1,60 @@ +using OrdinaryDiffEq +using Trixi +using ForwardDiff + +############################################################################### +# Semidiscretization of the quasi 1d compressible Euler equations +# See Chan et al. https://doi.org/10.48550/arXiv.2307.12089 for details + +equations = CompressibleEulerEquationsQuasi1D(1.4) + +initial_condition = initial_condition_convergence_test + +surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal) +volume_flux = surface_flux +solver = DGSEM(polydeg = 4, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index b8110cf5bdd..e7b849e2642 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -139,6 +139,7 @@ export AcousticPerturbationEquations2D, CompressibleEulerEquations3D, CompressibleEulerMulticomponentEquations1D, CompressibleEulerMulticomponentEquations2D, + CompressibleEulerEquationsQuasi1D, IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl new file mode 100644 index 00000000000..0a543277ee4 --- /dev/null +++ b/src/equations/compressible_euler_quasi_1d.jl @@ -0,0 +1,328 @@ +# 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 + +@doc raw""" + CompressibleEulerEquationsQuasi1D(gamma) + +The quasi-1d compressible Euler equations (see Chan et al. [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089) for details) +```math +\frac{\partial}{\partial t} +\begin{pmatrix} +a \rho \\ a \rho v_1 \\ a e +\end{pmatrix} ++ +\frac{\partial}{\partial x} +\begin{pmatrix} +a \rho v_1 \\ a \rho v_1^2 \\ a v_1 (e +p) +\end{pmatrix} ++ +a \frac{\partial}{\partial x} +\begin{pmatrix} +0 \\ p \\ 0 +\end{pmatrix} += +\begin{pmatrix} +0 \\ 0 \\ 0 +\end{pmatrix} +``` +for an ideal gas with ratio of specific heats `gamma` in one space dimension. +Here, ``\rho`` is the density, ``v_1`` the velocity, ``e`` the specific total energy **rather than** specific internal energy, +``a`` the (possibly) variable nozzle width, and +```math +p = (\gamma - 1) \left( e - \frac{1}{2} \rho v_1^2 \right) +``` +the pressure. + +The nozzle width function ``a(x)`` is set inside the initial condition routine +for a particular problem setup. To test the conservative form of the compressible Euler equations one can set the +nozzle width variable ``a`` to one. + +In addition to the unknowns, Trixi.jl currently stores the nozzle width values at the approximation points +despite being fixed in time. +This affects the implementation and use of these equations in various ways: +* The flux values corresponding to the nozzle width must be zero. +* The nozzle width values must be included when defining initial conditions, boundary conditions or + source terms. +* [`AnalysisCallback`](@ref) analyzes this variable. +* Trixi.jl's visualization tools will visualize the nozzle width by default. +""" +struct CompressibleEulerEquationsQuasi1D{RealT <: Real} <: + AbstractCompressibleEulerEquations{1, 4} + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + function CompressibleEulerEquationsQuasi1D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +have_nonconservative_terms(::CompressibleEulerEquationsQuasi1D) = True() +function varnames(::typeof(cons2cons), ::CompressibleEulerEquationsQuasi1D) + ("a_rho", "a_rho_v1", "a_e", "a") +end +function varnames(::typeof(cons2prim), ::CompressibleEulerEquationsQuasi1D) + ("rho", "v1", "p", "a") +end + +""" + initial_condition_convergence_test(x, t, equations::CompressibleEulerEquationsQuasi1D) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function initial_condition_convergence_test(x, t, + equations::CompressibleEulerEquationsQuasi1D) + c = 2 + A = 0.1 + L = 2 + f = 1 / L + ω = 2 * pi * f + ini = c + A * sin(ω * (x[1] - t)) + + rho = ini + v1 = 1.0 + e = ini^2 / rho + p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) + a = 1.5 - 0.5 * cos(x[1] * pi) + + return prim2cons(SVector(rho, v1, p, a), equations) +end + +""" + source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquationsQuasi1D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). + +This manufactured solution source term is specifically designed for the mozzle width 'a(x) = 1.5 - 0.5 * cos(x[1] * pi)' +as defined in [`initial_condition_convergence_test`](@ref). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::CompressibleEulerEquationsQuasi1D) + # Same settings as in `initial_condition_convergence_test`. + # Derivatives calculated with ForwardDiff.jl + c = 2 + A = 0.1 + L = 2 + f = 1 / L + ω = 2 * pi * f + x1, = x + ini(x1, t) = c + A * sin(ω * (x1 - t)) + + rho(x1, t) = ini(x1, t) + v1(x1, t) = 1.0 + e(x1, t) = ini(x1, t)^2 / rho(x1, t) + p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5 * rho(x1, t) * v1(x1, t)^2) + a(x1, t) = 1.5 - 0.5 * cos(x1 * pi) + + arho(x1, t) = a(x1, t) * rho(x1, t) + arhou(x1, t) = arho(x1, t) * v1(x1, t) + aE(x1, t) = a(x1, t) * e(x1, t) + + darho_dt(x1, t) = ForwardDiff.derivative(t -> arho(x1, t), t) + darhou_dx(x1, t) = ForwardDiff.derivative(x1 -> arhou(x1, t), x1) + + arhouu(x1, t) = arhou(x1, t) * v1(x1, t) + darhou_dt(x1, t) = ForwardDiff.derivative(t -> arhou(x1, t), t) + darhouu_dx(x1, t) = ForwardDiff.derivative(x1 -> arhouu(x1, t), x1) + dp1_dx(x1, t) = ForwardDiff.derivative(x1 -> p1(x1, t), x1) + + auEp(x1, t) = a(x1, t) * v1(x1, t) * (e(x1, t) + p1(x1, t)) + daE_dt(x1, t) = ForwardDiff.derivative(t -> aE(x1, t), t) + dauEp_dx(x1, t) = ForwardDiff.derivative(x1 -> auEp(x1, t), x1) + + du1 = darho_dt(x1, t) + darhou_dx(x1, t) + du2 = darhou_dt(x1, t) + darhouu_dx(x1, t) + a(x1, t) * dp1_dx(x1, t) + du3 = daE_dt(x1, t) + dauEp_dx(x1, t) + + return SVector(du1, du2, du3, 0.0) +end + +# Calculate 1D flux for a single point +@inline function flux(u, orientation::Integer, + equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + rho, v1, p, a = cons2prim(u, equations) + e = a_e / a + + # Ignore orientation since it is always "1" in 1D + f1 = a_rho_v1 + f2 = a_rho_v1 * v1 + f3 = a * v1 * (e + p) + + return SVector(f1, f2, f3, zero(eltype(u))) +end + +""" +@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsQuasi1D) + +Non-symmetric two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the pressure [`CompressibleEulerEquationsQuasi1D`](@ref) +and the nozzle width. + +Further details are available in the paper: +- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) + High order entropy stable schemes for the quasi-one-dimensional + shallow water and compressible Euler equations + [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089) +""" +@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsQuasi1D) + #Variables + _, _, p_ll, a_ll = cons2prim(u_ll, equations) + _, _, p_rr, _ = cons2prim(u_rr, equations) + + # For flux differencing using non-conservative terms, we return the + # non-conservative flux scaled by 2. This cancels with a factor of 0.5 + # in the arithmetic average of {p}. + p_avg = p_ll + p_rr + + z = zero(eltype(u_ll)) + + return SVector(z, a_ll * p_avg, z, z) +end + +""" +@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsQuasi1D) + +Conservative (symmetric) part of the entropy conservative flux for quasi 1D compressible Euler equations split form. +This flux is a generalization of [`flux_ranocha`](@ref) for [`CompressibleEulerEquations1D`](@ref). +Further details are available in the paper: +- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) + High order entropy stable schemes for the quasi-one-dimensional + shallow water and compressible Euler equations + [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089) +""" +@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsQuasi1D) + # Unpack left and right state + rho_ll, v1_ll, p_ll, a_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, p_rr, a_rr = cons2prim(u_rr, equations) + + # Compute the necessary mean values + rho_mean = ln_mean(rho_ll, rho_rr) + # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)` + # in exact arithmetic since + # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) + # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) + inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) + v1_avg = 0.5 * (v1_ll + v1_rr) + a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) + p_avg = 0.5 * (p_ll + p_rr) + velocity_square_avg = 0.5 * (v1_ll * v1_rr) + + # Calculate fluxes + # Ignore orientation since it is always "1" in 1D + f1 = rho_mean * a_v1_avg + f2 = rho_mean * a_v1_avg * v1_avg + f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + + 0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll) + + return SVector(f1, f2, f3, zero(eltype(u_ll))) +end + +# Calculate estimates for maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerEquationsQuasi1D) + a_rho_ll, a_rho_v1_ll, a_e_ll, a_ll = u_ll + a_rho_rr, a_rho_v1_rr, a_e_rr, a_rr = u_rr + + # Calculate primitive variables and speed of sound + rho_ll = a_rho_ll / a_ll + e_ll = a_e_ll / a_ll + v1_ll = a_rho_v1_ll / a_rho_ll + v_mag_ll = abs(v1_ll) + p_ll = (equations.gamma - 1) * (e_ll - 0.5 * rho_ll * v_mag_ll^2) + c_ll = sqrt(equations.gamma * p_ll / rho_ll) + rho_rr = a_rho_rr / a_rr + e_rr = a_e_rr / a_rr + v1_rr = a_rho_v1_rr / a_rho_rr + v_mag_rr = abs(v1_rr) + p_rr = (equations.gamma - 1) * (e_rr - 0.5 * rho_rr * v_mag_rr^2) + c_rr = sqrt(equations.gamma * p_rr / rho_rr) + + λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr) +end + +@inline function max_abs_speeds(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + rho = a_rho / a + v1 = a_rho_v1 / a_rho + e = a_e / a + p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) + c = sqrt(equations.gamma * p / rho) + + return (abs(v1) + c,) +end + +# Convert conservative variables to primitive. We use the convention that the primitive +# variables for the quasi-1D equations are `(rho, v1, p)` (i.e., the same as the primitive +# variables for `CompressibleEulerEquations1D`) +@inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + q = cons2prim(SVector(a_rho, a_rho_v1, a_e) / a, + CompressibleEulerEquations1D(equations.gamma)) + + return SVector(q[1], q[2], q[3], a) +end + +# The entropy for the quasi-1D compressible Euler equations is the entropy for the +# 1D compressible Euler equations scaled by the channel width `a`. +@inline function entropy(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + q = a * entropy(SVector(a_rho, a_rho_v1, a_e) / a, + CompressibleEulerEquations1D(equations.gamma)) + + return SVector(q[1], q[2], q[3], a) +end + +# Convert conservative variables to entropy. The entropy variables for the +# quasi-1D compressible Euler equations are identical to the entropy variables +# for the standard Euler equations for an appropriate definition of `entropy`. +@inline function cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + w = cons2entropy(SVector(a_rho, a_rho_v1, a_e) / a, + CompressibleEulerEquations1D(equations.gamma)) + + # we follow the convention for other spatially-varying equations such as + # `ShallowWaterEquations1D` and return the spatially varying coefficient + # `a` as the final entropy variable. + return SVector(w[1], w[2], w[3], a) +end + +# Convert primitive to conservative variables +@inline function prim2cons(u, equations::CompressibleEulerEquationsQuasi1D) + rho, v1, p, a = u + q = prim2cons(u, CompressibleEulerEquations1D(equations.gamma)) + + return SVector(a * q[1], a * q[2], a * q[3], a) +end + +@inline function density(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, _, _, a = u + rho = a_rho / a + return rho +end + +@inline function pressure(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + return pressure(SVector(a_rho, a_rho_v1, a_e) / a, + CompressibleEulerEquations1D(equations.gamma)) +end + +@inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D) + a_rho, a_rho_v1, a_e, a = u + return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a, + CompressibleEulerEquations1D(equations.gamma)) +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 582d672b756..7a3c326984d 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -409,6 +409,7 @@ abstract type AbstractCompressibleEulerEquations{NDIMS, NVARS} <: include("compressible_euler_1d.jl") include("compressible_euler_2d.jl") include("compressible_euler_3d.jl") +include("compressible_euler_quasi_1d.jl") # CompressibleEulerMulticomponentEquations abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 6cd7998ab02..39a1f6e30ba 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -393,6 +393,79 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_quasi_1d_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_source_terms.jl"), + l2=[ + 3.876288369618363e-7, + 2.2247043122302947e-7, + 2.964004224572679e-7, + 5.2716983399807875e-8, + ], + linf=[ + 2.3925118561862746e-6, + 1.3603693522767912e-6, + 1.821888865105592e-6, + 1.1166012159335992e-7, + ]) + # 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_quasi_1d_discontinuous.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_quasi_1d_discontinuous.jl"), + l2=[ + 0.045510421156346015, + 0.036750584788912195, + 0.2468985959132176, + 0.03684494180829024, + ], + linf=[ + 0.3313374853025697, + 0.11621933362158643, + 1.827403013568638, + 0.28045939999015723, + ]) + # 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_quasi_1d_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_ec.jl"), + l2=[ + 0.08889113985713998, + 0.16199235348889673, + 0.40316524365054346, + 2.9602775074723667e-16, + ], + linf=[ + 0.28891355898284043, + 0.3752709888964313, + 0.84477102402413, + 8.881784197001252e-16, + ]) + # 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 end end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index d2e744da62f..b3ed29d38e3 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -654,6 +654,19 @@ end end end +@timed_testset "Consistency check for flux_chan_etal: CEEQ" begin + + # Set up equations and dummy conservative variables state + equations = CompressibleEulerEquationsQuasi1D(1.4) + u = SVector(1.1, 2.34, 5.5, 2.73) + + orientations = [1] + for orientation in orientations + @test flux_chan_etal(u, u, orientation, equations) ≈ + flux(u, orientation, equations) + end +end + @timed_testset "Consistency check for HLL flux (naive): LEE" begin flux_hll = FluxHLL(min_max_speed_naive) From f19144435650e626c7c57d48060ea0a03e247895 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Thu, 14 Dec 2023 13:02:17 +0100 Subject: [PATCH 06/17] Enable save solution with time intervals for SimpleIntegratorSSP (#1677) * First attempt to enable save solution with time intervals for SimpleIntegratorSSP * Importing `add_tstop!` * First working version of SimpleIntegratorSSP with SaveSolutionCallback using time intervals * Improved formatting and updated reference solution for test * Modified initialization of tstops to ensure a stop at the end of the simulation * Added missing docstring * Removed OrdinaryDiffEq from Trixi's dependencies * Empty tstops BinaryHeap during the call to `terminate!(integrator::SimpleIntegratorSSP)` * Fixed bug and added explanatory comments * Updated Project.toml * format --------- Co-authored-by: Hendrik Ranocha --- Project.toml | 2 + .../elixir_euler_shockcapturing_subcell.jl | 2 +- src/Trixi.jl | 4 +- src/time_integration/methods_SSP.jl | 75 +++++++++++++++++-- test/test_tree_2d_euler.jl | 16 ++-- 5 files changed, 82 insertions(+), 17 deletions(-) diff --git a/Project.toml b/Project.toml index 539dafc3034..267a3aa7066 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.6.5-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" @@ -52,6 +53,7 @@ TrixiMakieExt = "Makie" [compat] CodeTracking = "1.0.5" ConstructionBase = "1.3" +DataStructures = "0.18.15" DiffEqCallbacks = "2.25" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl index 282805a0e03..44e63a0872e 100644 --- a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl +++ b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl @@ -67,7 +67,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval = 100, +save_solution = SaveSolutionCallback(dt = 0.1, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) diff --git a/src/Trixi.jl b/src/Trixi.jl index e7b849e2642..e18b2f6415c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -37,7 +37,8 @@ using SciMLBase: CallbackSet, DiscreteCallback, import SciMLBase: get_du, get_tmp_cache, u_modified!, AbstractODEIntegrator, init, step!, check_error, get_proposed_dt, set_proposed_dt!, - terminate!, remake + terminate!, remake, add_tstop!, has_tstop, first_tstop + using CodeTracking: CodeTracking using ConstructionBase: ConstructionBase using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect @@ -70,6 +71,7 @@ using TriplotBase: TriplotBase using TriplotRecipes: DGTriPseudocolor @reexport using SimpleUnPack: @unpack using SimpleUnPack: @pack! +using DataStructures: BinaryHeap, FasterForward, extract_all! # finite difference SBP operators using SummationByPartsOperators: AbstractDerivativeOperator, diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index dbb9e51121b..9d1e06488b4 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -55,17 +55,25 @@ struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1 -mutable struct SimpleIntegratorSSPOptions{Callback} +mutable struct SimpleIntegratorSSPOptions{Callback, TStops} callback::Callback # callbacks; used in Trixi adaptive::Bool # whether the algorithm is adaptive; ignored dtmax::Float64 # ignored maxiters::Int # maximal number of time steps - tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored + tstops::TStops # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored end function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...) - SimpleIntegratorSSPOptions{typeof(callback)}(callback, false, Inf, maxiters, - [last(tspan)]) + tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward()) + # We add last(tspan) to make sure that the time integration stops at the end time + push!(tstops_internal, last(tspan)) + # We add 2 * last(tspan) because add_tstop!(integrator, t) is only called by DiffEqCallbacks.jl if tstops contains a time that is larger than t + # (https://github.com/SciML/DiffEqCallbacks.jl/blob/025dfe99029bd0f30a2e027582744528eb92cd24/src/iterative_and_periodic.jl#L92) + push!(tstops_internal, 2 * last(tspan)) + SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback, + false, Inf, + maxiters, + tstops_internal) end # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77 @@ -78,6 +86,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, du::uType r0::uType t::RealT + tdir::RealT dt::RealT # current time step dtcache::RealT # ignored iter::Int # current number of time steps (iteration) @@ -87,8 +96,29 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg, alg::Alg opts::SimpleIntegratorSSPOptions finalstep::Bool # added for convenience + dtchangeable::Bool + force_stepfail::Bool end +""" + add_tstop!(integrator::SimpleIntegratorSSP, t) +Add a time stop during the time integration process. +This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution. +""" +function add_tstop!(integrator::SimpleIntegratorSSP, t) + integrator.tdir * (t - integrator.t) < zero(integrator.t) && + error("Tried to add a tstop that is behind the current time. This is strictly forbidden") + # We need to remove the first entry of tstops when a new entry is added. + # Otherwise, the simulation gets stuck at the previous tstop and dt is adjusted to zero. + if length(integrator.opts.tstops) > 1 + pop!(integrator.opts.tstops) + end + push!(integrator.opts.tstops, integrator.tdir * t) +end + +has_tstop(integrator::SimpleIntegratorSSP) = !isempty(integrator.opts.tstops) +first_tstop(integrator::SimpleIntegratorSSP) = first(integrator.opts.tstops) + # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771) function Base.getproperty(integrator::SimpleIntegratorSSP, field::Symbol) if field === :stats @@ -113,11 +143,13 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP; du = similar(u) r0 = similar(u) t = first(ode.tspan) + tdir = sign(ode.tspan[end] - ode.tspan[1]) iter = 0 - integrator = SimpleIntegratorSSP(u, du, r0, t, dt, zero(dt), iter, ode.p, + integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, zero(dt), iter, ode.p, (prob = ode,), ode.f, alg, SimpleIntegratorSSPOptions(callback, ode.tspan; - kwargs...), false) + kwargs...), + false, true, false) # resize container resize!(integrator.p, nelements(integrator.p.solver, integrator.p.cache)) @@ -160,6 +192,8 @@ function solve!(integrator::SimpleIntegratorSSP) terminate!(integrator) end + modify_dt_for_tstops!(integrator) + @. integrator.r0 = integrator.u for stage in eachindex(alg.c) t_stage = integrator.t + integrator.dt * alg.c[stage] @@ -198,6 +232,10 @@ function solve!(integrator::SimpleIntegratorSSP) end end + # Empty the tstops array. + # This cannot be done in terminate!(integrator::SimpleIntegratorSSP) because DiffEqCallbacks.PeriodicCallbackAffect would return at error. + extract_all!(integrator.opts.tstops) + for stage_callback in alg.stage_callbacks finalize_callback(stage_callback, integrator.p) end @@ -226,7 +264,30 @@ end # stop the time integration function terminate!(integrator::SimpleIntegratorSSP) integrator.finalstep = true - empty!(integrator.opts.tstops) +end + +""" + modify_dt_for_tstops!(integrator::SimpleIntegratorSSP) +Modify the time-step size to match the time stops specified in integrator.opts.tstops. +To avoid adding OrdinaryDiffEq to Trixi's dependencies, this routine is a copy of +https://github.com/SciML/OrdinaryDiffEq.jl/blob/d76335281c540ee5a6d1bd8bb634713e004f62ee/src/integrators/integrator_utils.jl#L38-L54 +""" +function modify_dt_for_tstops!(integrator::SimpleIntegratorSSP) + if has_tstop(integrator) + tdir_t = integrator.tdir * integrator.t + tdir_tstop = first_tstop(integrator) + if integrator.opts.adaptive + integrator.dt = integrator.tdir * + min(abs(integrator.dt), abs(tdir_tstop - tdir_t)) # step! to the end + elseif iszero(integrator.dtcache) && integrator.dtchangeable + integrator.dt = integrator.tdir * abs(tdir_tstop - tdir_t) + elseif integrator.dtchangeable && !integrator.force_stepfail + # always try to step! with dtcache, but lower if a tstop + # however, if force_stepfail then don't set to dtcache, and no tstop worry + integrator.dt = integrator.tdir * + min(abs(integrator.dtcache), abs(tdir_tstop - tdir_t)) # step! to the end + end + end end # used for AMR diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 04a295537a3..65899cd5263 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -214,16 +214,16 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_subcell.jl"), l2=[ - 0.08508147906199143, - 0.04510299017724501, - 0.045103019801950375, - 0.6930704343869766, + 0.08508152653623638, + 0.04510301725066843, + 0.04510304668512745, + 0.6930705064715306, ], linf=[ - 0.31123546471463326, - 0.5616274869594462, - 0.5619692712224448, - 2.88670199345138, + 0.31136518019691406, + 0.5617651935473419, + 0.5621200790240503, + 2.8866869108596056, ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 40e1af41e49747f22908fc07acfc6eac92aec204 Mon Sep 17 00:00:00 2001 From: Yik Haw Teoh <66682196+teohyikhaw@users.noreply.github.com> Date: Thu, 14 Dec 2023 09:08:08 -0800 Subject: [PATCH 07/17] Fixed `cons2entropy` files and implemented `entropy2cons` for `CompressibleEulerMulticomponent1D` and `CompressibleEulerMulticomponent2D` (#1767) * initial bug fix * formatting * Added test cases for entropy2cons and cons2entropy for compressible multicomponent euler 1d and 2d * Fixed total entropy function call * minor typo * Updated formatting Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Updated formatting Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> * Fixed formatting issues * Update test/test_tree_1d_eulermulti.jl * Update test/test_tree_2d_eulermulti.jl --------- Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> --- .../compressible_euler_multicomponent_1d.jl | 78 +++- .../compressible_euler_multicomponent_2d.jl | 80 +++- test/test_tree_1d_eulermulti.jl | 241 +++++----- test/test_tree_2d_eulermulti.jl | 441 +++++++++--------- 4 files changed, 505 insertions(+), 335 deletions(-) diff --git a/src/equations/compressible_euler_multicomponent_1d.jl b/src/equations/compressible_euler_multicomponent_1d.jl index 8ddb0dcd08f..6338e04c3ed 100644 --- a/src/equations/compressible_euler_multicomponent_1d.jl +++ b/src/equations/compressible_euler_multicomponent_1d.jl @@ -461,6 +461,7 @@ end # Convert conservative variables to entropy @inline function cons2entropy(u, equations::CompressibleEulerMulticomponentEquations1D) @unpack cv, gammas, gas_constants = equations + rho_v1, rho_e = u rho = density(u, equations) @@ -480,21 +481,86 @@ end s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p T = (rho_e - 0.5 * rho * v_square) / (help1) - entrop_rho = SVector{ncomponents(equations), real(equations)}(gas_constant * - ((gamma - s) / - (gamma - 1.0) - - (0.5 * v_square * - rho_p)) + + entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 2])) - + v1^2 / (2 * T)) for i in eachcomponent(equations)) w1 = gas_constant * v1 * rho_p - w2 = gas_constant * (-1.0 * rho_p) + w2 = gas_constant * (-rho_p) entrop_other = SVector{2, real(equations)}(w1, w2) return vcat(entrop_other, entrop_rho) end +# Convert entropy variables to conservative variables +@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations1D) + @unpack gammas, gas_constants, cv, cp = equations + T = -1 / w[2] + v1 = w[1] * T + cons_rho = SVector{ncomponents(equations), real(equations)}(exp(1 / + gas_constants[i] * + (-cv[i] * + log(-w[2]) - + cp[i] + w[i + 2] - + 0.5 * w[1]^2 / + w[2])) + for i in eachcomponent(equations)) + + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + gamma = help1 / help2 + u2 = p / (gamma - 1) + 0.5 * rho * v1^2 + cons_other = SVector{2, real(equations)}(u1, u2) + return vcat(cons_other, cons_rho) +end + +@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations1D) + @unpack cv, gammas, gas_constants = equations + rho_v1, rho_e = u + rho = density(u, equations) + T = temperature(u, equations) + + total_entropy = zero(u[1]) + for i in eachcomponent(equations) + total_entropy -= u[i + 2] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 2])) + end + + return total_entropy +end + +@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations1D) + @unpack cv, gammas, gas_constants = equations + + rho_v1, rho_e = u + + rho = density(u, equations) + help1 = zero(rho) + + for i in eachcomponent(equations) + help1 += u[i + 2] * cv[i] + end + + v1 = rho_v1 / rho + v_square = v1^2 + T = (rho_e - 0.5 * rho * v_square) / help1 + + return T +end + """ totalgamma(u, equations::CompressibleEulerMulticomponentEquations1D) diff --git a/src/equations/compressible_euler_multicomponent_2d.jl b/src/equations/compressible_euler_multicomponent_2d.jl index 940d88b1aa5..60fce222f21 100644 --- a/src/equations/compressible_euler_multicomponent_2d.jl +++ b/src/equations/compressible_euler_multicomponent_2d.jl @@ -665,22 +665,57 @@ end s = log(p) - gamma * log(rho) - log(gas_constant) rho_p = rho / p T = (rho_e - 0.5 * rho * v_square) / (help1) - entrop_rho = SVector{ncomponents(equations), real(equations)}(gas_constant * - ((gamma - s) / - (gamma - 1.0) - - (0.5 * v_square * - rho_p)) + + entrop_rho = SVector{ncomponents(equations), real(equations)}((cv[i] * + (1 - log(T)) + + gas_constants[i] * + (1 + log(u[i + 3])) - + v_square / (2 * T)) for i in eachcomponent(equations)) w1 = gas_constant * v1 * rho_p w2 = gas_constant * v2 * rho_p - w3 = gas_constant * rho_p * (-1) + w3 = gas_constant * (-rho_p) entrop_other = SVector{3, real(equations)}(w1, w2, w3) return vcat(entrop_other, entrop_rho) end +# Convert entropy variables to conservative variables +@inline function entropy2cons(w, equations::CompressibleEulerMulticomponentEquations2D) + @unpack gammas, gas_constants, cp, cv = equations + T = -1 / w[3] + v1 = w[1] * T + v2 = w[2] * T + v_squared = v1^2 + v2^2 + cons_rho = SVector{ncomponents(equations), real(equations)}(exp((w[i + 3] - + cv[i] * + (1 - log(T)) + + v_squared / + (2 * T)) / + gas_constants[i] - + 1) + for i in eachcomponent(equations)) + + rho = zero(cons_rho[1]) + help1 = zero(cons_rho[1]) + help2 = zero(cons_rho[1]) + p = zero(cons_rho[1]) + for i in eachcomponent(equations) + rho += cons_rho[i] + help1 += cons_rho[i] * cv[i] * gammas[i] + help2 += cons_rho[i] * cv[i] + p += cons_rho[i] * gas_constants[i] * T + end + u1 = rho * v1 + u2 = rho * v2 + gamma = help1 / help2 + u3 = p / (gamma - 1) + 0.5 * rho * v_squared + cons_other = SVector{3, real(equations)}(u1, u2, u3) + return vcat(cons_other, cons_rho) +end + # Convert primitive to conservative variables @inline function prim2cons(prim, equations::CompressibleEulerMulticomponentEquations2D) @unpack cv, gammas = equations @@ -700,6 +735,39 @@ end return vcat(cons_other, cons_rho) end +@inline function total_entropy(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations + rho = density(u, equations) + T = temperature(u, equations) + + total_entropy = zero(u[1]) + for i in eachcomponent(equations) + total_entropy -= u[i + 3] * (cv[i] * log(T) - gas_constants[i] * log(u[i + 3])) + end + + return total_entropy +end + +@inline function temperature(u, equations::CompressibleEulerMulticomponentEquations2D) + @unpack cv, gammas, gas_constants = equations + + rho_v1, rho_v2, rho_e = u + + rho = density(u, equations) + help1 = zero(rho) + + for i in eachcomponent(equations) + help1 += u[i + 3] * cv[i] + end + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + v_square = v1^2 + v2^2 + T = (rho_e - 0.5 * rho * v_square) / help1 + + return T +end + """ totalgamma(u, equations::CompressibleEulerMulticomponentEquations2D) diff --git a/test/test_tree_1d_eulermulti.jl b/test/test_tree_1d_eulermulti.jl index bd86de928e3..b6c79ce03d1 100644 --- a/test/test_tree_1d_eulermulti.jl +++ b/test/test_tree_1d_eulermulti.jl @@ -2,136 +2,153 @@ module TestExamples1DEulerMulti using Test using Trixi +using ForwardDiff include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Compressible Euler Multicomponent" begin -#! format: noindent + @trixi_testset "Testing entropy2cons and cons2entropy" begin + using ForwardDiff + gammas = (1.3272378792562836, 1.5269959187969864, 1.8362285750521512, + 1.0409061360276926, 1.4652015053812224, 1.3626493264184423) + gas_constants = (1.817636851910076, 6.760820475922636, 5.588953939749113, + 6.31574782981543, 3.362932038038397, 3.212779569399733) + equations = CompressibleEulerMulticomponentEquations1D(gammas = SVector{length(gammas)}(gammas...), + gas_constants = SVector{length(gas_constants)}(gas_constants...)) + u = [-1.4632513788889214, 0.9908786980927811, 0.2909066990257628, + 0.6256623915420473, 0.4905882754313441, 0.14481800501749112, + 1.0333532872771651, 0.6805599818745411] + w = cons2entropy(u, equations) + # test that the entropy variables match the gradients of the total entropy + @test w ≈ ForwardDiff.gradient(u -> Trixi.total_entropy(u, equations), u) + # test that `entropy2cons` is the inverse of `cons2entropy` + @test entropy2cons(w, equations) ≈ u + end -@trixi_testset "elixir_eulermulti_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), - l2=[0.15330089521538684, 0.4417674632047301, - 0.016888510510282385, 0.03377702102056477, - 0.06755404204112954], - linf=[0.29130548795961864, 0.8847009003152357, - 0.034686525099975274, 0.06937305019995055, - 0.1387461003999011]) - # 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 + @trixi_testset "elixir_eulermulti_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), + l2=[0.15330089521538684, 0.4417674632047301, + 0.016888510510282385, 0.03377702102056477, + 0.06755404204112954], + linf=[0.29130548795961864, 0.8847009003152357, + 0.034686525099975274, 0.06937305019995055, + 0.1387461003999011]) + # 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 -end -@trixi_testset "elixir_eulermulti_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), - l2=[ - 0.1522380497572071, - 0.43830846465313206, - 0.03907262116499431, - 0.07814524232998862, - ], - linf=[ - 0.24939193075537294, - 0.7139395740052739, - 0.06324208768391237, - 0.12648417536782475, - ]) - # 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 + @trixi_testset "elixir_eulermulti_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), + l2=[ + 0.1522380497572071, + 0.43830846465313206, + 0.03907262116499431, + 0.07814524232998862, + ], + linf=[ + 0.24939193075537294, + 0.7139395740052739, + 0.06324208768391237, + 0.12648417536782475, + ]) + # 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 -end -@trixi_testset "elixir_eulermulti_convergence_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), - l2=[ - 8.575236038539227e-5, - 0.00016387804318585358, - 1.9412699303977585e-5, - 3.882539860795517e-5, - ], - linf=[ - 0.00030593277277124464, - 0.0006244803933350696, - 7.253121435135679e-5, - 0.00014506242870271358, - ]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), + l2=[ + 8.575236038539227e-5, + 0.00016387804318585358, + 1.9412699303977585e-5, + 3.882539860795517e-5, + ], + linf=[ + 0.00030593277277124464, + 0.0006244803933350696, + 7.253121435135679e-5, + 0.00014506242870271358, + ]) + # 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 -end -@trixi_testset "elixir_eulermulti_convergence_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[1.8983933794407234e-5, 6.207744299844731e-5, - 1.5466205761868047e-6, 3.0932411523736094e-6, - 6.186482304747219e-6, 1.2372964609494437e-5], - linf=[0.00012014372605895218, 0.0003313207215800418, - 6.50836791016296e-6, 1.301673582032592e-5, - 2.603347164065184e-5, 5.206694328130368e-5]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[1.8983933794407234e-5, 6.207744299844731e-5, + 1.5466205761868047e-6, 3.0932411523736094e-6, + 6.186482304747219e-6, 1.2372964609494437e-5], + linf=[0.00012014372605895218, 0.0003313207215800418, + 6.50836791016296e-6, 1.301673582032592e-5, + 2.603347164065184e-5, 5.206694328130368e-5]) + # 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 -end -@trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[1.888450477353845e-5, 5.4910600482795386e-5, - 9.426737161533622e-7, 1.8853474323067245e-6, - 3.770694864613449e-6, 7.541389729226898e-6], - linf=[0.00011622351152063004, 0.0003079221967086099, - 3.2177423254231563e-6, 6.435484650846313e-6, - 1.2870969301692625e-5, 2.574193860338525e-5], - volume_flux=flux_chandrashekar) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[1.888450477353845e-5, 5.4910600482795386e-5, + 9.426737161533622e-7, 1.8853474323067245e-6, + 3.770694864613449e-6, 7.541389729226898e-6], + linf=[0.00011622351152063004, 0.0003079221967086099, + 3.2177423254231563e-6, 6.435484650846313e-6, + 1.2870969301692625e-5, 2.574193860338525e-5], + volume_flux=flux_chandrashekar) + # 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 -end -@trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_eulermulti_two_interacting_blast_waves.jl"), - l2=[1.288867611915533, 82.71335258388848, 0.00350680272313187, - 0.013698784353152794, - 0.019179518517518084], - linf=[29.6413044707026, 1322.5844802186496, 0.09191919374782143, - 0.31092970966717925, - 0.4417989757182038], - tspan=(0.0, 0.0001)) - # 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 + @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_two_interacting_blast_waves.jl"), + l2=[1.288867611915533, 82.71335258388848, 0.00350680272313187, + 0.013698784353152794, + 0.019179518517518084], + linf=[29.6413044707026, 1322.5844802186496, 0.09191919374782143, + 0.31092970966717925, + 0.4417989757182038], + tspan=(0.0, 0.0001)) + # 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 end -end end # module diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index 30d52b37b96..7c4a4e722e3 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -8,234 +8,253 @@ include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") @testset "Compressible Euler Multicomponent" begin -#! format: noindent + @trixi_testset "Testing entropy2cons and cons2entropy" begin + using ForwardDiff + gammas = (1.1546412974182538, 1.1171560258914812, 1.097107661471476, + 1.0587601652669245, 1.6209889683979308, 1.6732209755396386, + 1.2954303574165822) + gas_constants = (5.969461071171914, 3.6660802003290183, 6.639008614675539, + 8.116604827140456, 6.190706056680031, 1.6795013743693712, + 2.197737590916966) + equations = CompressibleEulerMulticomponentEquations2D(gammas = SVector{length(gammas)}(gammas...), + gas_constants = SVector{length(gas_constants)}(gas_constants...)) + u = [-1.7433292819144075, 0.8844413258376495, 0.6050737175812364, + 0.8261998359817043, 1.0801186290896465, 0.505654488367698, + 0.6364415555805734, 0.851669392285058, 0.31219606420306223, + 1.0930477805612038] + w = cons2entropy(u, equations) + # test that the entropy variables match the gradients of the total entropy + @test w ≈ ForwardDiff.gradient(u -> Trixi.total_entropy(u, equations), u) + # test that `entropy2cons` is the inverse of `cons2entropy` + @test entropy2cons(w, equations) ≈ u + end -# NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the -# simulations are set up with dimensional states. For example, the reference pressure in SI -# units is 101325 Pa, i.e., pressure has values of O(10^5) + # NOTE: Some of the L2/Linf errors are comparably large. This is due to the fact that some of the + # simulations are set up with dimensional states. For example, the reference pressure in SI + # units is 101325 Pa, i.e., pressure has values of O(10^5) -@trixi_testset "elixir_eulermulti_shock_bubble.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"), - l2=[ - 73.78467629094177, - 0.9174752929795251, - 57942.83587826468, - 0.1828847253029943, - 0.011127037850925347, - ], - linf=[ - 196.81051991521073, - 7.8456811648529605, - 158891.88930113698, - 0.811379581519794, - 0.08011973559187913, - ], - tspan=(0.0, 0.001)) - # 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 + @trixi_testset "elixir_eulermulti_shock_bubble.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_shock_bubble.jl"), + l2=[ + 73.78467629094177, + 0.9174752929795251, + 57942.83587826468, + 0.1828847253029943, + 0.011127037850925347, + ], + linf=[ + 196.81051991521073, + 7.8456811648529605, + 158891.88930113698, + 0.811379581519794, + 0.08011973559187913, + ], + tspan=(0.0, 0.001)) + # 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 -end -@trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin - rm("out/deviations.txt", force = true) - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"), - l2=[ - 81.52845664909304, - 2.5455678559421346, - 63229.190712645846, - 0.19929478404550321, - 0.011068604228443425, - ], - linf=[ - 249.21708417382013, - 40.33299887640794, - 174205.0118831558, - 0.6881458768113586, - 0.11274401158173972, - ], - initial_refinement_level=3, - tspan=(0.0, 0.001), - output_directory="out") - lines = readlines("out/deviations.txt") - @test lines[1] == "# iter, simu_time, rho1_min, rho2_min" - @test startswith(lines[end], "1") - # 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)) < 15000 + @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl" begin + rm("out/deviations.txt", force = true) + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl"), + l2=[ + 81.52845664909304, + 2.5455678559421346, + 63229.190712645846, + 0.19929478404550321, + 0.011068604228443425, + ], + linf=[ + 249.21708417382013, + 40.33299887640794, + 174205.0118831558, + 0.6881458768113586, + 0.11274401158173972, + ], + initial_refinement_level=3, + tspan=(0.0, 0.001), + output_directory="out") + lines = readlines("out/deviations.txt") + @test lines[1] == "# iter, simu_time, rho1_min, rho2_min" + @test startswith(lines[end], "1") + # 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)) < 15000 + end end -end -@trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"), - l2=[ - 73.10832638093902, - 1.4599215762968585, - 57176.014861335476, - 0.17812843581838675, - 0.010123079422717837, - ], - linf=[ - 214.50568817511956, - 25.40392579616452, - 152862.41011222568, - 0.564195553101797, - 0.0956331651771212, - ], - initial_refinement_level=3, - tspan=(0.0, 0.001)) - # 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)) < 15000 + @trixi_testset "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl"), + l2=[ + 73.10832638093902, + 1.4599215762968585, + 57176.014861335476, + 0.17812843581838675, + 0.010123079422717837, + ], + linf=[ + 214.50568817511956, + 25.40392579616452, + 152862.41011222568, + 0.564195553101797, + 0.0956331651771212, + ], + initial_refinement_level=3, + tspan=(0.0, 0.001)) + # 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)) < 15000 + end end -end -@trixi_testset "elixir_eulermulti_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), - l2=[ - 0.050182236154087095, - 0.050189894464434635, - 0.2258715597305131, - 0.06175171559771687, - ], - linf=[ - 0.3108124923284472, - 0.3107380389947733, - 1.054035804988521, - 0.29347582879608936, - ]) - # 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 + @trixi_testset "elixir_eulermulti_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"), + l2=[ + 0.050182236154087095, + 0.050189894464434635, + 0.2258715597305131, + 0.06175171559771687, + ], + linf=[ + 0.3108124923284472, + 0.3107380389947733, + 1.054035804988521, + 0.29347582879608936, + ]) + # 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 -end -@trixi_testset "elixir_eulermulti_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), - l2=[ - 0.0496546258404055, - 0.04965550099933263, - 0.22425206549856372, - 0.004087155041747821, - 0.008174310083495642, - 0.016348620166991283, - 0.032697240333982566, - ], - linf=[ - 0.2488251110766228, - 0.24832493304479406, - 0.9310354690058298, - 0.017452870465607374, - 0.03490574093121475, - 0.0698114818624295, - 0.139622963724859, - ]) - # 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 + @trixi_testset "elixir_eulermulti_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"), + l2=[ + 0.0496546258404055, + 0.04965550099933263, + 0.22425206549856372, + 0.004087155041747821, + 0.008174310083495642, + 0.016348620166991283, + 0.032697240333982566, + ], + linf=[ + 0.2488251110766228, + 0.24832493304479406, + 0.9310354690058298, + 0.017452870465607374, + 0.03490574093121475, + 0.0698114818624295, + 0.139622963724859, + ]) + # 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 -end -@trixi_testset "elixir_eulermulti_convergence_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), - l2=[ - 0.00012290225488326508, - 0.00012290225488321876, - 0.00018867397906337653, - 4.8542321753649044e-5, - 9.708464350729809e-5, - ], - linf=[ - 0.0006722819239133315, - 0.0006722819239128874, - 0.0012662292789555885, - 0.0002843844182700561, - 0.0005687688365401122, - ]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), + l2=[ + 0.00012290225488326508, + 0.00012290225488321876, + 0.00018867397906337653, + 4.8542321753649044e-5, + 9.708464350729809e-5, + ], + linf=[ + 0.0006722819239133315, + 0.0006722819239128874, + 0.0012662292789555885, + 0.0002843844182700561, + 0.0005687688365401122, + ]) + # 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 -end -@trixi_testset "elixir_eulermulti_convergence_es.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[ - 2.2661773867001696e-6, - 2.266177386666318e-6, - 6.593514692980009e-6, - 8.836308667348217e-7, - 1.7672617334696433e-6, - ], - linf=[ - 1.4713170997993075e-5, - 1.4713170997104896e-5, - 5.115618808515521e-5, - 5.3639516094383666e-6, - 1.0727903218876733e-5, - ]) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[ + 2.2661773867001696e-6, + 2.266177386666318e-6, + 6.593514692980009e-6, + 8.836308667348217e-7, + 1.7672617334696433e-6, + ], + linf=[ + 1.4713170997993075e-5, + 1.4713170997104896e-5, + 5.115618808515521e-5, + 5.3639516094383666e-6, + 1.0727903218876733e-5, + ]) + # 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 -end -@trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), - l2=[ - 1.8621737639352465e-6, - 1.862173764098385e-6, - 5.942585713809631e-6, - 6.216263279534722e-7, - 1.2432526559069443e-6, - ], - linf=[ - 1.6235495582606063e-5, - 1.6235495576388814e-5, - 5.854523678827661e-5, - 5.790274858807898e-6, - 1.1580549717615796e-5, - ], - volume_flux=flux_chandrashekar) - # 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 + @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"), + l2=[ + 1.8621737639352465e-6, + 1.862173764098385e-6, + 5.942585713809631e-6, + 6.216263279534722e-7, + 1.2432526559069443e-6, + ], + linf=[ + 1.6235495582606063e-5, + 1.6235495576388814e-5, + 5.854523678827661e-5, + 5.790274858807898e-6, + 1.1580549717615796e-5, + ], + volume_flux=flux_chandrashekar) + # 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 end -end end # module From 318e145e72110fe2acf0fd1f685bffe21e4df965 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Fri, 15 Dec 2023 17:09:16 +0100 Subject: [PATCH 08/17] CompatHelper: bump compat for T8code to 0.5, (keep existing compat) (#1775) Co-authored-by: CompatHelper Julia Co-authored-by: Hendrik Ranocha --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 267a3aa7066..a834ddeea73 100644 --- a/Project.toml +++ b/Project.toml @@ -86,7 +86,7 @@ StaticArrays = "1" StrideArrays = "0.1.18" StructArrays = "0.6" SummationByPartsOperators = "0.5.41" -T8code = "0.4.3" +T8code = "0.4.3, 0.5" TimerOutputs = "0.5.7" Triangulate = "2.0" TriplotBase = "0.1" From 6bef107cd322f7a31b50770219327de4345d93bc Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sat, 16 Dec 2023 15:06:09 +0100 Subject: [PATCH 09/17] Add get_proposed_dt to custom integrators (#1776) --- src/time_integration/methods_2N.jl | 5 +++++ src/time_integration/methods_3Sstar.jl | 5 +++++ 2 files changed, 10 insertions(+) diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index d2f22679c4f..f3b09b01e97 100644 --- a/src/time_integration/methods_2N.jl +++ b/src/time_integration/methods_2N.jl @@ -204,6 +204,11 @@ function set_proposed_dt!(integrator::SimpleIntegrator2N, dt) integrator.dt = dt end +# Required e.g. for `glm_speed_callback` +function get_proposed_dt(integrator::SimpleIntegrator2N) + return integrator.dt +end + # stop the time integration function terminate!(integrator::SimpleIntegrator2N) integrator.finalstep = true diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl index b0ce5930514..7b70466606c 100644 --- a/src/time_integration/methods_3Sstar.jl +++ b/src/time_integration/methods_3Sstar.jl @@ -282,6 +282,11 @@ function set_proposed_dt!(integrator::SimpleIntegrator3Sstar, dt) integrator.dt = dt end +# Required e.g. for `glm_speed_callback` +function get_proposed_dt(integrator::SimpleIntegrator3Sstar) + return integrator.dt +end + # stop the time integration function terminate!(integrator::SimpleIntegrator3Sstar) integrator.finalstep = true From b5d0a50214abb9b33dc4a4859a4a12c6e3149cb1 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Mon, 18 Dec 2023 08:36:23 -0600 Subject: [PATCH 10/17] Update NEWS.md (#1780) --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 265979c3508..abd9fd27882 100644 --- a/NEWS.md +++ b/NEWS.md @@ -38,7 +38,7 @@ for human readability. - Capability to set truly discontinuous initial conditions in 1D. - Wetting and drying feature and examples for 1D and 2D shallow water equations - Implementation of the polytropic Euler equations in 2D -- Implementation of the quasi-1D shallow water equations +- Implementation of the quasi-1D shallow water and compressible Euler equations - Subcell (positivity and local min/max) limiting support for conservative variables in 2D for `TreeMesh` - AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh` From 5a5424c1fcefe602b4854477bebf96dc77c2e9b7 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 20 Dec 2023 07:33:30 +0100 Subject: [PATCH 11/17] set version to v0.6.5 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a834ddeea73..1ef9ee13516 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.6.5-pre" +version = "0.6.5" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 622740a9db3c3e1b9bb1f5b0520b18b616d250e6 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 20 Dec 2023 07:33:42 +0100 Subject: [PATCH 12/17] set development version to v0.6.6-pre --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1ef9ee13516..6bedec18c78 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.6.5" +version = "0.6.6-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 24a365dbf2eb029b9dbd2be3e8404538bc54e89d Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 21 Dec 2023 09:13:10 +0100 Subject: [PATCH 13/17] hotfix: restrict DiffEqBase.jl to let CI pass (#1788) * hotfix: restrict DiffEqBase.jl to let CI pass * restrict DiffEqBase.jl in main Project.toml * Update Project.toml --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 6bedec18c78..faf9f82d335 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.6.6-pre" CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" @@ -54,6 +55,7 @@ TrixiMakieExt = "Makie" CodeTracking = "1.0.5" ConstructionBase = "1.3" DataStructures = "0.18.15" +DiffEqBase = "6 - 6.143" DiffEqCallbacks = "2.25" EllipsisNotation = "1.0" FillArrays = "0.13.2, 1" From 96c5bef1a68729ae61be5a6644d7b7b4756aac8e Mon Sep 17 00:00:00 2001 From: Ahmad Peyvan <115842305+apey236@users.noreply.github.com> Date: Sat, 23 Dec 2023 10:18:42 -0500 Subject: [PATCH 14/17] Parabolic Mortar for AMR `P4estMesh{3}` (#1765) * Clean branch * Un-Comment * un-comment * test coarsen * remove redundancy * Remove support for passive terms * expand resize * comments * format * Avoid code duplication * Update src/callbacks_step/amr_dg1d.jl Co-authored-by: Michael Schlottke-Lakemper * comment * comment & format * Try to increase coverage * Slightly more expressive names * Apply suggestions from code review * add specifier for 1d * Structs for resizing parabolic helpers * check if mortars are present * reuse `reinitialize_containers!` * resize calls for parabolic helpers * update analysis callbacks * Velocities for compr euler * Init container * correct copy-paste error * resize each dim * add dispatch * Add AMR for shear layer * USe only amr shear layer * first steps towards p4est parabolic amr * Add tests * remove plots * Format * remove redundant line * platform independent tests * No need for different flux_viscous comps after adding container_viscous to p4est * Laplace 3d * Longer times to allow converage to hit coarsen! * Increase testing of Laplace 3D * Add tests for velocities * remove comment * add elixir for amr testing * adding commented out mortar routines in 2D * Adding Mortar to 2d dg parabolic term * remove testing snippet * fix comments * add more arguments for dispatch * add some temporary todo notes * some updates for AP and KS * specialize mortar_fluxes_to_elements * BUGFIX: apply_jacobian_parabolic! was incorrect for P4estMesh * fixed rhs_parabolic! for mortars * more changes to elixir * indexing bug * comments * Adding the example for nonperiodic BCs with amr * hopefully this fixes AMR boundaries for parabolic terms * add elixir * Example with non periodic bopundary conditions * remove cruft * 3D parabolic amr * TGV elixir * Creating test for AMR 3D parabolic * Formatting * test formatting * Update src/Trixi.jl * Update src/equations/compressible_euler_1d.jl * Update src/equations/compressible_euler_2d.jl * Update src/equations/compressible_euler_3d.jl * Update src/solvers/dgsem_tree/container_viscous_2d.jl * Update src/solvers/dgsem_tree/container_viscous_2d.jl * Update src/solvers/dgsem_tree/container_viscous_2d.jl * Update src/solvers/dgsem_tree/container_viscous_3d.jl * Update src/solvers/dgsem_tree/container_viscous_3d.jl * Update src/solvers/dgsem_tree/container_viscous_3d.jl * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update src/solvers/dgsem_tree/dg_2d_parabolic.jl * Update test/test_parabolic_3d.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_tree/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_2d.jl * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_2d_parabolic.jl * Update test/test_parabolic_3d.jl * Update test/test_parabolic_3d.jl * Update test/test_parabolic_3d.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl * Update dg_3d_parabolic.jl * Update test_parabolic_3d.jl * Update elixir_navierstokes_3d_blast_wave_amr.jl * Update elixir_navierstokes_taylor_green_vortex_amr.jl * Update dg_3d_parabolic.jl * Update test_parabolic_3d.jl * Update test_parabolic_3d.jl * Update examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Update elixir_navierstokes_3d_blast_wave_amr.jl * Update elixir_navierstokes_taylor_green_vortex_amr.jl * Update dg_3d_parabolic.jl * Update test_parabolic_3d.jl * Delete examples/p4est_3d_dgsem/elixir_navierstokes_3d_blast_wave_amr.jl * Create elixir_navierstokes_blast_wave_amr.jl * Update test_parabolic_3d.jl * Update NEWS.md * Update NEWS.md * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl * Update src/solvers/dgsem_p4est/dg_3d_parabolic.jl --------- Co-authored-by: Daniel_Doehring Co-authored-by: Daniel Doehring Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com> Co-authored-by: Jesse Chan --- NEWS.md | 5 + .../elixir_navierstokes_blast_wave_amr.jl | 113 ++++++ ...ir_navierstokes_taylor_green_vortex_amr.jl | 106 ++++++ src/callbacks_step/amr_dg2d.jl | 4 +- src/equations/laplace_diffusion_3d.jl | 1 + src/solvers/dgsem_p4est/dg_3d_parabolic.jl | 347 +++++++++++++++++- test/test_parabolic_3d.jl | 67 ++++ 7 files changed, 637 insertions(+), 6 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl create mode 100644 examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl diff --git a/NEWS.md b/NEWS.md index abd9fd27882..e7bdd14eab2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,11 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes in the v0.6 lifecycle + +#### Added +- AMR for hyperbolic-parabolic equations on 3D `P4estMesh` + ## Changes when updating to v0.6 from v0.5.x #### Added diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl new file mode 100644 index 00000000000..5df89fbcdf2 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl @@ -0,0 +1,113 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 6.25e-4 # equivalent to Re = 1600 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), + Prandtl = prandtl_number()) + +function initial_condition_3d_blast_wave(x, t, equations::CompressibleEulerEquations3D) + rho_c = 1.0 + p_c = 1.0 + u_c = 0.0 + + rho_o = 0.125 + p_o = 0.1 + u_o = 0.0 + + rc = 0.5 + r = sqrt(x[1]^2 + x[2]^2 + x[3]^2) + if r < rc + rho = rho_c + v1 = u_c + v2 = u_c + v3 = u_c + p = p_c + else + rho = rho_o + v1 = u_o + v2 = u_o + v3 = u_o + p = p_o + end + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_3d_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +coordinates_min = (-1.0, -1.0, -1.0) .* pi +coordinates_max = (1.0, 1.0, 1.0) .* pi + +trees_per_dimension = (4, 4, 4) + +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = (true, true, true), initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.8) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +save_solution = SaveSolutionCallback(interval = analysis_interval, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.05, + max_level = 3, max_threshold = 0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval = 10, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + save_solution) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl new file mode 100644 index 00000000000..c15227a1c29 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl @@ -0,0 +1,106 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 6.25e-4 # equivalent to Re = 1600 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(), + Prandtl = prandtl_number()) + +""" + initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D) + +The classical Taylor-Green vortex. +""" +function initial_condition_taylor_green_vortex(x, t, + equations::CompressibleEulerEquations3D) + A = 1.0 # magnitude of speed + Ms = 0.1 # maximum Mach number + + rho = 1.0 + v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3]) + v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3]) + v3 = 0.0 + p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms + p = p + + 1.0 / 16.0 * A^2 * rho * + (cos(2 * x[1]) * cos(2 * x[3]) + 2 * cos(2 * x[2]) + 2 * cos(2 * x[1]) + + cos(2 * x[2]) * cos(2 * x[3])) + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_taylor_green_vortex + +@inline function vel_mag(u, equations::CompressibleEulerEquations3D) + rho, rho_v1, rho_v2, rho_v3, _ = u + return sqrt(rho_v1^2 + rho_v2^2 + rho_v3^2) / rho +end + +volume_flux = flux_ranocha +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0, -1.0) .* pi +coordinates_max = (1.0, 1.0, 1.0) .* pi + +trees_per_dimension = (2, 2, 2) + +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = (true, true, true), initial_refinement_level = 0) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 50 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (energy_kinetic, + energy_internal, + enstrophy)) +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +amr_indicator = IndicatorLöhner(semi, variable = vel_mag) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.1, + max_level = 3, max_threshold = 0.2) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + save_solution) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) +summary_callback() # print the timer summary diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 969f9c564f3..98e531295b7 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -137,7 +137,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM end function refine!(u_ode::AbstractVector, adaptor, - mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}}, + mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_refine) # Call `refine!` for the hyperbolic part, which does the heavy lifting of @@ -299,7 +299,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, end function coarsen!(u_ode::AbstractVector, adaptor, - mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}}, + mesh::Union{TreeMesh{2}, P4estMesh{2}, TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_remove) # Call `coarsen!` for the hyperbolic part, which does the heavy lifting of diff --git a/src/equations/laplace_diffusion_3d.jl b/src/equations/laplace_diffusion_3d.jl index 457e742430b..3988ce7144b 100644 --- a/src/equations/laplace_diffusion_3d.jl +++ b/src/equations/laplace_diffusion_3d.jl @@ -18,6 +18,7 @@ function varnames(variable_mapping, equations_parabolic::LaplaceDiffusion3D) varnames(variable_mapping, equations_parabolic.equations_hyperbolic) end +# no orientation specified since the flux is vector-valued function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion3D) dudx, dudy, dudz = gradients if orientation == 1 diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index b06cdd42127..0bb97c7af02 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -102,8 +102,20 @@ function rhs_parabolic!(du, u, t, mesh::P4estMesh{3}, dg.surface_integral, dg) end - # TODO: parabolic; extend to mortars - @assert nmortars(dg, cache) == 0 + # Prolong solution to mortars (specialized for AbstractEquationsParabolic) + # !!! NOTE: we reuse the hyperbolic cache here since it contains "mortars" and "u_threaded" + # !!! Is this OK? + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes (specialized for AbstractEquationsParabolic) + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values, + mesh, equations_parabolic, dg.mortar, + dg.surface_integral, dg, cache) + end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin @@ -230,8 +242,23 @@ function calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic, dg.surface_integral, dg) end - # TODO: parabolic; mortars - @assert nmortars(dg, cache) == 0 + # Prolong solution to mortars. These should reuse the hyperbolic version of `prolong2mortars` + # !!! NOTE: we reuse the hyperbolic cache here, since it contains both `mortars` and `u_threaded`. + # !!! should we have a separate mortars/u_threaded in cache_parabolic? + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes. These should reuse the hyperbolic version of `calc_mortar_flux`, + # along with a specialization on `calc_mortar_flux!` and `mortar_fluxes_to_elements!` for + # AbstractEquationsParabolic. + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, + mesh, False(), # False() = no nonconservative terms + equations_parabolic, + dg.mortar, dg.surface_integral, dg, cache) + end # Calculate surface integrals @trixi_timeit timer() "surface integral" begin @@ -324,6 +351,93 @@ function calc_gradient!(gradients, u_transformed, t, return nothing end +# This version is called during `calc_gradients!` and must be specialized because the flux +# in the gradient is {u} which doesn't depend on normals. Thus, you don't need to scale by +# 2 and flip the sign when storing the mortar fluxes back into surface_flux_values +@inline function mortar_fluxes_to_elements!(surface_flux_values, + mesh::P4estMesh{3}, + equations::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + dg::DGSEM, cache, mortar, fstar, u_buffer, + fstar_tmp) + @unpack neighbor_ids, node_indices = cache.mortars + index_range = eachnode(dg) + # Copy solution small to small + small_indices = node_indices[1, mortar] + small_direction = indices2direction(small_indices) + + for position in 1:4 # Loop over small elements + element = neighbor_ids[position, mortar] + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, j, small_direction, element] = fstar[v, i, j, + position] + end + end + end + + # Project small fluxes to large element. + multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_lower, mortar_l2.reverse_lower, + view(fstar, .., 1), + fstar_tmp) + add_multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_upper, mortar_l2.reverse_lower, + view(fstar, .., 2), + fstar_tmp) + add_multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_lower, mortar_l2.reverse_upper, + view(fstar, .., 3), + fstar_tmp) + add_multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_upper, mortar_l2.reverse_upper, + view(fstar, .., 4), + fstar_tmp) + + # The flux is calculated in the outward direction of the small elements, + # so the sign must be switched to get the flux in outward direction + # of the large element. + # The contravariant vectors of the large element (and therefore the normal + # vectors of the large element as well) are twice as large as the + # contravariant vectors of the small elements. Therefore, the flux needs + # to be scaled by a factor of 2 to obtain the flux of the large element. + # u_buffer .*= 0.5 + + # Copy interpolated flux values from buffer to large element face in the + # correct orientation. + # Note that the index of the small sides will always run forward but + # the index of the large side might need to run backwards for flipped sides. + large_element = neighbor_ids[5, mortar] + large_indices = node_indices[2, mortar] + large_direction = indices2direction(large_indices) + large_surface_indices = surface_indices(large_indices) + + i_large_start, i_large_step_i, i_large_step_j = index_to_start_step_3d(large_surface_indices[1], + index_range) + j_large_start, j_large_step_i, j_large_step_j = index_to_start_step_3d(large_surface_indices[2], + index_range) + + # Note that the indices of the small sides will always run forward but + # the large indices might need to run backwards for flipped sides. + i_large = i_large_start + j_large = j_large_start + for j in eachnode(dg) + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i_large, j_large, large_direction, large_element] = u_buffer[v, + i, + j] + end + i_large += i_large_step_i + j_large += j_large_step_i + end + i_large += i_large_step_j + j_large += j_large_step_j + end + + return nothing +end + # This version is used for parabolic gradient computations @inline function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{3}, nonconservative_terms::False, @@ -603,6 +717,231 @@ function calc_interface_flux!(surface_flux_values, return nothing end +function prolong2mortars_divergence!(cache, flux_viscous, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DGSEM) + @unpack neighbor_ids, node_indices = cache.mortars + @unpack fstar_tmp_threaded = cache + @unpack contravariant_vectors = cache.elements + index_range = eachnode(dg) + + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous + + @threaded for mortar in eachmortar(dg, cache) + # Copy solution data from the small elements using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + small_indices = node_indices[1, mortar] + direction_index = indices2direction(small_indices) + + i_small_start, i_small_step_i, i_small_step_j = index_to_start_step_3d(small_indices[1], + index_range) + j_small_start, j_small_step_i, j_small_step_j = index_to_start_step_3d(small_indices[2], + index_range) + k_small_start, k_small_step_i, k_small_step_j = index_to_start_step_3d(small_indices[3], + index_range) + + for position in 1:4 # Loop over small elements + i_small = i_small_start + j_small = j_small_start + k_small = k_small_start + element = neighbor_ids[position, mortar] + for j in eachnode(dg) + for i in eachnode(dg) + normal_direction = get_normal_direction(direction_index, + contravariant_vectors, + i_small, j_small, k_small, + element) + + for v in eachvariable(equations) + flux_viscous = SVector(flux_viscous_x[v, i_small, j_small, k_small, + element], + flux_viscous_y[v, i_small, j_small, k_small, + element], + flux_viscous_z[v, i_small, j_small, k_small, + element]) + + cache.mortars.u[1, v, position, i, j, mortar] = dot(flux_viscous, + normal_direction) + end + i_small += i_small_step_i + j_small += j_small_step_i + k_small += k_small_step_i + end + i_small += i_small_step_j + j_small += j_small_step_j + k_small += k_small_step_j + end + end + + # Buffer to copy solution values of the large element in the correct orientation + # before interpolating + u_buffer = cache.u_threaded[Threads.threadid()] + + # temporary buffer for projections + fstar_tmp = fstar_tmp_threaded[Threads.threadid()] + + # Copy solution of large element face to buffer in the + # correct orientation + large_indices = node_indices[2, mortar] + + i_large_start, i_large_step_i, i_large_step_j = index_to_start_step_3d(large_indices[1], + index_range) + j_large_start, j_large_step_i, j_large_step_j = index_to_start_step_3d(large_indices[2], + index_range) + k_large_start, k_large_step_i, k_large_step_j = index_to_start_step_3d(large_indices[3], + index_range) + + i_large = i_large_start + j_large = j_large_start + k_large = k_large_start + element = neighbor_ids[5, mortar] # Large element + for j in eachnode(dg) + for i in eachnode(dg) + normal_direction = get_normal_direction(direction_index, + contravariant_vectors, + i_large, j_large, k_large, element) + + for v in eachvariable(equations) + flux_viscous = SVector(flux_viscous_x[v, i_large, j_large, k_large, + element], + flux_viscous_y[v, i_large, j_large, k_large, + element], + flux_viscous_z[v, i_large, j_large, k_large, + element]) + + # We prolong the viscous flux dotted with respect the outward normal + # on the small element. We scale by -1/2 here because the normal + # direction on the large element is negative 2x that of the small + # element (these normal directions are "scaled" by the surface Jacobian) + u_buffer[v, i, j] = -0.5 * dot(flux_viscous, normal_direction) + end + i_large += i_large_step_i + j_large += j_large_step_i + k_large += k_large_step_i + end + i_large += i_large_step_j + j_large += j_large_step_j + k_large += k_large_step_j + end + + # Interpolate large element face data from buffer to small face locations + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, :, mortar), + mortar_l2.forward_lower, + mortar_l2.forward_lower, + u_buffer, + fstar_tmp) + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, :, mortar), + mortar_l2.forward_upper, + mortar_l2.forward_lower, + u_buffer, + fstar_tmp) + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 3, :, :, mortar), + mortar_l2.forward_lower, + mortar_l2.forward_upper, + u_buffer, + fstar_tmp) + multiply_dimensionwise!(view(cache.mortars.u, 2, :, 4, :, :, mortar), + mortar_l2.forward_upper, + mortar_l2.forward_upper, + u_buffer, + fstar_tmp) + end + + return nothing +end + +# We specialize `calc_mortar_flux!` for the divergence part of +# the parabolic terms. +function calc_mortar_flux_divergence!(surface_flux_values, + mesh::Union{P4estMesh{3}, T8codeMesh{3}}, + equations::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.mortars + @unpack contravariant_vectors = cache.elements + @unpack fstar_threaded, fstar_tmp_threaded = cache + index_range = eachnode(dg) + + @threaded for mortar in eachmortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar = fstar_threaded[Threads.threadid()] + fstar_tmp = fstar_tmp_threaded[Threads.threadid()] + + # Get index information on the small elements + small_indices = node_indices[1, mortar] + small_direction = indices2direction(small_indices) + + i_small_start, i_small_step_i, i_small_step_j = index_to_start_step_3d(small_indices[1], + index_range) + j_small_start, j_small_step_i, j_small_step_j = index_to_start_step_3d(small_indices[2], + index_range) + k_small_start, k_small_step_i, k_small_step_j = index_to_start_step_3d(small_indices[3], + index_range) + + for position in 1:4 # Loop over small elements + i_small = i_small_start + j_small = j_small_start + k_small = k_small_start + element = neighbor_ids[position, mortar] + for j in eachnode(dg) + for i in eachnode(dg) + for v in eachvariable(equations) + viscous_flux_normal_ll = cache.mortars.u[1, v, position, i, j, + mortar] + viscous_flux_normal_rr = cache.mortars.u[2, v, position, i, j, + mortar] + + # TODO: parabolic; only BR1 at the moment + fstar[v, i, j, position] = 0.5 * (viscous_flux_normal_ll + + viscous_flux_normal_rr) + end + + i_small += i_small_step_i + j_small += j_small_step_i + k_small += k_small_step_i + end + i_small += i_small_step_j + j_small += j_small_step_j + k_small += k_small_step_j + end + end + + # Buffer to interpolate flux values of the large element to before + # copying in the correct orientation + u_buffer = cache.u_threaded[Threads.threadid()] + + # this reuses the hyperbolic version of `mortar_fluxes_to_elements!` + mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar, u_buffer, fstar_tmp) + end + + return nothing +end + +# NOTE: Use analogy to "calc_mortar_flux!" for hyperbolic eqs with no nonconservative terms. +# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for +# hyperbolic terms with conserved terms only, i.e., no nonconservative terms. +@inline function calc_mortar_flux!(fstar, + mesh::P4estMesh{3}, + nonconservative_terms::False, + equations::AbstractEquationsParabolic, + surface_integral, dg::DG, cache, + mortar_index, position_index, normal_direction, + i_node_index, j_node_index) + @unpack u = cache.mortars + @unpack surface_flux = surface_integral + + u_ll, u_rr = get_surface_node_vars(u, equations, dg, position_index, i_node_index, + j_node_index, mortar_index) + + # TODO: parabolic; only BR1 at the moment + flux_ = 0.5 * (u_ll + u_rr) + # Copy flux to buffer + set_node_vars!(fstar, flux_, equations, dg, i_node_index, j_node_index, position_index) +end + # TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::P4estMesh{3}, diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index d6c720cf0d9..6fbfb8259d4 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -429,6 +429,14 @@ end "elixir_advection_diffusion_amr.jl"), l2=[0.000355780485397024], linf=[0.0010810770271614256]) + # 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 "TreeMesh3D: elixir_advection_diffusion_nonperiodic.jl" begin @@ -436,6 +444,65 @@ end "elixir_advection_diffusion_nonperiodic.jl"), l2=[0.0009808996243280868], linf=[0.01732621559135459]) + # 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 "P4estMesh3D: elixir_navierstokes_taylor_green_vortex_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_3d_dgsem", + "elixir_navierstokes_taylor_green_vortex_amr.jl"), + initial_refinement_level=0, tspan=(0.0, 0.5), + l2=[ + 0.0016588740573444188, + 0.03437058632045721, + 0.03437058632045671, + 0.041038898400430075, + 0.30978593009044153, + ], + linf=[ + 0.004173569912012121, + 0.09168674832979556, + 0.09168674832975021, + 0.12129218723807476, + 0.8433893297612087, + ]) + # 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 "P4estMesh3D: elixir_navierstokes_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_3d_dgsem", + "elixir_navierstokes_blast_wave_amr.jl"), + tspan=(0.0, 0.01), + l2=[ + 0.009472104410520866, 0.0017883742549557149, + 0.0017883742549557147, 0.0017883742549557196, + 0.024388540048562748, + ], + linf=[ + 0.6782397526873181, 0.17663702154066238, + 0.17663702154066266, 0.17663702154066238, 1.7327849844825238, + ]) + # 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 end From dc8e9272d4e2f4b8df7fc4e72e12a786666fde47 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 24 Dec 2023 13:12:13 +0100 Subject: [PATCH 15/17] reset the timer also on non-root MPI processes (#1787) Co-authored-by: Michael Schlottke-Lakemper --- src/callbacks_step/summary.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/summary.jl b/src/callbacks_step/summary.jl index 566f2c03418..21c7fc780a5 100644 --- a/src/callbacks_step/summary.jl +++ b/src/callbacks_step/summary.jl @@ -152,7 +152,13 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator; Polyester.reset_threads!() end - mpi_isroot() || return nothing + # The summary callback should only print information on the root process. + # However, all other MPI processes should also reset the timer so that + # it can be used to diagnose performance. + if !mpi_isroot() + reset_timer!(timer()) + return nothing + end print_startup_message() From a633f9c30b89bca886278a716cafd98c6046becd Mon Sep 17 00:00:00 2001 From: Benedict <135045760+bgeihe@users.noreply.github.com> Date: Sat, 30 Dec 2023 15:49:50 +0100 Subject: [PATCH 16/17] Add HLLC flux for non-cartesian meshes to CompressibleEulerEquations{2,3}D (#1790) * add HLLC flux for non-cartesian meshes * add tests for HLLC flux * Add 2D test with HLLC * Update test_p4est_3d.jl * Update test_p4est_2d.jl * Update test_p4est_3d.jl * Update src/equations/compressible_euler_3d.jl Co-authored-by: Hendrik Ranocha * Update src/equations/compressible_euler_2d.jl Co-authored-by: Hendrik Ranocha * Update compressible_euler_2d.jl * Update compressible_euler_3d.jl * Update test_p4est_2d.jl * Update test_p4est_3d.jl * Update compressible_euler_2d.jl * Update compressible_euler_2d.jl --------- Co-authored-by: Daniel Doehring Co-authored-by: Hendrik Ranocha --- NEWS.md | 1 + src/equations/compressible_euler_2d.jl | 102 +++++++++++++++++++-- src/equations/compressible_euler_3d.jl | 119 +++++++++++++++++++++++-- test/test_p4est_2d.jl | 26 ++++++ test/test_p4est_3d.jl | 28 ++++++ test/test_unit.jl | 46 +++++++++- 6 files changed, 305 insertions(+), 17 deletions(-) diff --git a/NEWS.md b/NEWS.md index e7bdd14eab2..cf695912ed7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,6 +8,7 @@ for human readability. #### Added - AMR for hyperbolic-parabolic equations on 3D `P4estMesh` +- `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D` ## Changes when updating to v0.6 from v0.5.x diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index a992f99eaf4..b0fd5c53f45 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1150,7 +1150,7 @@ end end """ - flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) + flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations2D) Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro [Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf) @@ -1185,18 +1185,18 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, if orientation == 1 # x-direction vel_L = v1_ll vel_R = v1_rr - ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 elseif orientation == 2 # y-direction vel_L = v2_ll vel_R = v2_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 end vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2) + v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr + v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr + vel_roe_mag = (v1_roe^2 + v2_roe^2) / sum_sqrt_rho^2 H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) Ssl = min(vel_L - c_ll, vel_roe - c_roe) Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L @@ -1252,6 +1252,98 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, return SVector(f1, f2, f3, f4) end +function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations2D) + # Calculate primitive variables and speed of sound + rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + norm_sq = norm_ * norm_ + inv_norm_sq = inv(norm_sq) + + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + # Obtain left and right fluxes + f_ll = flux(u_ll, normal_direction, equations) + f_rr = flux(u_rr, normal_direction, equations) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho + vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + vel_roe_mag = v1_roe^2 + v2_roe^2 + + e_ll = u_ll[4] / rho_ll + e_rr = u_rr[4] / rho_rr + + H_ll = (u_ll[4] + p_ll) / rho_ll + H_rr = (u_rr[4] + p_rr) / rho_rr + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_ + + Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe) + Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe) + sMu_L = Ssl - v_dot_n_ll + sMu_R = Ssr - v_dot_n_rr + + if Ssl >= 0.0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + elseif Ssr <= 0.0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + else + SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R + + (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R) + if Ssl <= 0.0 <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + enerStar = e_ll + + (SStar - v_dot_n_ll) * + (SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L)) + UStar1 = densStar + UStar2 = densStar * + (v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * enerStar + f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1]) + f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2]) + f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3]) + f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4]) + else + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + + (SStar - v_dot_n_rr) * + (SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R)) + UStar1 = densStar + UStar2 = densStar * + (v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * enerStar + f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1]) + f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2]) + f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3]) + f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4]) + end + end + return SVector(f1, f2, f3, f4) +end + """ min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations2D) diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index fc56f58025b..82c4a7efa32 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1192,7 +1192,7 @@ end end """ - flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D) + flux_hllc(u_ll, u_rr, orientation_or_normal_direction, equations::CompressibleEulerEquations3D) Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro [Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf) @@ -1231,25 +1231,22 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, if orientation == 1 # x-direction vel_L = v1_ll vel_R = v1_rr - ekin_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 + - (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2 elseif orientation == 2 # y-direction vel_L = v2_ll vel_R = v2_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 + - (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr)^2 else # z-direction vel_L = v3_ll vel_R = v3_rr - ekin_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr)^2 + - (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr)^2 end vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5 * (vel_roe^2 + ekin_roe / sum_sqrt_rho^2) + v1_roe = sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr + v2_roe = sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr + v3_roe = sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr + vel_roe_mag = (v1_roe^2 + v2_roe^2 + v3_roe^2) / sum_sqrt_rho^2 H_ll = (rho_e_ll + p_ll) / rho_ll H_rr = (rho_e_rr + p_rr) / rho_rr H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((equations.gamma - 1) * (H_roe - ekin_roe)) + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) Ssl = min(vel_L - c_ll, vel_roe - c_roe) Ssr = max(vel_R + c_rr, vel_roe + c_roe) sMu_L = Ssl - vel_L @@ -1321,6 +1318,110 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, return SVector(f1, f2, f3, f4, f5) end +function flux_hllc(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleEulerEquations3D) + # Calculate primitive variables and speed of sound + rho_ll, v1_ll, v2_ll, v3_ll, p_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, v3_rr, p_rr = cons2prim(u_rr, equations) + + v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + + v3_ll * normal_direction[3] + v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + v3_rr * normal_direction[3] + + norm_ = norm(normal_direction) + norm_sq = norm_ * norm_ + inv_norm_sq = inv(norm_sq) + + c_ll = sqrt(equations.gamma * p_ll / rho_ll) * norm_ + c_rr = sqrt(equations.gamma * p_rr / rho_rr) * norm_ + + # Obtain left and right fluxes + f_ll = flux(u_ll, normal_direction, equations) + f_rr = flux(u_rr, normal_direction, equations) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr + + v1_roe = (sqrt_rho_ll * v1_ll + sqrt_rho_rr * v1_rr) / sum_sqrt_rho + v2_roe = (sqrt_rho_ll * v2_ll + sqrt_rho_rr * v2_rr) / sum_sqrt_rho + v3_roe = (sqrt_rho_ll * v3_ll + sqrt_rho_rr * v3_rr) / sum_sqrt_rho + vel_roe = v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + + v3_roe * normal_direction[3] + vel_roe_mag = v1_roe^2 + v2_roe^2 + v3_roe^2 + + e_ll = u_ll[5] / rho_ll + e_rr = u_rr[5] / rho_rr + + H_ll = (u_ll[5] + p_ll) / rho_ll + H_rr = (u_rr[5] + p_rr) / rho_rr + + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho + c_roe = sqrt((equations.gamma - 1) * (H_roe - 0.5 * vel_roe_mag)) * norm_ + + Ssl = min(v_dot_n_ll - c_ll, vel_roe - c_roe) + Ssr = max(v_dot_n_rr + c_rr, vel_roe + c_roe) + sMu_L = Ssl - v_dot_n_ll + sMu_R = Ssr - v_dot_n_rr + + if Ssl >= 0.0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + f4 = f_ll[4] + f5 = f_ll[5] + elseif Ssr <= 0.0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + f4 = f_rr[4] + f5 = f_rr[5] + else + SStar = (rho_ll * v_dot_n_ll * sMu_L - rho_rr * v_dot_n_rr * sMu_R + + (p_rr - p_ll) * norm_sq) / (rho_ll * sMu_L - rho_rr * sMu_R) + if Ssl <= 0.0 <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + enerStar = e_ll + + (SStar - v_dot_n_ll) * + (SStar * inv_norm_sq + p_ll / (rho_ll * sMu_L)) + UStar1 = densStar + UStar2 = densStar * + (v1_ll + (SStar - v_dot_n_ll) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_ll + (SStar - v_dot_n_ll) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * + (v3_ll + (SStar - v_dot_n_ll) * normal_direction[3] * inv_norm_sq) + UStar5 = densStar * enerStar + f1 = f_ll[1] + Ssl * (UStar1 - u_ll[1]) + f2 = f_ll[2] + Ssl * (UStar2 - u_ll[2]) + f3 = f_ll[3] + Ssl * (UStar3 - u_ll[3]) + f4 = f_ll[4] + Ssl * (UStar4 - u_ll[4]) + f5 = f_ll[5] + Ssl * (UStar5 - u_ll[5]) + else + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + + (SStar - v_dot_n_rr) * + (SStar * inv_norm_sq + p_rr / (rho_rr * sMu_R)) + UStar1 = densStar + UStar2 = densStar * + (v1_rr + (SStar - v_dot_n_rr) * normal_direction[1] * inv_norm_sq) + UStar3 = densStar * + (v2_rr + (SStar - v_dot_n_rr) * normal_direction[2] * inv_norm_sq) + UStar4 = densStar * + (v3_rr + (SStar - v_dot_n_rr) * normal_direction[3] * inv_norm_sq) + UStar5 = densStar * enerStar + f1 = f_rr[1] + Ssr * (UStar1 - u_rr[1]) + f2 = f_rr[2] + Ssr * (UStar2 - u_rr[2]) + f3 = f_rr[3] + Ssr * (UStar3 - u_rr[3]) + f4 = f_rr[4] + Ssr * (UStar4 - u_rr[4]) + f5 = f_rr[5] + Ssr * (UStar5 - u_rr[5]) + end + end + return SVector(f1, f2, f3, f4, f5) +end + """ min_max_speed_einfeldt(u_ll, u_rr, orientation, equations::CompressibleEulerEquations3D) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index db34aecc168..cebc2917d52 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -203,6 +203,32 @@ end end end +@trixi_testset "elixir_euler_sedov.jl with HLLC Flux" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), + l2=[ + 0.4229948321239887, + 0.2559038337457483, + 0.2559038337457484, + 1.2990046683564136, + ], + linf=[ + 1.4989357969730492, + 1.325456585141623, + 1.3254565851416251, + 6.331283015053501, + ], + surface_flux=flux_hllc, + tspan=(0.0, 0.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_sedov.jl (HLLE)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index f2467f30204..4a2d2112c99 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -234,6 +234,34 @@ end end end +@trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"), + l2=[ + 8.444868392439035e-16, + 4.889826056731442e-15, + 2.2921260987087585e-15, + 4.268460455702414e-15, + 1.1356712092620279e-14, + ], + linf=[ + 7.749356711883593e-14, + 4.513472928735496e-13, + 2.9790059308254513e-13, + 1.057154364048074e-12, + 1.6271428648906294e-12, + ], + tspan=(0.0, 0.1), + surface_flux=flux_hllc) + # 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_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"), l2=[ diff --git a/test/test_unit.jl b/test/test_unit.jl index b3ed29d38e3..817b4cd550d 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -1076,6 +1076,46 @@ end end end +@timed_testset "Consistency check for HLLC flux: CEE" begin + # Set up equations and dummy conservative variables state + equations = CompressibleEulerEquations2D(1.4) + u = SVector(1.1, -0.5, 2.34, 5.5) + + orientations = [1, 2] + for orientation in orientations + @test flux_hllc(u, u, orientation, equations) ≈ flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for normal_direction in normal_directions + @test flux_hllc(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end + + equations = CompressibleEulerEquations3D(1.4) + u = SVector(1.1, -0.5, 2.34, 2.4, 5.5) + + orientations = [1, 2, 3] + for orientation in orientations + @test flux_hllc(u, u, orientation, equations) ≈ flux(u, orientation, equations) + end + + normal_directions = [SVector(1.0, 0.0, 0.0), + SVector(0.0, 1.0, 0.0), + SVector(0.0, 0.0, 1.0), + SVector(0.5, -0.5, 0.2), + SVector(-1.2, 0.3, 1.4)] + + for normal_direction in normal_directions + @test flux_hllc(u, u, normal_direction, equations) ≈ + flux(u, normal_direction, equations) + end +end + @timed_testset "Consistency check for Godunov flux" begin # Set up equations and dummy conservative variables state # Burgers' Equation @@ -1280,7 +1320,7 @@ end u_values = [SVector(1.0, 0.5, -0.7, 1.0), SVector(1.5, -0.2, 0.1, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] + flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc] for f_std in fluxes f_rot = FluxRotated(f_std) @@ -1303,8 +1343,8 @@ end u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0), SVector(1.5, -0.2, 0.1, 0.2, 5.0)] fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber, - FluxLMARS(340), - flux_hll, FluxHLL(min_max_speed_davis), flux_hlle] + FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc, + ] for f_std in fluxes f_rot = FluxRotated(f_std) From ad384c6530b2539175167350b6e56a9365d3521c Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 3 Jan 2024 07:24:11 +0100 Subject: [PATCH 17/17] Bump crate-ci/typos from 1.16.23 to 1.16.26 (#1793) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.16.23 to 1.16.26. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.16.23...v1.16.26) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 366cb1183a0..a780e975155 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.16.23 + uses: crate-ci/typos@v1.16.26