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/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/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