diff --git a/Project.toml b/Project.toml index 65074c5a1ca..267a3aa7066 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,12 @@ 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.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_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/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/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/Trixi.jl b/src/Trixi.jl index b8110cf5bdd..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, @@ -139,6 +141,7 @@ export AcousticPerturbationEquations2D, CompressibleEulerEquations3D, CompressibleEulerMulticomponentEquations1D, CompressibleEulerMulticomponentEquations2D, + CompressibleEulerEquationsQuasi1D, IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, 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/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 6bc789b116d..0c0b198800b 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -415,6 +415,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/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/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_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_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_euler.jl b/test/test_tree_2d_euler.jl index 9315ae4d281..06b2da681cd 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) 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 diff --git a/test/test_unit.jl b/test/test_unit.jl index 080f387e8fc..8b91f0875c4 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -655,6 +655,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) 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