Skip to content

Commit

Permalink
Add type parameter NDIMS_AMBIENT to P4estMesh for dimension of am…
Browse files Browse the repository at this point in the history 
…bient space (#2068)

* add ambient dimension as type parameter for P4estMesh

* add NDIMS_AMBIENT as parameter to serial/parallel aliases

* mesh.tree_node_coordinates should just be tree_node_coordinates in constructor.

* update Base.show to include NDIMS_AMBIENT type parameter

* fix typo

* format

* SerialP4estMesh{NDIMS} and ParallelP4estMesh{NDIMS} do not need NDIMS_AMBIENT

* add reference to PR in code comment

* add new type parameter to summary box

* describe use of NDIMS_AMBIENT in docstring and comment

* fix spelling

* remove reference to TrixiAtmo.jl and link to PR from comment

* add experimental feature note and fix NDIM -> NDIMS typo

---------

Co-authored-by: Hendrik Ranocha <[email protected]>
  • Loading branch information
@tristanmontoya @ranocha
tristanmontoya and ranocha authored Sep 11, 2024
1 parent 6bbcafc commit 3df1253
Show file tree
Hide file tree
Showing 5 changed files with 111 additions and 17 deletions.
33 changes: 31 additions & 2 deletions src/callbacks_step/analysis_dg2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,39 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2},
return (; u_local, u_tmp1, x_local, x_tmp1)
end

# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
function create_cache_analysis(analyzer, mesh::P4estMesh{2, NDIMS_AMBIENT},
equations, dg::DG, cache,
RealT, uEltype) where {NDIMS_AMBIENT}

# pre-allocate buffers
# We use `StrideArray`s here since these buffers are used in performance-critical
# places and the additional information passed to the compiler makes them faster
# than native `Array`s.
u_local = StrideArray(undef, uEltype,
StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)))
u_tmp1 = StrideArray(undef, uEltype,
StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(dg)))
x_local = StrideArray(undef, RealT,
StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)))
x_tmp1 = StrideArray(undef, RealT,
StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(dg)))
jacobian_local = StrideArray(undef, RealT,
StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)))
jacobian_tmp1 = StrideArray(undef, RealT,
StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))

return (; u_local, u_tmp1, x_local, x_tmp1, jacobian_local, jacobian_tmp1)
end

function create_cache_analysis(analyzer,
mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
UnstructuredMesh2D,
P4estMesh{2}, T8codeMesh{2}},
UnstructuredMesh2D, T8codeMesh{2}},
equations, dg::DG, cache,
RealT, uEltype)

Expand Down
46 changes: 44 additions & 2 deletions src/callbacks_step/analysis_dg3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,9 +35,51 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{3},
return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2)
end

# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
function create_cache_analysis(analyzer,
mesh::Union{StructuredMesh{3}, P4estMesh{3},
T8codeMesh{3}},
mesh::P4estMesh{3, NDIMS_AMBIENT},
equations, dg::DG, cache,
RealT, uEltype) where {NDIMS_AMBIENT}

# pre-allocate buffers
# We use `StrideArray`s here since these buffers are used in performance-critical
# places and the additional information passed to the compiler makes them faster
# than native `Array`s.
u_local = StrideArray(undef, uEltype,
StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
u_tmp1 = StrideArray(undef, uEltype,
StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
u_tmp2 = StrideArray(undef, uEltype,
StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
x_local = StrideArray(undef, RealT,
StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
x_tmp1 = StrideArray(undef, RealT,
StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
x_tmp2 = StrideArray(undef, RealT,
StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
jacobian_local = StrideArray(undef, RealT,
StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)))
jacobian_tmp1 = StrideArray(undef, RealT,
StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)),
StaticInt(nnodes(dg)))
jacobian_tmp2 = StrideArray(undef, RealT,
StaticInt(nnodes(analyzer)),
StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))

return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local,
jacobian_tmp1, jacobian_tmp2)
end

function create_cache_analysis(analyzer,
mesh::Union{StructuredMesh{3}, T8codeMesh{3}},
equations, dg::DG, cache,
RealT, uEltype)

Expand Down
39 changes: 30 additions & 9 deletions src/meshes/p4est_mesh.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,12 +6,23 @@
#! format: noindent

"""
P4estMesh{NDIMS} <: AbstractMesh{NDIMS}
P4estMesh{NDIMS, NDIMS_AMBIENT} <: AbstractMesh{NDIMS}
An unstructured curved mesh based on trees that uses the C library `p4est`
to manage trees and mesh refinement.
The parameter `NDIMS` denotes the dimension of the spatial domain or manifold represented
by the mesh itself, while `NDIMS_AMBIENT` denotes the dimension of the ambient space in
which the mesh is embedded. For example, the type `P4estMesh{3, 3}` corresponds to a
standard mesh for a three-dimensional volume, whereas `P4estMesh{2, 3}` corresponds to a
mesh for a two-dimensional surface or shell in three-dimensional space.
!!! warning "Experimental implementation"
The use of `NDIMS != NDIMS_AMBIENT` is an experimental feature and may change in future
releases.
"""
mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <:
mutable struct P4estMesh{NDIMS, NDIMS_AMBIENT, RealT <: Real, IsParallel, P, Ghost,
NDIMSP2, NNODES} <:
AbstractMesh{NDIMS}
p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t}
is_parallel :: IsParallel
Expand Down Expand Up @@ -48,7 +59,14 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
ghost = ghost_new_p4est(p4est)
ghost_pw = PointerWrapper(ghost)

mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel),
# To enable the treatment of a manifold of dimension NDIMS embedded within an
# ambient space of dimension NDIMS_AMBIENT, we store both as type parameters and
# allow them to differ in the general case. This functionality is used for
# constructing discretizations on spherical shell domains for applications in
# global atmospheric modelling. The ambient dimension NDIMS_AMBIENT is therefore
# set here in the inner constructor to size(tree_node_coordinates, 1).
mesh = new{NDIMS, size(tree_node_coordinates, 1),
eltype(tree_node_coordinates), typeof(is_parallel),
typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw,
is_parallel,
ghost_pw,
Expand All @@ -66,8 +84,8 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
end
end

const SerialP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:False}
const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:True}
const SerialP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Any, <:Real, <:False}
const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Any, <:Real, <:True}

@inline mpi_parallel(mesh::SerialP4estMesh) = False()
@inline mpi_parallel(mesh::ParallelP4estMesh) = True()
Expand All @@ -87,7 +105,8 @@ function destroy_mesh(mesh::P4estMesh{3})
end

@inline Base.ndims(::P4estMesh{NDIMS}) where {NDIMS} = NDIMS
@inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
@inline Base.real(::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} = RealT
@inline ndims_ambient(::P4estMesh{NDIMS, NDIMS_AMBIENT}) where {NDIMS, NDIMS_AMBIENT} = NDIMS_AMBIENT

@inline function ntrees(mesh::P4estMesh)
return mesh.p4est.trees.elem_count[]
Expand All @@ -97,7 +116,8 @@ end
@inline ncellsglobal(mesh::P4estMesh) = Int(mesh.p4est.global_num_quadrants[])

function Base.show(io::IO, mesh::P4estMesh)
print(io, "P4estMesh{", ndims(mesh), ", ", real(mesh), "}")
print(io, "P4estMesh{", ndims(mesh), ", ", ndims_ambient(mesh), ", ", real(mesh),
"}")
end

function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMesh)
Expand All @@ -110,8 +130,9 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMesh)
"polydeg" => length(mesh.nodes) - 1
]
summary_box(io,
"P4estMesh{" * string(ndims(mesh)) * ", " * string(real(mesh)) *
"}", setup)
"P4estMesh{" * string(ndims(mesh)) * ", " *
string(ndims_ambient(mesh)) *
", " * string(real(mesh)) * "}", setup)
end
end

Expand Down
3 changes: 2 additions & 1 deletion src/solvers/dgsem_p4est/containers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,8 @@ function Base.resize!(elements::P4estElementContainer, capacity)
end

# Create element container and initialize element data
function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, RealT}},
function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT},
T8codeMesh{NDIMS, RealT}},
equations,
basis,
::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real}
Expand Down
7 changes: 4 additions & 3 deletions src/solvers/dgsem_p4est/containers_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,14 @@ end

# Interpolate tree_node_coordinates to each quadrant at the specified nodes
function calc_node_coordinates!(node_coordinates,
mesh::P4estMesh{2},
nodes::AbstractVector)
mesh::P4estMesh{2, NDIMS_AMBIENT},
nodes::AbstractVector) where {NDIMS_AMBIENT}
# We use `StrideArray`s here since these buffers are used in performance-critical
# places and the additional information passed to the compiler makes them faster
# than native `Array`s.
tmp1 = StrideArray(undef, real(mesh),
StaticInt(2), static_length(nodes), static_length(mesh.nodes))
StaticInt(NDIMS_AMBIENT), static_length(nodes),
static_length(mesh.nodes))
matrix1 = StrideArray(undef, real(mesh),
static_length(nodes), static_length(mesh.nodes))
matrix2 = similar(matrix1)
Expand Down

0 comments on commit 3df1253

Please sign in to comment.