From 8a9fc7baeca9807de185592d5cb8d60040a24f09 Mon Sep 17 00:00:00 2001 From: Michael Schlottke-Lakemper Date: Fri, 10 May 2024 06:42:17 +0200 Subject: [PATCH] Add StructuredMeshView as proxy between solver and actual StructuredMesh (#1624) * Add StructuredMeshView as proxy between solver and actual StructuredMesh * Enable StructuredMeshView on submesh * Attempt to using coupled with meshviews. * Update structured_mesh_view.jl Format. * Applied autoformatter on smview elixir. * Applied autoformatter on meshview. * Corected minor typo with StructuredMeshView. * Corrected x-boundaries for smview example elixir. * Removed parent's periodicity for smview. * Removed redundant and problematic redifintion of StructuredMeshView function. * Applied auto formatter on files affected by the structured mesh view. * Applied autoforatter. * Added entries for meshview IO. * Added structuresd mesh view to parametric types. * Added unctionality of writing out mesh files for StructuredMeshView for different time steps. * Added temptative code for dynamically changing mesh view sizes. * Applied autoformatter. * Corrected the calculation of the coordinate mapping for mesh views. * Added cells_per_dimension to the StructuredMeshView. * Added StructuredMeshView to max_dt clculations. * Applied autoformatter. * Applied autoformatter. * Cleaned up meshview coupled example. * Added 2d structured mesh views to periodicity checks. * Cleaned up coupled strucutred mesh view example so that it can be used as test. * Added 2d structured mesh view to the tests. * Added analysis_interval variable. * Applied updated autoformatter. * Removed unused lines of code. Added comment about muladd. * Added explanatory comments. * Removed unused code. * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring * Update src/meshes/mesh_io.jl Co-authored-by: Daniel Doehring * Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring * Corrected comment on solver. * Changed order of imported meshes. * Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring * Added comment about coupling interface indexing. * Applied autoformatter. * Update src/meshes/structured_mesh_view.jl Co-authored-by: Daniel Doehring * Renamed index_min and index_max to indices_min and indices_max. * Removed relative tolerance from meshview test. * Removed further relative tolerance as it is not needed. * Update examples/structured_2d_dgsem/elixir_advection_smview.jl Co-authored-by: Daniel Doehring * Fixed typo on parent mesh generation. * Applied autoformatter. * Updated test results for elixir_advection_smview.jl. * Update src/meshes/mesh_io.jl Co-authored-by: Michael Schlottke-Lakemper * Update src/meshes/mesh_io.jl Co-authored-by: Michael Schlottke-Lakemper * Renamed elixir_advection_smview.jl to elixir_advection_meshview.jl since 'structured' is redundant for a file in this directory. * Renamed elixir_advection_smview.jl to elixir_advection_meshview.jl. * Moved save_mesh_file function for StructuredMeshView. * Update src/meshes/structured_mesh_view.jl Co-authored-by: Michael Schlottke-Lakemper * Removed redundant check for structured mesh views. --------- Co-authored-by: SimonCan Co-authored-by: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Co-authored-by: iomsn Co-authored-by: Daniel Doehring --- .../elixir_advection_meshview.jl | 122 ++++++++++++++++ src/Trixi.jl | 3 +- src/callbacks_step/analysis_dg2d.jl | 22 +-- src/callbacks_step/save_solution_dg.jl | 1 + src/callbacks_step/stepsize_dg2d.jl | 4 +- src/meshes/mesh_io.jl | 7 +- src/meshes/meshes.jl | 1 + src/meshes/structured_mesh_view.jl | 132 ++++++++++++++++++ .../semidiscretization_coupled.jl | 10 +- .../semidiscretization_hyperbolic.jl | 3 +- src/solvers/dg.jl | 5 +- src/solvers/dgsem_structured/containers.jl | 3 +- src/solvers/dgsem_structured/containers_2d.jl | 7 +- src/solvers/dgsem_structured/dg.jl | 11 +- src/solvers/dgsem_structured/dg_2d.jl | 23 +-- src/solvers/dgsem_tree/dg_2d.jl | 8 +- test/test_structured_2d.jl | 24 ++++ 17 files changed, 347 insertions(+), 39 deletions(-) create mode 100644 examples/structured_2d_dgsem/elixir_advection_meshview.jl create mode 100644 src/meshes/structured_mesh_view.jl diff --git a/examples/structured_2d_dgsem/elixir_advection_meshview.jl b/examples/structured_2d_dgsem/elixir_advection_meshview.jl new file mode 100644 index 00000000000..d8d27031090 --- /dev/null +++ b/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() diff --git a/src/Trixi.jl b/src/Trixi.jl index 5511f7e02e2..f3977f1f058 100644 --- a/src/Trixi.jl +++ b/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, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index a9e0cf87b0a..de6b9a2a4a6 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/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, diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 350aee7336a..7367886ca94 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/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, diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index c6d32c0f6dc..41251506a0d 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/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 diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index c1d78cbbf1e..d74a0c0cea1 100644 --- a/src/meshes/mesh_io.jl +++ b/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"]) diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index ed2158b169a..4d6016e5564 100644 --- a/src/meshes/meshes.jl +++ b/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") diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl new file mode 100644 index 00000000000..bd55115cc90 --- /dev/null +++ b/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 diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 33803752c2e..cc629d1a674 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -316,7 +316,8 @@ function save_mesh(semi::SemidiscretizationCoupled, output_directory, timestep = mesh, _, _, _ = mesh_equations_solver_cache(semi.semis[i]) if mesh.unsaved_changes - mesh.current_filename = save_mesh_file(mesh, output_directory, system = i) + mesh.current_filename = save_mesh_file(mesh, output_directory; system = i, + timestep = timestep) mesh.unsaved_changes = false end end @@ -630,7 +631,9 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionCoupled, - mesh::StructuredMesh, equations, + mesh::Union{StructuredMesh, + StructuredMeshView}, + equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @@ -662,7 +665,8 @@ end end end -function get_boundary_indices(element, orientation, mesh::StructuredMesh{2}) +function get_boundary_indices(element, orientation, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}) cartesian_indices = CartesianIndices(size(mesh)) if orientation == 1 # Get index of element in y-direction diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index f61378a7dca..dcd211671c8 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -259,7 +259,8 @@ function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{1}, end function check_periodicity_mesh_boundary_conditions(mesh::Union{TreeMesh{2}, - StructuredMesh{2}}, + StructuredMesh{2}, + StructuredMeshView{2}}, boundary_conditions::Union{NamedTuple, Tuple}) check_periodicity_mesh_boundary_conditions_x(mesh, boundary_conditions[1], diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 9b61df62cc3..0ab947e697a 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -436,8 +436,9 @@ function get_node_variables!(node_variables, mesh, equations, dg::DG, cache) get_node_variables!(node_variables, mesh, equations, dg.volume_integral, dg, cache) end -const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, - T8codeMesh} +const MeshesDGSEM = Union{TreeMesh, StructuredMesh, StructuredMeshView, + UnstructuredMesh2D, + P4estMesh, T8codeMesh} @inline function ndofs(mesh::MeshesDGSEM, dg::DG, cache) nelements(cache.elements) * nnodes(dg)^ndims(mesh) diff --git a/src/solvers/dgsem_structured/containers.jl b/src/solvers/dgsem_structured/containers.jl index 8adf005b782..7b0d275c5b5 100644 --- a/src/solvers/dgsem_structured/containers.jl +++ b/src/solvers/dgsem_structured/containers.jl @@ -23,7 +23,8 @@ struct ElementContainer{NDIMS, RealT <: Real, uEltype <: Real, NDIMSP1, NDIMSP2, end # Create element container and initialize element data -function init_elements(mesh::StructuredMesh{NDIMS, RealT}, +function init_elements(mesh::Union{StructuredMesh{NDIMS, RealT}, + StructuredMeshView{NDIMS, RealT}}, equations::AbstractEquations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index fb6db48e0a5..8a0722fc5d5 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -6,7 +6,8 @@ #! format: noindent # Initialize data structures in element container -function init_elements!(elements, mesh::StructuredMesh{2}, basis::LobattoLegendreBasis) +function init_elements!(elements, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, + basis::LobattoLegendreBasis) @unpack node_coordinates, left_neighbors, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -148,7 +149,9 @@ function calc_inverse_jacobian!(inverse_jacobian::AbstractArray{<:Any, 3}, eleme end # Save id of left neighbor of every element -function initialize_left_neighbor_connectivity!(left_neighbors, mesh::StructuredMesh{2}, +function initialize_left_neighbor_connectivity!(left_neighbors, + mesh::Union{StructuredMesh{2}, + StructuredMeshView{2}}, linear_indices) # Neighbors in x-direction for cell_y in 1:size(mesh, 2) diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl index 5cf4c4ef78c..00e321fba65 100644 --- a/src/solvers/dgsem_structured/dg.jl +++ b/src/solvers/dgsem_structured/dg.jl @@ -8,7 +8,8 @@ # 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::StructuredMesh, equations::AbstractEquations, dg::DG, ::Any, +function create_cache(mesh::Union{StructuredMesh, StructuredMeshView}, + equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype <: Real} elements = init_elements(mesh, equations, dg.basis, uEltype) @@ -30,7 +31,9 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition::BoundaryConditionPeriodic, - mesh::StructuredMesh, equations, + mesh::Union{StructuredMesh, + StructuredMeshView}, + equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) @@ -40,7 +43,9 @@ end @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation, boundary_condition, - mesh::StructuredMesh, equations, + mesh::Union{StructuredMesh, + StructuredMeshView}, + equations, surface_integral, dg::DG, cache, direction, node_indices, surface_node_indices, element) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 25a0eea096f..1eb244e06aa 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -6,7 +6,7 @@ #! format: noindent function rhs!(du, u, t, - mesh::StructuredMesh{2}, equations, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -58,8 +58,9 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 =# @inline function weak_form_kernel!(du, u, element, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -387,7 +388,7 @@ end end function calc_interface_flux!(cache, u, - mesh::StructuredMesh{2}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, nonconservative_terms, # can be True/False equations, surface_integral, dg::DG) @unpack elements = cache @@ -416,7 +417,8 @@ end @inline function calc_interface_flux!(surface_flux_values, left_element, right_element, orientation, u, - mesh::StructuredMesh{2}, + mesh::Union{StructuredMesh{2}, + StructuredMeshView{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache) # This is slow for LSA, but for some reason faster for Euler (see #519) @@ -551,13 +553,14 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, u, t, boundary_condition::BoundaryConditionPeriodic, - mesh::StructuredMesh{2}, equations, surface_integral, - dg::DG) + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, + equations, surface_integral, dg::DG) @assert isperiodic(mesh) end function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, - mesh::StructuredMesh{2}, equations, surface_integral, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}}, + equations, surface_integral, dg::DG) @unpack surface_flux_values = cache.elements linear_indices = LinearIndices(size(mesh)) @@ -616,8 +619,8 @@ function calc_boundary_flux!(cache, u, t, boundary_conditions::NamedTuple, end function apply_jacobian!(du, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, equations, dg::DG, cache) @unpack inverse_jacobian = cache.elements diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 547ed352ef3..6b66c2d4bfa 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -180,8 +180,8 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + StructuredMeshView{2}, UnstructuredMesh2D, + P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -1085,7 +1085,9 @@ end return nothing end -function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, +function calc_surface_integral!(du, u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, + StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DG, cache) @unpack boundary_interpolation = dg.basis diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index c3a792a0ffd..f095c97b19e 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -70,6 +70,30 @@ end end end +@trixi_testset "elixir_advection_meshview.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_meshview.jl"), + l2=[ + 8.311947673083206e-6, + 8.311947673068427e-6, + ], + linf=[ + 6.627000273318195e-5, + 6.62700027264096e-5, + ], + coverage_override=(maxiters = 10^5,)) + + @testset "analysis_callback(sol) for AnalysisCallbackCoupled" begin + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + @trixi_testset "elixir_advection_extended.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"), l2=[4.220397559713772e-6],