Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove ndims where not necessary #2081

Closed
wants to merge 3 commits into from
Closed

Remove ndims where not necessary #2081

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f7a7128
Remove `ndims` where not necessary
DanielDoehring Sep 17, 2024
70bc058
remove eval
DanielDoehring Sep 17, 2024
f17e398
Apply suggestions from code review
DanielDoehring Sep 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/callbacks_stage/positivity_zhang_shu_dg1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable,
u_mean += u_node * weights[i]
end
# note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
u_mean = u_mean / 2^ndims(mesh)
u_mean = u_mean / 2

# We compute the value directly with the mean values, as we assume that
# Jensen's inequality holds (e.g. pressure for compressible Euler equations).
Expand Down
14 changes: 8 additions & 6 deletions src/callbacks_step/amr_dg1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
@assert nelements(dg, cache) > old_n_elements

resize!(u_ode,
nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
nvariables(equations) * nnodes(dg) * 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
Expand All @@ -43,7 +43,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
# 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)
element_id += 2
else
# Copy old element data to new element container
@views u[:, .., element_id] .= old_u[:, .., old_element_id]
Expand All @@ -55,7 +55,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
# 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))"
1||element_id == nelements(dg, cache) + 2 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
end # GC.@preserve old_u_ode

# re-initialize interfaces container
Expand Down Expand Up @@ -170,7 +170,7 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
@assert nelements(dg, cache) < old_n_elements

resize!(u_ode,
nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
nvariables(equations) * nnodes(dg) * 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
Expand All @@ -187,13 +187,15 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
# 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"
@assert all(to_be_removed[old_element_id:(old_element_id + 2 - 1)]) "bad cell/element order"
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One could move skip = 1 already up here


# 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` = 1 in 1D
# because 2 children elements become 1 parent element
element_id += 1
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
skip = 2^ndims(mesh) - 1
skip = 1
else
# Copy old element data to new element container
@views u[:, .., element_id] .= old_u[:, .., old_element_id]
Expand Down
26 changes: 13 additions & 13 deletions src/callbacks_step/amr_dg2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations,
end

resize!(u_ode,
nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
nvariables(equations) * nnodes(dg)^2 * nelements(dg, cache))
u = wrap_array_native(u_ode, mesh, equations, dg, cache)

# Get new list of leaf cells
Expand Down Expand Up @@ -119,7 +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)^2 * 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
Expand All @@ -144,7 +144,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM

# Increment `element_id` on the refined mesh with the number
# of children, i.e., 4 in 2D
element_id += 2^ndims(mesh)
element_id += 4
else
if mesh isa P4estMesh
# Copy old element data to new element container and remove Jacobian scaling
Expand All @@ -165,13 +165,13 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM
# 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))"
1||element_id == nelements(dg, cache) + 2^2 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
end # GC.@preserve old_u_ode old_inverse_jacobian

# 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")
@assert ninterfaces(cache.interfaces)==2 * nelements(dg, cache) ("For 2D and periodic domains and conforming elements, the number of interfaces must be 2 times the number of elements")
end

return nothing
Expand All @@ -193,8 +193,8 @@ function refine!(u_ode::AbstractVector, adaptor,
# Sanity check
if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 &&
!mpi_isparallel()
@assert ninterfaces(cache_parabolic.interfaces)==ndims(mesh) *
nelements(dg, cache_parabolic) ("For $(ndims(mesh))D and periodic domains and conforming elements, the number of interfaces must be $(ndims(mesh)) times the number of elements")
@assert ninterfaces(cache_parabolic.interfaces)==2 *
nelements(dg, cache_parabolic) ("For 2D and periodic domains and conforming elements, the number of interfaces must be 2 times the number of elements")
end

return nothing
Expand Down Expand Up @@ -313,7 +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)^2 * nelements(dg, cache))
u = wrap_array(u_ode, mesh, equations, dg, cache)
Comment on lines 315 to 317
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think parts like these become much harder to understand with this PR. What do you think, @sloede?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I have to agree with Hendrik.

The advantage of explicitly adding ndims(mesh) is that it makes the factor self-explanatory. Also in the other parts of the code, one avoids having to explain why a factor is what it is and how it was computed, or even how to change it when switching to another dimension.

I see the argument for more "elegant" code (avoiding the impression of having a varying dimension). However, I've seen too many codes where "magic" constants appear out of nowhere, and everyone reading it after 3 years has to re-understand what the intent/thought process of the original author was.

IMHO, the current version is concise and self-explanatory in terms of what is being done, thus I'd prefer to keep it as it is.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To throw in my two cents. I want to keep the ndims calls in the code. I think it is much more informative to have terms like 2^ndims(mesh) - 1 rather than values like 3 and 7. This is coming from the perspective of someone who recently became very familiar with the AMR implementation as I debugged and tracked down the loss of conservation. I think the current code is perfectly understandable.

Copy link
Contributor Author

@DanielDoehring DanielDoehring Sep 24, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I will close this PR then. Opinions about the status of the issue?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would close it in light of the comments above


# Loop over all elements in old container and either copy them or coarsen them
Expand All @@ -330,7 +330,7 @@ function coarsen!(u_ode::AbstractVector, adaptor,
# 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"
@assert all(to_be_removed[old_element_id:(old_element_id + 2^2 - 1)]) "bad cell/element order"
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One could move skip = 3 already up here


# Coarsen elements and store solution directly in new data structure
coarsen_elements!(u, element_id, old_u, old_element_id,
Expand All @@ -349,7 +349,7 @@ function coarsen!(u_ode::AbstractVector, adaptor,
# 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
skip = 3
else
if mesh isa P4estMesh
# Copy old element data to new element container and remove Jacobian scaling
Expand All @@ -372,7 +372,7 @@ function coarsen!(u_ode::AbstractVector, adaptor,
# 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")
@assert ninterfaces(cache.interfaces)==2 * nelements(dg, cache) ("For 2D and periodic domains and conforming elements, the number of interfaces must be 2 times the number of elements")
end

return nothing
Expand All @@ -394,8 +394,8 @@ function coarsen!(u_ode::AbstractVector, adaptor,
# Sanity check
if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 &&
!mpi_isparallel()
@assert ninterfaces(cache_parabolic.interfaces)==ndims(mesh) *
nelements(dg, cache_parabolic) ("For $(ndims(mesh))D and periodic domains and conforming elements, the number of interfaces must be $(ndims(mesh)) times the number of elements")
@assert ninterfaces(cache_parabolic.interfaces)==2 *
nelements(dg, cache_parabolic) ("For 2D and periodic domains and conforming elements, the number of interfaces must be 2 times the number of elements")
end

return nothing
Expand Down
16 changes: 8 additions & 8 deletions src/callbacks_step/amr_dg3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +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)^3 * 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
Expand Down Expand Up @@ -77,7 +77,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM

# Increment `element_id` on the refined mesh with the number
# of children, i.e., 8 in 3D
element_id += 2^ndims(mesh)
element_id += 8
else
if mesh isa P4estMesh
# Copy old element data to new element container and remove Jacobian scaling
Expand All @@ -98,12 +98,12 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estM
# 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))"
1||element_id == nelements(dg, cache) + 8 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
end # GC.@preserve old_u_ode old_inverse_jacobian

# 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")
@assert ninterfaces(cache.interfaces)==3 * nelements(dg, cache) ("For 3D and periodic domains and conforming elements, the number of interfaces must be 3 times the number of elements")
end

return nothing
Expand Down Expand Up @@ -229,7 +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)^3 * 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
Expand All @@ -250,7 +250,7 @@ function coarsen!(u_ode::AbstractVector, adaptor,
# 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"
@assert all(to_be_removed[old_element_id:(old_element_id + 2^3 - 1)]) "bad cell/element order"
Copy link
Contributor Author

@DanielDoehring DanielDoehring Sep 17, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One could move skip = 7 already up here


# Coarsen elements and store solution directly in new data structure
coarsen_elements!(u, element_id, old_u, old_element_id, adaptor,
Expand All @@ -269,7 +269,7 @@ function coarsen!(u_ode::AbstractVector, adaptor,
# 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
skip = 7
else
if mesh isa P4estMesh
# Copy old element data to new element container and remove Jacobian scaling
Expand All @@ -291,7 +291,7 @@ function coarsen!(u_ode::AbstractVector, adaptor,

# 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")
@assert ninterfaces(cache.interfaces)==3 * nelements(dg, cache) ("For 3D and periodic domains and conforming elements, the number of interfaces must be 3 times the number of elements")
end

return nothing
Expand Down
Loading