From b33cbc2c0e779209dc998cc8b08b52bec86022b9 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 14 Aug 2023 08:44:39 +0200 Subject: [PATCH 01/53] Clean branch --- ...ixir_navierstokes_convergence_walls_amr.jl | 172 ++++++++++++++ src/callbacks_step/amr.jl | 172 ++++++++++++++ src/callbacks_step/amr_dg1d.jl | 222 ++++++++++++++++++ src/solvers/dgsem_tree/cache_viscous.jl | 33 +++ src/solvers/dgsem_tree/dg.jl | 3 + src/solvers/dgsem_tree/dg_1d_parabolic.jl | 11 +- test/test_parabolic_1d.jl | 19 ++ 7 files changed, 627 insertions(+), 5 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl create mode 100644 src/solvers/dgsem_tree/cache_viscous.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..d742b17d262 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,168 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, return has_changed 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 +# and a better understanding of such a coupling. +# `passive_args` is expected to be an iterable of `Tuple`s of the form +# `(p_u_ode, p_mesh, p_equations, p_dg, p_cache)`. +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, + passive_args = ()) + @unpack controller, adaptor = amr_callback + + u = wrap_array(u_ode, mesh, equations, dg, cache) + # TODO: Keep indicator based on hyperbolic variables? + lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache, + t = t, iter = iter) + + if mpi_isparallel() + # 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) + for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args + @trixi_timeit timer() "passive solver" refine!(p_u_ode, adaptor, p_mesh, + p_equations, p_dg, p_cache, + elements_to_refine) + end + 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) + + for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args + @trixi_timeit timer() "passive solver" coarsen!(p_u_ode, adaptor, p_mesh, + p_equations, p_dg, p_cache, + elements_to_remove) + end + 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 + @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..e729c8eccc8 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -76,6 +76,114 @@ 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) + # Return early if there is nothing to do + if isempty(elements_to_refine) + return + end + + # Determine for each existing element whether it needs to be refined + needs_refinement = falses(nelements(dg, cache)) + needs_refinement[elements_to_refine] .= true + + # Retain current solution data + old_n_elements = nelements(dg, cache) + old_u_ode = copy(u_ode) + GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Get new list of leaf cells + leaf_cell_ids = local_leaf_cells(mesh.tree) + + # re-initialize elements container + @unpack elements = cache + resize!(elements, length(leaf_cell_ids)) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) + @assert nelements(dg, cache) > old_n_elements + + @unpack elements, cache_viscous = cache_parabolic + resize!(elements, length(leaf_cell_ids)) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) + @assert nelements(dg, cache_parabolic) > old_n_elements + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + # Resize parabolic helper variables + resize!(cache_viscous, + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) + cache_parabolic.cache_viscous.u_transformed = unsafe_wrap(Array, + pointer(cache_parabolic.cache_viscous._u_transformed), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + cache_parabolic.cache_viscous.gradients = unsafe_wrap(Array, + pointer(cache_parabolic.cache_viscous._gradients), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + cache_parabolic.cache_viscous.flux_viscous = unsafe_wrap(Array, + pointer(cache_parabolic.cache_viscous._flux_viscous), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + # Loop over all elements in old container and either copy them or refine them + element_id = 1 + for old_element_id in 1:old_n_elements + if needs_refinement[old_element_id] + # Refine element and store solution directly in new data structure + refine_element!(u, element_id, old_u, old_element_id, + adaptor, equations, dg) + element_id += 2^ndims(mesh) + else + # Copy old element data to new element container + @views u[:, .., element_id] .= old_u[:, .., old_element_id] + element_id += 1 + end + end + # If everything is correct, we should have processed all elements. + # Depending on whether the last element processed above had to be refined or not, + # the counter `element_id` can have two different values at the end. + @assert element_id == + nelements(dg, cache) + + 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" + #= + @assert element_id == + nelements(dg, cache_parabolic) + + 1||element_id == nelements(dg, cache_parabolic) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache_parabolic) = $(nelements(dg, cache_parabolic))" + =# + end # GC.@preserve old_u_ode + + # re-initialize interfaces container + @unpack interfaces = cache + resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids)) + init_interfaces!(interfaces, elements, mesh) + + @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 + resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids)) + init_boundaries!(boundaries, elements, mesh) + + @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) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") + #@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 +309,120 @@ 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) + # Return early if there is nothing to do + if isempty(elements_to_remove) + return + end + + # Determine for each old element whether it needs to be removed + to_be_removed = falses(nelements(dg, cache)) + to_be_removed[elements_to_remove] .= true + + # Retain current solution data + old_n_elements = nelements(dg, cache) + old_u_ode = copy(u_ode) + GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Get new list of leaf cells + leaf_cell_ids = local_leaf_cells(mesh.tree) + + # re-initialize elements container + @unpack elements = cache + resize!(elements, length(leaf_cell_ids)) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) + @assert nelements(dg, cache) < old_n_elements + + @unpack elements, cache_viscous = cache_parabolic + resize!(elements, length(leaf_cell_ids)) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) + @assert nelements(dg, cache_parabolic) < old_n_elements + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + # Resize parabolic helper variables + resize!(cache_viscous, + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) + cache_parabolic.cache_viscous.u_transformed = unsafe_wrap(Array, + pointer(cache_parabolic.cache_viscous._u_transformed), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + cache_parabolic.cache_viscous.gradients = unsafe_wrap(Array, + pointer(cache_parabolic.cache_viscous._gradients), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + cache_parabolic.cache_viscous.flux_viscous = unsafe_wrap(Array, + pointer(cache_parabolic.cache_viscous._flux_viscous), + (nvariables(equations), + nnodes(dg), + nelements(dg, cache))) + + # Loop over all elements in old container and either copy them or coarsen them + skip = 0 + element_id = 1 + for old_element_id in 1:old_n_elements + # If skip is non-zero, we just coarsened 2^ndims elements and need to omit the following elements + if skip > 0 + skip -= 1 + continue + end + + if to_be_removed[old_element_id] + # If an element is to be removed, sanity check if the following elements + # are also marked - otherwise there would be an error in the way the + # cells/elements are sorted + @assert all(to_be_removed[old_element_id:(old_element_id + 2^ndims(mesh) - 1)]) "bad cell/element order" + + # Coarsen elements and store solution directly in new data structure + coarsen_elements!(u, element_id, old_u, old_element_id, + adaptor, equations, dg) + element_id += 1 + skip = 2^ndims(mesh) - 1 + else + # Copy old element data to new element container + @views u[:, .., element_id] .= old_u[:, .., old_element_id] + element_id += 1 + end + end + # If everything is correct, we should have processed all elements. + @assert element_id==nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" + @assert element_id==nelements(dg, cache_parabolic) + 1 "element_id = $element_id, nelements(dg, cache_parabolic) = $(nelements(dg, cache_parabolic))" + end # GC.@preserve old_u_ode + + # re-initialize interfaces container + @unpack interfaces = cache + resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids)) + init_interfaces!(interfaces, elements, mesh) + + @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 + resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids)) + init_boundaries!(boundaries, elements, mesh) + + @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) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") + @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/cache_viscous.jl b/src/solvers/dgsem_tree/cache_viscous.jl new file mode 100644 index 00000000000..17a2774db06 --- /dev/null +++ b/src/solvers/dgsem_tree/cache_viscous.jl @@ -0,0 +1,33 @@ +mutable struct CacheViscous1D{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 CacheViscous1D{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 + +# 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!(cache_viscous::CacheViscous1D, capacity) + resize!(cache_viscous._u_transformed, capacity) + resize!(cache_viscous._gradients, capacity) + resize!(cache_viscous._flux_viscous, capacity) + + return nothing +end diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index cb28dad968c..9b0874e5ae0 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("cache_viscous.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 c2aa75388c8..596a6a8acbf 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 cache_viscous = cache_parabolic + @unpack u_transformed, gradients, flux_viscous = cache_viscous # Convert conservative variables to a form more suitable for viscous flux calculations @trixi_timeit timer() "transform variables" begin @@ -530,18 +531,18 @@ function create_cache_parabolic(mesh::TreeMesh{1}, elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) + # Additions for parabolic 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) + + cache_viscous = CacheViscous1D{uEltype}(n_vars, n_nodes, n_elements) 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, cache_viscous) return cache end diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 06a55100d62..ee8c9ea9d5e 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -53,6 +53,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 From 7e32e24cedb655f7f13b126be28f10fa06d790e4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 14 Aug 2023 08:58:32 +0200 Subject: [PATCH 02/53] Un-Comment --- src/callbacks_step/amr_dg1d.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index e729c8eccc8..a04eef031f9 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -150,11 +150,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, @assert element_id == nelements(dg, cache) + 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - #= @assert element_id == nelements(dg, cache_parabolic) + 1||element_id == nelements(dg, cache_parabolic) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache_parabolic) = $(nelements(dg, cache_parabolic))" - =# end # GC.@preserve old_u_ode # re-initialize interfaces container From 5e1997b5b252c207b5e34fd5cceb37c00203bfae Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 14 Aug 2023 09:00:59 +0200 Subject: [PATCH 03/53] un-comment --- src/callbacks_step/amr_dg1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index a04eef031f9..ed1341ce3c3 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -176,7 +176,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # Sanity check if isperiodic(mesh.tree) @assert ninterfaces(interfaces)==1 * nelements(dg, cache) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") - #@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") + @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 From 54e328e378b301aa64db897f2cc2f2babdb1dd35 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 14 Aug 2023 09:20:27 +0200 Subject: [PATCH 04/53] test coarsen --- test/test_parabolic_1d.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index ee8c9ea9d5e..45210c6795e 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -20,6 +20,29 @@ 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 = 6) + tspan=(0.0, 1.0) + ode = semidiscretize(semi, tspan) + initial_refinement_level = 6 + 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) + ac_sol = analysis_callback(sol) + @test ac_sol.l2[1] ≈ 5.888985009128617e-6 + @test ac_sol.linf[1] ≈ 2.885052015366707e-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], From a9e4cb796816509923efcbe42e9569251c2e14e3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 14 Aug 2023 14:27:50 +0200 Subject: [PATCH 05/53] remove redundancy --- test/test_parabolic_1d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 45210c6795e..ccf7b0f5c7c 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -25,7 +25,6 @@ isdir(outdir) && rm(outdir, recursive=true) tspan=(0.0, 0.0), initial_refinement_level = 6) tspan=(0.0, 1.0) ode = semidiscretize(semi, tspan) - initial_refinement_level = 6 amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), base_level=4, med_level=5, med_threshold=0.1, From ebc9cc8be602140867202a32dbf5b09c86e37561 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 18 Aug 2023 09:51:55 +0200 Subject: [PATCH 06/53] Remove support for passive terms --- src/callbacks_step/amr.jl | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index d742b17d262..9ea76157ae1 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -356,19 +356,12 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, return has_changed 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 -# and a better understanding of such a coupling. -# `passive_args` is expected to be an iterable of `Tuple`s of the form -# `(p_u_ode, p_mesh, p_equations, p_dg, p_cache)`. 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, - passive_args = ()) + only_refine = false, only_coarsen = false) @unpack controller, adaptor = amr_callback u = wrap_array(u_ode, mesh, equations, dg, cache) @@ -416,11 +409,6 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @trixi_timeit timer() "solver" refine!(u_ode, adaptor, mesh, equations, dg, cache, cache_parabolic, elements_to_refine) - for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args - @trixi_timeit timer() "passive solver" refine!(p_u_ode, adaptor, p_mesh, - p_equations, p_dg, p_cache, - elements_to_refine) - end else # If there is nothing to refine, create empty array for later use refined_original_cells = Int[] @@ -485,12 +473,6 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @trixi_timeit timer() "solver" coarsen!(u_ode, adaptor, mesh, equations, dg, cache, cache_parabolic, elements_to_remove) - - for (p_u_ode, p_mesh, p_equations, p_dg, p_cache) in passive_args - @trixi_timeit timer() "passive solver" coarsen!(p_u_ode, adaptor, p_mesh, - p_equations, p_dg, p_cache, - elements_to_remove) - end else # If there is nothing to coarsen, create empty array for later use coarsened_original_cells = Int[] From 9c766b93306680af3a31b4fc2866f7c93fa6d1b4 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 18 Aug 2023 10:28:52 +0200 Subject: [PATCH 07/53] expand resize --- src/callbacks_step/amr_dg1d.jl | 41 ++----------- src/solvers/dgsem_tree/cache_viscous.jl | 33 ----------- .../dgsem_tree/container_viscous_1d.jl | 58 +++++++++++++++++++ src/solvers/dgsem_tree/dg.jl | 2 +- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 15 ++--- 5 files changed, 69 insertions(+), 80 deletions(-) delete mode 100644 src/solvers/dgsem_tree/cache_viscous.jl create mode 100644 src/solvers/dgsem_tree/container_viscous_1d.jl diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index ed1341ce3c3..dbbdf1ca397 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -103,7 +103,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache) > old_n_elements - @unpack elements, cache_viscous = cache_parabolic + @unpack elements, viscous_container = cache_parabolic resize!(elements, length(leaf_cell_ids)) init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache_parabolic) > old_n_elements @@ -113,23 +113,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, u = wrap_array(u_ode, mesh, equations, dg, cache) # Resize parabolic helper variables - resize!(cache_viscous, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - cache_parabolic.cache_viscous.u_transformed = unsafe_wrap(Array, - pointer(cache_parabolic.cache_viscous._u_transformed), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - cache_parabolic.cache_viscous.gradients = unsafe_wrap(Array, - pointer(cache_parabolic.cache_viscous._gradients), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - cache_parabolic.cache_viscous.flux_viscous = unsafe_wrap(Array, - pointer(cache_parabolic.cache_viscous._flux_viscous), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) + resize!(viscous_container, equations, dg, cache) + # Loop over all elements in old container and either copy them or refine them element_id = 1 for old_element_id in 1:old_n_elements @@ -334,7 +319,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache) < old_n_elements - @unpack elements, cache_viscous = cache_parabolic + @unpack elements, viscous_container = cache_parabolic resize!(elements, length(leaf_cell_ids)) init_elements!(elements, leaf_cell_ids, mesh, dg.basis) @assert nelements(dg, cache_parabolic) < old_n_elements @@ -344,23 +329,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, u = wrap_array(u_ode, mesh, equations, dg, cache) # Resize parabolic helper variables - resize!(cache_viscous, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - cache_parabolic.cache_viscous.u_transformed = unsafe_wrap(Array, - pointer(cache_parabolic.cache_viscous._u_transformed), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - cache_parabolic.cache_viscous.gradients = unsafe_wrap(Array, - pointer(cache_parabolic.cache_viscous._gradients), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) - cache_parabolic.cache_viscous.flux_viscous = unsafe_wrap(Array, - pointer(cache_parabolic.cache_viscous._flux_viscous), - (nvariables(equations), - nnodes(dg), - nelements(dg, cache))) + resize!(viscous_container, equations, dg, cache) # Loop over all elements in old container and either copy them or coarsen them skip = 0 diff --git a/src/solvers/dgsem_tree/cache_viscous.jl b/src/solvers/dgsem_tree/cache_viscous.jl deleted file mode 100644 index 17a2774db06..00000000000 --- a/src/solvers/dgsem_tree/cache_viscous.jl +++ /dev/null @@ -1,33 +0,0 @@ -mutable struct CacheViscous1D{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 CacheViscous1D{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 - -# 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!(cache_viscous::CacheViscous1D, capacity) - resize!(cache_viscous._u_transformed, capacity) - resize!(cache_viscous._gradients, capacity) - resize!(cache_viscous._flux_viscous, capacity) - - return nothing -end 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..47de34473a8 --- /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 9b0874e5ae0..33132fee774 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -55,7 +55,7 @@ include("containers.jl") include("dg_parallel.jl") # Helper struct for parabolic AMR -include("cache_viscous.jl") +include("container_viscous_1d.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 596a6a8acbf..64620edfb05 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -17,8 +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 cache_viscous = cache_parabolic - @unpack u_transformed, gradients, flux_viscous = cache_viscous + @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 @@ -530,19 +530,14 @@ function create_cache_parabolic(mesh::TreeMesh{1}, elements = init_elements(leaf_cell_ids, mesh, equations_hyperbolic, dg.basis, RealT, uEltype) - - # Additions for parabolic - n_vars = nvariables(equations_hyperbolic) - n_nodes = nnodes(elements) - n_elements = nelements(elements) - - cache_viscous = CacheViscous1D{uEltype}(n_vars, n_nodes, n_elements) + + 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, cache_viscous) + cache = (; elements, interfaces, boundaries, viscous_container) return cache end From 82894d71086f98abd389166892d0fde3512b1a13 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 18 Aug 2023 10:35:43 +0200 Subject: [PATCH 08/53] comments --- src/callbacks_step/amr.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 9ea76157ae1..0ea2500f6aa 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -370,6 +370,8 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, 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 @@ -401,7 +403,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, 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! + # Note: This assumes same indices for hyperbolic and parabolic part. elements_to_refine = findall(in(refined_original_cells), cache.elements.cell_ids) @@ -466,7 +468,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, end # Find all indices of elements whose cell ids are in removed_child_cells - # NOTE: This assumes same indices for hyperbolic and parabolic part! + # Note: This assumes same indices for hyperbolic and parabolic part. elements_to_remove = findall(in(removed_child_cells), cache.elements.cell_ids) # coarsen solver @@ -487,6 +489,8 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, # 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) From 6ef88ca6343fb173f939660011ff1e1a896124ce Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Fri, 18 Aug 2023 10:37:38 +0200 Subject: [PATCH 09/53] format --- src/solvers/dgsem_tree/container_viscous_1d.jl | 4 ++-- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 6 ++++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_1d.jl b/src/solvers/dgsem_tree/container_viscous_1d.jl index 47de34473a8..a4919f75396 100644 --- a/src/solvers/dgsem_tree/container_viscous_1d.jl +++ b/src/solvers/dgsem_tree/container_viscous_1d.jl @@ -20,7 +20,8 @@ mutable struct ViscousContainer1D{uEltype <: Real} end function init_viscous_container(n_vars::Integer, n_nodes::Integer, - n_elements::Integer, ::Type{uEltype}) where {uEltype <: Real} + n_elements::Integer, + ::Type{uEltype}) where {uEltype <: Real} return ViscousContainer1D{uEltype}(n_vars, n_nodes, n_elements) end @@ -30,7 +31,6 @@ end # `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) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 64620edfb05..561bd714032 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -530,8 +530,10 @@ 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) + + viscous_container = init_viscous_container(nvariables(equations_hyperbolic), + nnodes(elements), nelements(elements), + uEltype) interfaces = init_interfaces(leaf_cell_ids, mesh, elements) From 53f599145e001384c952843b0e540e1721f4a3b1 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 20 Aug 2023 14:45:48 +0200 Subject: [PATCH 10/53] Avoid code duplication --- src/callbacks_step/amr_dg1d.jl | 155 ++++----------------------------- 1 file changed, 16 insertions(+), 139 deletions(-) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index dbbdf1ca397..d5f300ecb04 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -79,88 +79,30 @@ end function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_refine) - # Return early if there is nothing to do - if isempty(elements_to_refine) - return - end - - # Determine for each existing element whether it needs to be refined - needs_refinement = falses(nelements(dg, cache)) - needs_refinement[elements_to_refine] .= true - - # Retain current solution data - old_n_elements = nelements(dg, cache) - old_u_ode = copy(u_ode) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - - # Get new list of leaf cells - leaf_cell_ids = local_leaf_cells(mesh.tree) - - # re-initialize elements container - @unpack elements = cache - resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - @assert nelements(dg, cache) > old_n_elements - - @unpack elements, viscous_container = cache_parabolic - resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - @assert nelements(dg, cache_parabolic) > old_n_elements + refine!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_refine) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + # Get new list of leaf cells + leaf_cell_ids = local_leaf_cells(mesh.tree) - # Resize parabolic helper variables - resize!(viscous_container, equations, dg, cache) + @unpack elements, viscous_container = cache_parabolic + resize!(elements, length(leaf_cell_ids)) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - # Loop over all elements in old container and either copy them or refine them - element_id = 1 - for old_element_id in 1:old_n_elements - if needs_refinement[old_element_id] - # Refine element and store solution directly in new data structure - refine_element!(u, element_id, old_u, old_element_id, - adaptor, equations, dg) - element_id += 2^ndims(mesh) - else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] - element_id += 1 - end - end - # If everything is correct, we should have processed all elements. - # Depending on whether the last element processed above had to be refined or not, - # the counter `element_id` can have two different values at the end. - @assert element_id == - nelements(dg, cache) + - 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - @assert element_id == - nelements(dg, cache_parabolic) + - 1||element_id == nelements(dg, cache_parabolic) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache_parabolic) = $(nelements(dg, cache_parabolic))" - end # GC.@preserve old_u_ode + # Resize parabolic helper variables + resize!(viscous_container, equations, dg, cache) # re-initialize interfaces container - @unpack interfaces = cache - resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids)) - init_interfaces!(interfaces, elements, mesh) - @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 - resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids)) - init_boundaries!(boundaries, elements, mesh) - @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) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") @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 @@ -295,95 +237,30 @@ end function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, equations, dg::DGSEM, cache, cache_parabolic, elements_to_remove) - # Return early if there is nothing to do - if isempty(elements_to_remove) - return - end + coarsen!(u_ode, adaptor, mesh, equations, dg, cache, elements_to_remove) - # Determine for each old element whether it needs to be removed - to_be_removed = falses(nelements(dg, cache)) - to_be_removed[elements_to_remove] .= true - - # Retain current solution data - old_n_elements = nelements(dg, cache) - old_u_ode = copy(u_ode) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - - # Get new list of leaf cells - leaf_cell_ids = local_leaf_cells(mesh.tree) - - # re-initialize elements container - @unpack elements = cache - resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - @assert nelements(dg, cache) < old_n_elements - - @unpack elements, viscous_container = cache_parabolic - resize!(elements, length(leaf_cell_ids)) - init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - @assert nelements(dg, cache_parabolic) < old_n_elements - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - # Resize parabolic helper variables - resize!(viscous_container, equations, dg, cache) - - # Loop over all elements in old container and either copy them or coarsen them - skip = 0 - element_id = 1 - for old_element_id in 1:old_n_elements - # If skip is non-zero, we just coarsened 2^ndims elements and need to omit the following elements - if skip > 0 - skip -= 1 - continue - end + # Get new list of leaf cells + leaf_cell_ids = local_leaf_cells(mesh.tree) - if to_be_removed[old_element_id] - # If an element is to be removed, sanity check if the following elements - # are also marked - otherwise there would be an error in the way the - # cells/elements are sorted - @assert all(to_be_removed[old_element_id:(old_element_id + 2^ndims(mesh) - 1)]) "bad cell/element order" + @unpack elements, viscous_container = cache_parabolic + resize!(elements, length(leaf_cell_ids)) + init_elements!(elements, leaf_cell_ids, mesh, dg.basis) - # Coarsen elements and store solution directly in new data structure - coarsen_elements!(u, element_id, old_u, old_element_id, - adaptor, equations, dg) - element_id += 1 - skip = 2^ndims(mesh) - 1 - else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] - element_id += 1 - end - end - # If everything is correct, we should have processed all elements. - @assert element_id==nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))" - @assert element_id==nelements(dg, cache_parabolic) + 1 "element_id = $element_id, nelements(dg, cache_parabolic) = $(nelements(dg, cache_parabolic))" - end # GC.@preserve old_u_ode + # Resize parabolic helper variables + resize!(viscous_container, equations, dg, cache) # re-initialize interfaces container - @unpack interfaces = cache - resize!(interfaces, count_required_interfaces(mesh, leaf_cell_ids)) - init_interfaces!(interfaces, elements, mesh) - @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 - resize!(boundaries, count_required_boundaries(mesh, leaf_cell_ids)) - init_boundaries!(boundaries, elements, mesh) - @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) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements") @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 From cdfe93b092c28e363e9dc8bbaf5c64a938728f8f Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Tue, 22 Aug 2023 13:16:22 +0200 Subject: [PATCH 11/53] Update src/callbacks_step/amr_dg1d.jl Co-authored-by: Michael Schlottke-Lakemper --- src/callbacks_step/amr_dg1d.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index d5f300ecb04..1b273b05ef5 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -79,7 +79,12 @@ 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) From 0f2b7790176c648cff5b3014b019b29f06ce644f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 22 Aug 2023 13:22:52 +0200 Subject: [PATCH 12/53] comment --- src/callbacks_step/amr_dg1d.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 1b273b05ef5..602d7645009 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -242,6 +242,8 @@ 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 From 376f99e39e20670397a302a4bed09cc2b3c08aa6 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 22 Aug 2023 14:28:02 +0200 Subject: [PATCH 13/53] comment & format --- src/callbacks_step/amr.jl | 2 +- src/callbacks_step/amr_dg1d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl index 0ea2500f6aa..ba840ff9675 100644 --- a/src/callbacks_step/amr.jl +++ b/src/callbacks_step/amr.jl @@ -365,7 +365,7 @@ function (amr_callback::AMRCallback)(u_ode::AbstractVector, mesh::TreeMesh, @unpack controller, adaptor = amr_callback u = wrap_array(u_ode, mesh, equations, dg, cache) - # TODO: Keep indicator based on hyperbolic variables? + # Indicator kept based on hyperbolic variables lambda = @trixi_timeit timer() "indicator" controller(u, mesh, equations, dg, cache, t = t, iter = iter) diff --git a/src/callbacks_step/amr_dg1d.jl b/src/callbacks_step/amr_dg1d.jl index 602d7645009..e721ccc61cb 100644 --- a/src/callbacks_step/amr_dg1d.jl +++ b/src/callbacks_step/amr_dg1d.jl @@ -82,7 +82,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1}, # 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 From 826553f118bb6843de913b6ec5c4e694f72614d3 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 28 Aug 2023 14:35:39 +0200 Subject: [PATCH 14/53] Try to increase coverage --- test/test_parabolic_1d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index ccf7b0f5c7c..3c4b34d23cb 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -22,7 +22,7 @@ isdir(outdir) && rm(outdir, recursive=true) @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 = 6) + 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), @@ -38,8 +38,8 @@ isdir(outdir) && rm(outdir, recursive=true) sol = solve(ode, KenCarp4(autodiff=false), abstol=time_abs_tol, reltol=time_int_tol, save_everystep=false, callback=callbacks) ac_sol = analysis_callback(sol) - @test ac_sol.l2[1] ≈ 5.888985009128617e-6 - @test ac_sol.linf[1] ≈ 2.885052015366707e-5 + @test ac_sol.l2[1] ≈ 6.4878111416468355e-6 + @test ac_sol.linf[1] ≈ 3.258075790424364e-5 end @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl" begin From d20962904a0308c8e7085b7d2a1083ea62cdc228 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 29 Aug 2023 11:50:58 +0200 Subject: [PATCH 15/53] Slightly more expressive names --- test/test_parabolic_1d.jl | 6 +++--- test/test_parabolic_2d.jl | 12 ++++++------ test/test_parabolic_3d.jl | 12 ++++++------ 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 3c4b34d23cb..7962151c1f3 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -37,9 +37,9 @@ isdir(outdir) && rm(outdir, recursive=true) 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) - ac_sol = analysis_callback(sol) - @test ac_sol.l2[1] ≈ 6.4878111416468355e-6 - @test ac_sol.linf[1] ≈ 3.258075790424364e-5 + analysis_callback_sol = analysis_callback(sol) + @test analysis_callback_sol.l2[1] ≈ 6.4878111416468355e-6 + @test analysis_callback_sol.linf[1] ≈ 3.258075790424364e-5 end @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl" begin diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 1564a33dc41..ca0b8b198a2 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 + analysis_callback_sol = analysis_callback(sol) + @test analysis_callback_sol.l2[1] ≈ 1.67452550744728e-6 + @test analysis_callback_sol.linf[1] ≈ 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] + analysis_callback_sol = analysis_callback(sol) + @test analysis_callback_sol.l2 ≈ [0.00024296959173852447; 0.0002093263158670915; 0.0005390572390977262; 0.00026753561392341537] + @test analysis_callback_sol.linf ≈ [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..14eb9471b41 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] + analysis_callback_sol = analysis_callback(sol) + @test analysis_callback_sol.l2 ≈ [0.0003991794175622818; 0.0008853745163670504; 0.0010658655552066817; 0.0008785559918324284; 0.001403163458422815] + @test analysis_callback_sol.linf ≈ [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] + analysis_callback_sol = analysis_callback(sol) + @test analysis_callback_sol.l2 ≈ [0.0013666103707729502; 0.2313581629543744; 0.2308164306264533; 0.17460246787819503; 0.28121914446544005] + @test analysis_callback_sol.linf ≈ [0.006938093883741336; 1.028235074139312; 1.0345438209717241; 1.0821111605203542; 1.2669636522564645] # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From a7d56a1ae56a020420bd1bfb46207d215b7732f2 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 1 Sep 2023 21:55:38 +0200 Subject: [PATCH 16/53] Apply suggestions from code review --- test/test_parabolic_1d.jl | 6 +++--- test/test_parabolic_2d.jl | 12 ++++++------ test/test_parabolic_3d.jl | 12 ++++++------ 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 7962151c1f3..3c2b8855ce8 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -37,9 +37,9 @@ isdir(outdir) && rm(outdir, recursive=true) 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) - analysis_callback_sol = analysis_callback(sol) - @test analysis_callback_sol.l2[1] ≈ 6.4878111416468355e-6 - @test analysis_callback_sol.linf[1] ≈ 3.258075790424364e-5 + 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 diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index ca0b8b198a2..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) - analysis_callback_sol = analysis_callback(sol) - @test analysis_callback_sol.l2[1] ≈ 1.67452550744728e-6 - @test analysis_callback_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) - analysis_callback_sol = analysis_callback(sol) - @test analysis_callback_sol.l2 ≈ [0.00024296959173852447; 0.0002093263158670915; 0.0005390572390977262; 0.00026753561392341537] - @test analysis_callback_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 14eb9471b41..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) - analysis_callback_sol = analysis_callback(sol) - @test analysis_callback_sol.l2 ≈ [0.0003991794175622818; 0.0008853745163670504; 0.0010658655552066817; 0.0008785559918324284; 0.001403163458422815] - @test analysis_callback_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); - analysis_callback_sol = analysis_callback(sol) - @test analysis_callback_sol.l2 ≈ [0.0013666103707729502; 0.2313581629543744; 0.2308164306264533; 0.17460246787819503; 0.28121914446544005] - @test analysis_callback_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) From abb6dfb81392dc5d4a35cb09b6cb83a4a163c67f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 09:55:25 +0200 Subject: [PATCH 17/53] add specifier for 1d --- src/solvers/dgsem_tree/container_viscous_1d.jl | 2 +- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_1d.jl b/src/solvers/dgsem_tree/container_viscous_1d.jl index a4919f75396..3f9943aea1c 100644 --- a/src/solvers/dgsem_tree/container_viscous_1d.jl +++ b/src/solvers/dgsem_tree/container_viscous_1d.jl @@ -19,7 +19,7 @@ mutable struct ViscousContainer1D{uEltype <: Real} end end -function init_viscous_container(n_vars::Integer, n_nodes::Integer, +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) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 561bd714032..20626471941 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -531,7 +531,7 @@ 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), + viscous_container = init_viscous_container_1d(nvariables(equations_hyperbolic), nnodes(elements), nelements(elements), uEltype) From bfa3a2448e1c9045e478da7f2dc1cb3b08353576 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 09:56:22 +0200 Subject: [PATCH 18/53] Structs for resizing parabolic helpers --- .../dgsem_tree/container_viscous_2d.jl | 60 ++++++++++++++++++ .../dgsem_tree/container_viscous_3d.jl | 61 +++++++++++++++++++ src/solvers/dgsem_tree/containers_viscous.jl | 4 ++ src/solvers/dgsem_tree/dg.jl | 4 +- 4 files changed, 127 insertions(+), 2 deletions(-) create mode 100644 src/solvers/dgsem_tree/container_viscous_2d.jl create mode 100644 src/solvers/dgsem_tree/container_viscous_3d.jl create mode 100644 src/solvers/dgsem_tree/containers_viscous.jl 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..151011e9757 --- /dev/null +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -0,0 +1,60 @@ +mutable struct ViscousContainer2D{uEltype <: Real} + u_transformed::Array{uEltype, 4} + # IDEA: Use SVector for fixed sized vectors? + gradients::Vector{Array{uEltype, 4}} + flux_viscous::Vector{Array{uEltype, 4}} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + _gradients::Vector{Vector{uEltype}} + _flux_viscous::Vector{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) for _ in 1:2], + [Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements) for _ in 1:2]) + 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!(cache_viscous._u_transformed, capacity) + for dim in 1:2 + resize!(cache_viscous._gradients[dim], capacity) + resize!(cache_viscous._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))) + + viscous_container.gradients = unsafe_wrap(Array, + pointer(viscous_container._gradients), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + viscous_container.flux_viscous = unsafe_wrap(Array, + pointer(viscous_container._flux_viscous), + (nvariables(equations), + nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + return nothing +end \ No newline at end of file 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..9a7749d7a9b --- /dev/null +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -0,0 +1,61 @@ +mutable struct ViscousContainer3D{uEltype <: Real} + u_transformed::Array{uEltype, 5} + # IDEA: Use SVector for fixed sized vectors? + gradients::Vector{Array{uEltype, 5}} + flux_viscous::Vector{Array{uEltype, 5}} + + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + _gradients::Vector{Vector{uEltype}} + _flux_viscous::Vector{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) for _ in 1:3], + [Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements) for _ in 1:3]) + 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!(cache_viscous._u_transformed, capacity) + for dim in 1:3 + resize!(cache_viscous._gradients[dim], capacity) + resize!(cache_viscous._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))) + + viscous_container.gradients = unsafe_wrap(Array, + pointer(viscous_container._gradients), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + viscous_container.flux_viscous = unsafe_wrap(Array, + pointer(viscous_container._flux_viscous), + (nvariables(equations), + nnodes(dg), nnodes(dg), nnodes(dg), + nelements(dg, cache))) + + return nothing +end \ No newline at end of file diff --git a/src/solvers/dgsem_tree/containers_viscous.jl b/src/solvers/dgsem_tree/containers_viscous.jl new file mode 100644 index 00000000000..2a8989db55c --- /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") \ No newline at end of file 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") From 40ad2668cfc6844cb8fbcbde71a3925390c006f6 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 09:57:18 +0200 Subject: [PATCH 19/53] check if mortars are present --- src/solvers/dgsem_tree/containers.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) 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 From 4d1914bcb257d441eb2e963a68061f046cb60d7f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 10:17:49 +0200 Subject: [PATCH 20/53] reuse `reinitialize_containers!` --- src/callbacks_step/amr_dg1d.jl | 43 +++++----------------------------- 1 file changed, 6 insertions(+), 37 deletions(-) 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 From 70c8e6641e078a55fefd0d01398ff63b7531893c Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 10:18:04 +0200 Subject: [PATCH 21/53] resize calls for parabolic helpers --- src/callbacks_step/amr_dg2d.jl | 44 ++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 1d37dfce034..5b36f22c82e 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -136,6 +136,28 @@ 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 +297,28 @@ 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, From a5db7d0b7d28ea0bbbe1591ef7c78a9c45b55f87 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 10:20:31 +0200 Subject: [PATCH 22/53] update analysis callbacks --- src/callbacks_step/analysis_dg2d.jl | 2 +- src/callbacks_step/analysis_dg3d.jl | 2 +- src/callbacks_step/analysis_dgmulti.jl | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) 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/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index dc294de9e7b..82a091250a8 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -140,7 +140,7 @@ function integrate(func::typeof(enstrophy), u, equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, dg::DGMulti, cache, cache_parabolic; normalize = true) - gradients_x, gradients_y, gradients_z = cache_parabolic.gradients + gradients_x, gradients_y, gradients_z = cache_parabolic.viscous_container.gradients # allocate local storage for gradients. # TODO: can we avoid allocating here? From c9d98a21ca831ce93ffb896c3f3ebb4ca135e26f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 10:24:43 +0200 Subject: [PATCH 23/53] Velocities for compr euler --- src/Trixi.jl | 2 +- src/equations/compressible_euler_1d.jl | 5 +++++ src/equations/compressible_euler_2d.jl | 10 ++++++++++ src/equations/compressible_euler_3d.jl | 15 +++++++++++++++ 4 files changed, 31 insertions(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index ec4d20558e5..b95ede67163 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -205,7 +205,7 @@ export initial_condition_eoc_test_coupled_euler_gravity, export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, cons2entropy, entropy2cons -export density, pressure, density_pressure, velocity, global_mean_vars, +export density, pressure, density_pressure, velocity, v1, v2, v3, global_mean_vars, equilibrium_distribution, waterheight_pressure export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 9204989e8be..6d9ad996973 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -973,6 +973,11 @@ end return rho_times_p end +@inline function v1(u, equations::CompressibleEulerEquations1D) + rho, rho_v1, _ = u + return rho_v1/rho +end + # Calculate thermodynamic entropy for a conservative state `cons` @inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations1D) # Pressure diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index 27b92f41953..ccc208d911d 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1428,6 +1428,16 @@ end return rho_times_p end +@inline function v1(u, equations::CompressibleEulerEquations2D) + rho, rho_v1, _, _ = u + return rho_v1/rho +end + +@inline function v2(u, equations::CompressibleEulerEquations2D) + rho, _, rho_v2, _ = u + return rho_v2/rho +end + # Calculates the entropy flux in direction "orientation" and the entropy variables for a state cons # NOTE: This method seems to work currently (b82534e) but is never used anywhere. Thus it is # commented here until someone uses it or writes a test for it. diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 7f25bde31fd..357851a18d7 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1514,6 +1514,21 @@ end return rho_times_p end +@inline function v1(u, equations::CompressibleEulerEquations3D) + rho, rho_v1, _, _, _ = u + return rho_v1/rho +end + +@inline function v2(u, equations::CompressibleEulerEquations3D) + rho, _, rho_v2, _, _ = u + return rho_v2/rho +end + +@inline function v3(u, equations::CompressibleEulerEquations3D) + rho, _, _, rho_v3, _ = u + return rho_v3/rho +end + # Calculate thermodynamic entropy for a conservative state `u` @inline function entropy_thermodynamic(u, equations::CompressibleEulerEquations3D) rho, _ = u From 694e6bc8a1379e4e76fe6cc400a5ba083770505d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 10:39:50 +0200 Subject: [PATCH 24/53] Init container --- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 8 ++++---- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 16 +++++++--------- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 16 +++++++--------- 3 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 20626471941..1d0a7e247ee 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -531,14 +531,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_1d(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 0da25230380..3ac41c278ee 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 @@ -906,21 +907,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 2745d312b37..b5719247fce 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 @@ -1096,21 +1097,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)...) From 07655a47b43a4a3446ca3728be73a032371a9ce0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 10:40:17 +0200 Subject: [PATCH 25/53] correct copy-paste error --- src/solvers/dgsem_tree/container_viscous_2d.jl | 6 +++--- src/solvers/dgsem_tree/container_viscous_3d.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl index 151011e9757..2a44d8f44aa 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -32,10 +32,10 @@ end # internal storage. function Base.resize!(viscous_container::ViscousContainer2D, equations, dg, cache) capacity = nvariables(equations) * nnodes(dg) * nnodes(dg) * nelements(dg, cache) - resize!(cache_viscous._u_transformed, capacity) + resize!(viscous_container._u_transformed, capacity) for dim in 1:2 - resize!(cache_viscous._gradients[dim], capacity) - resize!(cache_viscous._flux_viscous[dim], capacity) + resize!(viscous_container._gradients[dim], capacity) + resize!(viscous_container._flux_viscous[dim], capacity) end viscous_container.u_transformed = unsafe_wrap(Array, diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl index 9a7749d7a9b..778a1df737a 100644 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -33,10 +33,10 @@ end # internal storage. function Base.resize!(viscous_container::ViscousContainer3D, equations, dg, cache) capacity = nvariables(equations) * nnodes(dg) * nnodes(dg) * nnodes(dg) * nelements(dg, cache) - resize!(cache_viscous._u_transformed, capacity) + resize!(viscous_container._u_transformed, capacity) for dim in 1:3 - resize!(cache_viscous._gradients[dim], capacity) - resize!(cache_viscous._flux_viscous[dim], capacity) + resize!(viscous_container._gradients[dim], capacity) + resize!(viscous_container._flux_viscous[dim], capacity) end viscous_container.u_transformed = unsafe_wrap(Array, From edd82ce95b02b5a59478710164d8df4fdd037898 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 11:37:00 +0200 Subject: [PATCH 26/53] resize each dim --- .../dgsem_tree/container_viscous_2d.jl | 25 ++++++++++--------- .../dgsem_tree/container_viscous_3d.jl | 25 ++++++++++--------- 2 files changed, 26 insertions(+), 24 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl index 2a44d8f44aa..80554365aa5 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -44,17 +44,18 @@ function Base.resize!(viscous_container::ViscousContainer2D, equations, dg, cach nnodes(dg), nnodes(dg), nelements(dg, cache))) - viscous_container.gradients = unsafe_wrap(Array, - pointer(viscous_container._gradients), - (nvariables(equations), - nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - viscous_container.flux_viscous = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous), - (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 \ No newline at end of file diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl index 778a1df737a..30fa864825e 100644 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -45,17 +45,18 @@ function Base.resize!(viscous_container::ViscousContainer3D, equations, dg, cach nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache))) - viscous_container.gradients = unsafe_wrap(Array, - pointer(viscous_container._gradients), - (nvariables(equations), - nnodes(dg), nnodes(dg), nnodes(dg), - nelements(dg, cache))) - - viscous_container.flux_viscous = unsafe_wrap(Array, - pointer(viscous_container._flux_viscous), - (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 \ No newline at end of file From ba1ef1bedeff57d340803bb72f947ab3813d7ae0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 11:41:03 +0200 Subject: [PATCH 27/53] add dispatch --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 20 ++++++++++---------- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 21 ++++++++++----------- 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 3ac41c278ee..c4584a2cda5 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -132,10 +132,10 @@ function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{2}, P4estMe end # This is the version used when calculating the divergence of the viscous fluxes -function calc_volume_integral!(du, flux_viscous, +function calc_volume_integral!(du, flux_viscous::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache) + dg::DGSEM, cache) where {uEltype <: Real} @unpack derivative_dhat = dg.basis flux_viscous_x, flux_viscous_y = flux_viscous @@ -164,10 +164,10 @@ end # This is the version used when calculating the divergence of the viscous fluxes # We pass the `surface_integral` argument solely for dispatch -function prolong2interfaces!(cache_parabolic, flux_viscous, +function prolong2interfaces!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache) where {uEltype <: Real} @unpack interfaces = cache_parabolic @unpack orientations = interfaces @@ -240,10 +240,10 @@ function calc_interface_flux!(surface_flux_values, end # This is the version used when calculating the divergence of the viscous fluxes -function prolong2boundaries!(cache_parabolic, flux_viscous, +function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache) where {uEltype <: Real} @unpack boundaries = cache_parabolic @unpack orientations, neighbor_sides = boundaries flux_viscous_x, flux_viscous_y = flux_viscous @@ -288,10 +288,10 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, +function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, gradients, u_transformed, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic) + dg::DG, cache, cache_parabolic) where {uEltype <: Real} gradients_x, gradients_y = gradients flux_viscous_x, flux_viscous_y = flux_viscous # output arrays @@ -511,11 +511,11 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra return nothing end -function prolong2mortars!(cache, flux_viscous::Tuple{AbstractArray, AbstractArray}, +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] diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index b5719247fce..efade560303 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -132,10 +132,10 @@ function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{3}, P4estMe end # This is the version used when calculating the divergence of the viscous fluxes -function calc_volume_integral!(du, flux_viscous, +function calc_volume_integral!(du, flux_viscous::Vector{Array{uEltype, 5}}, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache) + dg::DGSEM, cache) where {uEltype <: Real} @unpack derivative_dhat = dg.basis flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @@ -171,10 +171,10 @@ end # This is the version used when calculating the divergence of the viscous fluxes # We pass the `surface_integral` argument solely for dispatch -function prolong2interfaces!(cache_parabolic, flux_viscous, +function prolong2interfaces!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 5}}, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache) where {uEltype <: Real} @unpack interfaces = cache_parabolic @unpack orientations = interfaces @@ -261,10 +261,10 @@ function calc_interface_flux!(surface_flux_values, end # This is the version used when calculating the divergence of the viscous fluxes -function prolong2boundaries!(cache_parabolic, flux_viscous, +function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 5}}, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) + surface_integral, dg::DG, cache) where {uEltype <: Real} @unpack boundaries = cache_parabolic @unpack orientations, neighbor_sides = boundaries flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @@ -336,10 +336,10 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, +function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, gradients, u_transformed, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic) + dg::DG, cache, cache_parabolic) where {uEltype <: Real} gradients_x, gradients_y, gradients_z = gradients flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays @@ -595,12 +595,11 @@ function calc_boundary_flux_by_direction_divergence!(surface_flux_values::Abstra end 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 From 7d351e512f967fdbea03733123f147f5dbf96adf Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 11:47:19 +0200 Subject: [PATCH 28/53] Add AMR for shear layer --- ...r.jl => elixir_navierstokes_shearlayer.jl} | 4 +- .../elixir_navierstokes_shearlayer_amr.jl | 85 +++++++++++++++++++ 2 files changed, 87 insertions(+), 2 deletions(-) rename examples/tree_2d_dgsem/{elixir_navierstokes_shear_layer.jl => elixir_navierstokes_shearlayer.jl} (94%) create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl similarity index 94% rename from examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl rename to examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl index dd26fd8097b..6b111659273 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shear_layer.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl @@ -22,8 +22,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) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl new file mode 100644 index 00000000000..10132f4ec1a --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -0,0 +1,85 @@ + +using OrdinaryDiffEq, Plots +using Trixi + +############################################################################### +# semidiscretization of the compressible Navier-Stokes equations + +# TODO: parabolic; unify names of these accessor functions +prandtl_number() = 0.72 +mu() = 1.0/3.0 * 10^(-5) + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), + Prandtl=prandtl_number()) + +function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) + # Shear layer parameters + k = 80 + delta = 0.05 + u0 = 1.0 + + Ms = 0.1 # maximum Mach number + + rho = 1.0 + 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) +end +initial_condition = initial_condition_shear_layer + +volume_flux = flux_ranocha +solver = DGSEM(polydeg=3, surface_flux=flux_hllc, + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + n_cells_max=100_000) + + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.7) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval,) + +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=10, + adapt_initial_condition=true, + adapt_initial_condition_only_refine=true) + +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, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary + +pd = PlotData2D(sol) +Plots.plot(pd["v1"]) +Plots.plot(getmesh(pd)) \ No newline at end of file From e60817454b7637133fa25ee2c5158a69d7090e74 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 12:01:28 +0200 Subject: [PATCH 29/53] USe only amr shear layer --- .../elixir_navierstokes_shearlayer.jl | 73 ------------------- .../elixir_navierstokes_shearlayer_amr.jl | 6 +- test/test_parabolic_2d.jl | 7 ++ 3 files changed, 10 insertions(+), 76 deletions(-) delete mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl deleted file mode 100644 index 6b111659273..00000000000 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer.jl +++ /dev/null @@ -1,73 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the compressible Navier-Stokes equations - -# TODO: parabolic; unify names of these accessor functions -prandtl_number() = 0.72 -mu() = 1.0/3.0 * 10^(-3) # equivalent to Re = 3000 - -equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), - Prandtl=prandtl_number()) - -function initial_condition_shear_layer(x, t, equations::CompressibleEulerEquations2D) - # Shear layer parameters - k = 80 - delta = 0.05 - u0 = 1.0 - - Ms = 0.1 # maximum Mach number - - rho = 1.0 - 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) -end -initial_condition = initial_condition_shear_layer - -volume_flux = flux_ranocha -solver = DGSEM(polydeg=3, surface_flux=flux_hllc, - volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) - -coordinates_min = (0.0, 0.0) -coordinates_max = (1.0, 1.0) -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=4, - n_cells_max=100_000) - - -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), - initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 2.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)) - -alive_callback = AliveCallback(analysis_interval=analysis_interval,) - -callbacks = CallbackSet(summary_callback, - analysis_callback, - alive_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) -summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 10132f4ec1a..1da63588bcf 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -52,7 +52,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 500 +analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval,) @@ -63,7 +63,7 @@ amr_controller = ControllerThreeLevel(semi, amr_indicator, med_level = 5, med_threshold=0.2, max_level = 7, max_threshold=0.5) amr_callback = AMRCallback(semi, amr_controller, - interval=10, + interval=50, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) @@ -75,7 +75,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation -time_int_tol = 1e-8 +time_int_tol = 1e-7 sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, ode_default_options()..., callback=callbacks) summary_callback() # print the timer summary diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 3fff4382cd1..a6cfd089c47 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -249,6 +249,13 @@ isdir(outdir) && rm(outdir, recursive=true) ) 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.005352707453147317, 0.4077270907626271, 0.43758381263181434, 1.1758133763041454], + linf = [0.034784713829601244, 1.1650897299747531, 1.484016014155391, 8.71087741820071] + ) + 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), From 6574bf59536e797b169137c8a9fbb24197c27473 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 13:26:19 +0200 Subject: [PATCH 30/53] first steps towards p4est parabolic amr --- src/callbacks_step/analysis_dgmulti.jl | 2 +- src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 14 ++++---- src/solvers/dgsem_p4est/dg_3d_parabolic.jl | 14 ++++---- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 33 ++++++++++++++++- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 42 +++++++++++++++++++++- 5 files changed, 86 insertions(+), 19 deletions(-) diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index 82a091250a8..dc294de9e7b 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/src/callbacks_step/analysis_dgmulti.jl @@ -140,7 +140,7 @@ function integrate(func::typeof(enstrophy), u, equations, equations_parabolic::CompressibleNavierStokesDiffusion3D, dg::DGMulti, cache, cache_parabolic; normalize = true) - gradients_x, gradients_y, gradients_z = cache_parabolic.viscous_container.gradients + gradients_x, gradients_y, gradients_z = cache_parabolic.gradients # allocate local storage for gradients. # TODO: can we avoid allocating here? 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/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index c4584a2cda5..3d3e42d1a40 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -288,9 +288,40 @@ function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype return nothing end -function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, 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) + gradients_x, gradients_y = gradients + flux_viscous_x, flux_viscous_y = flux_viscous # output arrays + + @threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + # Get solution and gradients + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, + element) + gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, + element) + gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, + element) + + # Calculate viscous flux and store each component for later use + flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node), 1, + equations_parabolic) + flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node), 2, + equations_parabolic) + set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, + i, j, element) + set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, + i, j, element) + end + end +end + +function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, + gradients::Vector{Array{uEltype, 4}}, u_transformed, + mesh::TreeMesh{2}, + equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) where {uEltype <: Real} gradients_x, gradients_y = gradients flux_viscous_x, flux_viscous_y = flux_viscous # output arrays diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index efade560303..1f625b74b3c 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -336,9 +336,49 @@ function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype return nothing end -function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, 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) + gradients_x, gradients_y, gradients_z = gradients + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays + + @threaded for element in eachelement(dg, cache) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + # Get solution and gradients + u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, k, + element) + gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, + k, element) + gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, + k, element) + gradients_3_node = get_node_vars(gradients_z, equations_parabolic, dg, i, j, + k, element) + + # Calculate viscous flux and store each component for later use + flux_viscous_node_x = flux(u_node, + (gradients_1_node, gradients_2_node, + gradients_3_node), 1, equations_parabolic) + flux_viscous_node_y = flux(u_node, + (gradients_1_node, gradients_2_node, + gradients_3_node), 2, equations_parabolic) + flux_viscous_node_z = flux(u_node, + (gradients_1_node, gradients_2_node, + gradients_3_node), 3, equations_parabolic) + set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, + i, j, k, element) + set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, + i, j, k, element) + set_node_vars!(flux_viscous_z, flux_viscous_node_z, equations_parabolic, dg, + i, j, k, element) + end + end +end + +function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, + gradients::Vector{Array{uEltype, 5}}, u_transformed, + mesh::TreeMesh{3}, + equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) where {uEltype <: Real} gradients_x, gradients_y, gradients_z = gradients flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays From 4c354301d6d6467f8fe1451c4b2a3f1bebd813b0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 13:28:03 +0200 Subject: [PATCH 31/53] Add tests --- .../elixir_navierstokes_shearlayer_amr.jl | 10 +- .../elixir_advection_diffusion_amr.jl | 94 +++++++++++++++++++ test/test_parabolic_2d.jl | 4 +- test/test_parabolic_3d.jl | 7 ++ 4 files changed, 112 insertions(+), 3 deletions(-) create mode 100644 examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 1da63588bcf..86a278f2c60 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -13,6 +13,14 @@ 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 @@ -37,7 +45,7 @@ solver = DGSEM(polydeg=3, surface_flux=flux_hllc, coordinates_min = (0.0, 0.0) coordinates_max = (1.0, 1.0) mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=3, + initial_refinement_level=4, n_cells_max=100_000) 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..08f678138ed --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl @@ -0,0 +1,94 @@ + +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 = LaplaceDiffusion2D(diffusivity(), equations) + +initial_condition = initial_condition_gauss +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.1) +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/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index a6cfd089c47..24de7f13679 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -251,8 +251,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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.005352707453147317, 0.4077270907626271, 0.43758381263181434, 1.1758133763041454], - linf = [0.034784713829601244, 1.1650897299747531, 1.484016014155391, 8.71087741820071] + l2 = [0.005353306916161318, 0.4077336028873465, 0.43756829240941464, 1.1759497303478266], + linf = [0.034791663521213545, 1.1650592394066237, 1.4841924833734612, 8.713274087336629] ) end diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index ded052fb9d3..05a173ec5d2 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -157,6 +157,13 @@ isdir(outdir) && rm(outdir, recursive=true) linf = [0.0006696415247340326, 0.03442565722527785, 0.03442565722577423, 0.06295407168705314, 0.032857472756916195] ) 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.00017912610056132068], + linf = [0.0007821751390747478] + ) + end end # Clean up afterwards: delete Trixi.jl output directory From 21d29f84d5332a971ee8e4f616ae0883e43ff0b6 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 13:30:19 +0200 Subject: [PATCH 32/53] remove plots --- .../tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 86a278f2c60..7661cc5b44b 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -1,5 +1,5 @@ -using OrdinaryDiffEq, Plots +using OrdinaryDiffEq using Trixi ############################################################################### @@ -86,8 +86,4 @@ callbacks = CallbackSet(summary_callback, time_int_tol = 1e-7 sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, ode_default_options()..., callback=callbacks) -summary_callback() # print the timer summary - -pd = PlotData2D(sol) -Plots.plot(pd["v1"]) -Plots.plot(getmesh(pd)) \ No newline at end of file +summary_callback() # print the timer summary \ No newline at end of file From 29bd2132c4b38e3208d9962893a80f65eb2dcb10 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 13:31:24 +0200 Subject: [PATCH 33/53] Format --- src/callbacks_step/amr_dg2d.jl | 10 ++- src/equations/compressible_euler_1d.jl | 2 +- src/equations/compressible_euler_2d.jl | 4 +- src/equations/compressible_euler_3d.jl | 6 +- .../dgsem_tree/container_viscous_1d.jl | 4 +- .../dgsem_tree/container_viscous_2d.jl | 55 ++++++------ .../dgsem_tree/container_viscous_3d.jl | 86 ++++++++++--------- src/solvers/dgsem_tree/containers_viscous.jl | 2 +- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 2 +- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 2 +- 10 files changed, 90 insertions(+), 83 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 5b36f22c82e..6395a9f348f 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -151,8 +151,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeMe # 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") + !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 @@ -312,8 +313,9 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, TreeM # 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") + !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 diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 6d9ad996973..cef5a268c82 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -975,7 +975,7 @@ end @inline function v1(u, equations::CompressibleEulerEquations1D) rho, rho_v1, _ = u - return rho_v1/rho + return rho_v1 / rho end # Calculate thermodynamic entropy for a conservative state `cons` diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index ccc208d911d..d61bba8e6a6 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1430,12 +1430,12 @@ end @inline function v1(u, equations::CompressibleEulerEquations2D) rho, rho_v1, _, _ = u - return rho_v1/rho + return rho_v1 / rho end @inline function v2(u, equations::CompressibleEulerEquations2D) rho, _, rho_v2, _ = u - return rho_v2/rho + return rho_v2 / rho end # Calculates the entropy flux in direction "orientation" and the entropy variables for a state cons diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index 357851a18d7..af746c37025 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1516,17 +1516,17 @@ end @inline function v1(u, equations::CompressibleEulerEquations3D) rho, rho_v1, _, _, _ = u - return rho_v1/rho + return rho_v1 / rho end @inline function v2(u, equations::CompressibleEulerEquations3D) rho, _, rho_v2, _, _ = u - return rho_v2/rho + return rho_v2 / rho end @inline function v3(u, equations::CompressibleEulerEquations3D) rho, _, _, rho_v3, _ = u - return rho_v3/rho + return rho_v3 / rho end # Calculate thermodynamic entropy for a conservative state `u` diff --git a/src/solvers/dgsem_tree/container_viscous_1d.jl b/src/solvers/dgsem_tree/container_viscous_1d.jl index 3f9943aea1c..71c68dfc6df 100644 --- a/src/solvers/dgsem_tree/container_viscous_1d.jl +++ b/src/solvers/dgsem_tree/container_viscous_1d.jl @@ -20,8 +20,8 @@ mutable struct ViscousContainer1D{uEltype <: Real} end function init_viscous_container_1d(n_vars::Integer, n_nodes::Integer, - n_elements::Integer, - ::Type{uEltype}) where {uEltype <: Real} + 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 index 80554365aa5..eb62178b66d 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -9,14 +9,15 @@ mutable struct ViscousContainer2D{uEltype <: Real} _gradients::Vector{Vector{uEltype}} _flux_viscous::Vector{Vector{uEltype}} - function ViscousContainer2D{uEltype}(n_vars::Integer, n_nodes::Integer, n_elements::Integer) where {uEltype <: Real} + 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) for _ in 1:2], [Vector{uEltype}(undef, n_vars * n_nodes^2 * n_elements) for _ in 1:2]) - end + end end function init_viscous_container_2d(n_vars::Integer, n_nodes::Integer, @@ -31,31 +32,31 @@ end # `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 + 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))) + 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))) + 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 \ No newline at end of file + 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 index 30fa864825e..0b44cc6cc18 100644 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -1,25 +1,27 @@ mutable struct ViscousContainer3D{uEltype <: Real} - u_transformed::Array{uEltype, 5} - # IDEA: Use SVector for fixed sized vectors? - gradients::Vector{Array{uEltype, 5}} - flux_viscous::Vector{Array{uEltype, 5}} + u_transformed::Array{uEltype, 5} + # IDEA: Use SVector for fixed sized vectors? + gradients::Vector{Array{uEltype, 5}} + flux_viscous::Vector{Array{uEltype, 5}} - # internal `resize!`able storage - _u_transformed::Vector{uEltype} - _gradients::Vector{Vector{uEltype}} - _flux_viscous::Vector{Vector{uEltype}} + # internal `resize!`able storage + _u_transformed::Vector{uEltype} + _gradients::Vector{Vector{uEltype}} + _flux_viscous::Vector{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) for _ in 1:3], - [Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements) for _ in 1:3]) - end + 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) for _ in 1:3], + [Vector{uEltype}(undef, n_vars * n_nodes^3 * n_elements) for _ in 1:3]) + end end - function init_viscous_container_3d(n_vars::Integer, n_nodes::Integer, n_elements::Integer, ::Type{uEltype}) where {uEltype <: Real} @@ -32,31 +34,33 @@ end # `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))) + 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 - 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.u_transformed = unsafe_wrap(Array, + pointer(viscous_container._u_transformed), + (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), + 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))) - end - return nothing -end \ No newline at end of file + + 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_viscous.jl b/src/solvers/dgsem_tree/containers_viscous.jl index 2a8989db55c..444f2cb7303 100644 --- a/src/solvers/dgsem_tree/containers_viscous.jl +++ b/src/solvers/dgsem_tree/containers_viscous.jl @@ -1,4 +1,4 @@ # Dimension-specific implementations include("container_viscous_1d.jl") include("container_viscous_2d.jl") -include("container_viscous_3d.jl") \ No newline at end of file +include("container_viscous_3d.jl") diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 3d3e42d1a40..582b624a199 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -318,7 +318,7 @@ function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, end end -function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, +function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, gradients::Vector{Array{uEltype, 4}}, u_transformed, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 1f625b74b3c..efed90d83aa 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -375,7 +375,7 @@ function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, end end -function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, +function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, gradients::Vector{Array{uEltype, 5}}, u_transformed, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, From cb3eac8e2c5371c345485265de7ab539af41c402 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 14:22:25 +0200 Subject: [PATCH 34/53] remove redundant line --- examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl index 08f678138ed..ded984bff2d 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl @@ -11,7 +11,6 @@ equations = LinearScalarAdvectionEquation3D(advection_velocity) diffusivity() = 5.0e-4 equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) -initial_condition = initial_condition_gauss solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) coordinates_min = (-1.0, -1.0, -1.0) From 7e68d945c2ec64fa27cce7df2c85b408b988299b Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 11 Sep 2023 15:13:50 +0200 Subject: [PATCH 35/53] platform independent tests --- .../elixir_navierstokes_shearlayer_amr.jl | 13 ++++++++----- test/test_parabolic_2d.jl | 4 ++-- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 7661cc5b44b..82c85aa16ba 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -60,7 +60,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 1000 +analysis_interval = 2000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval,) @@ -75,15 +75,18 @@ amr_callback = AMRCallback(semi, amr_controller, adapt_initial_condition=true, adapt_initial_condition_only_refine=true) +stepsize_callback = StepsizeCallback(cfl=1.3) + callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - amr_callback) + amr_callback, + stepsize_callback) ############################################################################### # run the simulation -time_int_tol = 1e-7 -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/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 24de7f13679..97447c0f438 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -251,8 +251,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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.005353306916161318, 0.4077336028873465, 0.43756829240941464, 1.1759497303478266], - linf = [0.034791663521213545, 1.1650592394066237, 1.4841924833734612, 8.713274087336629] + l2 = [0.00526017743452336, 0.4130430692895672, 0.4310996183791349, 1.1544344171604635], + linf = [0.03492185879198495, 1.392635891671335, 1.357551616406459, 8.713760873018146] ) end From 5c20c4fb0a7da37a4178ac08dad2243ee0f60eef Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 12 Sep 2023 15:44:41 +0200 Subject: [PATCH 36/53] No need for different flux_viscous comps after adding container_viscous to p4est --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 32 +----------------- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 41 +---------------------- 2 files changed, 2 insertions(+), 71 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 582b624a199..42d992ab52e 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -288,39 +288,9 @@ function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, - equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic) - gradients_x, gradients_y = gradients - flux_viscous_x, flux_viscous_y = flux_viscous # output arrays - - @threaded for element in eachelement(dg, cache) - for j in eachnode(dg), i in eachnode(dg) - # Get solution and gradients - u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, - element) - gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, - element) - gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, - element) - - # Calculate viscous flux and store each component for later use - flux_viscous_node_x = flux(u_node, (gradients_1_node, gradients_2_node), 1, - equations_parabolic) - flux_viscous_node_y = flux(u_node, (gradients_1_node, gradients_2_node), 2, - equations_parabolic) - set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, - i, j, element) - set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, - i, j, element) - end - end -end - function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, gradients::Vector{Array{uEltype, 4}}, u_transformed, - mesh::TreeMesh{2}, + mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) where {uEltype <: Real} gradients_x, gradients_y = gradients diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index efed90d83aa..5cb2448a1bd 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -336,48 +336,9 @@ function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, - equations_parabolic::AbstractEquationsParabolic, - dg::DG, cache, cache_parabolic) - gradients_x, gradients_y, gradients_z = gradients - flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays - - @threaded for element in eachelement(dg, cache) - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - # Get solution and gradients - u_node = get_node_vars(u_transformed, equations_parabolic, dg, i, j, k, - element) - gradients_1_node = get_node_vars(gradients_x, equations_parabolic, dg, i, j, - k, element) - gradients_2_node = get_node_vars(gradients_y, equations_parabolic, dg, i, j, - k, element) - gradients_3_node = get_node_vars(gradients_z, equations_parabolic, dg, i, j, - k, element) - - # Calculate viscous flux and store each component for later use - flux_viscous_node_x = flux(u_node, - (gradients_1_node, gradients_2_node, - gradients_3_node), 1, equations_parabolic) - flux_viscous_node_y = flux(u_node, - (gradients_1_node, gradients_2_node, - gradients_3_node), 2, equations_parabolic) - flux_viscous_node_z = flux(u_node, - (gradients_1_node, gradients_2_node, - gradients_3_node), 3, equations_parabolic) - set_node_vars!(flux_viscous_x, flux_viscous_node_x, equations_parabolic, dg, - i, j, k, element) - set_node_vars!(flux_viscous_y, flux_viscous_node_y, equations_parabolic, dg, - i, j, k, element) - set_node_vars!(flux_viscous_z, flux_viscous_node_z, equations_parabolic, dg, - i, j, k, element) - end - end -end - function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, gradients::Vector{Array{uEltype, 5}}, u_transformed, - mesh::TreeMesh{3}, + mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) where {uEltype <: Real} gradients_x, gradients_y, gradients_z = gradients From 21ccff481c8918f3732791ec4fddeb1e797ab660 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 12 Sep 2023 16:23:59 +0200 Subject: [PATCH 37/53] Laplace 3d --- .../elixir_advection_diffusion_amr.jl | 2 +- src/Trixi.jl | 2 +- src/equations/equations_parabolic.jl | 1 + src/equations/laplace_diffusion_3d.jl | 72 +++++++++++++++++++ test/test_parabolic_3d.jl | 4 +- 5 files changed, 77 insertions(+), 4 deletions(-) create mode 100644 src/equations/laplace_diffusion_3d.jl diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl index ded984bff2d..4f498af5421 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl @@ -9,7 +9,7 @@ advection_velocity = (0.2, -0.7, 0.5) equations = LinearScalarAdvectionEquation3D(advection_velocity) diffusivity() = 5.0e-4 -equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) +equations_parabolic = LaplaceDiffusion3D(diffusivity(), equations) solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs) diff --git a/src/Trixi.jl b/src/Trixi.jl index b95ede67163..57e756031c2 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -151,7 +151,7 @@ export AcousticPerturbationEquations2D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, LinearizedEulerEquations2D -export LaplaceDiffusion1D, LaplaceDiffusion2D, +export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D, CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D 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_3d.jl b/src/equations/laplace_diffusion_3d.jl new file mode 100644 index 00000000000..3988ce7144b --- /dev/null +++ b/src/equations/laplace_diffusion_3d.jl @@ -0,0 +1,72 @@ +@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 + +# no orientation specified since the flux is vector-valued +function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion3D) + dudx, dudy, dudz = gradients + if orientation == 1 + 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/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 05a173ec5d2..95cb6052383 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -160,8 +160,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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.00017912610056132068], - linf = [0.0007821751390747478] + l2 = [0.0001791112986201992], + linf = [0.0007775704284449514] ) end end From 8df117b00be1a71c89ab9cdcd5c2493370f97ec0 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 20 Sep 2023 11:45:30 +0200 Subject: [PATCH 38/53] Longer times to allow converage to hit coarsen! --- examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl | 2 +- test/test_parabolic_3d.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl index 4f498af5421..0fab685b642 100644 --- a/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl +++ b/examples/tree_3d_dgsem/elixir_advection_diffusion_amr.jl @@ -50,7 +50,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 0.1) +tspan = (0.0, 0.2) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index 923511860a5..ef91550436d 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -160,8 +160,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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.0001791112986201992], - linf = [0.0007775704284449514] + l2 = [0.000355780485397024], + linf = [0.0010810770271614256] ) end end From cdaa865da7ad10de528e53ce19956b6f95eb908e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 20 Sep 2023 11:57:51 +0200 Subject: [PATCH 39/53] Increase testing of Laplace 3D --- .../elixir_advection_diffusion_nonperiodic.jl | 91 +++++++++++++++++++ test/test_parabolic_3d.jl | 7 ++ 2 files changed, 98 insertions(+) create mode 100644 examples/tree_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl 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/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index ef91550436d..3014ee9c6cf 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -164,6 +164,13 @@ isdir(outdir) && rm(outdir, recursive=true) 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 From 4699a1003fcb40d6b756be38a0259c6951d8225d Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 20 Sep 2023 14:14:24 +0200 Subject: [PATCH 40/53] Add tests for velocities --- test/test_unit.jl | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/test/test_unit.jl b/test/test_unit.jl index 5c5291c2430..30bc9b6d97b 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -528,6 +528,28 @@ isdir(outdir) && rm(outdir, recursive=true) end end + @timed_testset "Euler test velocity functions" begin + rho, v1, v2, v3, p = 1.0, 0.1, 0.2, 0.3, 2.0 + + let equations = CompressibleEulerEquations1D(1.4) + cons_vars = prim2cons(SVector(rho,v1,p),equations) + @test v1 ≈ Trixi.v1(cons_vars, equations) + end + + let equations = CompressibleEulerEquations2D(1.4) + cons_vars = prim2cons(SVector(rho,v1,v2,p),equations) + @test v1 ≈ Trixi.v1(cons_vars, equations) + @test v2 ≈ Trixi.v2(cons_vars, equations) + end + + let equations = CompressibleEulerEquations3D(1.4) + cons_vars = prim2cons(SVector(rho,v1,v2,v3,p),equations) + @test v1 ≈ Trixi.v1(cons_vars, equations) + @test v2 ≈ Trixi.v2(cons_vars, equations) + @test v3 ≈ Trixi.v3(cons_vars, equations) + end + end + @timed_testset "Shallow water conversion between conservative/entropy variables" begin H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 From 389d89c7ec3624c4a306dc26213c6e49c1d1d200 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 5 Oct 2023 08:55:37 +0200 Subject: [PATCH 41/53] remove comment --- src/solvers/dgsem_tree/container_viscous_2d.jl | 1 - src/solvers/dgsem_tree/container_viscous_3d.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl index eb62178b66d..fd28e12a99a 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -1,6 +1,5 @@ mutable struct ViscousContainer2D{uEltype <: Real} u_transformed::Array{uEltype, 4} - # IDEA: Use SVector for fixed sized vectors? gradients::Vector{Array{uEltype, 4}} flux_viscous::Vector{Array{uEltype, 4}} diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl index 0b44cc6cc18..6eec0c61279 100644 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -1,6 +1,5 @@ mutable struct ViscousContainer3D{uEltype <: Real} u_transformed::Array{uEltype, 5} - # IDEA: Use SVector for fixed sized vectors? gradients::Vector{Array{uEltype, 5}} flux_viscous::Vector{Array{uEltype, 5}} From a5ddabe6a99e96f5df73e109d0600034b8564f82 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 11 Oct 2023 08:50:47 +0200 Subject: [PATCH 42/53] Add comments --- examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl | 1 + src/solvers/dgsem_tree/dg_2d_parabolic.jl | 2 ++ src/solvers/dgsem_tree/dg_3d_parabolic.jl | 2 ++ 3 files changed, 5 insertions(+) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 82c85aa16ba..2df05fbf210 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -65,6 +65,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval) alive_callback = AliveCallback(analysis_interval=analysis_interval,) +# This uses velocity-based AMR amr_indicator = IndicatorLöhner(semi, variable=v1) amr_controller = ControllerThreeLevel(semi, amr_indicator, base_level = 3, diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 220ad75fbe6..e1c2daf8f08 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -515,6 +515,8 @@ 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`. function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 4bcbdb0b8c5..34c75070d82 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -598,6 +598,8 @@ 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`. function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 5}}, mesh::TreeMesh{3}, From 2d960fc58cfe96cddf1a41d31a4d0d8c273520ce Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 11 Oct 2023 17:39:03 +0200 Subject: [PATCH 43/53] Remove some specializations --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 18 +++++++++--------- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 18 +++++++++--------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index e1c2daf8f08..a249e373fcb 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -133,10 +133,10 @@ function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{2}, P4estMe end # This is the version used when calculating the divergence of the viscous fluxes -function calc_volume_integral!(du, flux_viscous::Vector{Array{uEltype, 4}}, +function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache) where {uEltype <: Real} + dg::DGSEM, cache) @unpack derivative_dhat = dg.basis flux_viscous_x, flux_viscous_y = flux_viscous @@ -165,10 +165,10 @@ end # This is the version used when calculating the divergence of the viscous fluxes # We pass the `surface_integral` argument solely for dispatch -function prolong2interfaces!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 4}}, +function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) where {uEltype <: Real} + surface_integral, dg::DG, cache) @unpack interfaces = cache_parabolic @unpack orientations, neighbor_ids = interfaces interfaces_u = interfaces.u @@ -242,10 +242,10 @@ function calc_interface_flux!(surface_flux_values, end # This is the version used when calculating the divergence of the viscous fluxes -function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 4}}, +function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) where {uEltype <: Real} + surface_integral, dg::DG, cache) @unpack boundaries = cache_parabolic @unpack orientations, neighbor_sides, neighbor_ids = boundaries boundaries_u = boundaries.u @@ -291,11 +291,11 @@ function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype return nothing end -function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 4}}, - gradients::Vector{Array{uEltype, 4}}, 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) where {uEltype <: Real} + dg::DG, cache, cache_parabolic) gradients_x, gradients_y = gradients flux_viscous_x, flux_viscous_y = flux_viscous # output arrays diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 34c75070d82..613ec0a0991 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -133,10 +133,10 @@ function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{3}, P4estMe end # This is the version used when calculating the divergence of the viscous fluxes -function calc_volume_integral!(du, flux_viscous::Vector{Array{uEltype, 5}}, +function calc_volume_integral!(du, flux_viscous, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, - dg::DGSEM, cache) where {uEltype <: Real} + dg::DGSEM, cache) @unpack derivative_dhat = dg.basis flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @@ -172,10 +172,10 @@ end # This is the version used when calculating the divergence of the viscous fluxes # We pass the `surface_integral` argument solely for dispatch -function prolong2interfaces!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 5}}, +function prolong2interfaces!(cache_parabolic, flux_viscous, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) where {uEltype <: Real} + surface_integral, dg::DG, cache) @unpack interfaces = cache_parabolic @unpack orientations, neighbor_ids = interfaces interfaces_u = interfaces.u @@ -263,10 +263,10 @@ function calc_interface_flux!(surface_flux_values, end # This is the version used when calculating the divergence of the viscous fluxes -function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype, 5}}, +function prolong2boundaries!(cache_parabolic, flux_viscous, mesh::TreeMesh{3}, equations_parabolic::AbstractEquationsParabolic, - surface_integral, dg::DG, cache) where {uEltype <: Real} + surface_integral, dg::DG, cache) @unpack boundaries = cache_parabolic @unpack orientations, neighbor_sides, neighbor_ids = boundaries boundaries_u = boundaries.u @@ -339,11 +339,11 @@ function prolong2boundaries!(cache_parabolic, flux_viscous::Vector{Array{uEltype return nothing end -function calc_viscous_fluxes!(flux_viscous::Vector{Array{uEltype, 5}}, - gradients::Vector{Array{uEltype, 5}}, 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) where {uEltype <: Real} + dg::DG, cache, cache_parabolic) gradients_x, gradients_y, gradients_z = gradients flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous # output arrays From 9435d4e9f3b65b8c710e00bd5c4482ddafd3dc76 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 11 Oct 2023 17:57:28 +0200 Subject: [PATCH 44/53] Add comments --- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 3 +++ src/solvers/dgsem_tree/dg_3d_parabolic.jl | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 36bd6302bd3..06abff5e85b 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -517,6 +517,9 @@ 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::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, equations_parabolic::AbstractEquationsParabolic, diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 86fafb0e0b9..2561c5fe5b0 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -600,6 +600,9 @@ 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::Vector{Array{uEltype, 5}}, mesh::TreeMesh{3}, From 77b3b02bc94ba5cb12bb2c69c8ac04ffc12ac84f Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Oct 2023 09:17:26 +0200 Subject: [PATCH 45/53] Use tuple for outer, fixed size datastruct for internal quantities --- src/solvers/dgsem_tree/container_viscous_2d.jl | 13 +++++++++---- src/solvers/dgsem_tree/container_viscous_3d.jl | 15 +++++++++++---- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl index fd28e12a99a..5381cd4dba5 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -1,12 +1,15 @@ mutable struct ViscousContainer2D{uEltype <: Real} u_transformed::Array{uEltype, 4} + # Using an outer fixed-size datastructure leads to nasty implementations, + # avoid this here. gradients::Vector{Array{uEltype, 4}} flux_viscous::Vector{Array{uEltype, 4}} # internal `resize!`able storage _u_transformed::Vector{uEltype} - _gradients::Vector{Vector{uEltype}} - _flux_viscous::Vector{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} @@ -14,8 +17,10 @@ mutable struct ViscousContainer2D{uEltype <: Real} [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) for _ in 1:2], - [Vector{uEltype}(undef, n_vars * n_nodes^2 * 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)),) end end diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl index 6eec0c61279..083bb5284a5 100644 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -1,12 +1,15 @@ mutable struct ViscousContainer3D{uEltype <: Real} u_transformed::Array{uEltype, 5} + # Using an outer fixed-size datastructure leads to nasty implementations, + # avoid this here. gradients::Vector{Array{uEltype, 5}} flux_viscous::Vector{Array{uEltype, 5}} # internal `resize!`able storage _u_transformed::Vector{uEltype} - _gradients::Vector{Vector{uEltype}} - _flux_viscous::Vector{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} @@ -16,8 +19,12 @@ mutable struct ViscousContainer3D{uEltype <: Real} [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) for _ in 1:3], - [Vector{uEltype}(undef, n_vars * n_nodes^3 * 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))) end end From 2297486bef1a9ce7cb9c04d6f35842ca486758fc Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Thu, 12 Oct 2023 09:27:26 +0200 Subject: [PATCH 46/53] Format --- src/solvers/dgsem_tree/container_viscous_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl index 5381cd4dba5..57cb3728a11 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -20,7 +20,7 @@ mutable struct ViscousContainer2D{uEltype <: Real} (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 From 307915377fe78ad57e1acd9cfb98b8a0dbee85ce Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 18 Oct 2023 17:29:22 +0200 Subject: [PATCH 47/53] Add comments --- src/solvers/dgsem_tree/container_viscous_2d.jl | 5 ++++- src/solvers/dgsem_tree/container_viscous_3d.jl | 5 ++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_tree/container_viscous_2d.jl b/src/solvers/dgsem_tree/container_viscous_2d.jl index 57cb3728a11..bd7ff413af5 100644 --- a/src/solvers/dgsem_tree/container_viscous_2d.jl +++ b/src/solvers/dgsem_tree/container_viscous_2d.jl @@ -1,7 +1,10 @@ mutable struct ViscousContainer2D{uEltype <: Real} u_transformed::Array{uEltype, 4} # Using an outer fixed-size datastructure leads to nasty implementations, - # avoid this here. + # 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}} diff --git a/src/solvers/dgsem_tree/container_viscous_3d.jl b/src/solvers/dgsem_tree/container_viscous_3d.jl index 083bb5284a5..64d283fe189 100644 --- a/src/solvers/dgsem_tree/container_viscous_3d.jl +++ b/src/solvers/dgsem_tree/container_viscous_3d.jl @@ -1,7 +1,10 @@ mutable struct ViscousContainer3D{uEltype <: Real} u_transformed::Array{uEltype, 5} # Using an outer fixed-size datastructure leads to nasty implementations, - # avoid this here. + # 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}} From 9c32ae26e774b690aa0101aa86a8e44549d4eca1 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sun, 22 Oct 2023 10:54:05 +0200 Subject: [PATCH 48/53] Update examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl Co-authored-by: Michael Schlottke-Lakemper --- examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 2df05fbf210..897c67588eb 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -7,7 +7,7 @@ using Trixi # TODO: parabolic; unify names of these accessor functions prandtl_number() = 0.72 -mu() = 1.0/3.0 * 10^(-5) +mu() = 1.0/3.0 * 10^(-5) # equivalent to Re = 30,000 equations = CompressibleEulerEquations2D(1.4) equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), From 67c1027b3033da64a59c37cf8d2dd203d8b5579e Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Sun, 22 Oct 2023 11:31:27 +0200 Subject: [PATCH 49/53] Update src/Trixi.jl Co-authored-by: Michael Schlottke-Lakemper --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 2374c5d9c51..5cb3cf0a9fe 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -208,7 +208,7 @@ export initial_condition_eoc_test_coupled_euler_gravity, export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean, cons2entropy, entropy2cons -export density, pressure, density_pressure, velocity, v1, v2, v3, global_mean_vars, +export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, From 09f1e8571e8d85fbc4ddd2c38826164437521e88 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Oct 2023 11:45:35 +0200 Subject: [PATCH 50/53] Move velocity into elixir --- .../elixir_navierstokes_shearlayer_amr.jl | 6 +++++- src/equations/compressible_euler_1d.jl | 5 ----- src/equations/compressible_euler_2d.jl | 10 ---------- src/equations/compressible_euler_3d.jl | 15 --------------- test/test_parabolic_2d.jl | 3 ++- 5 files changed, 7 insertions(+), 32 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl index 2df05fbf210..55bcbb75588 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_shearlayer_amr.jl @@ -55,7 +55,7 @@ semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabol ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 0.7) +tspan = (0.0, 1.0) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -66,6 +66,10 @@ 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, diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index cef5a268c82..9204989e8be 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/src/equations/compressible_euler_1d.jl @@ -973,11 +973,6 @@ end return rho_times_p end -@inline function v1(u, equations::CompressibleEulerEquations1D) - rho, rho_v1, _ = u - return rho_v1 / rho -end - # Calculate thermodynamic entropy for a conservative state `cons` @inline function entropy_thermodynamic(cons, equations::CompressibleEulerEquations1D) # Pressure diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl index d61bba8e6a6..27b92f41953 100644 --- a/src/equations/compressible_euler_2d.jl +++ b/src/equations/compressible_euler_2d.jl @@ -1428,16 +1428,6 @@ end return rho_times_p end -@inline function v1(u, equations::CompressibleEulerEquations2D) - rho, rho_v1, _, _ = u - return rho_v1 / rho -end - -@inline function v2(u, equations::CompressibleEulerEquations2D) - rho, _, rho_v2, _ = u - return rho_v2 / rho -end - # Calculates the entropy flux in direction "orientation" and the entropy variables for a state cons # NOTE: This method seems to work currently (b82534e) but is never used anywhere. Thus it is # commented here until someone uses it or writes a test for it. diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl index af746c37025..7f25bde31fd 100644 --- a/src/equations/compressible_euler_3d.jl +++ b/src/equations/compressible_euler_3d.jl @@ -1514,21 +1514,6 @@ end return rho_times_p end -@inline function v1(u, equations::CompressibleEulerEquations3D) - rho, rho_v1, _, _, _ = u - return rho_v1 / rho -end - -@inline function v2(u, equations::CompressibleEulerEquations3D) - rho, _, rho_v2, _, _ = u - return rho_v2 / rho -end - -@inline function v3(u, equations::CompressibleEulerEquations3D) - rho, _, _, rho_v3, _ = u - return rho_v3 / rho -end - # Calculate thermodynamic entropy for a conservative state `u` @inline function entropy_thermodynamic(u, equations::CompressibleEulerEquations3D) rho, _ = u diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 27db6c57376..a13fdf2f535 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -252,7 +252,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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] + linf = [0.03492185879198495, 1.392635891671335, 1.357551616406459, 8.713760873018146], + tspan = (0.0, 0.7) ) end From 00a827db542ed632e988e11f28848986ae2850db Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Sun, 22 Oct 2023 12:12:42 +0200 Subject: [PATCH 51/53] remove tests --- test/test_unit.jl | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 30bc9b6d97b..5c5291c2430 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -528,28 +528,6 @@ isdir(outdir) && rm(outdir, recursive=true) end end - @timed_testset "Euler test velocity functions" begin - rho, v1, v2, v3, p = 1.0, 0.1, 0.2, 0.3, 2.0 - - let equations = CompressibleEulerEquations1D(1.4) - cons_vars = prim2cons(SVector(rho,v1,p),equations) - @test v1 ≈ Trixi.v1(cons_vars, equations) - end - - let equations = CompressibleEulerEquations2D(1.4) - cons_vars = prim2cons(SVector(rho,v1,v2,p),equations) - @test v1 ≈ Trixi.v1(cons_vars, equations) - @test v2 ≈ Trixi.v2(cons_vars, equations) - end - - let equations = CompressibleEulerEquations3D(1.4) - cons_vars = prim2cons(SVector(rho,v1,v2,v3,p),equations) - @test v1 ≈ Trixi.v1(cons_vars, equations) - @test v2 ≈ Trixi.v2(cons_vars, equations) - @test v3 ≈ Trixi.v3(cons_vars, equations) - end - end - @timed_testset "Shallow water conversion between conservative/entropy variables" begin H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 From fa893cb2964313b82e7815dd9e456abf580a07ad Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Mon, 23 Oct 2023 08:30:43 +0200 Subject: [PATCH 52/53] Remove deprecated comments --- src/equations/laplace_diffusion_2d.jl | 1 - src/equations/laplace_diffusion_3d.jl | 1 - 2 files changed, 2 deletions(-) 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 index 3988ce7144b..457e742430b 100644 --- a/src/equations/laplace_diffusion_3d.jl +++ b/src/equations/laplace_diffusion_3d.jl @@ -18,7 +18,6 @@ function varnames(variable_mapping, equations_parabolic::LaplaceDiffusion3D) varnames(variable_mapping, equations_parabolic.equations_hyperbolic) end -# no orientation specified since the flux is vector-valued function flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion3D) dudx, dudy, dudz = gradients if orientation == 1 From b8a010dfd26460265de6c4026db4cf2508a98817 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 24 Oct 2023 14:13:41 +0200 Subject: [PATCH 53/53] Add news --- NEWS.md | 1 + 1 file changed, 1 insertion(+) 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