diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl new file mode 100644 index 00000000000..5e942f0ccd1 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -0,0 +1,80 @@ +# The same setup as tree_2d_dgsem/elixir_advection_basic.jl +# to verify the StructuredMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +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 P4estMesh with 8 x 8 trees and 16 x 16 elements and the mesh views +parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + initial_refinement_level = 1, periodicity = (false, true)) +mesh1 = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (4, 8), + coordinates_min = coordinates_min, coordinates_max = (0.0, 1.0), + periodicity = (false, true)) +mesh2 = P4estMeshView(parent_mesh; indices_min = (5, 1), indices_max = (8, 8), + coordinates_min = (0.0, -1.0), coordinates_max = coordinates_max, + periodicity = (false, true)) + +boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_convergence_test), + :x_pos => BoundaryConditionDirichlet(initial_condition_convergence_test)) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, + solver, boundary_conditions = boundary_conditions) +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, + solver, boundary_conditions = boundary_conditions) + +# Create a semidiscretization that bundles semi1 and semi2 +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) + +# 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, 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 = 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 d9b9245918e..7a29313530f 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -228,7 +228,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 4cd92ce7c5f..89ac43dba4f 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -61,7 +61,9 @@ end function create_cache_analysis(analyzer, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, T8codeMesh{2}}, + UnstructuredMesh2D, + P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -138,7 +140,8 @@ 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 @@ -208,6 +211,7 @@ end function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @@ -237,7 +241,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 @@ -247,7 +252,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} @@ -266,7 +271,8 @@ end function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, @@ -311,7 +317,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, @@ -379,7 +385,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 diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 57cdb3ef8a2..30176e8996d 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -9,7 +9,7 @@ function save_solution_file(u, time, dt, timestep, mesh::Union{SerialTreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, SerialP4estMesh, - SerialT8codeMesh}, + P4estMeshView, SerialT8codeMesh}, equations, dg::DG, cache, solution_callback, element_variables = Dict{Symbol, Any}(), diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index c7922cecc66..6f8f0ba7ec2 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -79,7 +79,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}, StructuredMeshView{2}}, + T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{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 @@ -115,7 +115,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}, StructuredMeshView{2}}, + T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}}, constant_speed::True, equations, dg::DG, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements 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.jl b/src/meshes/p4est_mesh.jl index 65c0a431b29..eac84ed2f01 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -107,6 +107,7 @@ end @inline Base.ndims(::P4estMesh{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} = RealT @inline ndims_ambient(::P4estMesh{NDIMS, NDIMS_AMBIENT}) where {NDIMS, NDIMS_AMBIENT} = NDIMS_AMBIENT +@inline Base.size(mesh::P4estMesh) = size(mesh.tree_node_coordinates)[end] @inline function ntrees(mesh::P4estMesh) return mesh.p4est.trees.elem_count[] diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl new file mode 100644 index 00000000000..e6110f20369 --- /dev/null +++ b/src/meshes/p4est_mesh_view.jl @@ -0,0 +1,200 @@ +@muladd begin +#! format: noindent + +""" + P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, + NNODES} <: + AbstractMesh{NDIMS} + +A view on a p4est mesh. +""" +mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, IsParallel, P, Ghost, + NDIMSP2, + NNODES} <: + AbstractMesh{NDIMS} + # Attributes from the original P4est mesh + p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t} + is_parallel :: IsParallel + ghost :: Ghost # Either PointerWrapper{p4est_ghost_t} or PointerWrapper{p8est_ghost_t} + # Coordinates at the nodes specified by the tensor product of `nodes` (NDIMS times). + # This specifies the geometry interpolation for each tree. + tree_node_coordinates::Array{RealT, NDIMSP2} # [dimension, i, j, k, tree] + nodes::SVector{NNODES, RealT} + boundary_names::Array{Symbol, 2} # [face direction, tree] + current_filename::String + unsaved_changes::Bool + p4est_partition_allow_for_coarsening::Bool + + # Attributes pertaining the views. + parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT} + indices_min::NTuple{NDIMS, Int} + indices_max::NTuple{NDIMS, Int} +end + +# function P4estMeshView(parent::P4estMesh{NDIMS, RealT}, view_cells::Array{Bool}) where {NDIMS, RealT} +function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS, RealT}; + indices_min = ntuple(_ -> 1, Val(NDIMS)), + indices_max = size(parent), + coordinates_min = nothing, coordinates_max = nothing, + periodicity = (true, true)) where {NDIMS, RealT} + trees_per_dimension = indices_max .- indices_min .+ (1, 1) + + @assert indices_min <= indices_max + @assert all(indices_min .> 0) + @assert prod(trees_per_dimension) <= size(parent) + + ghost = ghost_new_p4est(parent.p4est) + ghost_pw = PointerWrapper(ghost) + + # Extract mapping + if isnothing(coordinates_min) + coordinates_min = minimum(parent.tree_node_coordinates) + end + if isnothing(coordinates_max) + coordinates_max = maximum(parent.tree_node_coordinates) + end + + mapping = coordinates2mapping(coordinates_min, coordinates_max) + + tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, + ntuple(_ -> length(parent.nodes), + NDIMS)..., + prod(trees_per_dimension)) + calc_tree_node_coordinates!(tree_node_coordinates, parent.nodes, mapping, + trees_per_dimension) + + connectivity = connectivity_structured(trees_per_dimension..., periodicity) + + # TODO: The initial refinement level of 1 should not be hard-coded. + p4est = new_p4est(connectivity, 1) + p4est_pw = PointerWrapper(p4est) + + # Non-periodic boundaries + boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) + + structured_boundary_names!(boundary_names, trees_per_dimension, periodicity) + + return P4estMeshView{NDIMS, NDIMS, eltype(parent.tree_node_coordinates), + typeof(parent.is_parallel), + typeof(p4est_pw), typeof(parent.ghost), NDIMS + 2, + length(parent.nodes)}(PointerWrapper(p4est), + parent.is_parallel, + parent.ghost, + tree_node_coordinates, + parent.nodes, boundary_names, + parent.current_filename, + parent.unsaved_changes, + parent.p4est_partition_allow_for_coarsening, + parent, indices_min, indices_max) +end + +@inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS +@inline Base.real(::P4estMeshView{NDIMS, NDIMS, RealT}) where {NDIMS, RealT} = RealT +@inline ndims_ambient(::P4estMesh{NDIMS, NDIMS}) where {NDIMS} = NDIMS +@inline function ntrees(mesh::P4estMeshView) + return mesh.p4est.trees.elem_count[] +end +@inline ncellsglobal(mesh::P4estMeshView) = Int(mesh.p4est.global_num_quadrants[]) + +function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMeshView) + if get(io, :compact, false) + show(io, mesh) + else + setup = [ + "#trees" => ntrees(mesh), + "current #cells" => ncellsglobal(mesh), + "polydeg" => length(mesh.nodes) - 1 + ] + summary_box(io, + "P4estMeshView{" * string(ndims(mesh)) * ", " * string(real(mesh)) * + "}", setup) + end +end + +function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL) + p4est_balance(mesh.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.p4est, P4EST_CONNECT_FACE, init_fn) +end + +@inline ncells(mesh::P4estMeshView) = Int(mesh.p4est.local_num_quadrants[]) + +function calc_node_coordinates!(node_coordinates, + mesh::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 + # than native `Array`s. + tmp1 = StrideArray(undef, real(mesh), + StaticInt(2), static_length(nodes), static_length(mesh.nodes)) + matrix1 = StrideArray(undef, real(mesh), + static_length(nodes), static_length(mesh.nodes)) + matrix2 = similar(matrix1) + baryweights_in = barycentric_weights(mesh.nodes) + + # Macros from `p4est` + p4est_root_len = 1 << P4EST_MAXLEVEL + p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) + + trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) + + for tree in eachindex(trees) + offset = trees[tree].quadrants_offset + quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants) + + for i in eachindex(quadrants) + element = offset + i + quad = quadrants[i] + + quad_length = p4est_quadrant_len(quad.level) / p4est_root_len + + nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.x / p4est_root_len) .- 1 + nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.y / p4est_root_len) .- 1 + polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, element), + matrix1, matrix2, + view(mesh.tree_node_coordinates, :, :, :, tree), + tmp1) + end + end + + return node_coordinates +end + +function save_mesh_file(mesh::P4estMeshView, output_directory; timestep = 0, + system = "") + # Create output directory (if it does not exist) + mkpath(output_directory) + + filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) + p4est_filename = @sprintf("p4est_data_%s", system) + + p4est_file = joinpath(output_directory, p4est_filename) + + # Save the complete connectivity and `p4est` data to disk. + save_p4est!(p4est_file, mesh.p4est) + + # Open file (clobber existing content) + h5open(filename, "w") do file + # Add context information as attributes + attributes(file)["mesh_type"] = get_name(mesh.parent) + attributes(file)["ndims"] = ndims(mesh) + attributes(file)["p4est_file"] = p4est_filename + + file["tree_node_coordinates"] = mesh.tree_node_coordinates + file["nodes"] = Vector(mesh.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.boundary_names .|> String + end + + return filename +end +end # @muladd diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 745a8d3f6f8..d8d0afda5a2 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -514,10 +514,17 @@ end function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization) n_boundaries = 2 * ndims(semi) mesh, equations, solver, _ = mesh_equations_solver_cache(semi) + boundary_dictionary_names = [:x_neg, :x_pos, :y_neg, :y_pos] for direction in 1:n_boundaries - boundary_condition = semi.boundary_conditions[direction] - + boundary_condition = nothing + if typeof(semi.boundary_conditions) <: Trixi.UnstructuredSortedBoundaryTypes + if boundary_dictionary_names[direction] in keys(semi.boundary_conditions.boundary_dictionary) + boundary_condition = semi.boundary_conditions.boundary_dictionary[boundary_dictionary_names[direction]] + end + else + boundary_condition = semi.boundary_conditions[direction] + end allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, solver) diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index ffe74af14cf..3e8c37d4a0c 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 0e4d667fbc2..31300a44e46 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -435,7 +435,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..397ac54a386 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) @@ -514,7 +518,8 @@ function init_surfaces_iter_face_inner(info, user_data) return nothing end -function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) +function init_surfaces!(interfaces, mortars, boundaries, + mesh::Union{P4estMesh, P4estMeshView}) # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh) @@ -689,7 +694,7 @@ function cfunction(::typeof(count_surfaces_iter_face), ::Val{3}) (Ptr{p8est_iter_face_info_t}, Ptr{Cvoid})) end -function count_required_surfaces(mesh::P4estMesh) +function count_required_surfaces(mesh::Union{P4estMesh, P4estMeshView}) # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh))) 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..2d90b4a87ea 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/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::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 diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 17b9af04467..4d9594eab1d 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)}, @@ -58,7 +59,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 +115,8 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -182,7 +184,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}, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -206,7 +209,8 @@ 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, @@ -247,7 +251,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) @@ -276,7 +280,7 @@ function prolong2boundaries!(cache, u, end function calc_boundary_flux!(cache, t, boundary_condition::BC, boundary_indexing, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) where {BC} @unpack boundaries = cache @unpack surface_flux_values = cache.elements @@ -312,7 +316,8 @@ end # inlined version of the boundary flux calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -343,7 +348,8 @@ end # inlined version of the boundary flux with nonconservative terms calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -385,7 +391,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, surface_integral, dg::DGSEM) @unpack neighbor_ids, node_indices = cache.mortars @@ -452,7 +459,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) @@ -626,7 +633,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_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 39bece53ca9..22088a59278 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -60,6 +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}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) @@ -625,7 +626,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 03f28659f5e..176aa2fa007 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -38,14 +38,15 @@ end # and called from the basic `create_cache` method at the top. function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{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[] @@ -71,7 +72,8 @@ 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} @@ -93,7 +95,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 +113,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 +185,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) @@ -237,7 +241,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) @@ -334,7 +338,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) @@ -396,7 +400,7 @@ end @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 6ef65c4e9da..9b76da4ce06 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -98,7 +98,8 @@ end end # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha -function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, alpha, +function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, alpha, alpha_tmp, dg, cache) # Copy alpha values such that smoothing is indpedenent of the element access order diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 1f8c1ffa16a..3cad9fb9bbc 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -314,14 +314,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 @@ -335,7 +337,7 @@ end function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, BC_indices::NTuple{N, Vector{Int}}, mesh::Union{UnstructuredMesh2D, P4estMesh, - T8codeMesh}, + P4estMeshView, T8codeMesh}, equations, surface_integral, dg::DG) where {N} # Extract the boundary condition type and index vector boundary_condition = first(BCs) @@ -359,7 +361,7 @@ end # terminate the type-stable iteration over tuples function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, mesh::Union{UnstructuredMesh2D, P4estMesh, - T8codeMesh}, + P4estMeshView, T8codeMesh}, equations, surface_integral, dg::DG) nothing end diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index cba79b58b70..8e905dcf63f 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -125,6 +125,20 @@ end end end +@trixi_testset "elixir_advection_p4estview.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4estview.jl"), + l2=[8.286821350926049e-6, 8.286821351032878e-6], + linf=[6.674634841197236e-5, 6.674634841874472e-5]) + # 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"),