From 3df12530f0ed289f9b44149749af7ce60c2c323c Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Wed, 11 Sep 2024 15:40:49 +0200 Subject: [PATCH] Add type parameter `NDIMS_AMBIENT` to `P4estMesh` for dimension of ambient 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 --- src/callbacks_step/analysis_dg2d.jl | 33 +++++++++++++++-- src/callbacks_step/analysis_dg3d.jl | 46 ++++++++++++++++++++++-- src/meshes/p4est_mesh.jl | 39 +++++++++++++++----- src/solvers/dgsem_p4est/containers.jl | 3 +- src/solvers/dgsem_p4est/containers_2d.jl | 7 ++-- 5 files changed, 111 insertions(+), 17 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index de6b9a2a4a6..4cd92ce7c5f 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -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) diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 27e8a2b722f..fd501a7257c 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -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) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 768b24c31d4..65c0a431b29 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -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 @@ -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, @@ -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() @@ -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[] @@ -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) @@ -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 diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index f9830d0011c..3ef9cb2a421 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -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} diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 236d7d24c06..6af6fd6d90e 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -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)