Added p4estMeshView to various routines.

More to come, especially bug fixing.

examples/p4est_2d_dgsem/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

src/Trixi.jl

@@ -220,7 +220,8 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
 export lake_at_rest_error
 export ncomponents, eachcomponent
 
-export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh
+export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, P4estMeshView,
+       T8codeMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,

src/callbacks_step/analysis_dg2d.jl

@@ -31,7 +31,7 @@ end
 
 function create_cache_analysis(analyzer,
                                mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,
-                                           P4estMesh{2}, T8codeMesh{2}},
+                                           P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
 
@@ -108,7 +108,7 @@ end
 
 function calc_error_norms(func, u, t, analyzer,
                           mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,
-                                      P4estMesh{2}, T8codeMesh{2}}, equations,
+				      P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations,
                           initial_condition, dg::DGSEM, cache, cache_analysis)
     @unpack vandermonde, weights = analyzer
     @unpack node_coordinates, inverse_jacobian = cache.elements
@@ -176,7 +176,7 @@ end
 
 function integrate_via_indices(func::Func, u,
                                mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,
-                                           P4estMesh{2}, T8codeMesh{2}}, equations,
+					   P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations,
                                dg::DGSEM, cache, args...; normalize = true) where {Func}
     @unpack weights = dg.basis
 
@@ -204,7 +204,7 @@ end
 
 function integrate(func::Func, u,
                    mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                               P4estMesh{2}, T8codeMesh{2}},
+			       P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                    equations, dg::DG, cache; normalize = true) where {Func}
     integrate_via_indices(u, mesh, equations, dg, cache;
                           normalize = normalize) do u, i, j, element, equations, dg
@@ -214,7 +214,7 @@ function integrate(func::Func, u,
 end
 
 function integrate(func::Func, u,
-                   mesh::Union{TreeMesh{2}, P4estMesh{2}},
+   		mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}},
                    equations, equations_parabolic,
                    dg::DGSEM,
                    cache, cache_parabolic; normalize = true) where {Func}
@@ -233,7 +233,7 @@ end
 
 function analyze(::typeof(entropy_timederivative), du, u, t,
                  mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                             P4estMesh{2}, T8codeMesh{2}},
+			     P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                  equations, dg::DG, cache)
     # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ
     integrate_via_indices(u, mesh, equations, dg, cache,
@@ -278,7 +278,7 @@ end
 
 function analyze(::Val{:l2_divb}, du, u, t,
                  mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                             T8codeMesh{2}},
+			     P4estMeshView{2}, T8codeMesh{2}},
                  equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
     @unpack contravariant_vectors = cache.elements
     integrate_via_indices(u, mesh, equations, dg, cache, cache,
@@ -346,7 +346,7 @@ end
 
 function analyze(::Val{:linf_divb}, du, u, t,
                  mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                             T8codeMesh{2}},
+			     P4estMeshView{2}, T8codeMesh{2}},
                  equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
     @unpack derivative_matrix, weights = dg.basis
     @unpack contravariant_vectors = cache.elements

src/meshes/meshes.jl

@@ -12,6 +12,7 @@ include("unstructured_mesh.jl")
 include("face_interpolant.jl")
 include("transfinite_mappings_3d.jl")
 include("p4est_mesh.jl")
+include("p4est_mesh_view.jl")
 include("t8code_mesh.jl")
 include("mesh_io.jl")
 include("dgmulti_meshes.jl")

src/meshes/p4est_mesh_view.jl

@@ -3,20 +3,19 @@
 
 mutable struct P4estMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS}
     parent::P4estMesh{NDIMS, RealT}
-    mapping::Any # Not relevant for performance
-    index_min::NTuple{NDIMS, Int}
-    index_max::NTuple{NDIMS, Int}
+    nodes::SVector
+    index_min::Int
+    index_max::Int
 end
 
 function P4estMeshView(parent::P4estMesh{NDIMS, RealT};
-                        index_min = ntuple(_ -> 1, Val(NDIMS)),
-                        index_max = size(parent)) where {NDIMS, RealT}
+                        index_min = 1::Int,
+                        index_max = sizeof(unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees))::Int) where {NDIMS, RealT}
     @assert index_min <= index_max
     @assert all(index_min .> 0)
-    @assert index_max <= size(parent)
+    @assert index_max <= sizeof(unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees))
 
-    return P4estMeshView{NDIMS, RealT}(parent, parent.mapping, index_min,
-                                       index_max)
+    return P4estMeshView{NDIMS, RealT}(parent, parent.nodes, index_min, index_max)
 end
 
 # Check if mesh is periodic
@@ -45,6 +44,36 @@ end
 Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh))
 Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i))
 
+function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL)
+    p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn)
+    # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes
+    # See https://github.com/cburstedde/p4est/issues/112
+    p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn)
+end
+
+@inline ncells(mesh::P4estMeshView) = Int(mesh.parent.p4est.local_num_quadrants[])
+
+function count_required_surfaces(mesh::P4estMeshView)
+    # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face
+    iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh.parent)))
+
+    # interfaces, mortars, boundaries
+    user_data = [0, 0, 0]
+
+    iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c)
+
+    # Return counters
+    return (interfaces = user_data[1],
+            mortars = user_data[2],
+            boundaries = user_data[3])
+end
+
+function init_interfaces!(interfaces, mesh::P4estMeshView)
+    init_surfaces!(interfaces, nothing, nothing, mesh.parent)
+
+    return interfaces
+end
+
 # function calc_node_coordinates!(node_coordinates, element,
 #                                 cell_x, cell_y, mapping,
 #                                 mesh::P4estMeshView{2},
@@ -93,7 +122,6 @@ function calc_node_coordinates!(node_coordinates,
     p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l)
 
     trees = unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees)
-    @autoinfiltrate
 
     for tree in eachindex(trees)
         offset = trees[tree].quadrants_offset

src/solvers/dg.jl

@@ -436,7 +436,7 @@ function get_node_variables!(node_variables, mesh, equations, dg::DG, cache)
 end
 
 const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh,
-                          T8codeMesh}
+                          P4estMeshView, T8codeMesh}
 
 @inline function ndofs(mesh::MeshesDGSEM, dg::DG, cache)
     nelements(cache.elements) * nnodes(dg)^ndims(mesh)

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, RealT}, P4estMeshView{NDIMS, RealT},
+				   T8codeMesh{NDIMS, RealT}},
                        equations,
                        basis,
                        ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real}
@@ -166,7 +167,7 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity)
 end
 
 # Create interface container and initialize interface data.
-function init_interfaces(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elements)
+function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, basis, elements)
     NDIMS = ndims(elements)
     uEltype = eltype(elements)
 

src/solvers/dgsem_p4est/containers_2d.jl

@@ -6,7 +6,7 @@
 #! format: noindent
 
 # Initialize data structures in element container
-function init_elements!(elements, mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+function init_elements!(elements, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                         basis::LobattoLegendreBasis)
     @unpack node_coordinates, jacobian_matrix,
     contravariant_vectors, inverse_jacobian = elements
@@ -26,7 +26,7 @@ end
 
 # Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis
 function calc_node_coordinates!(node_coordinates,
-                                mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+		mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                                 basis::LobattoLegendreBasis)
     # Hanging nodes will cause holes in the mesh if its polydeg is higher
     # than the polydeg of the solver.
@@ -37,7 +37,7 @@ end
 
 # Interpolate tree_node_coordinates to each quadrant at the specified nodes
 function calc_node_coordinates!(node_coordinates,
-                                mesh::P4estMesh{2},
+		mesh::Union{P4estMesh{2}, P4estMeshView{2}},
                                 nodes::AbstractVector)
     # 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

src/solvers/dgsem_p4est/dg.jl

@@ -8,7 +8,7 @@
 # This method is called when a SemidiscretizationHyperbolic is constructed.
 # It constructs the basic `cache` used throughout the simulation to compute
 # the RHS etc.
-function create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::Any,
+function create_cache(mesh::Union{P4estMesh, P4estMeshView}, equations::AbstractEquations, dg::DG, ::Any,
                       ::Type{uEltype}) where {uEltype <: Real}
     # Make sure to balance the `p4est` before creating any containers
     # in case someone has tampered with the `p4est` after creating the mesh

src/solvers/dgsem_p4est/dg_2d.jl

@@ -7,7 +7,7 @@
 
 # The methods below are specialized on the mortar type
 # and called from the basic `create_cache` method at the top.
-function create_cache(mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations,
+function create_cache(mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations,
                       mortar_l2::LobattoLegendreMortarL2, uEltype)
     # TODO: Taal performance using different types
     MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)},
@@ -58,7 +58,7 @@ end
 
 # We pass the `surface_integral` argument solely for dispatch
 function prolong2interfaces!(cache, u,
-                             mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+		mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                              equations, surface_integral, dg::DG)
     @unpack interfaces = cache
     index_range = eachnode(dg)
@@ -114,7 +114,7 @@ function prolong2interfaces!(cache, u,
 end
 
 function calc_interface_flux!(surface_flux_values,
-                              mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+		mesh::Union{P4estMesh{2}, P4estMesh{2}, T8codeMesh{2}},
                               nonconservative_terms,
                               equations, surface_integral, dg::DG, cache)
     @unpack neighbor_ids, node_indices = cache.interfaces
@@ -182,7 +182,7 @@ end
 
 # Inlined version of the interface flux computation for conservation laws
 @inline function calc_interface_flux!(surface_flux_values,
-                                      mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+		mesh::Union{P4estMesh{2}, P4estMesh{2}, T8codeMesh{2}},
                                       nonconservative_terms::False, equations,
                                       surface_integral, dg::DG, cache,
                                       interface_index, normal_direction,
@@ -206,7 +206,7 @@ end
 
 # Inlined version of the interface flux computation for equations with conservative and nonconservative terms
 @inline function calc_interface_flux!(surface_flux_values,
-                                      mesh::Union{P4estMesh{2}, T8codeMesh{2}},
+		mesh::Union{P4estMesh{2}, P4estMesh{2}, T8codeMesh{2}},
                                       nonconservative_terms::True, equations,
                                       surface_integral, dg::DG, cache,
                                       interface_index, normal_direction,

src/solvers/dgsem_tree/dg_2d.jl

@@ -37,14 +37,14 @@ end
 # The methods below are specialized on the volume integral type
 # and called from the basic `create_cache` method at the top.
 function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                                  P4estMesh{2}, T8codeMesh{2}},
+				  P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                       equations, volume_integral::VolumeIntegralFluxDifferencing,
                       dg::DG, uEltype)
     NamedTuple()
 end
 
 function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                                  P4estMesh{2}, T8codeMesh{2}}, equations,
+				  P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations,
                       volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
     element_ids_dg = Int[]
     element_ids_dgfv = Int[]
@@ -70,7 +70,7 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMe
 end
 
 function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                                  P4estMesh{2}, T8codeMesh{2}}, equations,
+				  P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations,
                       volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG,
                       uEltype)
     A3dp1_x = Array{uEltype, 3}
@@ -92,7 +92,7 @@ end
 # The methods below are specialized on the mortar type
 # and called from the basic `create_cache` method at the top.
 function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                                  P4estMesh{2}, T8codeMesh{2}},
+				  P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
                       equations, mortar_l2::LobattoLegendreMortarL2, uEltype)
     # TODO: Taal performance using different types
     MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2,
@@ -110,7 +110,7 @@ end
 # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer?
 
 function rhs!(du, u, t,
-              mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, equations,
+		mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations,
               initial_condition, boundary_conditions, source_terms::Source,
               dg::DG, cache) where {Source}
     # Reset du
@@ -181,7 +181,7 @@ end
 function calc_volume_integral!(du, u,
                                mesh::Union{TreeMesh{2}, StructuredMesh{2},
                                            UnstructuredMesh2D, P4estMesh{2},
-                                           T8codeMesh{2}},
+					   P4estMeshView{2}, T8codeMesh{2}},
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralWeakForm,
                                dg::DGSEM, cache)
@@ -235,7 +235,7 @@ end
 function calc_volume_integral!(du, u,
                                mesh::Union{TreeMesh{2}, StructuredMesh{2},
                                            UnstructuredMesh2D, P4estMesh{2},
-                                           T8codeMesh{2}},
+					   P4estMeshView{2}, T8codeMesh{2}},
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralFluxDifferencing,
                                dg::DGSEM, cache)
@@ -332,7 +332,7 @@ end
 function calc_volume_integral!(du, u,
                                mesh::Union{TreeMesh{2}, StructuredMesh{2},
                                            UnstructuredMesh2D, P4estMesh{2},
-                                           T8codeMesh{2}},
+					   P4estMeshView{2}, T8codeMesh{2}},
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralShockCapturingHG,
                                dg::DGSEM, cache)
@@ -391,7 +391,8 @@ end
 
 @inline function fv_kernel!(du, u,
                             mesh::Union{TreeMesh{2}, StructuredMesh{2},
-                                        UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}
+                                        UnstructuredMesh2D, P4estMesh{2},
+					P4estMeshView{2}, T8codeMesh{2}
                                         },
                             nonconservative_terms, equations,
                             volume_flux_fv, dg::DGSEM, cache, element, alpha = true)