From 5370878953dbac65cb9a472368a76f915720b6a9 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 1 Aug 2024 19:54:02 +0200 Subject: [PATCH 01/31] small typo fixes --- src/meshes/structured_mesh.jl | 4 +++- src/meshes/surface_interpolant.jl | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 553aabbbc20..f816d0bda42 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -255,13 +255,15 @@ function correction_term_3d(x, y, z, faces) linear_interpolate(x, faces[1](1, z), faces[2](1, z)) + linear_interpolate(z, faces[5](x, 1), faces[6](x, 1))) - # Correction for x-terms + # Correction for z-terms c_z = linear_interpolate(z, linear_interpolate(x, faces[1](y, -1), faces[2](y, -1)) + linear_interpolate(y, faces[3](x, -1), faces[4](x, -1)), linear_interpolate(x, faces[1](y, 1), faces[2](y, 1)) + linear_interpolate(y, faces[3](x, 1), faces[4](x, 1))) + # Each of the 12 edges are counted twice above + # so we divide the correction term by one half return 0.5 * (c_x + c_y + c_z) end diff --git a/src/meshes/surface_interpolant.jl b/src/meshes/surface_interpolant.jl index 22d14e38c5c..6a7c5c7834b 100644 --- a/src/meshes/surface_interpolant.jl +++ b/src/meshes/surface_interpolant.jl @@ -52,7 +52,7 @@ function derivative_at(s, boundary_curve::CurvedSurface) y_coordinate_at_s_on_boundary_curve_prime end -# Chebysehv-Gauss-Lobatto nodes and weights for use with curved boundaries +# Chebyshev-Gauss-Lobatto nodes and weights for use with curved boundaries function chebyshev_gauss_lobatto_nodes_weights(n_nodes::Integer) # Initialize output From 47a3b9ec294d8e6424ac5f48862a02de2f8451cb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 1 Aug 2024 20:07:01 +0200 Subject: [PATCH 02/31] Interpolate and project Ju instead of u to retain discrete conservation --- src/callbacks_step/amr_dg2d.jl | 339 ++++++++++++++++++-------- src/callbacks_step/amr_dg3d.jl | 425 ++++++++++++++++++++++++--------- 2 files changed, 550 insertions(+), 214 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 94524b23a3a..d2ebe45bc21 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -77,6 +77,9 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations, end # Refine elements in the DG solver based on a list of cell_ids that should be refined +# If an element refines the solution scaled by the Jacobian `J*u` is interpolated from +# the parent element into the four children elements. The solution on each child +# element is then recovered by dividing by the new element Jacobians. function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do @@ -96,35 +99,79 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # Retain current solution data old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) 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) + GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / + old_inverse_jacobian[i, j, + old_element]) + end + end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, 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 - 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 + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, 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 + 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) + # Before `element_id` is incremented, divide off by the new Jacobians and save + # the result again in the appropriate places + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + element_id]) + u[v, i, j, element_id + 1] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + element_id + 1]) + u[v, i, j, element_id + 2] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + element_id + 2]) + u[v, i, j, element_id + 3] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + element_id + 3]) + end + end + element_id += 2^ndims(mesh) + else + # Copy old element data to new element container and remove Jacobian scaling + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * + old_inverse_jacobian[i, j, + old_element_id]) + end + end + element_id += 1 + end 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))" + # 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))" + end # GC.@preserve old_inverse_jacobian end # GC.@preserve old_u_ode # Sanity check @@ -229,6 +276,9 @@ function refine_element!(u::AbstractArray{<:Any, 4}, element_id, end # Coarsen elements in the DG solver based on a list of cell_ids that should be removed +# If an element coarsens the solution scaled by the Jacobian `J*u` is projected from +# the four children elements back onto the parent element. The solution on the parent +# element is then recovered by dividing by the new element Jacobian. function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DGSEM, cache, elements_to_remove) @@ -246,47 +296,79 @@ function coarsen!(u_ode::AbstractVector, adaptor, to_be_removed = falses(nelements(dg, cache)) to_be_removed[elements_to_remove] .= true - # Retain current solution data + # Retain current solution data and Jacobians old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) 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) + GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / + old_inverse_jacobian[i, j, + old_element]) + end + end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, 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 + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, 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 - 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 + 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) + # Before `element_id` is incremented, divide off by the new Jacobian and save + # the result again in the appropriate place + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] *= (4 * + cache.elements.inverse_jacobian[i, + j, + element_id]) + end + end + element_id += 1 + skip = 2^ndims(mesh) - 1 + else + # Copy old element data to new element container and remove Jacobian scaling + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * + old_inverse_jacobian[i, j, + old_element_id]) + end + end + element_id += 1 + end 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))" + # 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))" + end # GC.@preserve old_inverse_jacobian end # GC.@preserve old_u_ode # Sanity check @@ -404,50 +486,109 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # Note: This is true for `quads`. T8_CHILDREN = 4 - # Retain current solution data. + # Retain current solution and inverse Jacobian data. old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - while old_index <= old_nelems && new_index <= new_nelems - if difference[old_index] > 0 # Refine. - - # Refine element and store solution directly in new data structure. - refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg) - - old_index += 1 - new_index += T8_CHILDREN - - elseif difference[old_index] < 0 # Coarsen. - - # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< 0) "bad cell/element order" - - # Coarsen elements and store solution directly in new data structure. - coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, - dg) - - old_index += T8_CHILDREN - new_index += 1 - - else # No changes. + GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to interpolation or projection + for old_element in 1:old_nelems + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / + old_inverse_jacobian[i, j, + old_element]) + end + end + end - # Copy old element data to new element container. - @views u[:, .., new_index] .= old_u[:, .., old_index] + reinitialize_containers!(mesh, equations, dg, cache) - old_index += 1 - new_index += 1 - end - end # while + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + while old_index <= old_nelems && new_index <= new_nelems + if difference[old_index] > 0 # Refine. + + # Refine element and store solution directly in new data structure. + refine_element!(u, new_index, old_u, old_index, adaptor, equations, + dg) + + # Before indices are incremented divide off by the new Jacobians and save + # the result again in the appropriate places + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, new_index] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + new_index]) + u[v, i, j, new_index + 1] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + new_index + 1]) + u[v, i, j, new_index + 2] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + new_index + 2]) + u[v, i, j, new_index + 3] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + new_index + 3]) + end + end + + old_index += 1 + new_index += T8_CHILDREN + + elseif difference[old_index] < 0 # Coarsen. + + # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< + 0) "bad cell/element order" + + # Coarsen elements and store solution directly in new data structure. + coarsen_elements!(u, new_index, old_u, old_index, adaptor, + equations, + dg) + + # Before the indices are incremented divide off by the new Jacobian and save + # the result again in the appropriate place + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, new_index] *= (4 * + cache.elements.inverse_jacobian[i, + j, + new_index]) + end + end + + old_index += T8_CHILDREN + new_index += 1 + + else # No changes. + + # Copy old element data to new element container and remove Jacobian scaling + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, new_index] = (old_u[v, i, j, old_index] * + old_inverse_jacobian[i, j, + old_index]) + end + end + + old_index += 1 + new_index += 1 + end + end # while + end # GC.@preserve old_inverse_jacobian end # GC.@preserve old_u_ode return nothing diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 3f67951bafe..32b41319f48 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -6,6 +6,9 @@ #! format: noindent # Refine elements in the DG solver based on a list of cell_ids that should be refined +# If an element refines the solution scaled by the Jacobian `J*u` is interpolated from +# the parent element into the eight children elements. The solution on each child +# element is then recovered by dividing by the new element Jacobians. function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, elements_to_refine) @@ -26,39 +29,109 @@ function refine!(u_ode::AbstractVector, adaptor, # Retain current solution data old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) 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) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - # Loop over all elements in old container and either copy them or refine them - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - 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, u_tmp1, u_tmp2) - 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 + GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, + old_element] / + old_inverse_jacobian[i, j, k, + old_element]) + end + end 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))" + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and either copy them or refine them + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + 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, u_tmp1, u_tmp2) + # Before `element_id` is incremented, divide off by the new Jacobians and save + # the result again in the appropriate places + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id]) + u[v, i, j, k, element_id + 1] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 1]) + u[v, i, j, k, element_id + 2] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 2]) + u[v, i, j, k, element_id + 3] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 3]) + u[v, i, j, k, element_id + 4] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 4]) + u[v, i, j, k, element_id + 5] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 5]) + u[v, i, j, k, element_id + 6] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 6]) + u[v, i, j, k, element_id + 7] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + 7]) + end + end + element_id += 2^ndims(mesh) + else + # Copy old element data to new element container and remove Jacobian scaling + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] = (old_u[v, i, j, k, + old_element_id] * + old_inverse_jacobian[i, j, k, + old_element_id]) + end + end + 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))" + end # GC.@preserve old_inverse_jacobian end # GC.@preserve old_u_ode # Sanity check @@ -146,6 +219,9 @@ function refine_element!(u::AbstractArray{<:Any, 5}, element_id, end # Coarsen elements in the DG solver based on a list of cell_ids that should be removed +# If an element coarsens the solution scaled by the Jacobian `J*u` is projected from +# the eight children elements back onto the parent element. The solution on the parent +# element is then recovered by dividing by the new element Jacobian. function coarsen!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, elements_to_remove) @@ -166,48 +242,83 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Retain current solution data old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) 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) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - # Loop over all elements in old container and either copy them or coarsen them - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - 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 + GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, + old_element] / + old_inverse_jacobian[i, j, k, + old_element]) + end + end 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, u_tmp1, u_tmp2) - 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 + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and either copy them or coarsen them + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + 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, u_tmp1, u_tmp2) + # Before `element_id` is incremented, divide off by the new Jacobian and save + # the result again in the appropriate place + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] *= (8 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id]) + end + end + element_id += 1 + skip = 2^ndims(mesh) - 1 + else + # Copy old element data to new element container and remove Jacobian scaling + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] = (old_u[v, i, j, k, + old_element_id] * + old_inverse_jacobian[i, j, k, + old_element_id]) + end + end + element_id += 1 + end 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))" + # 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))" + end # GC.@preserve old_inverse_jacobian end # GC.@preserve old_u_ode # Sanity check @@ -335,56 +446,140 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, # Note: This is only true for `hexs`. T8_CHILDREN = 8 - # Retain current solution data. + # Retain current solution and inverse Jacobian data. old_u_ode = copy(u_ode) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * ndofs(mesh, dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - - while old_index <= old_nelems && new_index <= new_nelems - if difference[old_index] > 0 # Refine. - - # Refine element and store solution directly in new data structure. - refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg, - u_tmp1, u_tmp2) - - old_index += 1 - new_index += T8_CHILDREN - - elseif difference[old_index] < 0 # Coarsen. - - # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< 0) "bad cell/element order" - - # Coarsen elements and store solution directly in new data structure. - coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, - dg, u_tmp1, u_tmp2) - - old_index += T8_CHILDREN - new_index += 1 - - else # No changes. + GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed + GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to interpolation or projection + for old_element in 1:old_nelems + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, + old_element] / + old_inverse_jacobian[i, j, k, + old_element]) + end + end + end - # Copy old element data to new element container. - @views u[:, .., new_index] .= old_u[:, .., old_index] + reinitialize_containers!(mesh, equations, dg, cache) - old_index += 1 - new_index += 1 - end - end # while + resize!(u_ode, + nvariables(equations) * ndofs(mesh, dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + + while old_index <= old_nelems && new_index <= new_nelems + if difference[old_index] > 0 # Refine. + + # Refine element and store solution directly in new data structure. + refine_element!(u, new_index, old_u, old_index, adaptor, equations, + dg, + u_tmp1, u_tmp2) + + # Before indices are incremented divide off by the new Jacobians and save + # the result again in the appropriate places + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, new_index] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index]) + u[v, i, j, k, new_index + 1] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 1]) + u[v, i, j, k, new_index + 2] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 2]) + u[v, i, j, k, new_index + 3] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 3]) + u[v, i, j, k, new_index + 4] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 4]) + u[v, i, j, k, new_index + 5] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 5]) + u[v, i, j, k, new_index + 6] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 6]) + u[v, i, j, k, new_index + 7] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + 7]) + end + end + + old_index += 1 + new_index += T8_CHILDREN + + elseif difference[old_index] < 0 # Coarsen. + + # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< + 0) "bad cell/element order" + + # Coarsen elements and store solution directly in new data structure. + coarsen_elements!(u, new_index, old_u, old_index, adaptor, + equations, + dg, u_tmp1, u_tmp2) + + # Before the indices are incremented divide off by the new Jacobian and save + # the result again in the appropriate place + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, new_index] *= (8 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index]) + end + end + + old_index += T8_CHILDREN + new_index += 1 + + else # No changes. + + # Copy old element data to new element container and remove Jacobian scaling + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, new_index] = (old_u[v, i, j, k, old_index] * + old_inverse_jacobian[i, j, k, + old_index]) + end + end + + old_index += 1 + new_index += 1 + end + end # while + end # GC.@preserve old_inverse_jacobian end # GC.@preserve old_u_ode return nothing From 8f4b373f9adb5bf56641df64976d40e325328a6c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 1 Aug 2024 20:56:42 +0200 Subject: [PATCH 03/31] add specialized routines for TreeMesh and P4estMesh. Hopefully some tests pass now --- src/callbacks_step/amr_dg2d.jl | 137 ++++++++++++++++++++++++++++++- src/callbacks_step/amr_dg3d.jl | 143 ++++++++++++++++++++++++++++++++- 2 files changed, 275 insertions(+), 5 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index d2ebe45bc21..3bf0f13dbb5 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -76,11 +76,73 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations, end # GC.@preserve old_u_ode end +# Refine elements in the DG solver based on a list of cell_ids that should be refined +function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, + equations, dg::DGSEM, cache, elements_to_refine) + # Return early if there is nothing to do + if isempty(elements_to_refine) + if mpi_isparallel() + # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do + # locally (there still might be other MPI ranks that have refined elements) + reinitialize_containers!(mesh, equations, dg, cache) + end + 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) + old_inverse_jacobian = copy(cache.elements.inverse_jacobian) + 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) + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, 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 + 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))" + end # GC.@preserve old_u_ode + + # Sanity check + if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && + !mpi_isparallel() + @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 + # Refine elements in the DG solver based on a list of cell_ids that should be refined # If an element refines the solution scaled by the Jacobian `J*u` is interpolated from # the parent element into the four children elements. The solution on each child # element is then recovered by dividing by the new element Jacobians. -function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}}, +function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -275,12 +337,81 @@ function refine_element!(u::AbstractArray{<:Any, 4}, element_id, return nothing end +# Coarsen elements in the DG solver based on a list of cell_ids that should be removed +function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, + equations, dg::DGSEM, cache, elements_to_remove) + # Return early if there is nothing to do + if isempty(elements_to_remove) + if mpi_isparallel() + # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do + # locally (there still might be other MPI ranks that have coarsened elements) + reinitialize_containers!(mesh, equations, dg, cache) + end + 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 and Jacobians + 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) + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, 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 + + 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))" + end # GC.@preserve old_u_ode + + # Sanity check + if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && + !mpi_isparallel() + @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 + # Coarsen elements in the DG solver based on a list of cell_ids that should be removed # If an element coarsens the solution scaled by the Jacobian `J*u` is projected from # the four children elements back onto the parent element. The solution on the parent # element is then recovered by dividing by the new element Jacobian. -function coarsen!(u_ode::AbstractVector, adaptor, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, +function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, equations, dg::DGSEM, cache, elements_to_remove) # Return early if there is nothing to do if isempty(elements_to_remove) diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 32b41319f48..ebf887af884 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -9,8 +9,74 @@ # If an element refines the solution scaled by the Jacobian `J*u` is interpolated from # the parent element into the eight children elements. The solution on each child # element is then recovered by dividing by the new element Jacobians. -function refine!(u_ode::AbstractVector, adaptor, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, +function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, + equations, dg::DGSEM, cache, elements_to_refine) + # Return early if there is nothing to do + if isempty(elements_to_refine) + if mpi_isparallel() + # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do + # locally (there still might be other MPI ranks that have refined elements) + reinitialize_containers!(mesh, equations, dg, cache) + end + 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) + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and either copy them or refine them + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + 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, u_tmp1, u_tmp2) + 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))" + end # GC.@preserve old_u_ode + + # Sanity check + if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 + @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 + +# Refine elements in the DG solver based on a list of cell_ids that should be refined +# If an element refines the solution scaled by the Jacobian `J*u` is interpolated from +# the parent element into the eight children elements. The solution on each child +# element is then recovered by dividing by the new element Jacobians. +function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -218,6 +284,79 @@ function refine_element!(u::AbstractArray{<:Any, 5}, element_id, return nothing end +# Coarsen elements in the DG solver based on a list of cell_ids that should be removed +function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, + equations, dg::DGSEM, cache, elements_to_remove) + # Return early if there is nothing to do + if isempty(elements_to_remove) + if mpi_isparallel() + # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do + # locally (there still might be other MPI ranks that have coarsened elements) + reinitialize_containers!(mesh, equations, dg, cache) + end + 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) + + reinitialize_containers!(mesh, equations, dg, cache) + + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) + + # Loop over all elements in old container and either copy them or coarsen them + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + 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, u_tmp1, u_tmp2) + 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))" + end # GC.@preserve old_u_ode + + # Sanity check + if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 + @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 + # Coarsen elements in the DG solver based on a list of cell_ids that should be removed # If an element coarsens the solution scaled by the Jacobian `J*u` is projected from # the eight children elements back onto the parent element. The solution on the parent From 2ec952ae8efd36eb68c930505a9daef4d808a660 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 1 Aug 2024 20:58:47 +0200 Subject: [PATCH 04/31] remove unused copy --- src/callbacks_step/amr_dg2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 3bf0f13dbb5..ff60dc25f96 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -96,7 +96,6 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, # Retain current solution data old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) - old_inverse_jacobian = copy(cache.elements.inverse_jacobian) 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) From 109eda2e02f185a6af682925c74c77a18de148ee Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 1 Aug 2024 21:00:37 +0200 Subject: [PATCH 05/31] remove false comment --- src/callbacks_step/amr_dg2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index ff60dc25f96..b77d3b4b0b4 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -353,7 +353,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, to_be_removed = falses(nelements(dg, cache)) to_be_removed[elements_to_remove] .= true - # Retain current solution data and Jacobians + # 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 From c11224e805b29ce7e59f9fa36788e2d32c718811 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 09:34:16 +0200 Subject: [PATCH 06/31] simplify the method on the children --- src/callbacks_step/amr_dg2d.jl | 68 +++++++------------ src/callbacks_step/amr_dg3d.jl | 120 ++++++++------------------------- 2 files changed, 51 insertions(+), 137 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index b77d3b4b0b4..e19f418769c 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -191,26 +191,16 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg) - # Before `element_id` is incremented, divide off by the new Jacobians and save - # the result again in the appropriate places - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, element_id] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - element_id]) - u[v, i, j, element_id + 1] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - element_id + 1]) - u[v, i, j, element_id + 2] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - element_id + 2]) - u[v, i, j, element_id + 3] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - element_id + 3]) + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id + m] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + element_id + m]) + end end end element_id += 2^ndims(mesh) @@ -472,8 +462,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg) - # Before `element_id` is incremented, divide off by the new Jacobian and save - # the result again in the appropriate place + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) u[v, i, j, element_id] *= (4 * @@ -650,26 +640,16 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg) - # Before indices are incremented divide off by the new Jacobians and save - # the result again in the appropriate places - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, new_index] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - new_index]) - u[v, i, j, new_index + 1] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - new_index + 1]) - u[v, i, j, new_index + 2] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - new_index + 2]) - u[v, i, j, new_index + 3] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - new_index + 3]) + # Before indices are incremented divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, new_index + m] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + new_index + m]) + end end end @@ -689,8 +669,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, equations, dg) - # Before the indices are incremented divide off by the new Jacobian and save - # the result again in the appropriate place + # Before the indices are incremented divide by the new Jacobian and save + # the result again in the parent element for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) u[v, i, j, new_index] *= (4 * diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index ebf887af884..4d332847b12 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -131,50 +131,17 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) - # Before `element_id` is incremented, divide off by the new Jacobians and save - # the result again in the appropriate places - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, element_id] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id]) - u[v, i, j, k, element_id + 1] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 1]) - u[v, i, j, k, element_id + 2] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 2]) - u[v, i, j, k, element_id + 3] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 3]) - u[v, i, j, k, element_id + 4] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 4]) - u[v, i, j, k, element_id + 5] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 5]) - u[v, i, j, k, element_id + 6] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 6]) - u[v, i, j, k, element_id + 7] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + 7]) + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id + m] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + m]) + end end end element_id += 2^ndims(mesh) @@ -378,7 +345,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, to_be_removed = falses(nelements(dg, cache)) to_be_removed[elements_to_remove] .= true - # Retain current solution data + # Retain current solution data and Jacobians old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) @@ -429,8 +396,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) - # Before `element_id` is incremented, divide off by the new Jacobian and save - # the result again in the appropriate place + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) u[v, i, j, k, element_id] *= (8 * @@ -625,50 +592,17 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, dg, u_tmp1, u_tmp2) - # Before indices are incremented divide off by the new Jacobians and save - # the result again in the appropriate places - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, new_index] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index]) - u[v, i, j, k, new_index + 1] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 1]) - u[v, i, j, k, new_index + 2] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 2]) - u[v, i, j, k, new_index + 3] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 3]) - u[v, i, j, k, new_index + 4] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 4]) - u[v, i, j, k, new_index + 5] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 5]) - u[v, i, j, k, new_index + 6] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 6]) - u[v, i, j, k, new_index + 7] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + 7]) + # Before indices are incremented divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, new_index + m] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + new_index + m]) + end end end @@ -688,8 +622,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, equations, dg, u_tmp1, u_tmp2) - # Before the indices are incremented divide off by the new Jacobian and save - # the result again in the appropriate place + # Before the indices are incremented divide by the new Jacobian and save + # the result again in the parent element for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) u[v, i, j, k, new_index] *= (8 * From 101013a655bc5f2455ae0712f5e6616ac9ab1b3f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 12:05:25 +0200 Subject: [PATCH 07/31] update 2d p4est tests --- test/test_p4est_2d.jl | 92 +++++++++++++++++++++---------------------- 1 file changed, 46 insertions(+), 46 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ef5ac2c9de3..ee7247f95cc 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -78,8 +78,8 @@ end @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.0012766060609964525], - linf=[0.01750280631586159], + l2=[0.0012808538770535593], + linf=[0.01752690016659812], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -301,16 +301,16 @@ end @trixi_testset "elixir_euler_blast_wave_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"), l2=[ - 6.32183914e-01, - 3.86914231e-01, - 3.86869171e-01, - 1.06575688e+00, + 0.6321850210104147, + 0.38691446170269167, + 0.3868695626809587, + 1.0657553825683956, ], linf=[ - 2.76020890e+00, - 2.32659890e+00, - 2.32580837e+00, - 2.15778188e+00, + 2.7602280007469666, + 2.3265993814913672, + 2.3258078438689673, + 2.1577683028925416, ], tspan=(0.0, 0.3), coverage_override=(maxiters = 6,)) @@ -327,16 +327,16 @@ end @trixi_testset "elixir_euler_wall_bc_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_wall_bc_amr.jl"), l2=[ - 0.020291447969983396, - 0.017479614254319948, - 0.011387644425613437, - 0.0514420126021293, + 0.02026685991647352, + 0.017467584076280237, + 0.011378371604813321, + 0.05138942558296091, ], linf=[ - 0.3582779022370579, - 0.32073537890751663, - 0.221818049107692, - 0.9209559420400415, + 0.35924402060711524, + 0.32068389566068806, + 0.2361141752119986, + 0.9289840057748628, ], tspan=(0.0, 0.15)) # Ensure that we do not have excessive memory allocations @@ -352,16 +352,16 @@ end @trixi_testset "elixir_euler_forward_step_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_forward_step_amr.jl"), l2=[ - 0.004194875320833303, - 0.003785140699353966, - 0.0013696609105790351, - 0.03265268616046424, + 0.004191480950848891, + 0.003781298410569231, + 0.0013470418422981045, + 0.03262817609394949, ], linf=[ - 2.0585399781442852, - 2.213428805506876, - 3.862362410419163, - 17.75187237459251, + 2.0581500751947113, + 2.2051301367971288, + 3.8502467979250254, + 17.750333649853616, ], tspan=(0.0, 0.0001), rtol=1.0e-7, @@ -409,16 +409,16 @@ end @trixi_testset "elixir_euler_supersonic_cylinder.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_supersonic_cylinder.jl"), l2=[ - 0.026798021911954406, - 0.05118546368109259, - 0.03206703583774831, - 0.19680026567208672, + 0.02676082999794676, + 0.05110830068968181, + 0.03205164257040607, + 0.1965981012724311, ], linf=[ - 3.653905721692421, - 4.285035711361009, - 6.8544353186357645, - 31.748244912257533, + 3.6830683476364476, + 4.284442685012427, + 6.857777546171545, + 31.749285097390576, ], tspan=(0.0, 0.001), skip_coverage=true) @@ -529,16 +529,16 @@ end @trixi_testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.4552084651735862, 0.8918048264575757, 0.832471223081887, + l2=[0.4552094211937344, 0.8918052934760807, 0.8324715234680768, 0.0, - 0.9801264164951583, 0.10475690769435382, 0.1555132409149897, + 0.9801268321975978, 0.10475722739111007, 0.15551326369033164, 0.0, - 2.0597079362763556e-5], - linf=[10.194181233788775, 18.25472397868819, 10.031307436191334, + 2.0602990858239632e-5], + linf=[10.19421969147307, 18.254409357804683, 10.031954811332596, 0.0, - 19.647239392277378, 1.3938810140985936, 1.8724965294853084, + 19.646870938371492, 1.3938679692894465, 1.8725058401937984, 0.0, - 0.0016290067532561904], + 0.0016201762010257296], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -617,12 +617,12 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), l2=[ - 5.634402680811982e-7, 6.748066107517321e-6, - 1.091879472416885e-5, 0.0006686372064029146, + 5.56114097044427e-7, 6.62284247153255e-6, + 1.0823259724601275e-5, 0.000659804574787503, ], linf=[ - 0.0021456247890772823, 0.03957142889488085, - 0.03832024233032798, 2.6628739573358495, + 0.002157589754528455, 0.039163189253511164, + 0.038386804399707625, 2.6685831417913914, ], amr_interval=1, base_level=0, med_level=1, max_level=1, @@ -652,8 +652,8 @@ end lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, semi.cache, semi) - @test isapprox(lift, 0.029076443678087403, atol = 1e-13) - @test isapprox(drag, 0.13564720009197903, atol = 1e-13) + @test isapprox(lift, 0.029094009322876882, atol = 1e-13) + @test isapprox(drag, 0.13579200776643238, atol = 1e-13) end @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin From 31de8b9f57fbe753a18a26e21e302deeabd81577 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 12:26:04 +0200 Subject: [PATCH 08/31] update 3d p4est tests --- test/test_p4est_3d.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 7483cde2752..d8780a67ea9 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -78,8 +78,8 @@ end @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[1.6236411810065552e-5], - linf=[0.0010554006923731395], + l2=[1.6163120948209677e-5], + linf=[0.0010572201890564834], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, base_level = 0, med_level = 1, max_level = 2)) @@ -542,16 +542,16 @@ end @trixi_testset "elixir_mhd_shockcapturing_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"), - l2=[0.006298541670176575, 0.0064368506652601265, - 0.007108729762852636, 0.006530420607206385, - 0.02061185869237284, 0.005562033787605515, - 0.007571716276627825, 0.005571862660453231, - 3.909755063709152e-6], - linf=[0.20904054009050665, 0.18622917151105936, - 0.2347957890323218, 0.19432508025509926, - 0.6858860133405615, 0.15172116633332622, - 0.22432820727833747, 0.16805989780225183, - 0.000535219040687628], + l2=[0.006297229188299052, 0.0064363477630573936, + 0.007109134822960387, 0.0065295379843073945, + 0.02061487028361094, 0.005561406556868266, + 0.007570747563219415, 0.005571060186624124, + 3.910359570546058e-6], + linf=[0.20904050617411984, 0.18630026905465372, + 0.23476537952044518, 0.19430178061639747, + 0.6858488631108304, 0.15169972134884624, + 0.22431157069631724, 0.16823638724229162, + 0.0005352202836463904], tspan=(0.0, 0.04), coverage_override=(maxiters = 6, initial_refinement_level = 1, base_level = 1, max_level = 2)) From 95700e8d00c39f4eee1bbe62a9973ed352e36545 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 12:31:52 +0200 Subject: [PATCH 09/31] update 2d t8code tests --- test/test_t8code_2d.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index b63d2a105ac..a0cb65ebdd9 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -121,8 +121,8 @@ end @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.001993165013217687], - linf=[0.032891018571625796], + l2=[0.002019623611753929], + linf=[0.03542375961299987], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -305,16 +305,16 @@ end # This test is identical to the one in `test_p4est_2d.jl` besides minor # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.44211360369891683, 0.8805178316216257, 0.8262710688468049, + l2=[0.44207324634847545, 0.8804644301177857, 0.8262542320669426, 0.0, - 0.9616090460973586, 0.10386643568745411, - 0.15403457366543802, 0.0, - 2.8399715649715473e-5], - linf=[10.04369305341599, 17.995640564998403, 9.576041548174265, + 0.9615023124189027, 0.10386709616755131, 0.1540308191628843, 0.0, - 19.429658884314534, 1.3821395681242314, 1.818559351543182, + 2.8350276854372125e-5], + linf=[10.04548675437385, 17.998696852394836, 9.575802136190026, 0.0, - 0.002261930217575465], + 19.431290746184473, 1.3821685018474321, 1.8186235976551453, + 0.0, + 0.002309422702635547], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From 4d580e6437ee792641344a7f0ef15efad1353198 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 12:33:47 +0200 Subject: [PATCH 10/31] update 3d t8code tests --- test/test_t8code_3d.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl index 5c452444be0..8778b769291 100644 --- a/test/test_t8code_3d.jl +++ b/test/test_t8code_3d.jl @@ -61,7 +61,7 @@ mkdir(outdir) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"), l2=[0.00253595715323843], linf=[0.016486952252155795]) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -80,7 +80,7 @@ mkdir(outdir) linf=[0.0007889950196294793], coverage_override=(maxiters = 6, initial_refinement_level = 1, base_level = 1, med_level = 2, max_level = 3)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -95,12 +95,12 @@ mkdir(outdir) @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[2.0556575425846923e-5], - linf=[0.00105682693484822], + l2=[2.0535121347526814e-5], + linf=[0.0010586603797777504], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, base_level = 0, med_level = 1, max_level = 2)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -129,7 +129,7 @@ mkdir(outdir) 0.008526972236273522, ], tspan=(0.0, 0.01)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -158,7 +158,7 @@ mkdir(outdir) 0.01562861968368434, ], tspan=(0.0, 1.0)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -186,7 +186,7 @@ mkdir(outdir) 9.412914891981927e-12, ], tspan=(0.0, 0.03)) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -214,7 +214,7 @@ mkdir(outdir) 9.592326932761353e-13, ], tspan=(0.0, 0.1), atol=5.0e-13,) - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -243,7 +243,7 @@ mkdir(outdir) ], tspan=(0.0, 0.2), coverage_override=(polydeg = 3,)) # Prevent long compile time in CI - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] @@ -273,7 +273,7 @@ mkdir(outdir) ], tspan=(0.0, 0.3), coverage_override=(polydeg = 3,)) # Prevent long compile time in CI - # Ensure that we do not have excessive memory allocations + # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let t = sol.t[end] From 1b7056359a55549b7af86655117be05624f6629b Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 12:39:52 +0200 Subject: [PATCH 11/31] update mpi test values --- test/test_mpi_p4est_2d.jl | 4 ++-- test/test_mpi_p4est_3d.jl | 4 ++-- test/test_mpi_t8code_2d.jl | 4 ++-- test/test_mpi_t8code_3d.jl | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl index 6d66bc68a26..29de4efc6d0 100644 --- a/test/test_mpi_p4est_2d.jl +++ b/test/test_mpi_p4est_2d.jl @@ -97,8 +97,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.0012766060609964525], - linf=[0.01750280631586159], + l2=[0.0012808538770535593], + linf=[0.01752690016659812], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl index cca9093ec51..4f9465b85dc 100644 --- a/test/test_mpi_p4est_3d.jl +++ b/test/test_mpi_p4est_3d.jl @@ -69,8 +69,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[1.6236411810065552e-5], - linf=[0.0010554006923731395], + l2=[1.6163120948209677e-5], + linf=[0.0010572201890564834], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, diff --git a/test/test_mpi_t8code_2d.jl b/test/test_mpi_t8code_2d.jl index 7c7fc03898c..b92afdf2b6a 100644 --- a/test/test_mpi_t8code_2d.jl +++ b/test/test_mpi_t8code_2d.jl @@ -97,8 +97,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_flag.jl"), - l2=[0.001980652042312077], - linf=[0.0328882442132265], + l2=[0.002019623611753929], + linf=[0.03542375961299987], coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_mpi_t8code_3d.jl b/test/test_mpi_t8code_3d.jl index a15690a7629..2e837f79ad8 100644 --- a/test/test_mpi_t8code_3d.jl +++ b/test/test_mpi_t8code_3d.jl @@ -69,8 +69,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_dgsem") @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"), - l2=[2.0556575425846923e-5], - linf=[0.00105682693484822], + l2=[2.0535121347526814e-5], + linf=[0.0010586603797777504], tspan=(0.0, 1.0), coverage_override=(maxiters = 6, initial_refinement_level = 0, From 1f06cd6b220f1001a9259fcaf29a1704481e3f0f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 12:48:28 +0200 Subject: [PATCH 12/31] fix formatting --- test/test_p4est_2d.jl | 3 ++- test/test_p4est_3d.jl | 8 ++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ee7247f95cc..db1d256662d 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -531,7 +531,8 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), l2=[0.4552094211937344, 0.8918052934760807, 0.8324715234680768, 0.0, - 0.9801268321975978, 0.10475722739111007, 0.15551326369033164, + 0.9801268321975978, 0.10475722739111007, + 0.15551326369033164, 0.0, 2.0602990858239632e-5], linf=[10.19421969147307, 18.254409357804683, 10.031954811332596, diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index d8780a67ea9..d6d7682862f 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -548,10 +548,10 @@ end 0.007570747563219415, 0.005571060186624124, 3.910359570546058e-6], linf=[0.20904050617411984, 0.18630026905465372, - 0.23476537952044518, 0.19430178061639747, - 0.6858488631108304, 0.15169972134884624, - 0.22431157069631724, 0.16823638724229162, - 0.0005352202836463904], + 0.23476537952044518, 0.19430178061639747, + 0.6858488631108304, 0.15169972134884624, + 0.22431157069631724, 0.16823638724229162, + 0.0005352202836463904], tspan=(0.0, 0.04), coverage_override=(maxiters = 6, initial_refinement_level = 1, base_level = 1, max_level = 2)) From 93e01f192176f114981077172eb58d8273c7145a Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 2 Aug 2024 13:39:37 +0200 Subject: [PATCH 13/31] remove unnecessary comment --- src/callbacks_step/amr_dg3d.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 4d332847b12..04dd49f776d 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -6,9 +6,6 @@ #! format: noindent # Refine elements in the DG solver based on a list of cell_ids that should be refined -# If an element refines the solution scaled by the Jacobian `J*u` is interpolated from -# the parent element into the eight children elements. The solution on each child -# element is then recovered by dividing by the new element Jacobians. function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do From 9ed4f523171ce58b51708b5523ef8dcfdb68c137 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 5 Aug 2024 10:29:05 +0200 Subject: [PATCH 14/31] add new 2d tests including conservation --- .../elixir_euler_weak_blast_wave_amr.jl | 118 ++++++++++++++++++ .../elixir_euler_weak_blast_wave_amr.jl | 113 +++++++++++++++++ test/test_p4est_2d.jl | 34 +++++ test/test_t8code_2d.jl | 34 +++++ 4 files changed, 299 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl create mode 100644 examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..25f67521dca --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,118 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + r0 = 0.2 + E = 1 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +# Get the DG approximation space + +# Activate the shock capturing + flux differencing +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +############################################################################### + +# Affine type mapping to take the [-1,1]^2 domain +# and warp it as described in https://arxiv.org/abs/2012.12040 +# Warping with the coefficient 0.2 is even more extreme. +function mapping_twist(xi, eta) + y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y) + return SVector(x, y) +end + +# The mesh below can be made periodic +trees_per_dimension = (8, 8) + +# Create P4estMesh with 8 x 8 trees +mesh = P4estMesh(trees_per_dimension, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 0, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 400 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.2, + save_initial_solution = true, + save_final_solution = true) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.05, + max_level = 2, max_threshold = 0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +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);#, maxiters=4); +summary_callback() # print the timer summary diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..66e1a469660 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,113 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4) + +function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + r0 = 0.2 + E = 1 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +# Get the DG approximation space + +# Activate the shock capturing + flux differencing +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +############################################################################### + +# Affine type mapping to take the [-1,1]^2 domain +# and warp it as described in https://arxiv.org/abs/2012.12040 +# Warping with the coefficient 0.2 is even more extreme. +function mapping_twist(xi, eta) + y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y) + return SVector(x, y) +end + +# The mesh below can be made periodic +trees_per_dimension = (8, 8) + +# Create T8codeMesh with 8 x 8 trees +mesh = T8codeMesh(trees_per_dimension, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 0, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.2) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 400 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.05, + max_level = 2, max_threshold = 0.1) +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + 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);#, maxiters=4); +summary_callback() # print the timer summary diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index db1d256662d..f56dc9cd769 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -685,6 +685,40 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.11134260363848127, + 0.11752357091804219, + 0.11829112104640764, + 0.7557891142955036, + ], + linf=[ + 0.5728647031475109, + 0.8353132977670252, + 0.8266797080712205, + 3.9792506230548317, + ], + tspan=(0.0, 0.1), + coverage_override=(maxiters = 6,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index a0cb65ebdd9..c1fcc355218 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -325,6 +325,40 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.10823279736983638, + 0.1158152939803735, + 0.11633970342992006, + 0.751152651902375, + ], + linf=[ + 0.5581611332828653, + 0.8354026029724041, + 0.834485181423738, + 3.923553028014343, + ], + tspan=(0.0, 0.1), + coverage_override=(maxiters = 6,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory From 455fb3463a0220081efa1342549016ac62795853 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 5 Aug 2024 10:42:24 +0200 Subject: [PATCH 15/31] add new 3d tests including conservation --- .../elixir_euler_weak_blast_wave_amr.jl | 114 ++++++++++++++++++ .../elixir_euler_weak_blast_wave_amr.jl | 112 +++++++++++++++++ test/test_p4est_3d.jl | 37 ++++++ test/test_t8code_3d.jl | 37 ++++++ 4 files changed, 300 insertions(+) create mode 100644 examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl create mode 100644 examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl diff --git a/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..2debd80cc0c --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,114 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + r0 = 0.2 + E = 1.0 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements +trees_per_dimension = (4, 4, 4) + +# Affine type mapping to take the [-1,1]^3 domain +# and warp it as described in https://arxiv.org/abs/2012.12040 +function mapping_twist(xi, eta, zeta) + y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta)) + + x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta)) + + z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta)) + + return SVector(x, y, z) +end + +mesh = P4estMesh(trees_per_dimension, + polydeg = 2, + mapping = mapping_twist, + initial_refinement_level = 1, + periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 2, med_threshold = 0.05, + max_level = 3, max_threshold = 0.15) +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + 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 diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl new file mode 100644 index 00000000000..a846c896aae --- /dev/null +++ b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -0,0 +1,112 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + r0 = 0.2 + E = 1.0 + p0_inner = 3 + p0_outer = 1 + + # Calculate primitive variables + rho = 1.1 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0, + alpha_min = 0.001, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements +trees_per_dimension = (4, 4, 4) + +function mapping_twist(xi, eta, zeta) + y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta)) + + x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta)) + + z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta)) + + return SVector(x, y, z) +end + +mesh = T8codeMesh(trees_per_dimension, + polydeg = 2, + mapping = mapping_twist, + initial_refinement_level = 1, + periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 1, + med_level = 2, med_threshold = 0.05, + max_level = 3, max_threshold = 0.15) +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + 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 diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index d6d7682862f..7b69d1c0cf2 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -586,6 +586,43 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.011345993108796831, + 0.018525073963833696, + 0.019102348105917946, + 0.01920515438943838, + 0.15060493968460148, + ], + linf=[ + 0.2994949779783401, + 0.5530175050084679, + 0.5335803757792128, + 0.5647252867336123, + 3.6462732329242566, + ], + tspan=(0.0, 0.025), + coverage_override=(maxiters = 6,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) + @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl index 8778b769291..940d2c43372 100644 --- a/test/test_t8code_3d.jl +++ b/test/test_t8code_3d.jl @@ -316,6 +316,43 @@ mkdir(outdir) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"), + l2=[ + 0.010014531529951328, + 0.0176268986746271, + 0.01817514447099777, + 0.018271085903740675, + 0.15193033077438198, + ], + linf=[ + 0.2898958869606375, + 0.529717119064458, + 0.5567193302705906, + 0.570663236219957, + 3.5496520808512027, + ], + tspan=(0.0, 0.025), + coverage_override=(maxiters = 6,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + # Check for conservation + state_integrals = Trixi.integrate(sol.u[2], semi) + initial_state_integrals = analysis_callback.affect!.initial_state_integrals + + @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13) + @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13) + @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) + @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) + @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13) + end end # Clean up afterwards: delete Trixi.jl output directory From 7b251fe563c4a76a55879b40267da9d4d1b5960c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 7 Aug 2024 13:34:31 +0200 Subject: [PATCH 16/31] update t8code 2d mpi test --- .../t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl | 3 ++- test/test_mpi_t8code_2d.jl | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl index 9138586cccf..d8893854811 100644 --- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl +++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl @@ -61,7 +61,8 @@ amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first) amr_callback = AMRCallback(semi, amr_controller, interval = 5, adapt_initial_condition = true, - adapt_initial_condition_only_refine = true) + adapt_initial_condition_only_refine = true, + dynamic_load_balancing = true) stepsize_callback = StepsizeCallback(cfl = 0.7) diff --git a/test/test_mpi_t8code_2d.jl b/test/test_mpi_t8code_2d.jl index b92afdf2b6a..75e65c8c380 100644 --- a/test/test_mpi_t8code_2d.jl +++ b/test/test_mpi_t8code_2d.jl @@ -99,6 +99,7 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_dgsem") "elixir_advection_amr_unstructured_flag.jl"), l2=[0.002019623611753929], linf=[0.03542375961299987], + dynamic_load_balancing=false, coverage_override=(maxiters = 6,)) # Ensure that we do not have excessive memory allocations From 6eb2f64c0ff2eb41993543bd49f300dee33b950e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 8 Aug 2024 07:29:24 +0200 Subject: [PATCH 17/31] remove unnecessary Union --- src/callbacks_step/amr_dg3d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 04dd49f776d..7581db4f953 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -325,8 +325,7 @@ end # If an element coarsens the solution scaled by the Jacobian `J*u` is projected from # the eight children elements back onto the parent element. The solution on the parent # element is then recovered by dividing by the new element Jacobian. -function coarsen!(u_ode::AbstractVector, adaptor, - mesh::Union{TreeMesh{3}, P4estMesh{3}}, +function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, equations, dg::DGSEM, cache, elements_to_remove) # Return early if there is nothing to do if isempty(elements_to_remove) From 0eea3769018e18969dcff326e1152811a16e155e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 8 Aug 2024 08:22:43 +0200 Subject: [PATCH 18/31] add comments regarding how element ids are incremented during refine or coarsen --- src/callbacks_step/amr_dg2d.jl | 19 +++++++++++++++++++ src/callbacks_step/amr_dg3d.jl | 19 +++++++++++++++++++ 2 files changed, 38 insertions(+) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index e19f418769c..391a3dd5576 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -113,10 +113,13 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg) + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 4 in 2D element_id += 2^ndims(mesh) else # Copy old element data to new element container @views u[:, .., element_id] .= old_u[:, .., old_element_id] + # No refinement occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -203,6 +206,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, end end end + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 4 in 2D element_id += 2^ndims(mesh) else # Copy old element data to new element container and remove Jacobian scaling @@ -213,6 +218,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, old_element_id]) end end + # No refinement occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -375,11 +381,14 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg) + # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D + # because 4 children elements become 1 parent element 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] + # No coarsening occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -472,6 +481,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, element_id]) end end + # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D + # because 4 children elements become 1 parent element element_id += 1 skip = 2^ndims(mesh) - 1 else @@ -483,6 +494,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, old_element_id]) end end + # No coarsening occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -653,6 +665,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, end end + # Increment `old_index` on the original mesh and the `new_index` + # on the refined mesh with the number of children, i.e., T8_CHILDREN = 4 old_index += 1 new_index += T8_CHILDREN @@ -680,6 +694,9 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, end end + # Increment `old_index` on the original mesh with the number of children + # (T8_CHILDREN = 4 in 2D) and the `new_index` by one for the single + # coarsened element old_index += T8_CHILDREN new_index += 1 @@ -694,6 +711,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, end end + # No refinement / coarsening occured, so increment element index + # on each mesh by one old_index += 1 new_index += 1 end diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 7581db4f953..6d4f882c9c5 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -46,10 +46,13 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 8 in 2D element_id += 2^ndims(mesh) else # Copy old element data to new element container @views u[:, .., element_id] .= old_u[:, .., old_element_id] + # No refinement occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -141,6 +144,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, end end end + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 8 in 2D element_id += 2^ndims(mesh) else # Copy old element data to new element container and remove Jacobian scaling @@ -152,6 +157,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, old_element_id]) end end + # No refinement occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -301,11 +307,14 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) + # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D + # because 8 children elements become 1 parent element 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] + # No coarsening occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -403,6 +412,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, element_id]) end end + # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D + # because 8 children elements become 1 parent element element_id += 1 skip = 2^ndims(mesh) - 1 else @@ -415,6 +426,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, old_element_id]) end end + # No coarsening occured, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -602,6 +614,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, end end + # Increment `old_index` on the original mesh and the `new_index` + # on the refined mesh with the number of children, i.e., T8_CHILDREN = 8 old_index += 1 new_index += T8_CHILDREN @@ -630,6 +644,9 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, end end + # Increment `old_index` on the original mesh with the number of children + # (T8_CHILDREN = 8 in 3D) and the `new_index` by one for the single + # coarsened element old_index += T8_CHILDREN new_index += 1 @@ -644,6 +661,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, end end + # No refinement / coarsening occured, so increment element index + # on each mesh by one old_index += 1 new_index += 1 end From bde64ee8f0bb04f63cfa7c280f5aeed2e55e9a06 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 8 Aug 2024 08:25:03 +0200 Subject: [PATCH 19/31] typo fixes --- src/callbacks_step/amr_dg2d.jl | 10 +++++----- src/callbacks_step/amr_dg3d.jl | 10 +++++----- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 391a3dd5576..a53946b532a 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -119,7 +119,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, else # Copy old element data to new element container @views u[:, .., element_id] .= old_u[:, .., old_element_id] - # No refinement occured, so increment `element_id` on the new mesh by one + # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -218,7 +218,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, old_element_id]) end end - # No refinement occured, so increment `element_id` on the new mesh by one + # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -388,7 +388,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, else # Copy old element data to new element container @views u[:, .., element_id] .= old_u[:, .., old_element_id] - # No coarsening occured, so increment `element_id` on the new mesh by one + # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -494,7 +494,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, old_element_id]) end end - # No coarsening occured, so increment `element_id` on the new mesh by one + # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -711,7 +711,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, end end - # No refinement / coarsening occured, so increment element index + # No refinement / coarsening occurred, so increment element index # on each mesh by one old_index += 1 new_index += 1 diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 6d4f882c9c5..cb72300b094 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -52,7 +52,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, else # Copy old element data to new element container @views u[:, .., element_id] .= old_u[:, .., old_element_id] - # No refinement occured, so increment `element_id` on the new mesh by one + # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -157,7 +157,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, old_element_id]) end end - # No refinement occured, so increment `element_id` on the new mesh by one + # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -314,7 +314,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, else # Copy old element data to new element container @views u[:, .., element_id] .= old_u[:, .., old_element_id] - # No coarsening occured, so increment `element_id` on the new mesh by one + # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -426,7 +426,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, old_element_id]) end end - # No coarsening occured, so increment `element_id` on the new mesh by one + # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 end end @@ -661,7 +661,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, end end - # No refinement / coarsening occured, so increment element index + # No refinement / coarsening occurred, so increment element index # on each mesh by one old_index += 1 new_index += 1 From 4f37fcb9ab1d8b9aeb22cdc0638a92f7845f4ded Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 8 Aug 2024 08:27:22 +0200 Subject: [PATCH 20/31] another typo fix --- src/callbacks_step/amr_dg3d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index cb72300b094..c1e633bd472 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -47,7 +47,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) # Increment `element_id` on the refined mesh with the number - # of children, i.e., 8 in 2D + # of children, i.e., 8 in 3D element_id += 2^ndims(mesh) else # Copy old element data to new element container @@ -145,7 +145,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, end end # Increment `element_id` on the refined mesh with the number - # of children, i.e., 8 in 2D + # of children, i.e., 8 in 3D element_id += 2^ndims(mesh) else # Copy old element data to new element container and remove Jacobian scaling From 50aa47866b8992ea5e627fba44823d85e63ce2ca Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 8 Aug 2024 08:28:42 +0200 Subject: [PATCH 21/31] fix comment --- src/meshes/structured_mesh.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index f816d0bda42..9e5c55197d8 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -263,7 +263,7 @@ function correction_term_3d(x, y, z, faces) linear_interpolate(y, faces[3](x, 1), faces[4](x, 1))) # Each of the 12 edges are counted twice above - # so we divide the correction term by one half + # so we divide the correction term by two return 0.5 * (c_x + c_y + c_z) end From 7cf1bb84b381a91438bb7289a41895153e6e4ca4 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 12 Aug 2024 06:50:10 +0200 Subject: [PATCH 22/31] Apply suggestions from code review Co-authored-by: Daniel Doehring --- src/callbacks_step/amr_dg2d.jl | 4 ++-- src/callbacks_step/amr_dg3d.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index a53946b532a..151f197b9a5 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -140,7 +140,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, return nothing end -# Refine elements in the DG solver based on a list of cell_ids that should be refined +# Refine elements in the DG solver based on a list of cell_ids that should be refined. # If an element refines the solution scaled by the Jacobian `J*u` is interpolated from # the parent element into the four children elements. The solution on each child # element is then recovered by dividing by the new element Jacobians. @@ -405,7 +405,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, return nothing end -# Coarsen elements in the DG solver based on a list of cell_ids that should be removed +# Coarsen elements in the DG solver based on a list of cell_ids that should be removed. # If an element coarsens the solution scaled by the Jacobian `J*u` is projected from # the four children elements back onto the parent element. The solution on the parent # element is then recovered by dividing by the new element Jacobian. diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index c1e633bd472..175e3fc596a 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -72,7 +72,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, return nothing end -# Refine elements in the DG solver based on a list of cell_ids that should be refined +# Refine elements in the DG solver based on a list of cell_ids that should be refined. # If an element refines the solution scaled by the Jacobian `J*u` is interpolated from # the parent element into the eight children elements. The solution on each child # element is then recovered by dividing by the new element Jacobians. @@ -330,7 +330,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, return nothing end -# Coarsen elements in the DG solver based on a list of cell_ids that should be removed +# Coarsen elements in the DG solver based on a list of cell_ids that should be removed. # If an element coarsens the solution scaled by the Jacobian `J*u` is projected from # the eight children elements back onto the parent element. The solution on the parent # element is then recovered by dividing by the new element Jacobian. From aa428991de755d8021b4c61a6905b1a3fdd7e24c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 14 Aug 2024 11:12:24 +0200 Subject: [PATCH 23/31] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 3 +-- examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 2 +- examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 4 +--- examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 2 +- 4 files changed, 4 insertions(+), 7 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl index 25f67521dca..68680673712 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -61,9 +61,8 @@ function mapping_twist(xi, eta) end # The mesh below can be made periodic -trees_per_dimension = (8, 8) - # Create P4estMesh with 8 x 8 trees +trees_per_dimension = (8, 8) mesh = P4estMesh(trees_per_dimension, polydeg = 4, mapping = mapping_twist, initial_refinement_level = 0, diff --git a/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl index 2debd80cc0c..b5b56220004 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -70,7 +70,7 @@ mesh = P4estMesh(trees_per_dimension, initial_refinement_level = 1, periodicity = true) -# Create the semi discretization object +# Create the semidiscretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) ############################################################################### diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl index 66e1a469660..fd76f0a7f1a 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -49,7 +49,6 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) -############################################################################### # Affine type mapping to take the [-1,1]^2 domain # and warp it as described in https://arxiv.org/abs/2012.12040 @@ -61,9 +60,8 @@ function mapping_twist(xi, eta) end # The mesh below can be made periodic -trees_per_dimension = (8, 8) - # Create T8codeMesh with 8 x 8 trees +trees_per_dimension = (8, 8) mesh = T8codeMesh(trees_per_dimension, polydeg = 4, mapping = mapping_twist, initial_refinement_level = 0, diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl index a846c896aae..97e066b3763 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -68,7 +68,7 @@ mesh = T8codeMesh(trees_per_dimension, initial_refinement_level = 1, periodicity = true) -# Create the semi discretization object +# Create the semidiscretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) ############################################################################### From 627a1695ceaf1f3ca59c14b70df6d27401d57463 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 15 Aug 2024 14:20:45 +0200 Subject: [PATCH 24/31] combine GC.preserve blocks --- src/callbacks_step/amr_dg2d.jl | 226 +++++++++++++++-------------- src/callbacks_step/amr_dg3d.jl | 254 ++++++++++++++++----------------- 2 files changed, 238 insertions(+), 242 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 151f197b9a5..3b9dc459b2f 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -164,71 +164,70 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / - old_inverse_jacobian[i, j, - old_element]) - end + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / + old_inverse_jacobian[i, j, + old_element]) end end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, 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 - 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) - # Before `element_id` is incremented, divide by the new Jacobians on each - # child element and save the result - for m in 0:3 # loop over the children - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, element_id + m] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - element_id + m]) - end - end - end - # Increment `element_id` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id += 2^ndims(mesh) - else - # Copy old element data to new element container and remove Jacobian scaling + # 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) + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) - u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * - old_inverse_jacobian[i, j, - old_element_id]) + u[v, i, j, element_id + m] *= (0.25 * + cache.elements.inverse_jacobian[i, + j, + element_id + m]) end end - # No refinement occurred, so increment `element_id` on the new mesh by one - element_id += 1 end + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 4 in 2D + element_id += 2^ndims(mesh) + else + # Copy old element data to new element container and remove Jacobian scaling + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * + old_inverse_jacobian[i, j, + old_element_id]) + end + end + # No refinement occurred, so increment `element_id` on the new mesh by one + element_id += 1 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))" - end # GC.@preserve old_inverse_jacobian + 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))" end # GC.@preserve old_u_ode # Sanity check @@ -429,78 +428,77 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / - old_inverse_jacobian[i, j, - old_element]) - end + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / + old_inverse_jacobian[i, j, + old_element]) end end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, 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 + # 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) - # Before `element_id` is incremented, divide by the new Jacobian and save - # the result in the parent element - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, element_id] *= (4 * - cache.elements.inverse_jacobian[i, - j, - element_id]) - 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) + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] *= (4 * + cache.elements.inverse_jacobian[i, + j, + element_id]) end - # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D - # because 4 children elements become 1 parent element - element_id += 1 - skip = 2^ndims(mesh) - 1 - else - # Copy old element data to new element container and remove Jacobian scaling - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * - old_inverse_jacobian[i, j, - old_element_id]) - end + end + # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D + # because 4 children elements become 1 parent element + element_id += 1 + skip = 2^ndims(mesh) - 1 + else + # Copy old element data to new element container and remove Jacobian scaling + for j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * + old_inverse_jacobian[i, j, + old_element_id]) end - # No coarsening occurred, so increment `element_id` on the new mesh by one - element_id += 1 end + # No coarsening occurred, so increment `element_id` on the new mesh by one + element_id += 1 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))" - end # GC.@preserve old_inverse_jacobian + 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))" end # GC.@preserve old_u_ode # Sanity check diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 175e3fc596a..6a002437664 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -96,78 +96,77 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, - old_element] / - old_inverse_jacobian[i, j, k, - old_element]) - end + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, + old_element] / + old_inverse_jacobian[i, j, k, + old_element]) end end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and either copy them or refine them - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - 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, u_tmp1, u_tmp2) - # Before `element_id` is incremented, divide by the new Jacobians on each - # child element and save the result - for m in 0:7 # loop over the children - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, element_id + m] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + m]) - end - end - end - # Increment `element_id` on the refined mesh with the number - # of children, i.e., 8 in 3D - element_id += 2^ndims(mesh) - else - # Copy old element data to new element container and remove Jacobian scaling + # Loop over all elements in old container and either copy them or refine them + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + 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, u_tmp1, u_tmp2) + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) - u[v, i, j, k, element_id] = (old_u[v, i, j, k, - old_element_id] * - old_inverse_jacobian[i, j, k, - old_element_id]) + u[v, i, j, k, element_id + m] *= (0.125 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id + m]) end end - # No refinement occurred, so increment `element_id` on the new mesh by one - element_id += 1 end + # Increment `element_id` on the refined mesh with the number + # of children, i.e., 8 in 3D + element_id += 2^ndims(mesh) + else + # Copy old element data to new element container and remove Jacobian scaling + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] = (old_u[v, i, j, k, + old_element_id] * + old_inverse_jacobian[i, j, k, + old_element_id]) + end + end + # No refinement occurred, so increment `element_id` on the new mesh by one + element_id += 1 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))" - end # GC.@preserve old_inverse_jacobian + 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))" end # GC.@preserve old_u_ode # Sanity check @@ -354,85 +353,84 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, old_n_elements = nelements(dg, cache) old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode old_inverse_jacobian begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, - old_element] / - old_inverse_jacobian[i, j, k, - old_element]) - end + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element in 1:old_n_elements + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, + old_element] / + old_inverse_jacobian[i, j, k, + old_element]) end end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and either copy them or coarsen them - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - 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 + # Loop over all elements in old container and either copy them or coarsen them + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + 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, u_tmp1, u_tmp2) - # Before `element_id` is incremented, divide by the new Jacobian and save - # the result in the parent element - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, element_id] *= (8 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id]) - 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, u_tmp1, u_tmp2) + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] *= (8 * + cache.elements.inverse_jacobian[i, + j, + k, + element_id]) end - # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D - # because 8 children elements become 1 parent element - element_id += 1 - skip = 2^ndims(mesh) - 1 - else - # Copy old element data to new element container and remove Jacobian scaling - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, element_id] = (old_u[v, i, j, k, - old_element_id] * - old_inverse_jacobian[i, j, k, - old_element_id]) - end + end + # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D + # because 8 children elements become 1 parent element + element_id += 1 + skip = 2^ndims(mesh) - 1 + else + # Copy old element data to new element container and remove Jacobian scaling + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + u[v, i, j, k, element_id] = (old_u[v, i, j, k, + old_element_id] * + old_inverse_jacobian[i, j, k, + old_element_id]) end - # No coarsening occurred, so increment `element_id` on the new mesh by one - element_id += 1 end + # No coarsening occurred, so increment `element_id` on the new mesh by one + element_id += 1 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))" - end # GC.@preserve old_inverse_jacobian + 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))" end # GC.@preserve old_u_ode # Sanity check From b7f3911c957d4fd66a67354f4ad3887ecda24b3d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 15 Aug 2024 14:22:17 +0200 Subject: [PATCH 25/31] fix formatting in elixir --- examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl index fd76f0a7f1a..52da853c6ff 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -49,7 +49,6 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) - # Affine type mapping to take the [-1,1]^2 domain # and warp it as described in https://arxiv.org/abs/2012.12040 # Warping with the coefficient 0.2 is even more extreme. From 493bb16662fca4a5398bb4f7f9065c9ca639f261 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 15 Aug 2024 15:10:39 +0200 Subject: [PATCH 26/31] simplify implementation and add if blocks for P4estMesh specialities --- src/callbacks_step/amr_dg2d.jl | 234 +++++++------------------------ src/callbacks_step/amr_dg3d.jl | 246 +++++++-------------------------- 2 files changed, 104 insertions(+), 376 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 3b9dc459b2f..6b3817ceb20 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -76,75 +76,11 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations, end # GC.@preserve old_u_ode end -# Refine elements in the DG solver based on a list of cell_ids that should be refined -function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, - equations, dg::DGSEM, cache, elements_to_refine) - # Return early if there is nothing to do - if isempty(elements_to_refine) - if mpi_isparallel() - # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do - # locally (there still might be other MPI ranks that have refined elements) - reinitialize_containers!(mesh, equations, dg, cache) - end - 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) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, 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 - 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) - # Increment `element_id` on the refined mesh with the number - # of children, i.e., 4 in 2D - element_id += 2^ndims(mesh) - else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] - # No refinement occurred, so increment `element_id` on the new mesh by one - 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))" - end # GC.@preserve old_u_ode - - # Sanity check - if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && - !mpi_isparallel() - @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 - # Refine elements in the DG solver based on a list of cell_ids that should be refined. -# If an element refines the solution scaled by the Jacobian `J*u` is interpolated from -# the parent element into the four children elements. The solution on each child +# On `P4estMesh`, if an element refines the solution scaled by the Jacobian `J*u` is interpolated +# from the parent element into the four children elements. The solution on each child # element is then recovered by dividing by the new element Jacobians. -function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, +function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -168,14 +104,14 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements for v in eachvariable(equations) - old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / - old_inverse_jacobian[i, j, - old_element]) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) end end end @@ -194,29 +130,32 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg) - # Before `element_id` is incremented, divide by the new Jacobians on each - # child element and save the result - for m in 0:3 # loop over the children - for j in eachnode(dg), i in eachnode(dg) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children for v in eachvariable(equations) - u[v, i, j, element_id + m] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - element_id + m]) + u[v, .., element_id + m] .*= (0.25 .* + cache.elements.inverse_jacobian[.., + element_id + m]) end end end + # Increment `element_id` on the refined mesh with the number # of children, i.e., 4 in 2D element_id += 2^ndims(mesh) else - # Copy old element data to new element container and remove Jacobian scaling - for j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling for v in eachvariable(equations) - u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * - old_inverse_jacobian[i, j, - old_element_id]) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] end # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 @@ -331,84 +270,12 @@ function refine_element!(u::AbstractArray{<:Any, 4}, element_id, return nothing end -# Coarsen elements in the DG solver based on a list of cell_ids that should be removed -function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{2}, - equations, dg::DGSEM, cache, elements_to_remove) - # Return early if there is nothing to do - if isempty(elements_to_remove) - if mpi_isparallel() - # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do - # locally (there still might be other MPI ranks that have coarsened elements) - reinitialize_containers!(mesh, equations, dg, cache) - end - 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) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, 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 - - 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) - # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D - # because 4 children elements become 1 parent element - 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] - # No coarsening occurred, so increment `element_id` on the new mesh by one - 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))" - end # GC.@preserve old_u_ode - - # Sanity check - if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && - !mpi_isparallel() - @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 - # Coarsen elements in the DG solver based on a list of cell_ids that should be removed. -# If an element coarsens the solution scaled by the Jacobian `J*u` is projected from -# the four children elements back onto the parent element. The solution on the parent +# On `P4estMesh`, if an element coarsens the solution scaled by the Jacobian `J*u` is projected +# from the four children elements back onto the parent element. The solution on the parent # element is then recovered by dividing by the new element Jacobian. -function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, +function coarsen!(u_ode::AbstractVector, adaptor, + mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, dg::DGSEM, cache, elements_to_remove) # Return early if there is nothing to do if isempty(elements_to_remove) @@ -432,14 +299,14 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements for v in eachvariable(equations) - old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / - old_inverse_jacobian[i, j, - old_element]) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) end end end @@ -470,28 +337,31 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg) - # Before `element_id` is incremented, divide by the new Jacobian and save - # the result in the parent element - for j in eachnode(dg), i in eachnode(dg) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element for v in eachvariable(equations) - u[v, i, j, element_id] *= (4 * - cache.elements.inverse_jacobian[i, - j, - element_id]) + u[v, .., element_id] .*= (4 .* + cache.elements.inverse_jacobian[.., + element_id]) end end + # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D # because 4 children elements become 1 parent element element_id += 1 skip = 2^ndims(mesh) - 1 else - # Copy old element data to new element container and remove Jacobian scaling - for j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling for v in eachvariable(equations) - u[v, i, j, element_id] = (old_u[v, i, j, old_element_id] * - old_inverse_jacobian[i, j, - old_element_id]) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] end # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 6a002437664..39dc266e512 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -5,78 +5,11 @@ @muladd begin #! format: noindent -# Refine elements in the DG solver based on a list of cell_ids that should be refined -function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, - equations, dg::DGSEM, cache, elements_to_refine) - # Return early if there is nothing to do - if isempty(elements_to_refine) - if mpi_isparallel() - # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do - # locally (there still might be other MPI ranks that have refined elements) - reinitialize_containers!(mesh, equations, dg, cache) - end - 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) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - # Loop over all elements in old container and either copy them or refine them - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - 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, u_tmp1, u_tmp2) - # Increment `element_id` on the refined mesh with the number - # of children, i.e., 8 in 3D - element_id += 2^ndims(mesh) - else - # Copy old element data to new element container - @views u[:, .., element_id] .= old_u[:, .., old_element_id] - # No refinement occurred, so increment `element_id` on the new mesh by one - 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))" - end # GC.@preserve old_u_ode - - # Sanity check - if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 - @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 - # Refine elements in the DG solver based on a list of cell_ids that should be refined. -# If an element refines the solution scaled by the Jacobian `J*u` is interpolated from -# the parent element into the eight children elements. The solution on each child +# On `P4estMesh`, if an element refines the solution scaled by the Jacobian `J*u` is interpolated +# from the parent element into the eight children elements. The solution on each child # element is then recovered by dividing by the new element Jacobians. -function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, +function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -100,15 +33,14 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements for v in eachvariable(equations) - old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, - old_element] / - old_inverse_jacobian[i, j, k, - old_element]) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) end end end @@ -131,31 +63,32 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, # Refine element and store solution directly in new data structure refine_element!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) - # Before `element_id` is incremented, divide by the new Jacobians on each - # child element and save the result - for m in 0:7 # loop over the children - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children for v in eachvariable(equations) - u[v, i, j, k, element_id + m] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id + m]) + u[v, .., element_id + m] .*= (0.125 .* + cache.elements.inverse_jacobian[.., + element_id + m]) end end end + # Increment `element_id` on the refined mesh with the number # of children, i.e., 8 in 3D element_id += 2^ndims(mesh) else - # Copy old element data to new element container and remove Jacobian scaling - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling for v in eachvariable(equations) - u[v, i, j, k, element_id] = (old_u[v, i, j, k, - old_element_id] * - old_inverse_jacobian[i, j, k, - old_element_id]) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] end # No refinement occurred, so increment `element_id` on the new mesh by one element_id += 1 @@ -253,87 +186,12 @@ function refine_element!(u::AbstractArray{<:Any, 5}, element_id, return nothing end -# Coarsen elements in the DG solver based on a list of cell_ids that should be removed -function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{3}, - equations, dg::DGSEM, cache, elements_to_remove) - # Return early if there is nothing to do - if isempty(elements_to_remove) - if mpi_isparallel() - # MPICache init uses all-to-all communication -> reinitialize even if there is nothing to do - # locally (there still might be other MPI ranks that have coarsened elements) - reinitialize_containers!(mesh, equations, dg, cache) - end - 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) - - reinitialize_containers!(mesh, equations, dg, cache) - - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) - - # Loop over all elements in old container and either copy them or coarsen them - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - 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, u_tmp1, u_tmp2) - # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D - # because 8 children elements become 1 parent element - 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] - # No coarsening occurred, so increment `element_id` on the new mesh by one - 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))" - end # GC.@preserve old_u_ode - - # Sanity check - if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 - @assert ninterfaces(cache.interfaces)==ndims(mesh) * nelements(dg, cache) ("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 - # Coarsen elements in the DG solver based on a list of cell_ids that should be removed. -# If an element coarsens the solution scaled by the Jacobian `J*u` is projected from -# the eight children elements back onto the parent element. The solution on the parent +# On `P4estMesh`, if an element coarsens the solution scaled by the Jacobian `J*u` is projected +# from the eight children elements back onto the parent element. The solution on the parent # element is then recovered by dividing by the new element Jacobian. -function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, +function coarsen!(u_ode::AbstractVector, adaptor, + mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg::DGSEM, cache, elements_to_remove) # Return early if there is nothing to do if isempty(elements_to_remove) @@ -357,15 +215,14 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, GC.@preserve old_u_ode old_inverse_jacobian begin old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to projection - for old_element in 1:old_n_elements - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to projection + for old_element_id in 1:old_n_elements for v in eachvariable(equations) - old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, - old_element] / - old_inverse_jacobian[i, j, k, - old_element]) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) end end end @@ -400,30 +257,31 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{3}, # Coarsen elements and store solution directly in new data structure coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, equations, dg, u_tmp1, u_tmp2) - # Before `element_id` is incremented, divide by the new Jacobian and save - # the result in the parent element - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + + if mesh isa P4estMesh + # Before `element_id` is incremented, divide by the new Jacobian and save + # the result in the parent element for v in eachvariable(equations) - u[v, i, j, k, element_id] *= (8 * - cache.elements.inverse_jacobian[i, - j, - k, - element_id]) + u[v, .., element_id] .*= (8 .* + cache.elements.inverse_jacobian[.., + element_id]) end end + # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D # because 8 children elements become 1 parent element element_id += 1 skip = 2^ndims(mesh) - 1 else - # Copy old element data to new element container and remove Jacobian scaling - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + if mesh isa P4estMesh + # Copy old element data to new element container and remove Jacobian scaling for v in eachvariable(equations) - u[v, i, j, k, element_id] = (old_u[v, i, j, k, - old_element_id] * - old_inverse_jacobian[i, j, k, - old_element_id]) + u[v, .., element_id] .= (old_u[v, .., old_element_id] .* + old_inverse_jacobian[.., + old_element_id]) end + else # isa TreeMesh + @views u[:, .., element_id] .= old_u[:, .., old_element_id] end # No coarsening occurred, so increment `element_id` on the new mesh by one element_id += 1 From c58e5677faea433e51207be1279e83344b03bc1c Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 15 Aug 2024 20:11:11 +0200 Subject: [PATCH 27/31] unify structure across P4estMesh and T8codeMesh --- src/callbacks_step/amr_dg2d.jl | 150 ++++++++++++++---------------- src/callbacks_step/amr_dg3d.jl | 161 +++++++++++++++------------------ 2 files changed, 140 insertions(+), 171 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 6b3817ceb20..2b7944aac94 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -167,7 +167,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM @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))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && @@ -369,7 +369,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, 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))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 && @@ -490,103 +490,89 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to interpolation or projection - for old_element in 1:old_nelems - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - old_u[v, i, j, old_element] = (old_u[v, i, j, old_element] / - old_inverse_jacobian[i, j, - old_element]) - end - end + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to interpolation or projection + for old_element_id in 1:old_nelems + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + resize!(u_ode, + nvariables(equations) * nnodes(dg)^ndims(mesh) * + nelements(dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) - while old_index <= old_nelems && new_index <= new_nelems - if difference[old_index] > 0 # Refine. + while old_index <= old_nelems && new_index <= new_nelems + if difference[old_index] > 0 # Refine. - # Refine element and store solution directly in new data structure. - refine_element!(u, new_index, old_u, old_index, adaptor, equations, - dg) + # Refine element and store solution directly in new data structure. + refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg) - # Before indices are incremented divide by the new Jacobians on each - # child element and save the result - for m in 0:3 # loop over the children - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, new_index + m] *= (0.25 * - cache.elements.inverse_jacobian[i, - j, - new_index + m]) - end - end + # Before indices are incremented divide by the new Jacobians on each + # child element and save the result + for m in 0:3 # loop over the children + for v in eachvariable(equations) + u[v, .., new_index + m] .*= (0.25 .* + cache.elements.inverse_jacobian[.., + new_index + m]) end + end - # Increment `old_index` on the original mesh and the `new_index` - # on the refined mesh with the number of children, i.e., T8_CHILDREN = 4 - old_index += 1 - new_index += T8_CHILDREN - - elseif difference[old_index] < 0 # Coarsen. + # Increment `old_index` on the original mesh and the `new_index` + # on the refined mesh with the number of children, i.e., T8_CHILDREN = 4 + old_index += 1 + new_index += T8_CHILDREN - # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< - 0) "bad cell/element order" + elseif difference[old_index] < 0 # Coarsen. - # Coarsen elements and store solution directly in new data structure. - coarsen_elements!(u, new_index, old_u, old_index, adaptor, - equations, - dg) + # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< + 0) "bad cell/element order" - # Before the indices are incremented divide by the new Jacobian and save - # the result again in the parent element - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, new_index] *= (4 * - cache.elements.inverse_jacobian[i, - j, - new_index]) - end - end + # Coarsen elements and store solution directly in new data structure. + coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, + dg) - # Increment `old_index` on the original mesh with the number of children - # (T8_CHILDREN = 4 in 2D) and the `new_index` by one for the single - # coarsened element - old_index += T8_CHILDREN - new_index += 1 + # Before the indices are incremented divide by the new Jacobian and save + # the result again in the parent element + for v in eachvariable(equations) + u[v, .., new_index] .*= (4 .* + cache.elements.inverse_jacobian[.., + new_index]) + end - else # No changes. + # Increment `old_index` on the original mesh with the number of children + # (T8_CHILDREN = 4 in 2D) and the `new_index` by one for the single + # coarsened element + old_index += T8_CHILDREN + new_index += 1 - # Copy old element data to new element container and remove Jacobian scaling - for j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, new_index] = (old_u[v, i, j, old_index] * - old_inverse_jacobian[i, j, - old_index]) - end - end + else # No changes. - # No refinement / coarsening occurred, so increment element index - # on each mesh by one - old_index += 1 - new_index += 1 + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., new_index] .= (old_u[v, .., old_index] .* + old_inverse_jacobian[.., old_index]) end - end # while - end # GC.@preserve old_inverse_jacobian - end # GC.@preserve old_u_ode + + # No refinement / coarsening occurred, so increment element index + # on each mesh by one + old_index += 1 + new_index += 1 + end + end # while + end # GC.@preserve old_u_ode old_inverse_jacobian return nothing end diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 39dc266e512..3e68eebccb7 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -100,7 +100,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM @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))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 @@ -289,7 +289,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, 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))" - end # GC.@preserve old_u_ode + end # GC.@preserve old_u_ode old_inverse_jacobian # Sanity check if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 @@ -420,111 +420,94 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, old_u_ode = copy(u_ode) old_inverse_jacobian = copy(cache.elements.inverse_jacobian) - GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed - GC.@preserve old_inverse_jacobian begin # OBS! If we don't GC.@preserve old_inverse_jacobian, it might be GC'ed - old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) + # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed + GC.@preserve old_u_ode begin + old_u = wrap_array(old_u_ode, mesh, equations, dg, cache) - # Loop over all elements in old container and scale the old solution by the Jacobian - # prior to interpolation or projection - for old_element in 1:old_nelems - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - old_u[v, i, j, k, old_element] = (old_u[v, i, j, k, - old_element] / - old_inverse_jacobian[i, j, k, - old_element]) - end - end + # Loop over all elements in old container and scale the old solution by the Jacobian + # prior to interpolation or projection + for old_element_id in 1:old_nelems + for v in eachvariable(equations) + old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./ + old_inverse_jacobian[.., + old_element_id]) end + end - reinitialize_containers!(mesh, equations, dg, cache) + reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * ndofs(mesh, dg, cache)) - u = wrap_array(u_ode, mesh, equations, dg, cache) + resize!(u_ode, + nvariables(equations) * ndofs(mesh, dg, cache)) + u = wrap_array(u_ode, mesh, equations, dg, cache) - u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) - u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), - nnodes(dg), nnodes(dg)) + u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) + u_tmp2 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), + nnodes(dg), nnodes(dg)) - while old_index <= old_nelems && new_index <= new_nelems - if difference[old_index] > 0 # Refine. + while old_index <= old_nelems && new_index <= new_nelems + if difference[old_index] > 0 # Refine. - # Refine element and store solution directly in new data structure. - refine_element!(u, new_index, old_u, old_index, adaptor, equations, - dg, - u_tmp1, u_tmp2) + # Refine element and store solution directly in new data structure. + refine_element!(u, new_index, old_u, old_index, adaptor, equations, + dg, u_tmp1, u_tmp2) - # Before indices are incremented divide by the new Jacobians on each - # child element and save the result - for m in 0:7 # loop over the children - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, new_index + m] *= (0.125 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index + m]) - end - end + # Before indices are incremented divide by the new Jacobians on each + # child element and save the result + for m in 0:7 # loop over the children + for v in eachvariable(equations) + u[v, .., new_index + m] .*= (0.125 .* + cache.elements.inverse_jacobian[.., + new_index + m]) end + end - # Increment `old_index` on the original mesh and the `new_index` - # on the refined mesh with the number of children, i.e., T8_CHILDREN = 8 - old_index += 1 - new_index += T8_CHILDREN - - elseif difference[old_index] < 0 # Coarsen. + # Increment `old_index` on the original mesh and the `new_index` + # on the refined mesh with the number of children, i.e., T8_CHILDREN = 8 + old_index += 1 + new_index += T8_CHILDREN - # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< - 0) "bad cell/element order" + elseif difference[old_index] < 0 # Coarsen. - # Coarsen elements and store solution directly in new data structure. - coarsen_elements!(u, new_index, old_u, old_index, adaptor, - equations, - dg, u_tmp1, u_tmp2) + # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< + 0) "bad cell/element order" - # Before the indices are incremented divide by the new Jacobian and save - # the result again in the parent element - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, new_index] *= (8 * - cache.elements.inverse_jacobian[i, - j, - k, - new_index]) - end - end + # Coarsen elements and store solution directly in new data structure. + coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, + dg, u_tmp1, u_tmp2) - # Increment `old_index` on the original mesh with the number of children - # (T8_CHILDREN = 8 in 3D) and the `new_index` by one for the single - # coarsened element - old_index += T8_CHILDREN - new_index += 1 + # Before the indices are incremented divide by the new Jacobian and save + # the result again in the parent element + for v in eachvariable(equations) + u[v, .., new_index] .*= (8 .* + cache.elements.inverse_jacobian[.., + new_index]) + end - else # No changes. + # Increment `old_index` on the original mesh with the number of children + # (T8_CHILDREN = 8 in 3D) and the `new_index` by one for the single + # coarsened element + old_index += T8_CHILDREN + new_index += 1 - # Copy old element data to new element container and remove Jacobian scaling - for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) - for v in eachvariable(equations) - u[v, i, j, k, new_index] = (old_u[v, i, j, k, old_index] * - old_inverse_jacobian[i, j, k, - old_index]) - end - end + else # No changes. - # No refinement / coarsening occurred, so increment element index - # on each mesh by one - old_index += 1 - new_index += 1 + # Copy old element data to new element container and remove Jacobian scaling + for v in eachvariable(equations) + u[v, .., new_index] .= (old_u[v, .., old_index] .* + old_inverse_jacobian[.., old_index]) end - end # while - end # GC.@preserve old_inverse_jacobian - end # GC.@preserve old_u_ode + + # No refinement / coarsening occurred, so increment element index + # on each mesh by one + old_index += 1 + new_index += 1 + end + end # while + end # GC.@preserve old_u_ode old_inverse_jacobian return nothing end From 7761bfdf7787091015c5afca522d7d616df9edc3 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 Aug 2024 08:59:04 +0200 Subject: [PATCH 28/31] fix some of the bad formatting --- src/callbacks_step/amr_dg2d.jl | 16 +++++----------- src/callbacks_step/amr_dg3d.jl | 27 +++++++++++---------------- 2 files changed, 16 insertions(+), 27 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 2b7944aac94..6fc02e12b21 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -119,8 +119,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) # Loop over all elements in old container and either copy them or refine them @@ -314,8 +313,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) # Loop over all elements in old container and either copy them or coarsen them @@ -506,9 +504,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) + resize!(u_ode, nvariables(equations) * ndofs(mesh, dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) while old_index <= old_nelems && new_index <= new_nelems @@ -537,8 +533,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< - 0) "bad cell/element order" + @assert all(difference[old_index:(old_index + T8_CHILDREN - 1)] .< 0) "bad cell/element order" # Coarsen elements and store solution directly in new data structure. coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, @@ -547,8 +542,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # Before the indices are incremented divide by the new Jacobian and save # the result again in the parent element for v in eachvariable(equations) - u[v, .., new_index] .*= (4 .* - cache.elements.inverse_jacobian[.., + u[v, .., new_index] .*= (4 .* cache.elements.inverse_jacobian[.., new_index]) end diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 3e68eebccb7..8b37be3fa5d 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -48,8 +48,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) # Loop over all elements in old container and either copy them or refine them @@ -61,8 +60,8 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM 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, u_tmp1, u_tmp2) + refine_element!(u, element_id, old_u, old_element_id, adaptor, + equations, dg, u_tmp1, u_tmp2) if mesh isa P4estMesh # Before `element_id` is incremented, divide by the new Jacobians on each @@ -230,8 +229,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, reinitialize_containers!(mesh, equations, dg, cache) resize!(u_ode, - nvariables(equations) * nnodes(dg)^ndims(mesh) * - nelements(dg, cache)) + nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) # Loop over all elements in old container and either copy them or coarsen them @@ -255,8 +253,8 @@ function coarsen!(u_ode::AbstractVector, adaptor, @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, u_tmp1, u_tmp2) + coarsen_elements!(u, element_id, old_u, old_element_id, adaptor, + equations, dg, u_tmp1, u_tmp2) if mesh isa P4estMesh # Before `element_id` is incremented, divide by the new Jacobian and save @@ -436,8 +434,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, reinitialize_containers!(mesh, equations, dg, cache) - resize!(u_ode, - nvariables(equations) * ndofs(mesh, dg, cache)) + resize!(u_ode, nvariables(equations) * ndofs(mesh, dg, cache)) u = wrap_array(u_ode, mesh, equations, dg, cache) u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg), @@ -449,8 +446,8 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, if difference[old_index] > 0 # Refine. # Refine element and store solution directly in new data structure. - refine_element!(u, new_index, old_u, old_index, adaptor, equations, - dg, u_tmp1, u_tmp2) + refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg, + u_tmp1, u_tmp2) # Before indices are incremented divide by the new Jacobians on each # child element and save the result @@ -472,8 +469,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, # 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(difference[old_index:(old_index + T8_CHILDREN - 1)] .< - 0) "bad cell/element order" + @assert all(difference[old_index:(old_index + T8_CHILDREN - 1)] .< 0) "bad cell/element order" # Coarsen elements and store solution directly in new data structure. coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations, @@ -482,8 +478,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, # Before the indices are incremented divide by the new Jacobian and save # the result again in the parent element for v in eachvariable(equations) - u[v, .., new_index] .*= (8 .* - cache.elements.inverse_jacobian[.., + u[v, .., new_index] .*= (8 .* cache.elements.inverse_jacobian[.., new_index]) end From 570cf75c9f4eaebcef6e003d9322bdb434f8be38 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 Aug 2024 11:26:22 +0200 Subject: [PATCH 29/31] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 4 ++++ examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl | 4 ++++ src/callbacks_step/amr_dg2d.jl | 4 ++-- src/callbacks_step/amr_dg3d.jl | 4 ++-- 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl index 52da853c6ff..a3366caa317 100644 --- a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -108,3 +108,7 @@ 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);#, maxiters=4); summary_callback() # print the timer summary + +# Finalize `T8codeMesh` to make sure MPI related objects in t8code are +# released before `MPI` finalizes. +!isinteractive() && finalize(mesh) diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl index 97e066b3763..106b4821144 100644 --- a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl +++ b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl @@ -110,3 +110,7 @@ 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 + +# Finalize `T8codeMesh` to make sure MPI related objects in t8code are +# released before `MPI` finalizes. +!isinteractive() && finalize(mesh) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 6fc02e12b21..062020b6d42 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -135,7 +135,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM # child element and save the result for m in 0:3 # loop over the children for v in eachvariable(equations) - u[v, .., element_id + m] .*= (0.25 .* + u[v, .., element_id + m] .*= (0.25f0 .* cache.elements.inverse_jacobian[.., element_id + m]) end @@ -517,7 +517,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations, # child element and save the result for m in 0:3 # loop over the children for v in eachvariable(equations) - u[v, .., new_index + m] .*= (0.25 .* + u[v, .., new_index + m] .*= (0.25f0 .* cache.elements.inverse_jacobian[.., new_index + m]) end diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl index 8b37be3fa5d..594c30dcca5 100644 --- a/src/callbacks_step/amr_dg3d.jl +++ b/src/callbacks_step/amr_dg3d.jl @@ -68,7 +68,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM # child element and save the result for m in 0:7 # loop over the children for v in eachvariable(equations) - u[v, .., element_id + m] .*= (0.125 .* + u[v, .., element_id + m] .*= (0.125f0 .* cache.elements.inverse_jacobian[.., element_id + m]) end @@ -453,7 +453,7 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations, # child element and save the result for m in 0:7 # loop over the children for v in eachvariable(equations) - u[v, .., new_index + m] .*= (0.125 .* + u[v, .., new_index + m] .*= (0.125f0 .* cache.elements.inverse_jacobian[.., new_index + m]) end From cc9f1d9c8e56ef31087f54c70afeefc61bd751b3 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Aug 2024 11:32:25 +0200 Subject: [PATCH 30/31] add news --- NEWS.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/NEWS.md b/NEWS.md index bca7d03e954..43bb5d8c355 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes in the v0.8 lifecycle + +#### Changed + +- The AMR routing for `P4estMesh` and `T8codeMesh` was changed to work on the product + of the Jacobian and the conserved variables instead of the conserved variables only + to make AMR fully conservative ([#2028]). This may change AMR results slightly. + ## Changes when updating to v0.8 from v0.7.x #### Added From b650fff3acb9fd7057613f9afb69b3d0225a4464 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Mon, 19 Aug 2024 12:10:28 +0200 Subject: [PATCH 31/31] Update NEWS.md Co-authored-by: Andrew Winters --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 43bb5d8c355..9371cafa07f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -8,7 +8,7 @@ for human readability. #### Changed -- The AMR routing for `P4estMesh` and `T8codeMesh` was changed to work on the product +- The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product of the Jacobian and the conserved variables instead of the conserved variables only to make AMR fully conservative ([#2028]). This may change AMR results slightly.