Add StructuredMeshView as proxy between solver and actual StructuredMesh

examples/structured_2d_dgsem/elixir_advection_meshview.jl

@@ -0,0 +1,122 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Coupled semidiscretization of two linear advection systems using converter functions
+# and mesh views for the semidiscretizations. First we define a parent mesh
+# for the entire physical domain, then we define the two mesh views on this parent.
+#
+# In this elixir, we have a square domain that is divided into left and right subdomains.
+# On each half of the domain, a completely independent `SemidiscretizationHyperbolic`
+# is created for the linear advection equations. The two systems are coupled in the
+# x-direction.
+# For a high-level overview, see also the figure below:
+#
+# (-1,  1)                                   ( 1,  1)
+#     ┌────────────────────┬────────────────────┐
+#     │    ↑ periodic ↑    │    ↑ periodic ↑    │
+#     │                    │                    │
+#     │     =========      │     =========      │
+#     │     system #1      │     system #2      │
+#     │     =========      │     =========      │
+#     │                    │                    │
+#     │<-- coupled         │<-- coupled         │
+#     │         coupled -->│         coupled -->│
+#     │                    │                    │
+#     │    ↓ periodic ↓    │    ↓ periodic ↓    │
+#     └────────────────────┴────────────────────┘
+# (-1, -1)                                   ( 1, -1)
+
+advection_velocity = (0.2, -0.7)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+# Domain size of the parent mesh.
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+# Cell dimensions of the parent mesh.
+cells_per_dimension = (16, 16)
+
+# Create parent mesh with 16 x 16 elements
+parent_mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+# Create the two mesh views, each of which takes half of the parent mesh.
+mesh1 = StructuredMeshView(parent_mesh; indices_min = (1, 1), indices_max = (8, 16))
+mesh2 = StructuredMeshView(parent_mesh; indices_min = (9, 1), indices_max = (16, 16))
+
+# The coupling function is simply the identity, as we are dealing with two identical systems.
+coupling_function = (x, u, equations_other, equations_own) -> u
+
+# Define the coupled boundary conditions
+# The indices (:end, :i_forward) and (:begin, :i_forward) denote the interface indexing.
+# For a system with coupling in x and y see examples/structured_2d_dgsem/elixir_advection_coupled.jl.
+boundary_conditions1 = (
+                        # Connect left boundary with right boundary of left mesh
+                        x_neg = BoundaryConditionCoupled(2, (:end, :i_forward), Float64,
+                                                         coupling_function),
+                        x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64,
+                                                         coupling_function),
+                        y_neg = boundary_condition_periodic,
+                        y_pos = boundary_condition_periodic)
+boundary_conditions2 = (
+                        # Connect left boundary with right boundary of left mesh
+                        x_neg = BoundaryConditionCoupled(1, (:end, :i_forward), Float64,
+                                                         coupling_function),
+                        x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64,
+                                                         coupling_function),
+                        y_neg = boundary_condition_periodic,
+                        y_pos = boundary_condition_periodic)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = boundary_conditions1)
+semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test,
+                                     solver,
+                                     boundary_conditions = boundary_conditions2)
+semi = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# 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_callback1 = AnalysisCallback(semi1, interval = 100)
+analysis_callback2 = AnalysisCallback(semi2, interval = 100)
+analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2)
+
+analysis_interval = 100
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     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, analysis_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 = 5.0e-2, # 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

@@ -225,7 +225,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, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
+       T8codeMesh
 
 export DG,
        DGSEM, LobattoLegendreBasis,

src/callbacks_step/analysis_dg2d.jl

@@ -30,7 +30,8 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2},
 end
 
 function create_cache_analysis(analyzer,
-                               mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,
+                               mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
+                                           UnstructuredMesh2D,
                                            P4estMesh{2}, T8codeMesh{2}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
@@ -107,8 +108,9 @@ function calc_error_norms(func, u, t, analyzer,
 end
 
 function calc_error_norms(func, u, t, analyzer,
-                          mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,
-                                      P4estMesh{2}, T8codeMesh{2}}, equations,
+                          mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
+                                      UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
+                          equations,
                           initial_condition, dg::DGSEM, cache, cache_analysis)
     @unpack vandermonde, weights = analyzer
     @unpack node_coordinates, inverse_jacobian = cache.elements
@@ -175,8 +177,10 @@ function integrate_via_indices(func::Func, u,
 end
 
 function integrate_via_indices(func::Func, u,
-                               mesh::Union{StructuredMesh{2}, UnstructuredMesh2D,
-                                           P4estMesh{2}, T8codeMesh{2}}, equations,
+                               mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
+                                           UnstructuredMesh2D, P4estMesh{2},
+                                           T8codeMesh{2}},
+                               equations,
                                dg::DGSEM, cache, args...; normalize = true) where {Func}
     @unpack weights = dg.basis
 
@@ -203,8 +207,8 @@ function integrate_via_indices(func::Func, u,
 end
 
 function integrate(func::Func, u,
-                   mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                               P4estMesh{2}, T8codeMesh{2}},
+                   mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2},
+                               UnstructuredMesh2D, P4estMesh{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
@@ -232,8 +236,8 @@ function integrate(func::Func, u,
 end
 
 function analyze(::typeof(entropy_timederivative), du, u, t,
-                 mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
-                             P4estMesh{2}, T8codeMesh{2}},
+                 mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2},
+                             UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}},
                  equations, dg::DG, cache)
     # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ
     integrate_via_indices(u, mesh, equations, dg, cache,

src/callbacks_step/save_solution_dg.jl

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

src/callbacks_step/stepsize_dg2d.jl

@@ -77,7 +77,7 @@ end
 
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                            T8codeMesh{2}},
+                            T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::False, equations, dg::DG, cache)
     # to avoid a division by zero if the speed vanishes everywhere,
     # e.g. for steady-state linear advection
@@ -113,7 +113,7 @@ end
 
 function max_dt(u, t,
                 mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
-                            T8codeMesh{2}},
+                            T8codeMesh{2}, StructuredMeshView{2}},
                 constant_speed::True, equations, dg::DG, cache)
     @unpack contravariant_vectors, inverse_jacobian = cache.elements
 

src/meshes/mesh_io.jl

@@ -97,7 +97,10 @@ end
 # of the mesh, like its size and the type of boundary mapping function.
 # Then, within Trixi2Vtk, the StructuredMesh and its node coordinates are reconstructured from
 # these attributes for plotting purposes
-function save_mesh_file(mesh::StructuredMesh, output_directory; system = "")
+# Note: the `timestep` argument is needed for compatibility with the method for
+# `StructuredMeshView`
+function save_mesh_file(mesh::StructuredMesh, output_directory; system = "",
+                        timestep = 0)
     # Create output directory (if it does not exist)
     mkpath(output_directory)
 
@@ -256,7 +259,7 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
         end
         mesh = TreeMesh(SerialTree{ndims}, max(n_cells_max, capacity))
         load_mesh!(mesh, mesh_file)
-    elseif mesh_type == "StructuredMesh"
+    elseif mesh_type in ("StructuredMesh", "StructuredMeshView")
         size_, mapping_as_string = h5open(mesh_file, "r") do file
             return read(attributes(file)["size"]),
                    read(attributes(file)["mapping"])

src/meshes/meshes.jl

@@ -7,6 +7,7 @@
 
 include("tree_mesh.jl")
 include("structured_mesh.jl")
+include("structured_mesh_view.jl")
 include("surface_interpolant.jl")
 include("unstructured_mesh.jl")
 include("face_interpolant.jl")

src/meshes/structured_mesh_view.jl

@@ -0,0 +1,132 @@
+# 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
+
+"""
+    StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS}
+
+A view on a structured curved mesh.
+"""
+mutable struct StructuredMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS}
+    parent::StructuredMesh{NDIMS, RealT}
+    cells_per_dimension::NTuple{NDIMS, Int}
+    mapping::Any # Not relevant for performance
+    mapping_as_string::String
+    current_filename::String
+    indices_min::NTuple{NDIMS, Int}
+    indices_max::NTuple{NDIMS, Int}
+    unsaved_changes::Bool
+end
+
+"""
+    StructuredMeshView(parent; indices_min, indices_max)
+
+Create a StructuredMeshView on a StructuredMesh parent.
+
+# Arguments
+- `parent`: the parent StructuredMesh.
+- `indices_min`: starting indices of the parent mesh.
+- `indices_max`: ending indices of the parent mesh.
+"""
+function StructuredMeshView(parent::StructuredMesh{NDIMS, RealT};
+                            indices_min = ntuple(_ -> 1, Val(NDIMS)),
+                            indices_max = size(parent)) where {NDIMS, RealT}
+    @assert indices_min <= indices_max
+    @assert all(indices_min .> 0)
+    @assert indices_max <= size(parent)
+
+    cells_per_dimension = indices_max .- indices_min .+ 1
+
+    # Compute cell sizes `deltas`
+    deltas = (parent.mapping.coordinates_max .- parent.mapping.coordinates_min) ./
+             parent.cells_per_dimension
+    # Calculate the domain boundaries.
+    coordinates_min = parent.mapping.coordinates_min .+ deltas .* (indices_min .- 1)
+    coordinates_max = parent.mapping.coordinates_min .+ deltas .* indices_max
+    mapping = coordinates2mapping(coordinates_min, coordinates_max)
+    mapping_as_string = """
+        coordinates_min = $coordinates_min
+        coordinates_max = $coordinates_max
+        mapping = coordinates2mapping(coordinates_min, coordinates_max)
+        """
+
+    return StructuredMeshView{NDIMS, RealT}(parent, cells_per_dimension, mapping,
+                                            mapping_as_string,
+                                            parent.current_filename,
+                                            indices_min, indices_max,
+                                            parent.unsaved_changes)
+end
+
+# Check if mesh is periodic
+function isperiodic(mesh::StructuredMeshView)
+    @unpack parent = mesh
+    return isperiodic(parent) && size(parent) == size(mesh)
+end
+
+function isperiodic(mesh::StructuredMeshView, dimension)
+    @unpack parent, indices_min, indices_max = mesh
+    return (isperiodic(parent, dimension) &&
+            indices_min[dimension] == 1 &&
+            indices_max[dimension] == size(parent, dimension))
+end
+
+@inline Base.ndims(::StructuredMeshView{NDIMS}) where {NDIMS} = NDIMS
+@inline Base.real(::StructuredMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT
+function Base.size(mesh::StructuredMeshView)
+    @unpack indices_min, indices_max = mesh
+    return indices_max .- indices_min .+ 1
+end
+function Base.size(mesh::StructuredMeshView, i)
+    @unpack indices_min, indices_max = mesh
+    return indices_max[i] - indices_min[i] + 1
+end
+Base.axes(mesh::StructuredMeshView) = map(Base.OneTo, size(mesh))
+Base.axes(mesh::StructuredMeshView, i) = Base.OneTo(size(mesh, i))
+
+function calc_node_coordinates!(node_coordinates, element,
+                                cell_x, cell_y, mapping,
+                                mesh::StructuredMeshView{2},
+                                basis)
+    @unpack nodes = basis
+
+    # Get cell length in reference mesh
+    dx = 2 / size(mesh, 1)
+    dy = 2 / size(mesh, 2)
+
+    # Calculate node coordinates of reference mesh
+    cell_x_offset = -1 + (cell_x - 1) * dx + dx / 2
+    cell_y_offset = -1 + (cell_y - 1) * dy + dy / 2
+
+    for j in eachnode(basis), i in eachnode(basis)
+        # node_coordinates are the mapped reference node_coordinates
+        node_coordinates[:, i, j, element] .= mapping(cell_x_offset + dx / 2 * nodes[i],
+                                                      cell_y_offset + dy / 2 * nodes[j])
+    end
+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 StructuredMesh and its node coordinates are reconstructured from
+# these attributes for plotting purposes.
+function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "",
+                        timestep = 0)
+    # Create output directory (if it does not exist)
+    mkpath(output_directory)
+
+    filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep))
+
+    # 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)["size"] = collect(size(mesh))
+        attributes(file)["mapping"] = mesh.mapping_as_string
+    end
+
+    return filename
+end
+end # @muladd