WIP: Sc/Enhanced p4est mesh views

examples/p4est_2d_dgsem/elixir_advection_meshview.jl

@@ -0,0 +1,65 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Most basic p4est mesh view setup where the entire domain
+# is part of the single mesh view.
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
+
+trees_per_dimension = (8, 8)
+
+# Create parent P4estMesh with 8 x 8 trees and 8 x 8 elements
+parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                        coordinates_min = coordinates_min,
+                        coordinates_max = coordinates_max,
+                        initial_refinement_level = 0)
+
+# Define the mesh view covering the whole parent mesh.
+cell_ids = Vector(1:prod(trees_per_dimension))
+mesh = P4estMeshView(parent_mesh, cell_ids)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+ode = semidiscretize(semi, (0.0, 1.0))
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, save_solution, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

src/Trixi.jl

@@ -235,7 +235,7 @@ export lake_at_rest_error
 export ncomponents, eachcomponent
 
 export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
-       T8codeMesh
+       P4estMeshView, T8codeMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,

src/callbacks_step/analysis_dg2d.jl

@@ -30,7 +30,9 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2},
 end
 
 # Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
-function create_cache_analysis(analyzer, mesh::P4estMesh{2, NDIMS_AMBIENT},
+function create_cache_analysis(analyzer,
+                               mesh::Union{P4estMesh{2, NDIMS_AMBIENT},
+                                           P4estMeshView{2, NDIMS_AMBIENT}},
                                equations, dg::DG, cache,
                                RealT, uEltype) where {NDIMS_AMBIENT}
 
@@ -138,7 +140,9 @@ end
 
 function calc_error_norms(func, u, t, analyzer,
                           mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
-                                      UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
+                                      UnstructuredMesh2D, P4estMesh{2},
+                                      P4estMeshView{2},
+                                      T8codeMesh{2}},
                           equations,
                           initial_condition, dg::DGSEM, cache, cache_analysis)
     @unpack vandermonde, weights = analyzer
@@ -239,7 +243,8 @@ end
 
 function integrate(func::Func, u,
                    mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2},
-                               UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
+                               UnstructuredMesh2D, 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

src/callbacks_step/save_solution_dg.jl

@@ -9,6 +9,7 @@ function save_solution_file(u, time, dt, timestep,
                             mesh::Union{SerialTreeMesh, StructuredMesh,
                                         StructuredMeshView,
                                         UnstructuredMesh2D, SerialP4estMesh,
+                                        P4estMeshView,
                                         SerialT8codeMesh},
                             equations, dg::DG, cache,
                             solution_callback,

src/callbacks_step/stepsize_dg2d.jl

@@ -115,7 +115,7 @@ end
 
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                            T8codeMesh{2}, StructuredMeshView{2}},
+                            P4estMeshView{2}, T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::True, equations, dg::DG, cache)
     @unpack contravariant_vectors, inverse_jacobian = cache.elements
 

src/meshes/mesh_io.jl

@@ -6,7 +6,8 @@
 #! format: noindent
 
 # Save current mesh with some context information as an HDF5 file.
-function save_mesh_file(mesh::Union{TreeMesh, P4estMesh, T8codeMesh}, output_directory,
+function save_mesh_file(mesh::Union{TreeMesh, P4estMesh, P4estMeshView, T8codeMesh},
+                        output_directory,
                         timestep = 0)
     save_mesh_file(mesh, output_directory, timestep, mpi_parallel(mesh))
 end
@@ -382,7 +383,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
                                   periodicity = periodicity_,
                                   unsaved_changes = false)
         mesh.current_filename = mesh_file
-    elseif mesh_type == "P4estMesh"
+    elseif mesh_type == "P4estMesh" || mesh_type == "P4estMeshView"
         p4est_filename, tree_node_coordinates,
         nodes, boundary_names_ = h5open(mesh_file, "r") do file
             return read(attributes(file)["p4est_file"]),

src/meshes/meshes.jl

@@ -13,6 +13,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

@@ -0,0 +1,143 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+"""
+    P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: AbstractMesh{NDIMS}
+
+A view on a p4est mesh.
+"""
+mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <:
+               AbstractMesh{NDIMS}
+    parent::Parent
+    cell_ids::Vector{Int}
+    unsaved_changes::Bool
+    current_filename::String
+end
+
+"""
+    P4estMeshView(parent; cell_ids)
+
+Create a P4estMeshView on a P4estMesh parent.
+
+# Arguments
+- `parent`: the parent P4estMesh.
+- `cell_ids`: array of cell ids that are part of this view.
+"""
+function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT},
+                       cell_ids::Vector) where {NDIMS, NDIMS_AMBIENT, RealT}
+    return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids,
+                                                                      parent.unsaved_changes,
+                                                                      parent.current_filename)
+end
+
+@inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS
+@inline Base.real(::P4estMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT
+@inline ndims_ambient(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS
+
+@inline balance!(::P4estMeshView) = nothing
+@inline ncells(mesh::P4estMeshView) = length(mesh.cell_ids)
+
+function extract_p4est_mesh_view(elements_parent,
+                                 interfaces_parent,
+                                 boundaries_parent,
+                                 mortars_parent,
+                                 mesh)
+    elements = elements_parent
+    elements.inverse_jacobian = elements_parent.inverse_jacobian[.., mesh.cell_ids]
+    elements.jacobian_matrix = elements_parent.jacobian_matrix[.., mesh.cell_ids]
+    elements.node_coordinates = elements_parent.node_coordinates[.., mesh.cell_ids]
+    elements.contravariant_vectors = elements_parent.contravariant_vectors[..,
+                                                                           mesh.cell_ids]
+    elements.surface_flux_values = elements_parent.surface_flux_values[..,
+                                                                       mesh.cell_ids]
+    elements._inverse_jacobian = vec(elements.inverse_jacobian)
+    elements._jacobian_matrix = vec(elements.jacobian_matrix)
+    elements._node_coordinates = vec(elements.node_coordinates)
+    elements._node_coordinates = vec(elements.node_coordinates)
+    elements._surface_flux_values = vec(elements.surface_flux_values)
+    interfaces = extract_interfaces(mesh, interfaces_parent)
+
+    return elements, interfaces, boundaries_parent, mortars_parent
+end
+
+function extract_interfaces(mesh::P4estMeshView, interfaces_parent)
+    """
+    Remove all interfaces that have a tuple of neighbor_ids where at least one is
+    not part of this meshview, i.e. mesh.cell_ids.
+    """
+
+    mask = BitArray(undef, size(interfaces_parent.neighbor_ids)[2])
+    for interface in 1:size(interfaces_parent.neighbor_ids)[2]
+        mask[interface] = (interfaces_parent.neighbor_ids[1, interface] in mesh.cell_ids) &&
+                          (interfaces_parent.neighbor_ids[2, interface] in mesh.cell_ids)
+    end
+    interfaces = interfaces_parent
+    interfaces.u = interfaces_parent.u[.., mask]
+    interfaces.node_indices = interfaces_parent.node_indices[.., mask]
+    neighbor_ids = interfaces_parent.neighbor_ids[.., mask]
+    # Transform the global (parent) indices into local (view) indices.
+    interfaces.neighbor_ids = zeros(Int, size(neighbor_ids))
+    for interface in 1:size(neighbor_ids)[2]
+        interfaces.neighbor_ids[1, interface] = findall(id -> id ==
+                                                              neighbor_ids[1, interface],
+                                                        mesh.cell_ids)[1]
+        interfaces.neighbor_ids[2, interface] = findall(id -> id ==
+                                                              neighbor_ids[2, interface],
+                                                        mesh.cell_ids)[1]
+    end
+
+    # Flatten the arrays.
+    interfaces._u = vec(interfaces.u)
+    interfaces._node_indices = vec(interfaces.node_indices)
+    interfaces._neighbor_ids = vec(interfaces.neighbor_ids)
+
+    return interfaces
+end
+
+# Does not save the mesh itself to an HDF5 file. Instead saves important attributes
+# of the mesh, like its size and the type of boundary mapping function.
+# Then, within Trixi2Vtk, the P4estMeshView and its node coordinates are reconstructured from
+# these attributes for plotting purposes
+function save_mesh_file(mesh::P4estMeshView, output_directory, timestep,
+                        mpi_parallel::False)
+    # Create output directory (if it does not exist)
+    mkpath(output_directory)
+
+    # Determine file name based on existence of meaningful time step
+    if timestep > 0
+        filename = joinpath(output_directory, @sprintf("mesh_%09d.h5", timestep))
+        p4est_filename = @sprintf("p4est_data_%09d", timestep)
+    else
+        filename = joinpath(output_directory, "mesh.h5")
+        p4est_filename = "p4est_data"
+    end
+
+    p4est_file = joinpath(output_directory, p4est_filename)
+
+    # Save the complete connectivity and `p4est` data to disk.
+    save_p4est!(p4est_file, mesh.parent.p4est)
+
+    # Open file (clobber existing content)
+    h5open(filename, "w") do file
+        # Add context information as attributes
+        attributes(file)["mesh_type"] = get_name(mesh)
+        attributes(file)["ndims"] = ndims(mesh)
+        attributes(file)["p4est_file"] = p4est_filename
+
+        file["tree_node_coordinates"] = mesh.parent.tree_node_coordinates[..,
+                                                                          mesh.cell_ids]
+        file["nodes"] = Vector(mesh.parent.nodes) # the mesh uses `SVector`s for the nodes
+        # to increase the runtime performance
+        # but HDF5 can only handle plain arrays
+        file["boundary_names"] = mesh.parent.boundary_names .|> String
+    end
+
+    return filename
+end
+
+@inline mpi_parallel(mesh::P4estMeshView) = False()
+end # @muladd

src/semidiscretization/semidiscretization_hyperbolic.jl

@@ -214,6 +214,7 @@ end
 
 # No checks for these meshes yet available
 function check_periodicity_mesh_boundary_conditions(mesh::Union{P4estMesh,
+                                                                P4estMeshView,
                                                                 UnstructuredMesh2D,
                                                                 T8codeMesh,
                                                                 DGMultiMesh},

src/solvers/dg.jl

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