From 3523c49120d7c282518769a5b3d40ce7c9cc5882 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 12 Sep 2023 15:00:58 +0200 Subject: [PATCH] AMR for 1D Parabolic Eqs (Clean branch) (#1605) * Clean branch * Un-Comment * un-comment * test coarsen * remove redundancy * Remove support for passive terms * expand resize * comments * format * Avoid code duplication * Update src/callbacks_step/amr_dg1d.jl Co-authored-by: Michael Schlottke-Lakemper * comment * comment & format * Try to increase coverage * Slightly more expressive names * Apply suggestions from code review --------- Co-authored-by: Michael Schlottke-Lakemper --- ...ixir_navierstokes_convergence_walls_amr.jl | 172 ++++++++++++++++++ src/callbacks_step/amr.jl | 158 ++++++++++++++++ src/callbacks_step/amr_dg1d.jl | 73 ++++++++ .../dgsem_tree/container_viscous_1d.jl | 58 ++++++ src/solvers/dgsem_tree/dg.jl | 3 + src/solvers/dgsem_tree/dg_1d_parabolic.jl | 14 +- test/test_parabolic_1d.jl | 41 +++++ test/test_parabolic_2d.jl | 12 +- test/test_parabolic_3d.jl | 12 +- 9 files changed, 523 insertions(+), 20 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl create mode 100644 src/solvers/dgsem_tree/container_viscous_1d.jl diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl new file mode 100644 index 00000000000..1daeab04a71 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl @@ -0,0 +1,172 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations1D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesEntropy()) + +# 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, + volume_integral=VolumeIntegralWeakForm()) + +coordinates_min = -1.0 +coordinates_max = 1.0 + +# 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 + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion1D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion1D`) +# and by the initial condition (which passes in `CompressibleEulerEquations1D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x[1] + pi_t = pi * t + + rho = c + A * cos(pi_x) * cos(pi_t) + v1 = log(x[1] + 2.0) * (1.0 - exp(-A * (x[1] - 1.0)) ) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + x = x[1] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5 + c = 2.0 + + # convenience values for trig. functions + pi_x = pi * x + pi_t = pi * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * cos(pi_x) * cos(pi_t) + rho_t = -pi * A * cos(pi_x) * sin(pi_t) + rho_x = -pi * A * sin(pi_x) * cos(pi_t) + rho_xx = -pi * pi * A * cos(pi_x) * cos(pi_t) + + v1 = log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * cos(pi_t) + v1_t = -pi * log(x + 2.0) * (1.0 - exp(-A * (x - 1.0))) * sin(pi_t) + v1_x = (A * log(x + 2.0) * exp(-A * (x - 1.0)) + (1.0 - exp(-A * (x - 1.0))) / (x + 2.0)) * cos(pi_t) + v1_xx = (( 2.0 * A * exp(-A * (x - 1.0)) / (x + 2.0) + - A * A * log(x + 2.0) * exp(-A * (x - 1.0)) + - (1.0 - exp(-A * (x - 1.0))) / ((x + 2.0) * (x + 2.0))) * cos(pi_t)) + + p = rho * rho + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 0.5 * rho_t * v1^2 + rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 0.5 * rho_x * v1^2 + rho * v1 * v1_x + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1.0 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + + # y-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + + p_x + rho_x * v1^2 + 2.0 * rho * v1 * v1_x + # stress tensor from y-direction + - v1_xx * mu_) + + # total energy equation + du3 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + # stress tensor and temperature gradient terms from x-direction + - v1_xx * v1 * mu_ + - v1_x * v1_x * mu_ + - T_const * inv_rho_cubed * ( p_xx * rho * rho + - 2.0 * p_x * rho * rho_x + + 2.0 * p * rho_x * rho_x + - p * rho * rho_xx ) * mu_ ) + + return SVector(du1, du2, du3) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2]) + +heat_bc_left = Isothermal((x, t, equations) -> + Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, t, equations), + equations_parabolic)) +heat_bc_right = Adiabatic((x, t, equations) -> 0.0) + +boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_left) +boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, heat_bc_right) + +# define inviscid boundary conditions +boundary_conditions = (; x_neg = boundary_condition_slip_wall, + x_pos = boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_left, + x_pos = boundary_condition_right) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +amr_controller = ControllerThreeLevel(semi, + IndicatorLöhner(semi, variable=Trixi.density), + base_level=3, + med_level=4, med_threshold=0.005, + max_level=5, max_threshold=0.01) + +amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true) + +# 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, amr_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary \ No newline at end of file diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 4d80e6e1139..ba840ff9675 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -192,6 +192,16 @@ end amr_callback(u_ode, mesh_equations_solver_cache(semi)..., semi, t, iter; kwargs...) end +@inline function (amr_callback::AMRCallback)(u_ode::AbstractVector, + semi::SemidiscretizationHyperbolicParabolic, + t, iter; + kwargs...) + # Note that we don't `wrap_array` the vector `u_ode` to be able to `resize!` + # it when doing AMR while still dispatching on the `mesh` etc. + amr_callback(u_ode, mesh_equations_solver_cache(semi)..., semi.cache_parabolic, + semi, t, iter; kwargs...) +end + # `passive_args` is currently used for Euler with self-gravity to adapt the gravity solver # passively without querying its indicator, based on the assumption that both solvers use # the same mesh. That's a hack and should be improved in the future once we have more examples @@ -346,6 +356,154 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, return has_changed end +function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, + equations, dg::DG, + cache, cache_parabolic, + semi::SemidiscretizationHyperbolicParabolic, + t, iter; + only_refine = false, only_coarsen = false) + @unpack controller, adaptor = amr_callback + + u = wrap_array(u_ode, mesh, equations, dg, cache) + # Indicator kept based on hyperbolic variables + lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache, + t = t, iter = iter) + + if mpi_isparallel() + error("MPI has not been verified yet for parabolic AMR") + + # Collect lambda for all elements + lambda_global = Vector{eltype(lambda)}(undef, nelementsglobal(dg, cache)) + # Use parent because n_elements_by_rank is an OffsetArray + recvbuf = MPI.VBuffer(lambda_global, parent(cache.mpi_cache.n_elements_by_rank)) + MPI.Allgatherv!(lambda, recvbuf, mpi_comm()) + lambda = lambda_global + end + + leaf_cell_ids = leaf_cells(mesh.tree) + @boundscheck begin + @assert axes(lambda)==axes(leaf_cell_ids) ("Indicator (axes = $(axes(lambda))) and leaf cell (axes = $(axes(leaf_cell_ids))) arrays have different axes") + end + + @unpack to_refine, to_coarsen = amr_callback.amr_cache + empty!(to_refine) + empty!(to_coarsen) + for element in 1:length(lambda) + controller_value = lambda[element] + if controller_value > 0 + push!(to_refine, leaf_cell_ids[element]) + elseif controller_value < 0 + push!(to_coarsen, leaf_cell_ids[element]) + end + end + + @trixi_timeit timer() "refine" if !only_coarsen && !isempty(to_refine) + # refine mesh + refined_original_cells = @trixi_timeit timer() "mesh" refine!(mesh.tree, + to_refine) + + # Find all indices of elements whose cell ids are in refined_original_cells + # Note: This assumes same indices for hyperbolic and parabolic part. + elements_to_refine = findall(in(refined_original_cells), + cache.elements.cell_ids) + + # refine solver + @trixi_timeit timer() "solver" refine!(u_ode, adaptor, mesh, equations, dg, + cache, cache_parabolic, + elements_to_refine) + else + # If there is nothing to refine, create empty array for later use + refined_original_cells = Int[] + end + + @trixi_timeit timer() "coarsen" if !only_refine && !isempty(to_coarsen) + # Since the cells may have been shifted due to refinement, first we need to + # translate the old cell ids to the new cell ids + if !isempty(to_coarsen) + to_coarsen = original2refined(to_coarsen, refined_original_cells, mesh) + end + + # Next, determine the parent cells from which the fine cells are to be + # removed, since these are needed for the coarsen! function. However, since + # we only want to coarsen if *all* child cells are marked for coarsening, + # we count the coarsening indicators for each parent cell and only coarsen + # if all children are marked as such (i.e., where the count is 2^ndims). At + # the same time, check if a cell is marked for coarsening even though it is + # *not* a leaf cell -> this can only happen if it was refined due to 2:1 + # smoothing during the preceding refinement operation. + parents_to_coarsen = zeros(Int, length(mesh.tree)) + for cell_id in to_coarsen + # If cell has no parent, it cannot be coarsened + if !has_parent(mesh.tree, cell_id) + continue + end + + # If cell is not leaf (anymore), it cannot be coarsened + if !is_leaf(mesh.tree, cell_id) + continue + end + + # Increase count for parent cell + parent_id = mesh.tree.parent_ids[cell_id] + parents_to_coarsen[parent_id] += 1 + end + + # Extract only those parent cells for which all children should be coarsened + to_coarsen = collect(1:length(parents_to_coarsen))[parents_to_coarsen .== 2^ndims(mesh)] + + # Finally, coarsen mesh + coarsened_original_cells = @trixi_timeit timer() "mesh" coarsen!(mesh.tree, + to_coarsen) + + # Convert coarsened parent cell ids to the list of child cell ids that have + # been removed, since this is the information that is expected by the solver + removed_child_cells = zeros(Int, + n_children_per_cell(mesh.tree) * + length(coarsened_original_cells)) + for (index, coarse_cell_id) in enumerate(coarsened_original_cells) + for child in 1:n_children_per_cell(mesh.tree) + removed_child_cells[n_children_per_cell(mesh.tree) * (index - 1) + child] = coarse_cell_id + + child + end + end + + # Find all indices of elements whose cell ids are in removed_child_cells + # Note: This assumes same indices for hyperbolic and parabolic part. + elements_to_remove = findall(in(removed_child_cells), cache.elements.cell_ids) + + # coarsen solver + @trixi_timeit timer() "solver" coarsen!(u_ode, adaptor, mesh, equations, dg, + cache, cache_parabolic, + elements_to_remove) + else + # If there is nothing to coarsen, create empty array for later use + coarsened_original_cells = Int[] + end + + # Store whether there were any cells coarsened or refined + has_changed = !isempty(refined_original_cells) || !isempty(coarsened_original_cells) + if has_changed # TODO: Taal decide, where shall we set this? + # don't set it to has_changed since there can be changes from earlier calls + mesh.unsaved_changes = true + end + + # Dynamically balance computational load by first repartitioning the mesh and then redistributing the cells/elements + if has_changed && mpi_isparallel() && amr_callback.dynamic_load_balancing + error("MPI has not been verified yet for parabolic AMR") + + @trixi_timeit timer() "dynamic load balancing" begin + old_mpi_ranks_per_cell = copy(mesh.tree.mpi_ranks) + + partition!(mesh) + + rebalance_solver!(u_ode, mesh, equations, dg, cache, old_mpi_ranks_per_cell) + end + end + + # Return true if there were any cells coarsened or refined, otherwise false + return has_changed +end + # Copy controller values to quad user data storage, will be called below function copy_to_quad_iter_volume(info, user_data) info_pw = PointerWrapper(info) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index e31a74730ea..e721ccc61cb 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -76,6 +76,44 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, return nothing end +function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, + 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) + + # 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 + 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) + + # Sanity check + 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 + + 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, 3}, element_id, @@ -201,6 +239,41 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, return nothing end +function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, + 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) + + # 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 + 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) + + # Sanity check + 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 + + return nothing +end + # TODO: Taal compare performance of different implementations # Coarsen solution data u for two elements, using L2 projection function coarsen_elements!(u::AbstractArray{<:Any, 3}, element_id, diff --git a/src/solvers/dgsem_tree/container_viscous_1d.jl b/src/solvers/dgsem_tree/container_viscous_1d.jl new file mode 100644 index 00000000000..a4919f75396 --- /dev/null +++ b/src/solvers/dgsem_tree/container_viscous_1d.jl @@ -0,0 +1,58 @@ +mutable struct ViscousContainer1D{uEltype <: Real} + u_transformed::Array{uEltype, 3} + gradients::Array{uEltype, 3} + flux_viscous::Array{uEltype, 3} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + _gradients::Vector{uEltype} + _flux_viscous::Vector{uEltype} + + function ViscousContainer1D{uEltype}(n_vars::Integer, n_nodes::Integer, + n_elements::Integer) where {uEltype <: Real} + new(Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), + Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), + Array{uEltype, 3}(undef, n_vars, n_nodes, n_elements), + Vector{uEltype}(undef, n_vars * n_nodes * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes * n_elements), + Vector{uEltype}(undef, n_vars * n_nodes * n_elements)) + end +end + +function init_viscous_container(n_vars::Integer, n_nodes::Integer, + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} + return ViscousContainer1D{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::ViscousContainer1D, equations, dg, cache) + capacity = nvariables(equations) * nnodes(dg) * nelements(dg, cache) + resize!(viscous_container._u_transformed, capacity) + resize!(viscous_container._gradients, capacity) + resize!(viscous_container._flux_viscous, capacity) + + viscous_container.u_transformed = unsafe_wrap(Array, + pointer(viscous_container._u_transformed), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + viscous_container.gradients = unsafe_wrap(Array, + pointer(viscous_container._gradients), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + viscous_container.flux_viscous = unsafe_wrap(Array, + pointer(viscous_container._flux_viscous), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + return nothing +end diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 6e02bc1d94a..ff37bad3b3a 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -54,6 +54,9 @@ include("containers.jl") # Dimension-agnostic parallel setup include("dg_parallel.jl") +# Helper struct for parabolic AMR +include("container_viscous_1d.jl") + # 1D DG implementation include("dg_1d.jl") include("dg_1d_parabolic.jl") diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 7602331d7c8..97e31e0e22b 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -17,7 +17,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1}, 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 @@ -534,18 +535,15 @@ function create_cache_parabolic(mesh::TreeMesh{1}, 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_elements) - gradients = similar(u_transformed) - flux_viscous = similar(u_transformed) + 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) - cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + cache = (; elements, interfaces, boundaries, viscous_container) return cache end diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 06a55100d62..3c2b8855ce8 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -20,6 +20,28 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "TreeMesh1D: elixir_advection_diffusion.jl (AMR)" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion.jl"), + tspan=(0.0, 0.0), initial_refinement_level = 5) + tspan=(0.0, 1.0) + ode = semidiscretize(semi, tspan) + amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), + base_level=4, + med_level=5, med_threshold=0.1, + max_level=6, max_threshold=0.6) + amr_callback = AMRCallback(semi, amr_controller, + interval=5, + adapt_initial_condition=true) + + # 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, amr_callback) + sol = solve(ode, KenCarp4(autodiff=false), abstol=time_abs_tol, reltol=time_int_tol, + save_everystep=false, callback=callbacks) + l2_error, linf_error = analysis_callback(sol) + @test l2_error ≈ [6.4878111416468355e-6] + @test linf_error ≈ [3.258075790424364e-5] + end + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_periodic.jl"), l2 = [0.0001133835907077494, 6.226282245610444e-5, 0.0002820171699999139], @@ -53,6 +75,25 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [0.002754803146635787, 0.0028567714697580906, 0.012941794048176192] ) end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls_amr.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls_amr.jl"), + equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number()), + l2 = [2.527877257772131e-5, 2.5539911566937718e-5, 0.0001211860451244785], + linf = [0.00014663867588948776, 0.00019422448348348196, 0.0009556439394007299] + ) + end + + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_walls_amr.jl: GradientVariablesEntropy" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_navierstokes_convergence_walls_amr.jl"), + equations_parabolic = CompressibleNavierStokesDiffusion1D(equations, mu=mu(), + Prandtl=prandtl_number(), + gradient_variables = GradientVariablesEntropy()), + l2 = [2.4593699163175966e-5, 2.392863645712634e-5, 0.00011252526651714956], + linf = [0.00011850555445525046, 0.0001898777490968537, 0.0009597561467877824] + ) + end end # Clean up afterwards: delete Trixi output directory diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 1564a33dc41..3fff4382cd1 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -143,9 +143,9 @@ isdir(outdir) && rm(outdir, recursive=true) callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, ode_default_options()..., callback=callbacks) - ac_sol = analysis_callback(sol) - @test ac_sol.l2[1] ≈ 1.67452550744728e-6 - @test ac_sol.linf[1] ≈ 7.905059166368744e-6 + l2_error, linf_error = analysis_callback(sol) + @test l2_error ≈ [1.67452550744728e-6] + @test linf_error ≈ [7.905059166368744e-6] # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -229,9 +229,9 @@ isdir(outdir) && rm(outdir, recursive=true) callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, ode_default_options()..., callback=callbacks) - ac_sol = analysis_callback(sol) - @test ac_sol.l2 ≈ [0.00024296959173852447; 0.0002093263158670915; 0.0005390572390977262; 0.00026753561392341537] - @test ac_sol.linf ≈ [0.0016210102053424436; 0.002593287648655501; 0.002953907343823712; 0.002077119120180271] + l2_error, linf_error = analysis_callback(sol) + @test l2_error ≈ [0.00024296959173852447; 0.0002093263158670915; 0.0005390572390977262; 0.00026753561392341537] + @test linf_error ≈ [0.0016210102053424436; 0.002593287648655501; 0.002953907343823712; 0.002077119120180271] end @trixi_testset "TreeMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index d607962afa0..ded052fb9d3 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -94,9 +94,9 @@ isdir(outdir) && rm(outdir, recursive=true) callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, ode_default_options()..., callback=callbacks) - ac_sol = analysis_callback(sol) - @test ac_sol.l2 ≈ [0.0003991794175622818; 0.0008853745163670504; 0.0010658655552066817; 0.0008785559918324284; 0.001403163458422815] - @test ac_sol.linf ≈ [0.0035306410538458177; 0.01505692306169911; 0.008862444161110705; 0.015065647972869856; 0.030402714743065218] + l2_error, linf_error = analysis_callback(sol) + @test l2_error ≈ [0.0003991794175622818; 0.0008853745163670504; 0.0010658655552066817; 0.0008785559918324284; 0.001403163458422815] + @test linf_error ≈ [0.0035306410538458177; 0.01505692306169911; 0.008862444161110705; 0.015065647972869856; 0.030402714743065218] end @trixi_testset "TreeMesh3D: elixir_navierstokes_taylor_green_vortex.jl" begin @@ -127,9 +127,9 @@ isdir(outdir) && rm(outdir, recursive=true) sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), dt=5e-3, save_everystep=false, callback=callbacks); - ac_sol = analysis_callback(sol) - @test ac_sol.l2 ≈ [0.0013666103707729502; 0.2313581629543744; 0.2308164306264533; 0.17460246787819503; 0.28121914446544005] - @test ac_sol.linf ≈ [0.006938093883741336; 1.028235074139312; 1.0345438209717241; 1.0821111605203542; 1.2669636522564645] + l2_error, linf_error = analysis_callback(sol) + @test l2_error ≈ [0.0013666103707729502; 0.2313581629543744; 0.2308164306264533; 0.17460246787819503; 0.28121914446544005] + @test linf_error ≈ [0.006938093883741336; 1.028235074139312; 1.0345438209717241; 1.0821111605203542; 1.2669636522564645] # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)