From 47a3b9ec294d8e6424ac5f48862a02de2f8451cb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 1 Aug 2024 20:07:01 +0200 Subject: [PATCH] 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