diff --git a/examples/p4est_2d_dgsem/elixir_advection_meshview.jl b/examples/p4est_2d_dgsem/elixir_advection_meshview.jl new file mode 100644 index 00000000000..9435930b538 --- /dev/null +++ b/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() diff --git a/src/Trixi.jl b/src/Trixi.jl index f8cc3a1a4d3..2fbf48bd6ba 100644 --- a/src/Trixi.jl +++ b/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, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 658778828b9..c538d444f5a 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/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 diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index ee18ee1acbd..091318b74b9 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/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, diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index b895dcd5b5d..1ce1385dd6a 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/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 diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index a9292c0b2ff..6c8a743efd3 100644 --- a/src/meshes/mesh_io.jl +++ b/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"]), diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 4d6016e5564..2a4c806ef94 100644 --- a/src/meshes/meshes.jl +++ b/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") diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl new file mode 100644 index 00000000000..4c1546de06e --- /dev/null +++ b/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 diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 02cff56a1d0..c909196b5db 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/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}, diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index a0fb0d95079..20b989da334 100644 --- a/src/solvers/dg.jl +++ b/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) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 3ef9cb2a421..a070db6b701 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -82,6 +82,7 @@ end # Create element container and initialize element data function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT}, + P4estMeshView{NDIMS, NDIMS, RealT}, T8codeMesh{NDIMS, RealT}}, equations, basis, @@ -167,7 +168,8 @@ 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) @@ -197,7 +199,7 @@ function init_interfaces(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, e return interfaces end -function init_interfaces!(interfaces, mesh::P4estMesh) +function init_interfaces!(interfaces, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(interfaces, nothing, nothing, mesh) return interfaces @@ -242,7 +244,8 @@ function Base.resize!(boundaries::P4estBoundaryContainer, capacity) end # Create interface container and initialize interface data in `elements`. -function init_boundaries(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elements) +function init_boundaries(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, + basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -271,7 +274,7 @@ function init_boundaries(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, e return boundaries end -function init_boundaries!(boundaries, mesh::P4estMesh) +function init_boundaries!(boundaries, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(nothing, nothing, boundaries, mesh) return boundaries @@ -373,7 +376,8 @@ function Base.resize!(mortars::P4estMortarContainer, capacity) end # Create mortar container and initialize mortar data. -function init_mortars(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elements) +function init_mortars(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, + basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -408,7 +412,7 @@ function init_mortars(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elem return mortars end -function init_mortars!(mortars, mesh::P4estMesh) +function init_mortars!(mortars, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(nothing, mortars, nothing, mesh) return mortars diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 6af6fd6d90e..04a2b9f9ecd 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -6,7 +6,8 @@ #! 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 +27,8 @@ 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. diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 8197ad4a2d0..d4a1b2d3a50 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -29,6 +29,39 @@ function create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::A return cache end +# 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::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 + balance!(mesh.parent) + + elements_parent = init_elements(mesh.parent, equations, dg.basis, uEltype) + interfaces_parent = init_interfaces(mesh.parent, equations, dg.basis, + elements_parent) + boundaries_parent = init_boundaries(mesh.parent, equations, dg.basis, + elements_parent) + mortars_parent = init_mortars(mesh.parent, equations, dg.basis, elements_parent) + + # Extract data for views. + elements, interfaces, boundaries, mortars = extract_p4est_mesh_view(elements_parent, + interfaces_parent, + boundaries_parent, + mortars_parent, + mesh) + + cache = (; elements, interfaces, boundaries, mortars) + + # Add specialized parts of the cache required to compute the volume integral etc. + cache = (; cache..., + create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) + cache = (; cache..., create_cache(mesh, equations, dg.mortar, uEltype)...) + + return cache +end + # Extract outward-pointing normal direction # (contravariant vector ±Ja^i, i = index) # Note that this vector is not normalized diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 6a7bbbb71d7..b65fd9b56d9 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -7,7 +7,8 @@ # 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)}, @@ -62,7 +63,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) @@ -118,7 +119,8 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -186,7 +188,8 @@ 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}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -245,7 +248,7 @@ end end function prolong2boundaries!(cache, u, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack boundaries = cache index_range = eachnode(dg) @@ -382,7 +385,8 @@ end end function prolong2mortars!(cache, u, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM) @unpack neighbor_ids, node_indices = cache.mortars @@ -449,7 +453,7 @@ function prolong2mortars!(cache, u, end function calc_mortar_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -635,7 +639,8 @@ end end function calc_surface_integral!(du, u, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index 9f8306031c1..2f84510df6a 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -520,7 +520,7 @@ end # This is the version used when calculating the divergence of the viscous fluxes # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache_parabolic, flux_viscous, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) (; interfaces) = cache_parabolic diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index cb4e75149ea..2c56d0bceb5 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -60,7 +60,7 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 element, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] @@ -619,7 +619,8 @@ end function apply_jacobian!(du, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{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 98bd0889206..be5d5b5600a 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -90,7 +90,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, @@ -111,7 +111,8 @@ 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, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -182,7 +183,8 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 48d4fe153c6..d3ee2f019c2 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -310,14 +310,16 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{UnstructuredMesh2D, P4estMesh, T8codeMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, + T8codeMesh}, equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end # Function barrier for type stability function calc_boundary_flux!(cache, t, boundary_conditions, - mesh::Union{UnstructuredMesh2D, P4estMesh, T8codeMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, + T8codeMesh}, equations, surface_integral, dg::DG) @unpack boundary_condition_types, boundary_indices = boundary_conditions diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 0e70282e77c..9541c36c7a9 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -125,6 +125,20 @@ end end end +@trixi_testset "elixir_advection_meshview.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_meshview.jl"), + l2=[0.00013773915040249946], + linf=[0.0010140184322192658]) + # 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 + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),