diff --git a/NEWS.md b/NEWS.md index 0c78484a782..bc213fcea55 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,7 @@ for human readability. - Wetting and drying feature and examples for 1D and 2D shallow water equations - Implementation of the polytropic Euler equations in 2D - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh` +- AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh` #### Changed diff --git a/Project.toml b/Project.toml index 48bb6b30182..f19f7fdecc3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.5.47-pre" +version = "0.5.48-pre" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl similarity index 53% rename from examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl rename to examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index dd26fd8097b..06e8f06d3ca 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -7,12 +7,20 @@ using Trixi # TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 1.0/3.0 * 10^(-3) # equivalent to Re = 3000 +mu() = 1.0/3.0 * 10^(-5) # equivalent to Re = 30,000 equations = CompressibleEulerEquations2D(1.4) equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number()) +""" +A compressible version of the double shear layer initial condition. Adapted from +Brown and Minion (1995). + +- David L. Brown and Michael L. Minion (1995) + Performance of Under-resolved Two-Dimensional Incompressible Flow Simulations. + [DOI: 10.1006/jcph.1995.1205](https://doi.org/10.1006/jcph.1995.1205) +""" function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) # Shear layer parameters k = 80 @@ -22,8 +30,8 @@ function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquatio Ms = 0.1 # maximum Mach number rho = 1.0 - v1 = x[2] <= 0.5 ? u0 * tanh(k*(x[2]*0.5 - 0.25)) : u0 * tanh(k*(0.75 -x[2]*0.5)) - v2 = u0 * delta * sin(2*pi*(x[1]*0.5 + 0.25)) + v1 = x[2] <= 0.5 ? u0 * tanh(k*(x[2] - 0.25)) : u0 * tanh(k*(0.75 -x[2])) + v2 = u0 * delta * sin(2*pi*(x[1]+ 0.25)) p = (u0 / Ms)^2 * rho / equations.gamma # scaling to get Ms return prim2cons(SVector(rho, v1, v2, p), equations) @@ -47,27 +55,43 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 2.0) +tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 50 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true, - extra_analysis_integrals=(energy_kinetic, - energy_internal, - enstrophy)) +analysis_interval = 2000 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval,) +# This uses velocity-based AMR +@inline function v1(u, equations::CompressibleEulerEquations2D) + rho, rho_v1, _, _ = u + return rho_v1 / rho +end +amr_indicator = IndicatorLöhner(semi, variable=v1) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 3, + med_level = 5, med_threshold=0.2, + max_level = 7, max_threshold=0.5) +amr_callback = AMRCallback(semi, amr_controller, + interval=50, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + +stepsize_callback = StepsizeCallback(cfl=1.3) + callbacks = CallbackSet(summary_callback, analysis_callback, - alive_callback) + alive_callback, + amr_callback, + stepsize_callback) ############################################################################### # run the simulation -time_int_tol = 1e-8 -sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, - ode_default_options()..., callback=callbacks) +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 \ No newline at end of file diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl new file mode 100644 index 00000000000..0fab685b642 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl @@ -0,0 +1,93 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +diffusivity() = 5.0e-4 +equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations) + +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0, -1.0) +coordinates_max = ( 1.0, 1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=80_000) + +# Define initial condition +function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation3D) + # Store translated coordinate for easy use of exact solution + x_trans = x - equation.advection_velocity * t + + nu = diffusivity() + c = 1.0 + A = 0.5 + L = 2 + f = 1/L + omega = 2 * pi * f + scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t) + return SVector(scalar) +end +initial_condition = initial_condition_diffusive_convergence_test + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_periodic +boundary_conditions_parabolic = boundary_condition_periodic + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_integrals=(entropy,)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=100, + save_initial_solution=true, + save_final_solution=true, + solution_variables=cons2prim) + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=3, + med_level=4, med_threshold=1.2, + max_level=5, max_threshold=1.45) +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + +stepsize_callback = StepsizeCallback(cfl=1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + amr_callback, + 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 \ No newline at end of file diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl new file mode 100644 index 00000000000..5dc6e6338a7 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl @@ -0,0 +1,91 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +diffusivity() = 5.0e-2 +advection_velocity = (1.0, 0.0, 0.0) +equations = LinearScalarAdvectionEquation3D(advection_velocity) +equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) + +coordinates_min = (-1.0, -0.5, -0.25) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = ( 0.0, 0.5, 0.25) # maximum coordinates (max(x), max(y), max(z)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + periodicity=false, + n_cells_max=30_000) # set maximum capacity of tree data structure + +# Example setup taken from +# - Truman Ellis, Jesse Chan, and Leszek Demkowicz (2016). +# Robust DPG methods for transient convection-diffusion. +# In: Building bridges: connections and challenges in modern approaches +# to numerical partial differential equations. +# [DOI](https://doi.org/10.1007/978-3-319-41640-3_6). +function initial_condition_eriksson_johnson(x, t, equations) + l = 4 + epsilon = diffusivity() # TODO: this requires epsilon < .6 due to sqrt + lambda_1 = (-1 + sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + lambda_2 = (-1 - sqrt(1 - 4 * epsilon * l)) / (-2 * epsilon) + r1 = (1 + sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + s1 = (1 - sqrt(1 + 4 * pi^2 * epsilon^2)) / (2 * epsilon) + u = exp(-l * t) * (exp(lambda_1 * x[1]) - exp(lambda_2 * x[1])) + + cos(pi * x[2]) * (exp(s1 * x[1]) - exp(r1 * x[1])) / (exp(-s1) - exp(-r1)) + return SVector{1}(u) +end +initial_condition = initial_condition_eriksson_johnson + +boundary_conditions = (; x_neg = BoundaryConditionDirichlet(initial_condition), + y_neg = BoundaryConditionDirichlet(initial_condition), + z_neg = boundary_condition_do_nothing, + y_pos = BoundaryConditionDirichlet(initial_condition), + x_pos = boundary_condition_do_nothing, + z_pos = boundary_condition_do_nothing) + +boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions=(boundary_conditions, + boundary_conditions_parabolic)) + + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan); + +# 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_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +time_int_tol = 1.0e-11 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, + ode_default_options()..., callback=callbacks) + +# Print the timer summary +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index 457d9dc336d..5cb3cf0a9fe 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -153,7 +153,7 @@ export AcousticPerturbationEquations2D, LinearizedEulerEquations2D, PolytropicEulerEquations2D -export LaplaceDiffusion1D, LaplaceDiffusion2D, +export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index e721ccc61cb..b4cd6a00271 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -83,30 +83,13 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # actually transferring the solution to the refined cells refine!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_refine) - # The remaining function only handles the necessary adaptation of the data structures - # for the parabolic part of the semidiscretization - - # Get new list of leaf cells - leaf_cell_ids = local_leaf_cells(mesh.tree) - - @unpack elements, viscous_container = cache_parabolic - resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - # Resize parabolic helper variables + @unpack viscous_container = cache_parabolic resize!(viscous_container, equations, dg, cache) - - # re-initialize interfaces container - @unpack interfaces = cache_parabolic - resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids)) - init_interfaces!(interfaces, elements, mesh) - - # re-initialize boundaries container - @unpack boundaries = cache_parabolic - resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids)) - init_boundaries!(boundaries, elements, mesh) + reinitialize_containers!(mesh, equations, dg, cache_parabolic) # Sanity check + @unpack interfaces = cache_parabolic if isperiodic(mesh.tree) @assert ninterfaces(interfaces)==1 * nelements(dg, cache_parabolic) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") end @@ -246,27 +229,13 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # actually transferring the solution to the coarsened cells coarsen!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_remove) - # Get new list of leaf cells - leaf_cell_ids = local_leaf_cells(mesh.tree) - - @unpack elements, viscous_container = cache_parabolic - resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - # Resize parabolic helper variables + @unpack viscous_container = cache_parabolic resize!(viscous_container, equations, dg, cache) - - # re-initialize interfaces container - @unpack interfaces = cache_parabolic - resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids)) - init_interfaces!(interfaces, elements, mesh) - - # re-initialize boundaries container - @unpack boundaries = cache_parabolic - resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids)) - init_boundaries!(boundaries, elements, mesh) + reinitialize_containers!(mesh, equations, dg, cache_parabolic) # Sanity check + @unpack interfaces = cache_parabolic if isperiodic(mesh.tree) @assert ninterfaces(interfaces)==1 * nelements(dg, cache_parabolic) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") end diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 1d37dfce034..6395a9f348f 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -136,6 +136,29 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM return nothing end +# AMR for hyperbolic-parabolic equations currently only supported on TreeMeshes +function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMesh{3}}, + equations, dg::DGSEM, cache, cache_parabolic, + elements_to_refine) + # Call `refine!` for the hyperbolic part, which does the heavy lifting of + # actually transferring the solution to the refined cells + refine!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_refine) + + # Resize parabolic helper variables + @unpack viscous_container = cache_parabolic + resize!(viscous_container, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache_parabolic) + + # Sanity check + if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && + !mpi_isparallel() + @assert ninterfaces(cache_parabolic.interfaces)==ndims(mesh) * + nelements(dg, cache_parabolic) ("For $(ndims(mesh))D and periodic domains and conforming elements, the number of interfaces must be $(ndims(mesh)) times the number of elements") + end + + return nothing +end + # TODO: Taal compare performance of different implementations # Refine solution data u for an element, using L2 projection (interpolation) function refine_element!(u::AbstractArray{<:Any, 4}, element_id, @@ -275,6 +298,29 @@ function coarsen!(u_ode::AbstractVector, adaptor, return nothing end +# AMR for hyperbolic-parabolic equations currently only supported on TreeMeshes +function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMesh{3}}, + equations, dg::DGSEM, cache, cache_parabolic, + elements_to_remove) + # Call `coarsen!` for the hyperbolic part, which does the heavy lifting of + # actually transferring the solution to the coarsened cells + coarsen!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_remove) + + # Resize parabolic helper variables + @unpack viscous_container = cache_parabolic + resize!(viscous_container, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache_parabolic) + + # Sanity check + if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && + !mpi_isparallel() + @assert ninterfaces(cache_parabolic.interfaces)==ndims(mesh) * + nelements(dg, cache_parabolic) ("For $(ndims(mesh))D and periodic domains and conforming elements, the number of interfaces must be $(ndims(mesh)) times the number of elements") + end + + return nothing +end + # TODO: Taal compare performance of different implementations # Coarsen solution data u for four elements, using L2 projection function coarsen_elements!(u::AbstractArray{<:Any, 4}, element_id, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index aecabf0e4b7..a9e0cf87b0a 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -218,7 +218,7 @@ function integrate(func::Func, u, equations, equations_parabolic, dg::DGSEM, cache, cache_parabolic; normalize = true) where {Func} - gradients_x, gradients_y = cache_parabolic.gradients + gradients_x, gradients_y = cache_parabolic.viscous_container.gradients integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 3d9b38fd2a5..81d0795a159 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -232,7 +232,7 @@ function integrate(func::Func, u, equations, equations_parabolic, dg::DGSEM, cache, cache_parabolic; normalize = true) where {Func} - gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + gradients_x, gradients_y, gradients_z = cache_parabolic.viscous_container.gradients integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, k, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, k, element) diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index 66214025044..47a76174cb1 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -7,6 +7,7 @@ abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end include("laplace_diffusion_1d.jl") include("laplace_diffusion_2d.jl") +include("laplace_diffusion_3d.jl") # Compressible Navier-Stokes equations abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS} <: diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl index 3443e9c097b..b848633fbcb 100644 --- a/src/equations/laplace_diffusion_2d.jl +++ b/src/equations/laplace_diffusion_2d.jl @@ -18,7 +18,6 @@ function varnames(variable_mapping, equations_parabolic::LaplaceDiffusion2D) varnames(variable_mapping, equations_parabolic.equations_hyperbolic) end -# no orientation specified since the flux is vector-valued function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion2D) dudx, dudy = gradients if orientation == 1 diff --git a/src/equations/laplace_diffusion_3d.jl b/src/equations/laplace_diffusion_3d.jl new file mode 100644 index 00000000000..457e742430b --- /dev/null +++ b/src/equations/laplace_diffusion_3d.jl @@ -0,0 +1,71 @@ +@doc raw""" + LaplaceDiffusion3D(diffusivity, equations) + +`LaplaceDiffusion3D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))`` +with diffusivity ``\kappa`` applied to each solution component defined by `equations`. +""" +struct LaplaceDiffusion3D{E, N, T} <: AbstractLaplaceDiffusion{3, N} + diffusivity::T + equations_hyperbolic::E +end + +function LaplaceDiffusion3D(diffusivity, equations_hyperbolic) + LaplaceDiffusion3D{typeof(equations_hyperbolic), nvariables(equations_hyperbolic), + typeof(diffusivity)}(diffusivity, equations_hyperbolic) +end + +function varnames(variable_mapping, equations_parabolic::LaplaceDiffusion3D) + varnames(variable_mapping, equations_parabolic.equations_hyperbolic) +end + +function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion3D) + dudx, dudy, dudz = gradients + if orientation == 1 + return SVector(equations_parabolic.diffusivity * dudx) + elseif orientation == 2 + return SVector(equations_parabolic.diffusivity * dudy) + else # if orientation == 3 + return SVector(equations_parabolic.diffusivity * dudz) + end +end + +# TODO: parabolic; should this remain in the equations file, be moved to solvers, or live in the elixir? +# The penalization depends on the solver, but also depends explicitly on physical parameters, +# and would probably need to be specialized for every different equation. +function penalty(u_outer, u_inner, inv_h, equations_parabolic::LaplaceDiffusion3D, + dg::ViscousFormulationLocalDG) + return dg.penalty_parameter * (u_outer - u_inner) * equations_parabolic.diffusivity +end + +# Dirichlet-type boundary condition for use with a parabolic solver in weak form +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations_parabolic::LaplaceDiffusion3D) + return boundary_condition.boundary_value_function(x, t, equations_parabolic) +end + +@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations_parabolic::LaplaceDiffusion3D) + return flux_inner +end + +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, + normal::AbstractVector, + x, t, + operator_type::Divergence, + equations_parabolic::LaplaceDiffusion3D) + return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic) +end + +@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner, + normal::AbstractVector, + x, t, + operator_type::Gradient, + equations_parabolic::LaplaceDiffusion3D) + return flux_inner +end diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index a04523d2fb4..cf07645b949 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -10,14 +10,11 @@ function create_cache_parabolic(mesh::P4estMesh{2}, equations_hyperbolic::Abstra interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements) boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements) - n_vars = nvariables(equations_hyperbolic) - n_elements = nelements(elements) - n_nodes = nnodes(dg.basis) # nodes in one direction - u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_elements) - gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) - flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) + viscous_container = init_viscous_container_2d(nvariables(equations_hyperbolic), + nnodes(dg.basis), nelements(elements), + uEltype) - cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + cache = (; elements, interfaces, boundaries, viscous_container) return cache end @@ -28,7 +25,8 @@ function rhs_parabolic!(du, u, t, mesh::P4estMesh{2}, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions_parabolic, source_terms, dg::DG, parabolic_scheme, cache, cache_parabolic) - (; u_transformed, gradients, flux_viscous) = cache_parabolic + @unpack viscous_container = cache_parabolic + @unpack u_transformed, gradients, flux_viscous = viscous_container # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 2d26c1aff50..b06cdd42127 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -10,14 +10,11 @@ function create_cache_parabolic(mesh::P4estMesh{3}, equations_hyperbolic::Abstra interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements) boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements) - n_vars = nvariables(equations_hyperbolic) - n_elements = nelements(elements) - n_nodes = nnodes(dg.basis) # nodes in one direction - u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) - gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) - flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) + viscous_container = init_viscous_container_3d(nvariables(equations_hyperbolic), + nnodes(dg.basis), nelements(elements), + uEltype) - cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + cache = (; elements, interfaces, boundaries, viscous_container) return cache end @@ -36,7 +33,8 @@ function rhs_parabolic!(du, u, t, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions_parabolic, source_terms, dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack u_transformed, gradients, flux_viscous = cache_parabolic + @unpack viscous_container = cache_parabolic + @unpack u_transformed, gradients, flux_viscous = viscous_container # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin diff --git a/src/solvers/dgsem_tree/container_viscous_1d.jl b/src/solvers/dgsem_tree/container_viscous_1d.jl index a4919f75396..71c68dfc6df 100644 --- a/src/solvers/dgsem_tree/container_viscous_1d.jl +++ b/src/solvers/dgsem_tree/container_viscous_1d.jl @@ -19,9 +19,9 @@ mutable struct ViscousContainer1D{uEltype <: Real} end end -function init_viscous_container(n_vars::Integer, n_nodes::Integer, - n_elements::Integer, - ::Type{uEltype}) where {uEltype <: Real} +function init_viscous_container_1d(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} return ViscousContainer1D{uEltype}(n_vars, n_nodes, n_elements) end diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl new file mode 100644 index 00000000000..bd7ff413af5 --- /dev/null +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -0,0 +1,69 @@ +mutable struct ViscousContainer2D{uEltype <: Real} + u_transformed::Array{uEltype, 4} + # Using an outer fixed-size datastructure leads to nasty implementations, + # see https://github.com/trixi-framework/Trixi.jl/pull/1629#discussion_r1355293953. + # Also: This does not result in speed up compared to using tuples for the internal + # datastructures, see + # https://github.com/trixi-framework/Trixi.jl/pull/1629#discussion_r1363352188. + gradients::Vector{Array{uEltype, 4}} + flux_viscous::Vector{Array{uEltype, 4}} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + # Use Tuple for outer, fixed-size datastructure + _gradients::Tuple{Vector{uEltype}, Vector{uEltype}} + _flux_viscous::Tuple{Vector{uEltype}, Vector{uEltype}} + + function ViscousContainer2D{uEltype}(n_vars::Integer, n_nodes::Integer, + n_elements::Integer) where {uEltype <: Real} + new(Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements), + [Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements) for _ in 1:2], + [Array{uEltype, 4}(undef, n_vars, n_nodes, n_nodes, n_elements) for _ in 1:2], + Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), + (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements)), + (Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements))) + end +end + +function init_viscous_container_2d(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} + return ViscousContainer2D{uEltype}(n_vars, n_nodes, n_elements) +end + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(viscous_container::ViscousContainer2D, equations, dg, cache) + capacity = nvariables(equations) * nnodes(dg) * nnodes(dg) * nelements(dg, cache) + resize!(viscous_container._u_transformed, capacity) + for dim in 1:2 + resize!(viscous_container._gradients[dim], capacity) + resize!(viscous_container._flux_viscous[dim], capacity) + end + + viscous_container.u_transformed = unsafe_wrap(Array, + pointer(viscous_container._u_transformed), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + for dim in 1:2 + viscous_container.gradients[dim] = unsafe_wrap(Array, + pointer(viscous_container._gradients[dim]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + viscous_container.flux_viscous[dim] = unsafe_wrap(Array, + pointer(viscous_container._flux_viscous[dim]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + end + return nothing +end diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl new file mode 100644 index 00000000000..64d283fe189 --- /dev/null +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -0,0 +1,75 @@ +mutable struct ViscousContainer3D{uEltype <: Real} + u_transformed::Array{uEltype, 5} + # Using an outer fixed-size datastructure leads to nasty implementations, + # see https://github.com/trixi-framework/Trixi.jl/pull/1629#discussion_r1355293953. + # Also: This does not result in speed up compared to using tuples for the internal + # datastructures, see + # https://github.com/trixi-framework/Trixi.jl/pull/1629#discussion_r1363352188. + gradients::Vector{Array{uEltype, 5}} + flux_viscous::Vector{Array{uEltype, 5}} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + # Use Tuple for outer, fixed-size datastructure + _gradients::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}} + _flux_viscous::Tuple{Vector{uEltype}, Vector{uEltype}, Vector{uEltype}} + + function ViscousContainer3D{uEltype}(n_vars::Integer, n_nodes::Integer, + n_elements::Integer) where {uEltype <: Real} + new(Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements), + [Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) + for _ in 1:3], + [Array{uEltype, 5}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) + for _ in 1:3], + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + (Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements)), + (Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements))) + end +end + +function init_viscous_container_3d(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} + return ViscousContainer3D{uEltype}(n_vars, n_nodes, n_elements) +end + +# Only one-dimensional `Array`s are `resize!`able in Julia. +# Hence, we use `Vector`s as internal storage and `resize!` +# them whenever needed. Then, we reuse the same memory by +# `unsafe_wrap`ping multi-dimensional `Array`s around the +# internal storage. +function Base.resize!(viscous_container::ViscousContainer3D, equations, dg, cache) + capacity = nvariables(equations) * nnodes(dg) * nnodes(dg) * nnodes(dg) * + nelements(dg, cache) + resize!(viscous_container._u_transformed, capacity) + for dim in 1:3 + resize!(viscous_container._gradients[dim], capacity) + resize!(viscous_container._flux_viscous[dim], capacity) + end + + viscous_container.u_transformed = unsafe_wrap(Array, + pointer(viscous_container._u_transformed), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + for dim in 1:3 + viscous_container.gradients[dim] = unsafe_wrap(Array, + pointer(viscous_container._gradients[dim]), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + viscous_container.flux_viscous[dim] = unsafe_wrap(Array, + pointer(viscous_container._flux_viscous[dim]), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nnodes(dg), + nelements(dg, cache))) + end + return nothing +end diff --git a/src/solvers/dgsem_tree/containers.jl b/src/solvers/dgsem_tree/containers.jl index bba8b83b23a..3f05daf81d8 100644 --- a/src/solvers/dgsem_tree/containers.jl +++ b/src/solvers/dgsem_tree/containers.jl @@ -28,9 +28,11 @@ function reinitialize_containers!(mesh::TreeMesh, equations, dg::DGSEM, cache) init_boundaries!(boundaries, elements, mesh) # re-initialize mortars container - @unpack mortars = cache - resize!(mortars, count_required_mortars(mesh, leaf_cell_ids)) - init_mortars!(mortars, elements, mesh) + if hasproperty(cache, :mortars) # cache_parabolic does not carry mortars + @unpack mortars = cache + resize!(mortars, count_required_mortars(mesh, leaf_cell_ids)) + init_mortars!(mortars, elements, mesh) + end if mpi_isparallel() # re-initialize mpi_interfaces container diff --git a/src/solvers/dgsem_tree/containers_viscous.jl b/src/solvers/dgsem_tree/containers_viscous.jl new file mode 100644 index 00000000000..444f2cb7303 --- /dev/null +++ b/src/solvers/dgsem_tree/containers_viscous.jl @@ -0,0 +1,4 @@ +# Dimension-specific implementations +include("container_viscous_1d.jl") +include("container_viscous_2d.jl") +include("container_viscous_3d.jl") diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index ff37bad3b3a..ef9a42b4c1a 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -54,8 +54,8 @@ include("containers.jl") # Dimension-agnostic parallel setup include("dg_parallel.jl") -# Helper struct for parabolic AMR -include("container_viscous_1d.jl") +# Helper structs for parabolic AMR +include("containers_viscous.jl") # 1D DG implementation include("dg_1d.jl") diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 97e31e0e22b..90007b05b3d 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -535,14 +535,14 @@ function create_cache_parabolic(mesh::TreeMesh{1}, elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) - viscous_container = init_viscous_container(nvariables(equations_hyperbolic), - nnodes(elements), nelements(elements), - uEltype) - interfaces = init_interfaces(leaf_cell_ids, mesh, elements) boundaries = init_boundaries(leaf_cell_ids, mesh, elements) + viscous_container = init_viscous_container_1d(nvariables(equations_hyperbolic), + nnodes(elements), nelements(elements), + uEltype) + cache = (; elements, interfaces, boundaries, viscous_container) return cache diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 1c32703c7c3..06abff5e85b 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -17,7 +17,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions_parabolic, source_terms, dg::DG, parabolic_scheme, cache, cache_parabolic) - (; u_transformed, gradients, flux_viscous) = cache_parabolic + @unpack viscous_container = cache_parabolic + @unpack u_transformed, gradients, flux_viscous = viscous_container # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin @@ -290,7 +291,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, +function calc_viscous_fluxes!(flux_viscous, + gradients, u_transformed, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) @@ -513,11 +515,16 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra return nothing end -function prolong2mortars!(cache, flux_viscous::Tuple{AbstractArray, AbstractArray}, +# `cache` is the hyperbolic cache, i.e., in particular not `cache_parabolic`. +# This is because mortar handling is done in the (hyperbolic) `cache`. +# Specialization `flux_viscous::Vector{Array{uEltype, 4}}` needed since +#`prolong2mortars!` in dg_2d.jl is used for both purely hyperbolic and +# hyperbolic-parabolic systems. +function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, mortar_l2::LobattoLegendreMortarL2, surface_integral, - dg::DGSEM) + dg::DGSEM) where {uEltype <: Real} flux_viscous_x, flux_viscous_y = flux_viscous @threaded for mortar in eachmortar(dg, cache) large_element = cache.mortars.neighbor_ids[3, mortar] @@ -909,21 +916,18 @@ function create_cache_parabolic(mesh::TreeMesh{2}, elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) - n_vars = nvariables(equations_hyperbolic) - n_nodes = nnodes(elements) - n_elements = nelements(elements) - u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_elements) - gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) - flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) - interfaces = init_interfaces(leaf_cell_ids, mesh, elements) boundaries = init_boundaries(leaf_cell_ids, mesh, elements) # mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar) + viscous_container = init_viscous_container_2d(nvariables(equations_hyperbolic), + nnodes(elements), nelements(elements), + uEltype) + # cache = (; elements, interfaces, boundaries, mortars) - cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + cache = (; elements, interfaces, boundaries, viscous_container) # Add specialized parts of the cache required to compute the mortars etc. # cache = (;cache..., create_cache(mesh, equations_parabolic, dg.mortar, uEltype)...) diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 37492dbcb91..2561c5fe5b0 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -17,7 +17,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions_parabolic, source_terms, dg::DG, parabolic_scheme, cache, cache_parabolic) - @unpack u_transformed, gradients, flux_viscous = cache_parabolic + @unpack viscous_container = cache_parabolic + @unpack u_transformed, gradients, flux_viscous = viscous_container # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin @@ -338,7 +339,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, +function calc_viscous_fluxes!(flux_viscous, + gradients, u_transformed, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) @@ -596,13 +598,17 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra return nothing end +# `cache` is the hyperbolic cache, i.e., in particular not `cache_parabolic`. +# This is because mortar handling is done in the (hyperbolic) `cache`. +# Specialization `flux_viscous::Vector{Array{uEltype, 4}}` needed since +#`prolong2mortars!` in dg_2d.jl is used for both purely hyperbolic and +# hyperbolic-parabolic systems. function prolong2mortars!(cache, - flux_viscous::Tuple{AbstractArray, AbstractArray, - AbstractArray}, + flux_viscous::Vector{Array{uEltype, 5}}, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, mortar_l2::LobattoLegendreMortarL2, - surface_integral, dg::DGSEM) + surface_integral, dg::DGSEM) where {uEltype <: Real} # temporary buffer for projections @unpack fstar_tmp1_threaded = cache @@ -1099,21 +1105,18 @@ function create_cache_parabolic(mesh::TreeMesh{3}, elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) - n_vars = nvariables(equations_hyperbolic) - n_nodes = nnodes(elements) - n_elements = nelements(elements) - u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) - gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) - flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) - interfaces = init_interfaces(leaf_cell_ids, mesh, elements) boundaries = init_boundaries(leaf_cell_ids, mesh, elements) # mortars = init_mortars(leaf_cell_ids, mesh, elements, dg.mortar) + viscous_container = init_viscous_container_3d(nvariables(equations_hyperbolic), + nnodes(elements), nelements(elements), + uEltype) + # cache = (; elements, interfaces, boundaries, mortars) - cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + cache = (; elements, interfaces, boundaries, viscous_container) # Add specialized parts of the cache required to compute the mortars etc. # cache = (;cache..., create_cache(mesh, equations_parabolic, dg.mortar, uEltype)...) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 6424dde7d31..733f89c2158 100644 --- a/src/time_integration/methods_SSP.jl +++ b/src/time_integration/methods_SSP.jl @@ -148,7 +148,7 @@ function solve!(integrator::SimpleIntegratorSSP) callbacks = integrator.opts.callback integrator.finalstep = false - while !integrator.finalstep + @trixi_timeit timer() "main loop" while !integrator.finalstep if isnan(integrator.dt) error("time step size `dt` is NaN") end diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 31dfe1d35a5..5baa090c8d7 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -17,27 +17,42 @@ isdir(outdir) && rm(outdir, recursive=true) # Expected errors are exactly the same as with TreeMesh! l2 = [8.311947673061856e-6], linf = [6.627000273229378e-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 @trixi_testset "elixir_advection_nonconforming_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming_flag.jl"), l2 = [3.198940059144588e-5], linf = [0.00030636069494005547]) - - # 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 + # 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_advection_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"), l2 = [0.0005379687442422346], linf = [0.007438525029884735]) + # 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_advection_amr_solution_independent.jl" begin @@ -46,6 +61,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [4.949660644033807e-5], linf = [0.0004867846262313763], coverage_override = (maxiters=6,)) + # 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_advection_amr_unstructured_flag.jl" begin @@ -53,18 +76,42 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.0012766060609964525], linf = [0.01750280631586159], coverage_override = (maxiters=6,)) + # 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_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2 = [4.507575525876275e-6], linf = [6.21489667023134e-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 @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), l2 = [0.0034516244508588046, 0.0023420334036925493, 0.0024261923964557187, 0.004731710454271893], linf = [0.04155789011775046, 0.024772109862748914, 0.03759938693042297, 0.08039824959535657]) + # 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_free_stream.jl" begin @@ -73,6 +120,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [1.9539925233402755e-14, 2e-12, 4.8e-12, 4e-12], atol = 2.0e-12, # required to make CI tests pass on macOS ) + # 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_shockcapturing_ec.jl" begin @@ -80,6 +135,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [9.53984675e-02, 1.05633455e-01, 1.05636158e-01, 3.50747237e-01], linf = [2.94357464e-01, 4.07893014e-01, 3.97334516e-01, 1.08142520e+00], tspan = (0.0, 1.0)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_sedov.jl" begin @@ -87,6 +150,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [3.76149952e-01, 2.46970327e-01, 2.46970327e-01, 1.28889042e+00], linf = [1.22139001e+00, 1.17742626e+00, 1.17742626e+00, 6.20638482e+00], tspan = (0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_blast_wave_amr.jl" begin @@ -95,6 +166,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [2.76020890e+00, 2.32659890e+00, 2.32580837e+00, 2.15778188e+00], tspan = (0.0, 0.3), coverage_override = (maxiters=6,)) + # 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_wall_bc_amr.jl" begin @@ -102,15 +181,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.020291447969983396, 0.017479614254319948, 0.011387644425613437, 0.0514420126021293], linf = [0.3582779022370579, 0.32073537890751663, 0.221818049107692, 0.9209559420400415], tspan = (0.0, 0.15)) - - # 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 + # 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_forward_step_amr.jl" begin @@ -120,6 +198,16 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.0001), rtol = 1.0e-7, skip_coverage=true) + if @isdefined sol # Skipped in coverage run + # 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_euler_double_mach_amr.jl" begin @@ -128,6 +216,16 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [6.902000373057003, 53.95714139820832, 24.241610279839758, 561.0630401858057], tspan = (0.0, 0.0001), skip_coverage=true) + if @isdefined sol # Skipped in coverage run + # 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_euler_supersonic_cylinder.jl" begin @@ -136,6 +234,16 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [3.653905721692421, 4.285035711361009, 6.8544353186357645, 31.748244912257533], tspan = (0.0, 0.001), skip_coverage=true) + if @isdefined sol # Skipped in coverage run + # 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_eulergravity_convergence.jl" begin @@ -143,6 +251,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697], linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758], tspan = (0.0, 0.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)) < 1000 + end end @trixi_testset "elixir_shallowwater_source_terms.jl" begin @@ -150,6 +266,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [9.168126407325352e-5, 0.0009795410115453788, 0.002546408320320785, 3.941189812642317e-6], linf = [0.0009903782521019089, 0.0059752684687262025, 0.010941106525454103, 1.2129488214718265e-5], tspan = (0.0, 0.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)) < 1000 + end end @trixi_testset "elixir_mhd_alfven_wave.jl" begin @@ -160,6 +284,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [4.255466285174592e-5, 1.0029706745823264e-5, 1.0029706747467781e-5, 1.2122265939010224e-5, 5.4791097160444835e-6, 5.18922042269665e-6, 5.189220422141538e-6, 9.552667261422676e-6, 1.4237578427628152e-6]) + # 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_mhd_rotor.jl" begin @@ -171,12 +303,28 @@ isdir(outdir) && rm(outdir, recursive=true) 19.647239392277378, 1.3938810140985936, 1.8724965294853084, 0.0, 0.0016290067532561904], tspan = (0.0, 0.02)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_linearizedeuler_gaussian_source.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_linearizedeuler_gaussian_source.jl"), l2 = [0.006047938590548741, 0.0040953286019907035, 0.004222698522497298, 0.006269492499336128], linf = [0.06386175207349379, 0.0378926444850457, 0.041759728067967065, 0.06430136016259067]) + # 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 diff --git a/test/test_paper_self_gravitating_gas_dynamics.jl b/test/test_paper_self_gravitating_gas_dynamics.jl index 6a35543fe47..68aa2992601 100644 --- a/test/test_paper_self_gravitating_gas_dynamics.jl +++ b/test/test_paper_self_gravitating_gas_dynamics.jl @@ -17,6 +17,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence.jl"), l2 = [0.0001740977055972079, 0.0003369355182519592, 0.0003369355182518708, 0.0006099171220334989], linf = [0.001079347149189669, 0.0018836938381321389, 0.001883693838132583, 0.003971575376718217]) + # 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_convergence.jl with polydeg=4" begin @@ -24,6 +32,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [1.7187201161597772e-5, 2.678065111772951e-5, 2.678065111783027e-5, 4.952504160091526e-5], linf = [0.0001501749544159381, 0.00016549482504535362, 0.00016549482504601976, 0.0004372960291432193], polydeg = 4) + # 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 @@ -31,6 +47,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_convergence.jl"), l2 = [0.003154024896093942, 0.012394432074951856, 0.02185973823794725], linf = [0.01731850928579215, 0.07843510773347553, 0.11242300176349201]) + # 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_hypdiff_convergence.jl with polydeg=4" begin @@ -38,6 +62,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [0.0002511283012128458, 0.0008808243846610255, 0.0016313343228567005], linf = [0.0017290715087938668, 0.003129184465704738, 0.01000728849316701], polydeg = 4) + # 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 @@ -46,6 +78,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697], linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758], tspan = (0.0, 0.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)) < 1000 + end end @trixi_testset "elixir_eulergravity_convergence.jl with polydeg=4" begin @@ -53,6 +93,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [1.9537712148648045e-5, 2.7564396197947587e-5, 2.7564396197967635e-5, 5.688838772067586e-5], linf = [0.00012335710672761735, 0.00020086268350816283, 0.00020086268350727465, 0.0004962155455632278], tspan = (0.0, 0.1), polydeg = 4) + # 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_eulergravity_convergence.jl with 1st order RK3S*" begin @@ -60,6 +108,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [0.00024871265138959434, 0.000337007710281087, 0.0003370077102811394, 0.0007231525515231289], linf = [0.0015813032941613958, 0.002049428843978518, 0.0020494288439798503, 0.004793821198143977], tspan = (0.0, 0.1), timestep_gravity=Trixi.timestep_gravity_erk51_3Sstar!) + # 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_eulergravity_convergence.jl with 3rd order RK3S*" begin @@ -67,6 +123,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [0.0002487126513894034, 0.00033700771023049785, 0.00033700771023048245, 0.0007231525514158737], linf = [0.0015813032943847727, 0.002049428842844314, 0.0020494288428452023, 0.004793821195971937], tspan = (0.0, 0.1), timestep_gravity=Trixi.timestep_gravity_erk53_3Sstar!) + # 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 @@ -77,6 +141,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam tspan = (0.0, 0.1), atol = 4.0e-6 # the background field is reatively large, so this corresponds to our usual atol ) + # 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_eulergravity_jeans_instability.jl with RK3S*" begin @@ -91,6 +163,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam resid_tol=1.0e-4, n_iterations_max=1000, timestep_gravity=timestep_gravity_erk52_3Sstar!)) + # 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 "Printing" begin @@ -108,6 +188,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam show(stdout, semi_euler.boundary_conditions) show(stdout, TrivialCallback()) show(stdout, equations_euler) + # 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 @@ -117,6 +205,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam linf = [2.3874843337593776, 4.07876384374792, 4.07876384374792, 16.23914384809855], tspan = (0.0, 0.05), coverage_override = (maxiters=2,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_eulergravity_sedov_blast_wave.jl with ref-level=8 and no AMR" begin @@ -124,6 +220,14 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "paper_self_gravitating_gas_dynam l2 = [0.00289222135995042, 0.013724813590853825, 0.013724813590853832, 0.05822904710548214], linf = [0.26361780693997594, 1.3908873830688688, 1.3908873830688688, 4.066701303607613], tspan = (0.0, 0.005), initial_refinement_level=8, amr_callback=TrivialCallback()) + # 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 diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 046494e1000..a57462ef7ea 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -377,6 +377,14 @@ isdir(outdir) && rm(outdir, recursive=true) end end + @trixi_testset "TreeMesh2D: elixir_navierstokes_shearlayer_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_navierstokes_shearlayer_amr.jl"), + l2 = [0.00526017743452336, 0.4130430692895672, 0.4310996183791349, 1.1544344171604635], + linf = [0.03492185879198495, 1.392635891671335, 1.357551616406459, 8.713760873018146], + tspan = (0.0, 0.7) + ) + end + @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_advection_diffusion_periodic.jl"), trees_per_dimension = (1, 1), initial_refinement_level = 2, tspan=(0.0, 0.5), diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index f74546d0146..276e37518ee 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -252,6 +252,20 @@ isdir(outdir) && rm(outdir, recursive=true) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "TreeMesh3D: elixir_advection_diffusion_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_advection_diffusion_amr.jl"), + l2 = [0.000355780485397024], + linf = [0.0010810770271614256] + ) + end + + @trixi_testset "TreeMesh3D: elixir_advection_diffusion_nonperiodic.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_3d_dgsem", "elixir_advection_diffusion_nonperiodic.jl"), + l2 = [0.0009808996243280868], + linf = [0.01732621559135459] + ) + end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 6d528abc7af..c3192f52b2c 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -19,6 +19,15 @@ isdir(outdir) && rm(outdir, recursive=true) # Expected errors are exactly the same as with TreeMesh! l2 = [8.311947673061856e-6], linf = [6.627000273229378e-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_advection_coupled.jl" begin @@ -31,6 +40,14 @@ isdir(outdir) && rm(outdir, recursive=true) errors = analysis_callback(sol) @test errors.l2 ≈ [7.816742843181738e-6, 7.816742843196112e-6] rtol=1.0e-4 @test errors.linf ≈ [6.314906965543265e-5, 6.314906965410039e-5] rtol=1.0e-4 + # 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 @@ -38,6 +55,14 @@ isdir(outdir) && rm(outdir, recursive=true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"), l2 = [4.220397559713772e-6], linf = [3.477948874874848e-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 @trixi_testset "elixir_advection_extended.jl with polydeg=4" begin @@ -48,6 +73,14 @@ isdir(outdir) && rm(outdir, recursive=true) cells_per_dimension = (16, 23), polydeg = 4, cfl = 1.4) + # 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 @testset "elixir_advection_rotated.jl" begin @@ -57,6 +90,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [8.311947673061856e-6], linf = [6.627000273229378e-5], alpha = 0.0) + # 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_advection_rotated.jl with α = 0.1" begin @@ -65,6 +106,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [8.3122750550501e-6], linf = [6.626802581322089e-5], alpha = 0.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)) < 1000 + end end @trixi_testset "elixir_advection_rotated.jl with α = 0.5 * pi" begin @@ -73,6 +122,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [8.311947673061856e-6], linf = [6.627000273229378e-5], alpha = 0.5 * pi) + # 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 @@ -81,30 +138,70 @@ isdir(outdir) && rm(outdir, recursive=true) # Expected errors are exactly the same as in elixir_advection_basic! l2 = [8.311947673061856e-6], linf = [6.627000273229378e-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 @trixi_testset "elixir_advection_waving_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_waving_flag.jl"), l2 = [0.00018553859900545866], linf = [0.0016167719118129753]) + # 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_advection_free_stream.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_free_stream.jl"), l2 = [6.8925194184204476e-15], linf = [9.903189379656396e-14]) + # 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_advection_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"), l2 = [0.00025552740731641223], linf = [0.007252625722805939]) + # 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_advection_restart.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"), l2 = [4.219208035582454e-6], linf = [3.438434404412494e-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 @trixi_testset "elixir_advection_restart.jl with waving flag mesh" begin @@ -114,6 +211,14 @@ isdir(outdir) && rm(outdir, recursive=true) rtol = 5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS) elixir_file="elixir_advection_waving_flag.jl", restart_file="restart_000021.h5") + # 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_advection_restart.jl with free stream mesh" begin @@ -122,6 +227,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [1.0857981180834031e-13], elixir_file="elixir_advection_free_stream.jl", restart_file="restart_000036.h5") + # 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_source_terms.jl" begin @@ -129,6 +242,14 @@ isdir(outdir) && rm(outdir, recursive=true) # Expected errors are exactly the same as with TreeMesh! l2 = [9.321181253186009e-7, 1.4181210743438511e-6, 1.4181210743487851e-6, 4.824553091276693e-6], linf = [9.577246529612893e-6, 1.1707525976012434e-5, 1.1707525976456523e-5, 4.8869615580926506e-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 @testset "elixir_euler_source_terms_rotated.jl" begin @@ -138,7 +259,15 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [9.321181253186009e-7, 1.4181210743438511e-6, 1.4181210743487851e-6, 4.824553091276693e-6], linf = [9.577246529612893e-6, 1.1707525976012434e-5, 1.1707525976456523e-5, 4.8869615580926506e-5], alpha = 0.0) - end + # 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_source_terms_rotated.jl with α = 0.1" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"), @@ -146,7 +275,15 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [9.321188057029291e-7, 1.3195106906473365e-6, 1.510307360354032e-6, 4.82455408101712e-6], linf = [9.57723626271445e-6, 1.0480225511866337e-5, 1.2817828088262928e-5, 4.886962393513272e-5], alpha = 0.1) - end + # 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_source_terms_rotated.jl with α = 0.2 * pi" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"), @@ -154,7 +291,15 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [9.32127973957391e-7, 8.477824799744325e-7, 1.8175286311402784e-6, 4.824562453521076e-6], linf = [9.576898420737834e-6, 5.057704352218195e-6, 1.635260719945464e-5, 4.886978754825577e-5], alpha = 0.2 * pi) - end + # 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_source_terms_rotated.jl with α = 0.5 * pi" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"), @@ -162,19 +307,42 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [9.321181253186009e-7, 1.4181210743438511e-6, 1.4181210743487851e-6, 4.824553091276693e-6], linf = [9.577246529612893e-6, 1.1707525976012434e-5, 1.1707525976456523e-5, 4.8869615580926506e-5], alpha = 0.5 * pi) - end + # 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_source_terms_parallelogram.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_parallelogram.jl"), l2 = [1.1167802955144833e-5, 1.0805775514153104e-5, 1.953188337010932e-5, 5.5033856574857146e-5], linf = [8.297006495561199e-5, 8.663281475951301e-5, 0.00012264160606778596, 0.00041818802502024965]) + # 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_source_terms_waving_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_waving_flag.jl"), l2 = [2.991891317562739e-5, 3.6063177168283174e-5, 2.7082941743640572e-5, 0.00011414695350996946], linf = [0.0002437454930492855, 0.0003438936171968887, 0.00024217622945688078, 0.001266380414757684]) + # 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_free_stream.jl" begin @@ -182,6 +350,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14], linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 6.502853810985698e-13, 2.7000623958883807e-13], atol = 7.0e-13) + # 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_free_stream.jl with FluxRotated(flux_lax_friedrichs)" begin @@ -190,12 +366,28 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14], linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 6.502853810985698e-13, 2.7000623958883807e-13], atol = 7.0e-13) + # 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_source_terms_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"), l2 = [2.259440511901724e-6, 2.3188881559075347e-6, 2.3188881559568146e-6, 6.332786324137878e-6], linf = [1.4987382622067003e-5, 1.918201192063762e-5, 1.918201192019353e-5, 6.052671713430158e-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 @trixi_testset "elixir_euler_ec.jl" begin @@ -203,6 +395,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.03774907669925568, 0.02845190575242045, 0.028262802829412605, 0.13785915638851698], linf = [0.3368296929764073, 0.27644083771519773, 0.27990039685141377, 1.1971436487402016], tspan = (0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_sedov.jl" begin @@ -210,6 +410,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [3.69856202e-01, 2.35242180e-01, 2.41444928e-01, 1.28807120e+00], linf = [1.82786223e+00, 1.30452904e+00, 1.40347257e+00, 6.21791658e+00], tspan = (0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_euler_rayleigh_taylor_instability.jl" begin @@ -218,30 +426,70 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [0.4799214336153155, 0.024595483032220266, 0.02059808120543466, 0.03190756362943725], cells_per_dimension = (8,8), tspan = (0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"), l2 = [0.0016688820596537988, 0.0025921681885685425, 0.003280950351435014], linf = [0.010994679664394269, 0.01331197845637, 0.020080117011346488]) + # 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_eulerpolytropic_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_ec.jl"), l2 = [0.03647890611450939, 0.025284915444045052, 0.025340697771609126], linf = [0.32516731565355583, 0.37509762516540046, 0.29812843284727336]) + # 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_eulerpolytropic_isothermal_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_isothermal_wave.jl"), l2 = [0.004998778491726366, 0.004998916000294425, 9.259136963058664e-17], linf = [0.010001103673834888, 0.010051165098399503, 7.623942913643681e-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 @trixi_testset "elixir_eulerpolytropic_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_wave.jl"), l2 = [0.23642682112204072, 0.20904264390331334, 8.174982691297391e-17], linf = [0.4848250368349989, 0.253350873815695, 4.984552457753618e-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 @trixi_testset "elixir_hypdiff_nonperiodic.jl" begin @@ -250,6 +498,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [1.0771947577311836, 1.9143913544309838, 2.149549109115789], tspan = (0.0, 0.1), coverage_override = (polydeg=3,)) # Prevent long compile time in CI + # 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 @trixi_testset "elixir_hypdiff_harmonic_nonperiodic.jl" begin @@ -258,6 +514,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [0.35026352556630114, 0.8344372248051408, 0.8344372248051408], tspan = (0.0, 0.1), coverage_override = (polydeg=3,)) # Prevent long compile time in CI + # 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_mhd_ec.jl" begin @@ -269,6 +533,14 @@ isdir(outdir) && rm(outdir, recursive=true) 0.9757376320946505, 0.12123736788315098, 0.12837436699267113, 0.17793825293524734, 0.03460761690059514], tspan = (0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end @trixi_testset "elixir_mhd_alfven_wave.jl" begin @@ -280,6 +552,14 @@ isdir(outdir) && rm(outdir, recursive=true) 0.03817506476831778, 0.042847170999492534, 0.03761563456810613, 0.048184237474911844, 0.04114666955364693], tspan = (0.0, 1.0)) + # 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_shallowwater_source_terms.jl" begin @@ -287,6 +567,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.0017285599436729316, 0.025584610912606776, 0.028373834961180594, 6.274146767730866e-5], linf = [0.012972309788264802, 0.108283714215621, 0.15831585777928936, 0.00018196759554722775], 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 "elixir_shallowwater_well_balanced.jl" begin @@ -294,6 +582,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.7920927046419308, 9.92129670988898e-15, 1.0118635033124588e-14, 0.7920927046419308], linf = [2.408429868800133, 5.5835419986809516e-14, 5.448874313931364e-14, 2.4084298688001335], tspan = (0.0, 0.25)) + # 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_shallowwater_well_balanced_wet_dry.jl" begin @@ -301,6 +597,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.019731646454942086, 1.0694532773278277e-14, 1.1969913383405568e-14, 0.0771517260037954], linf = [0.4999999999998892, 6.067153702623552e-14, 4.4849667259339357e-14, 1.9999999999999993], tspan = (0.0, 0.25)) + # 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_shallowwater_conical_island.jl" begin @@ -308,6 +612,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.04593154164306353, 0.1644534881916908, 0.16445348819169076, 0.0011537702354532122], linf = [0.21100717610846442, 0.9501592344310412, 0.950159234431041, 0.021790250683516296], tspan = (0.0, 0.025)) + # 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_shallowwater_parabolic_bowl.jl" begin @@ -315,6 +627,14 @@ isdir(outdir) && rm(outdir, recursive=true) l2 = [0.00015285369980313484, 1.9536806395943226e-5, 9.936906607758672e-5, 5.0686313334616055e-15], linf = [0.003316119030459211, 0.0005075409427972817, 0.001986721761060583, 4.701794509287538e-14], tspan = (0.0, 0.025), cells_per_dimension = (40, 40)) + # 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_mhd_ec_shockcapturing.jl" begin @@ -325,6 +645,14 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [0.25144373906033013, 0.32881947152723745, 0.3053266801502693, 0.20989755319972866, 0.9927517314507455, 0.1105172121361323, 0.1257708104676617, 0.1628334844841588, 0.02624301627479052]) + # 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 diff --git a/test/test_tree_1d_hypdiff.jl b/test/test_tree_1d_hypdiff.jl index 560f77b2a13..19200b26892 100644 --- a/test/test_tree_1d_hypdiff.jl +++ b/test/test_tree_1d_hypdiff.jl @@ -14,6 +14,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") l2 = [1.3655114954641076e-7, 1.0200345025539218e-6], linf = [7.173286515893551e-7, 4.507116363683394e-6], atol = 2.5e-13) + # 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_hypdiff_harmonic_nonperiodic.jl" begin @@ -21,6 +29,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") l2 = [3.0130941075207524e-12, 2.6240829677090014e-12], linf = [4.054534485931072e-12, 3.8826719617190975e-12], atol = 2.5e-13) + # 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)) < 10000 + end end end diff --git a/test/test_tree_2d_advection.jl b/test/test_tree_2d_advection.jl index 36cb1e882cc..365011eef53 100644 --- a/test/test_tree_2d_advection.jl +++ b/test/test_tree_2d_advection.jl @@ -15,6 +15,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") linf = [6.627000273229378e-5], # Let the small basic test run to the end coverage_override = (maxiters=10^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 @trixi_testset "elixir_advection_extended.jl with polydeg=1" begin @@ -22,6 +30,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") l2 = [0.02134571266411136], linf = [0.04347734797775926], polydeg=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)) < 1000 + end end @trixi_testset "elixir_advection_restart.jl" begin @@ -49,14 +65,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") l2 = [0.0015188466707237375], linf = [0.008446655719187679]) - # 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 + # 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_advection_amr.jl" begin @@ -66,6 +82,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") linf = [0.00045263895394385967], # Let this test run to the end to cover some AMR code coverage_override = (maxiters=10^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 @trixi_testset "elixir_advection_amr_nonperiodic.jl" begin @@ -74,6 +98,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") l2 = [3.2207388565869075e-5], linf = [0.0007508059772436404], coverage_override = (maxiters=6,)) + # 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_advection_amr_solution_independent.jl" begin @@ -81,6 +113,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") l2 = [4.949660644033807e-5], linf = [0.0004867846262313763], coverage_override = (maxiters=6,)) + # 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_advection_amr_visualization.jl" begin @@ -114,6 +154,15 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") # Let this test terminate by time instead of maxiters to cover some lines # in time_integration/methods_2N.jl coverage_override = (maxiters=10^5, tspan=(0.0, 0.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_advection_timeintegration.jl with carpenter_kennedy_erk43" begin @@ -122,6 +171,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") linf = [0.0005437136621948904], ode_algorithm=Trixi.CarpenterKennedy2N43(), cfl = 1.0) + # 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 @trixi_testset "elixir_advection_timeintegration.jl with carpenter_kennedy_erk43 with maxiters=1" begin @@ -131,6 +188,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") ode_algorithm=Trixi.CarpenterKennedy2N43(), cfl = 1.0, maxiters = 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 @trixi_testset "elixir_advection_timeintegration.jl with parsani_ketcheson_deconinck_erk94" begin @@ -138,6 +203,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") l2 = [2.4976673477385313e-5], linf = [0.0005534166916640881], ode_algorithm=Trixi.ParsaniKetchesonDeconinck3Sstar94()) + # 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 @trixi_testset "elixir_advection_timeintegration.jl with parsani_ketcheson_deconinck_erk32" begin @@ -146,6 +219,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") linf = [0.0005799465470165757], ode_algorithm=Trixi.ParsaniKetchesonDeconinck3Sstar32(), cfl = 1.0) + # 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 @trixi_testset "elixir_advection_timeintegration.jl with parsani_ketcheson_deconinck_erk32 with maxiters=1" begin @@ -155,14 +236,29 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") ode_algorithm=Trixi.ParsaniKetchesonDeconinck3Sstar32(), cfl = 1.0, maxiters = 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 @trixi_testset "elixir_advection_callbacks.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_callbacks.jl"), l2 = [8.311947673061856e-6], linf = [6.627000273229378e-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 # Coverage test for all initial conditions @testset "Linear scalar advection: Tests for initial conditions" begin @@ -173,6 +269,14 @@ end linf = [0.0007140570281718439], maxiters = 1, initial_condition = Trixi.initial_condition_sin_sin) + # 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_advection_extended.jl with initial_condition_constant" begin @@ -181,6 +285,14 @@ end linf = [1.3322676295501878e-15], maxiters = 1, initial_condition = initial_condition_constant) + # 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_advection_extended.jl with initial_condition_linear_x_y" begin @@ -191,6 +303,14 @@ end initial_condition = Trixi.initial_condition_linear_x_y, boundary_conditions = Trixi.boundary_condition_linear_x_y, periodicity=false) + # 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_advection_extended.jl with initial_condition_linear_x" begin @@ -201,6 +321,14 @@ end initial_condition = Trixi.initial_condition_linear_x, boundary_conditions = Trixi.boundary_condition_linear_x, periodicity=false) + # 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_advection_extended.jl with initial_condition_linear_y" begin @@ -211,6 +339,14 @@ end initial_condition = Trixi.initial_condition_linear_y, boundary_conditions = Trixi.boundary_condition_linear_y, periodicity=false) + # 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 diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index 660e4dcf1be..79a650ded8a 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -137,7 +137,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") 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 + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 end end diff --git a/test/test_tree_2d_eulermulti.jl b/test/test_tree_2d_eulermulti.jl index c38c78ba6fd..c454a6bcfbf 100644 --- a/test/test_tree_2d_eulermulti.jl +++ b/test/test_tree_2d_eulermulti.jl @@ -44,7 +44,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") 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 + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 end end diff --git a/test/test_tree_2d_hypdiff.jl b/test/test_tree_2d_hypdiff.jl index eb8f8b297b6..30481fe910a 100644 --- a/test/test_tree_2d_hypdiff.jl +++ b/test/test_tree_2d_hypdiff.jl @@ -12,18 +12,42 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_lax_friedrichs.jl"), l2 = [0.00015687751817403066, 0.001025986772216324, 0.0010259867722164071], linf = [0.001198695637957381, 0.006423873515531753, 0.006423873515533529]) + # 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 @trixi_testset "elixir_hypdiff_harmonic_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_harmonic_nonperiodic.jl"), l2 = [8.618132355121019e-8, 5.619399844384306e-7, 5.619399844844044e-7], linf = [1.1248618588430072e-6, 8.622436487026874e-6, 8.622436487915053e-6]) + # 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_hypdiff_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"), l2 = [8.523077653954864e-6, 2.8779323653020624e-5, 5.454942769125663e-5], linf = [5.522740952468297e-5, 0.00014544895978971679, 0.00032396328684924924]) + # 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_hypdiff_godunov.jl" begin @@ -31,6 +55,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem") l2 = [5.868147556427088e-6, 3.80517927324465e-5, 3.805179273249344e-5], linf = [3.701965498725812e-5, 0.0002122422943138247, 0.00021224229431116015], atol = 2.0e-12 #= required for CI on macOS =#) + # 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 diff --git a/test/test_tree_2d_kpp.jl b/test/test_tree_2d_kpp.jl index 16d240d4a9f..fb2a212ddf8 100644 --- a/test/test_tree_2d_kpp.jl +++ b/test/test_tree_2d_kpp.jl @@ -17,17 +17,16 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") atol = 1e-6, rtol = 1e-6, skip_coverage = true) - - #= Does fail in coverage run (since it is not executed) - # 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 + if @isdefined sol # Skipped in coverage run + # 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 diff --git a/test/test_tree_3d_hypdiff.jl b/test/test_tree_3d_hypdiff.jl index 3db0ce93186..42231a9aaf6 100644 --- a/test/test_tree_3d_hypdiff.jl +++ b/test/test_tree_3d_hypdiff.jl @@ -13,6 +13,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem") l2 = [0.001530331609036682, 0.011314177033289238, 0.011314177033289402, 0.011314177033289631], linf = [0.02263459033909354, 0.10139777904683545, 0.10139777904683545, 0.10139777904683545], initial_refinement_level=2) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 + end end @trixi_testset "elixir_hypdiff_lax_friedrichs.jl with surface_flux=flux_godunov)" begin @@ -20,12 +28,28 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_3d_dgsem") l2 = [0.0015377731806850128, 0.01137685274151801, 0.011376852741518175, 0.011376852741518494], linf = [0.022715420630041172, 0.10183745338964201, 0.10183745338964201, 0.1018374533896429], initial_refinement_level=2, surface_flux=flux_godunov) + # 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 @trixi_testset "elixir_hypdiff_nonperiodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"), l2 = [0.00022868320512754316, 0.0007974309948540525, 0.0015035143230654987, 0.0015035143230655293], linf = [0.0016405001653623241, 0.0029870057159104594, 0.009410031618285686, 0.009410031618287462]) + # 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