From 0a3a59f67c3bf887f36b4760eae59ec5dde993c2 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 14 Oct 2024 14:13:08 +0100 Subject: [PATCH 01/42] Sterted re-implementing p4est mesh views to be more powerful. --- src/Trixi.jl | 2 +- src/meshes/meshes.jl | 1 + src/solvers/dgsem_p4est/containers.jl | 1 + src/solvers/dgsem_p4est/dg.jl | 2 +- 4 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 0cedec782df..7777d2089e1 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/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/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 3ef9cb2a421..310026e6970 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,RealT}, T8codeMesh{NDIMS, RealT}}, equations, basis, diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 8197ad4a2d0..b6259ff0207 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -8,7 +8,7 @@ # 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 From 6e1043422ff72370cc2d1138454f03b5e3fe1541 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 14 Oct 2024 20:36:09 +0100 Subject: [PATCH 02/42] Added P4estMeshView specific methods for computing calc_node_coordinates. --- src/solvers/dgsem_p4est/containers.jl | 2 +- src/solvers/dgsem_p4est/containers_2d.jl | 66 +++++++++++++++++++++++- 2 files changed, 66 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 310026e6970..55859cf92d1 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -82,7 +82,7 @@ end # Create element container and initialize element data function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT}, - P4estMeshView{NDIMS,RealT}, + P4estMeshView{NDIMS, NDIMS, RealT}, T8codeMesh{NDIMS, RealT}}, equations, basis, diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 6af6fd6d90e..01a94a94094 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -6,7 +6,7 @@ #! 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 @@ -35,6 +35,17 @@ function calc_node_coordinates!(node_coordinates, calc_node_coordinates!(node_coordinates, mesh, basis.nodes) end +# Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis +function calc_node_coordinates!(node_coordinates, + mesh::P4estMeshView{2}, + basis::LobattoLegendreBasis) + # Hanging nodes will cause holes in the mesh if its polydeg is higher + # than the polydeg of the solver. + @assert length(basis.nodes)>=length(mesh.parent.nodes) "The solver can't have a lower polydeg than the mesh" + + calc_node_coordinates!(node_coordinates, mesh, basis.nodes) +end + # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, mesh::P4estMesh{2, NDIMS_AMBIENT}, @@ -85,6 +96,59 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +# Interpolate tree_node_coordinates to each quadrant at the specified nodes +function calc_node_coordinates!(node_coordinates, + mesh::P4estMeshView{2, NDIMS_AMBIENT}, + nodes::AbstractVector) where {NDIMS_AMBIENT} + # 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(NDIMS_AMBIENT), static_length(nodes), + static_length(mesh.parent.nodes)) + matrix1 = StrideArray(undef, real(mesh), + static_length(nodes), static_length(mesh.parent.nodes)) + matrix2 = similar(matrix1) + baryweights_in = barycentric_weights(mesh.parent.nodes) + + # Macros from `p4est` + p4est_root_len = 1 << P4EST_MAXLEVEL + p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) + + # SC: It seems that here we need to make some serious changes so it works + # with P4estMeshView without toching the p4est struct. + trees = unsafe_wrap_sc(p4est_tree_t, mesh.parent.p4est.trees)[mesh.cell_ids] + @autoinfiltrate + + 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.parent.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, mesh.parent.nodes, nodes_out_y, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, element), + matrix1, matrix2, + view(mesh.parent.tree_node_coordinates, :, :, :, tree), + tmp1) + end + end + + return node_coordinates +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{2}, faces, orientation, interface_id) From ac9b4175723999fb6326036b0cd2352204ee077f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 15 Oct 2024 09:07:50 +0100 Subject: [PATCH 03/42] Added p4est mesh views. --- src/meshes/p4est_mesh_view.jl | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 src/meshes/p4est_mesh_view.jl diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl new file mode 100644 index 00000000000..c408090523b --- /dev/null +++ b/src/meshes/p4est_mesh_view.jl @@ -0,0 +1,35 @@ +@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, Parent} <: AbstractMesh{NDIMS} + parent::Parent + cell_ids::Vector{Int} +end + +function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} + # SC: number of cells should be corrected. + @autoinfiltrate + cell_ids = Vector{Int}(undef, ncells(parent)) + # SC: do not populate this array. It needs to be given by the user. + for i in 1:ncells(parent) + cell_ids[i] = i + end + + return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids) +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) + +end # @muladd From 6f50db4ffa88a4766af954894f190e5f52d9f809 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 15 Oct 2024 11:23:24 +0100 Subject: [PATCH 04/42] Removed some debugging. --- src/meshes/p4est_mesh_view.jl | 1 - src/solvers/dgsem_p4est/containers_2d.jl | 1 - 2 files changed, 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index c408090523b..e82d0fa7e42 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -15,7 +15,6 @@ end function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} # SC: number of cells should be corrected. - @autoinfiltrate cell_ids = Vector{Int}(undef, ncells(parent)) # SC: do not populate this array. It needs to be given by the user. for i in 1:ncells(parent) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 01a94a94094..373645afe66 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -118,7 +118,6 @@ function calc_node_coordinates!(node_coordinates, # SC: It seems that here we need to make some serious changes so it works # with P4estMeshView without toching the p4est struct. trees = unsafe_wrap_sc(p4est_tree_t, mesh.parent.p4est.trees)[mesh.cell_ids] - @autoinfiltrate for tree in eachindex(trees) offset = trees[tree].quadrants_offset From 7056c4684a7cb1f792ac0a3821793a5ae7a9b449 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 15 Oct 2024 14:23:58 +0100 Subject: [PATCH 05/42] Corrected matrix calculation in calc_node_coordinates for p4est mesh view. --- src/solvers/dgsem_p4est/containers_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 373645afe66..065c92db852 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -106,7 +106,7 @@ function calc_node_coordinates!(node_coordinates, tmp1 = StrideArray(undef, real(mesh), StaticInt(NDIMS_AMBIENT), static_length(nodes), static_length(mesh.parent.nodes)) - matrix1 = StrideArray(undef, real(mesh), + matrix1 = StrideArray(undef, real(mesh.parent), static_length(nodes), static_length(mesh.parent.nodes)) matrix2 = similar(matrix1) baryweights_in = barycentric_weights(mesh.parent.nodes) @@ -116,7 +116,7 @@ function calc_node_coordinates!(node_coordinates, p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) # SC: It seems that here we need to make some serious changes so it works - # with P4estMeshView without toching the p4est struct. + # with P4estMeshView without touching the p4est struct. trees = unsafe_wrap_sc(p4est_tree_t, mesh.parent.p4est.trees)[mesh.cell_ids] for tree in eachindex(trees) From 805ea7196f883e4c7d06d2dde18822e2157af94e Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 16 Oct 2024 16:19:12 +0100 Subject: [PATCH 06/42] Plae holders for p4est mesh view. Thius version runs, but does not have all features. --- src/callbacks_step/analysis_dg2d.jl | 13 +- src/callbacks_step/stepsize_dg2d.jl | 2 +- .../semidiscretization_hyperbolic.jl | 2 +- src/solvers/dg.jl | 2 +- src/solvers/dgsem_p4est/containers.jl | 124 ++++++++++++++++++ src/solvers/dgsem_p4est/dg_2d.jl | 2 +- src/solvers/dgsem_tree/dg_2d.jl | 71 +++++++++- 7 files changed, 210 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 4cd92ce7c5f..06c23f3ceba 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 # 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} @@ -246,6 +247,16 @@ function integrate(func::Func, u, end end +function integrate(func::Func, u, + mesh::P4estMeshView{2}, + equations, dg::DG, cache; normalize = true) where {Func} + integrate_via_indices(u, mesh.parent, equations, dg, cache; + normalize = normalize) do u, i, j, element, equations, dg + u_local = get_node_vars(u, equations, dg, i, j, element) + return func(u_local, equations) + end +end + function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, equations_parabolic, diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index c7922cecc66..06eba6a11be 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/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index ffe74af14cf..54dc8761300 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -213,7 +213,7 @@ function digest_boundary_conditions(boundary_conditions::AbstractArray, mesh, so end # No checks for these meshes yet available -function check_periodicity_mesh_boundary_conditions(mesh::Union{P4estMesh, +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 55859cf92d1..eb7c57bb184 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -198,6 +198,37 @@ function init_interfaces(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, e return interfaces end +# Create interface container and initialize interface data. +function init_interfaces(mesh::P4estMeshView, equations, basis, elements) + NDIMS = ndims(elements) + uEltype = eltype(elements) + + # Initialize container + n_interfaces = count_required_surfaces(mesh).interfaces + + _u = Vector{uEltype}(undef, + 2 * nvariables(equations) * nnodes(basis)^(NDIMS - 1) * + n_interfaces) + u = unsafe_wrap(Array, pointer(_u), + (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., + n_interfaces)) + + _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces) + neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces)) + + _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces) + node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces)) + + interfaces = P4estInterfaceContainer{NDIMS, uEltype, NDIMS + 2}(u, neighbor_ids, + node_indices, + _u, _neighbor_ids, + _node_indices) + + init_interfaces!(interfaces, mesh.parent) + + return interfaces +end + function init_interfaces!(interfaces, mesh::P4estMesh) init_surfaces!(interfaces, nothing, nothing, mesh) @@ -272,12 +303,48 @@ function init_boundaries(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, e return boundaries end +# Create interface container and initialize interface data in `elements`. +function init_boundaries(mesh::P4estMeshView, equations, basis, elements) + NDIMS = ndims(elements) + uEltype = eltype(elements) + + # Initialize container + n_boundaries = count_required_surfaces(mesh.parent).boundaries + + _u = Vector{uEltype}(undef, + nvariables(equations) * nnodes(basis)^(NDIMS - 1) * + n_boundaries) + u = unsafe_wrap(Array, pointer(_u), + (nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., + n_boundaries)) + + neighbor_ids = Vector{Int}(undef, n_boundaries) + node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, n_boundaries) + names = Vector{Symbol}(undef, n_boundaries) + + boundaries = P4estBoundaryContainer{NDIMS, uEltype, NDIMS + 1}(u, neighbor_ids, + node_indices, names, + _u) + + if n_boundaries > 0 + init_boundaries!(boundaries, mesh) + end + + return boundaries +end + function init_boundaries!(boundaries, mesh::P4estMesh) init_surfaces!(nothing, nothing, boundaries, mesh) return boundaries end +function init_boundaries!(boundaries, mesh::P4estMeshView) + init_surfaces!(nothing, nothing, boundaries, mesh) + + return boundaries +end + # Function barrier for type stability function init_boundaries_iter_face_inner(info_pw, boundaries, boundary_id, mesh) # Extract boundary data @@ -409,12 +476,54 @@ function init_mortars(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elem return mortars end +# Create mortar container and initialize mortar data. +function init_mortars(mesh::P4estMeshView, equations, basis, elements) + NDIMS = ndims(elements) + uEltype = eltype(elements) + + # Initialize container + n_mortars = count_required_surfaces(mesh.parent).mortars + + _u = Vector{uEltype}(undef, + 2 * nvariables(equations) * 2^(NDIMS - 1) * + nnodes(basis)^(NDIMS - 1) * n_mortars) + u = unsafe_wrap(Array, pointer(_u), + (2, nvariables(equations), 2^(NDIMS - 1), + ntuple(_ -> nnodes(basis), NDIMS - 1)..., n_mortars)) + + _neighbor_ids = Vector{Int}(undef, (2^(NDIMS - 1) + 1) * n_mortars) + neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), + (2^(NDIMS - 1) + 1, n_mortars)) + + _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_mortars) + node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_mortars)) + + mortars = P4estMortarContainer{NDIMS, uEltype, NDIMS + 1, NDIMS + 3}(u, + neighbor_ids, + node_indices, + _u, + _neighbor_ids, + _node_indices) + + if n_mortars > 0 + init_mortars!(mortars, mesh) + end + + return mortars +end + function init_mortars!(mortars, mesh::P4estMesh) init_surfaces!(nothing, mortars, nothing, mesh) return mortars end +function init_mortars!(mortars, mesh::P4estMeshView) + init_surfaces!(nothing, mortars, nothing, mesh) + + return mortars +end + function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) # Re-initialize elements container @unpack elements = cache @@ -705,6 +814,21 @@ function count_required_surfaces(mesh::P4estMesh) boundaries = user_data[3]) end +function count_required_surfaces(mesh::P4estMeshView) + # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face + iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh))) + + # interfaces, mortars, boundaries + user_data = [0, 0, 0] + + iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) + + # Return counters + return (interfaces = user_data[1], + mortars = user_data[2], + boundaries = user_data[3]) +end + # Return direction of the face, which is indexed by node_indices @inline function indices2direction(indices) if indices[1] === :begin diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 3c868289181..4787cb3b113 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -7,7 +7,7 @@ # 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)}, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 03f28659f5e..ddfba31ed34 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -93,7 +93,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}}, + T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal performance using different types MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, @@ -179,6 +179,75 @@ function rhs!(du, u, t, return nothing end +function rhs!(du, u, t, + mesh::P4estMeshView{2}, equations, + boundary_conditions, source_terms::Source, + dg::DG, cache) where {Source} + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + + # Calculate volume integral + @trixi_timeit timer() "volume integral" begin + calc_volume_integral!(du, u, mesh.parent, + have_nonconservative_terms(equations), equations, + dg.volume_integral, dg, cache) + end + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(cache, u, mesh.parent, equations, + dg.surface_integral, dg) + end + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(cache.elements.surface_flux_values, mesh.parent, + have_nonconservative_terms(equations), equations, + dg.surface_integral, dg, cache) + end + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries" begin + prolong2boundaries!(cache, u, mesh.parent, equations, + dg.surface_integral, dg) + end + + # Calculate boundary fluxes + @trixi_timeit timer() "boundary flux" begin + calc_boundary_flux!(cache, t, boundary_conditions, mesh.parent, equations, + dg.surface_integral, dg) + end + + # Prolong solution to mortars + @trixi_timeit timer() "prolong2mortars" begin + prolong2mortars!(cache, u, mesh.parent, equations, + dg.mortar, dg.surface_integral, dg) + end + + # Calculate mortar fluxes + @trixi_timeit timer() "mortar flux" begin + calc_mortar_flux!(cache.elements.surface_flux_values, mesh.parent, + have_nonconservative_terms(equations), equations, + dg.mortar, dg.surface_integral, dg, cache) + end + + # Calculate surface integrals + @trixi_timeit timer() "surface integral" begin + calc_surface_integral!(du, u, mesh.parent, equations, + dg.surface_integral, dg, cache) + end + + # Apply Jacobian from mapping to reference element + @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh.parent, equations, dg, cache) + + # Calculate source terms + @trixi_timeit timer() "source terms" begin + calc_sources!(du, u, t, source_terms, equations, dg, cache) + end + + return nothing +end + function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, From c7664562360b392fa8fe7339a2121fd06ae13567 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 17 Oct 2024 08:52:51 +0100 Subject: [PATCH 07/42] Some refactoring. --- src/solvers/dgsem_p4est/containers.jl | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index eb7c57bb184..8ee7b50a39f 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -82,8 +82,8 @@ end # Create element container and initialize element data function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT}, - P4estMeshView{NDIMS, NDIMS, RealT}, - T8codeMesh{NDIMS, RealT}}, + P4estMeshView{NDIMS, NDIMS, RealT}, + T8codeMesh{NDIMS, RealT}}, equations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} @@ -333,13 +333,7 @@ function init_boundaries(mesh::P4estMeshView, equations, basis, elements) return boundaries end -function init_boundaries!(boundaries, mesh::P4estMesh) - init_surfaces!(nothing, nothing, boundaries, mesh) - - return boundaries -end - -function init_boundaries!(boundaries, mesh::P4estMeshView) +function init_boundaries!(boundaries, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(nothing, nothing, boundaries, mesh) return boundaries From ff167a74d5755f19959f06f91bc148b35886a76c Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 17 Oct 2024 13:37:46 +0100 Subject: [PATCH 08/42] Removed most references to parent mesh. This now goes through with mesh view being equal to the orginal mesh. --- src/callbacks_step/analysis_dg2d.jl | 12 +- src/solvers/dgsem_p4est/containers.jl | 123 +++------------------ src/solvers/dgsem_p4est/containers_2d.jl | 2 +- src/solvers/dgsem_p4est/dg_2d.jl | 14 +-- src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 2 +- src/solvers/dgsem_structured/dg_2d.jl | 5 +- src/solvers/dgsem_tree/dg_2d.jl | 73 +----------- src/solvers/dgsem_unstructured/dg_2d.jl | 2 +- 8 files changed, 31 insertions(+), 202 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 06c23f3ceba..03384151102 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -238,7 +238,7 @@ 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,16 +247,6 @@ function integrate(func::Func, u, end end -function integrate(func::Func, u, - mesh::P4estMeshView{2}, - equations, dg::DG, cache; normalize = true) where {Func} - integrate_via_indices(u, mesh.parent, equations, dg, cache; - normalize = normalize) do u, i, j, element, equations, dg - u_local = get_node_vars(u, equations, dg, i, j, element) - return func(u_local, equations) - end -end - function integrate(func::Func, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations, equations_parabolic, diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 8ee7b50a39f..2f534a2acf6 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -168,7 +168,7 @@ 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) @@ -198,38 +198,7 @@ function init_interfaces(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, e return interfaces end -# Create interface container and initialize interface data. -function init_interfaces(mesh::P4estMeshView, equations, basis, elements) - NDIMS = ndims(elements) - uEltype = eltype(elements) - - # Initialize container - n_interfaces = count_required_surfaces(mesh).interfaces - - _u = Vector{uEltype}(undef, - 2 * nvariables(equations) * nnodes(basis)^(NDIMS - 1) * - n_interfaces) - u = unsafe_wrap(Array, pointer(_u), - (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., - n_interfaces)) - - _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces) - neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces)) - - _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces) - node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces)) - - interfaces = P4estInterfaceContainer{NDIMS, uEltype, NDIMS + 2}(u, neighbor_ids, - node_indices, - _u, _neighbor_ids, - _node_indices) - - init_interfaces!(interfaces, mesh.parent) - - 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 @@ -274,7 +243,7 @@ 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) @@ -303,36 +272,6 @@ function init_boundaries(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, e return boundaries end -# Create interface container and initialize interface data in `elements`. -function init_boundaries(mesh::P4estMeshView, equations, basis, elements) - NDIMS = ndims(elements) - uEltype = eltype(elements) - - # Initialize container - n_boundaries = count_required_surfaces(mesh.parent).boundaries - - _u = Vector{uEltype}(undef, - nvariables(equations) * nnodes(basis)^(NDIMS - 1) * - n_boundaries) - u = unsafe_wrap(Array, pointer(_u), - (nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., - n_boundaries)) - - neighbor_ids = Vector{Int}(undef, n_boundaries) - node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, n_boundaries) - names = Vector{Symbol}(undef, n_boundaries) - - boundaries = P4estBoundaryContainer{NDIMS, uEltype, NDIMS + 1}(u, neighbor_ids, - node_indices, names, - _u) - - if n_boundaries > 0 - init_boundaries!(boundaries, mesh) - end - - return boundaries -end - function init_boundaries!(boundaries, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(nothing, nothing, boundaries, mesh) @@ -435,7 +374,7 @@ 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) @@ -470,49 +409,7 @@ function init_mortars(mesh::Union{P4estMesh, T8codeMesh}, equations, basis, elem return mortars end -# Create mortar container and initialize mortar data. -function init_mortars(mesh::P4estMeshView, equations, basis, elements) - NDIMS = ndims(elements) - uEltype = eltype(elements) - - # Initialize container - n_mortars = count_required_surfaces(mesh.parent).mortars - - _u = Vector{uEltype}(undef, - 2 * nvariables(equations) * 2^(NDIMS - 1) * - nnodes(basis)^(NDIMS - 1) * n_mortars) - u = unsafe_wrap(Array, pointer(_u), - (2, nvariables(equations), 2^(NDIMS - 1), - ntuple(_ -> nnodes(basis), NDIMS - 1)..., n_mortars)) - - _neighbor_ids = Vector{Int}(undef, (2^(NDIMS - 1) + 1) * n_mortars) - neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), - (2^(NDIMS - 1) + 1, n_mortars)) - - _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_mortars) - node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_mortars)) - - mortars = P4estMortarContainer{NDIMS, uEltype, NDIMS + 1, NDIMS + 3}(u, - neighbor_ids, - node_indices, - _u, - _neighbor_ids, - _node_indices) - - if n_mortars > 0 - init_mortars!(mortars, mesh) - end - - return mortars -end - -function init_mortars!(mortars, mesh::P4estMesh) - init_surfaces!(nothing, mortars, nothing, mesh) - - return mortars -end - -function init_mortars!(mortars, mesh::P4estMeshView) +function init_mortars!(mortars, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(nothing, mortars, nothing, mesh) return mortars @@ -628,6 +525,16 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) return interfaces end +function init_surfaces!(interfaces, mortars, boundaries, mesh::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.parent) + + iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) + + return interfaces +end + # Initialization of interfaces after the function barrier function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack interfaces, interface_id, mesh = user_data diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 065c92db852..a7a806fd621 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -26,7 +26,7 @@ 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_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 4787cb3b113..24a531c3724 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -58,7 +58,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 +114,7 @@ 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 @@ -182,7 +182,7 @@ 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, @@ -241,7 +241,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) @@ -378,7 +378,7 @@ 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 @@ -445,7 +445,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) @@ -618,7 +618,7 @@ 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 23968a2bfac..63fa50960ee 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -519,7 +519,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 ddfba31ed34..a8f302b1122 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -111,7 +111,7 @@ 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 @@ -179,79 +179,10 @@ function rhs!(du, u, t, return nothing end -function rhs!(du, u, t, - mesh::P4estMeshView{2}, equations, - boundary_conditions, source_terms::Source, - dg::DG, cache) where {Source} - # Reset du - @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) - - # Calculate volume integral - @trixi_timeit timer() "volume integral" begin - calc_volume_integral!(du, u, mesh.parent, - have_nonconservative_terms(equations), equations, - dg.volume_integral, dg, cache) - end - - # Prolong solution to interfaces - @trixi_timeit timer() "prolong2interfaces" begin - prolong2interfaces!(cache, u, mesh.parent, equations, - dg.surface_integral, dg) - end - - # Calculate interface fluxes - @trixi_timeit timer() "interface flux" begin - calc_interface_flux!(cache.elements.surface_flux_values, mesh.parent, - have_nonconservative_terms(equations), equations, - dg.surface_integral, dg, cache) - end - - # Prolong solution to boundaries - @trixi_timeit timer() "prolong2boundaries" begin - prolong2boundaries!(cache, u, mesh.parent, equations, - dg.surface_integral, dg) - end - - # Calculate boundary fluxes - @trixi_timeit timer() "boundary flux" begin - calc_boundary_flux!(cache, t, boundary_conditions, mesh.parent, equations, - dg.surface_integral, dg) - end - - # Prolong solution to mortars - @trixi_timeit timer() "prolong2mortars" begin - prolong2mortars!(cache, u, mesh.parent, equations, - dg.mortar, dg.surface_integral, dg) - end - - # Calculate mortar fluxes - @trixi_timeit timer() "mortar flux" begin - calc_mortar_flux!(cache.elements.surface_flux_values, mesh.parent, - have_nonconservative_terms(equations), equations, - dg.mortar, dg.surface_integral, dg, cache) - end - - # Calculate surface integrals - @trixi_timeit timer() "surface integral" begin - calc_surface_integral!(du, u, mesh.parent, equations, - dg.surface_integral, dg, cache) - end - - # Apply Jacobian from mapping to reference element - @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh.parent, equations, dg, cache) - - # Calculate source terms - @trixi_timeit timer() "source terms" begin - calc_sources!(du, u, t, source_terms, equations, dg, cache) - end - - return nothing -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..e5086a58db6 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -310,7 +310,7 @@ 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 From a5ce21045cd84f05e4a149db3470c46326ccc432 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 17 Oct 2024 14:11:20 +0100 Subject: [PATCH 09/42] Minor corrections. --- src/meshes/p4est_mesh_view.jl | 3 ++- src/solvers/dgsem_unstructured/dg_2d.jl | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index e82d0fa7e42..469f1c154e0 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -11,6 +11,7 @@ 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 end function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} @@ -21,7 +22,7 @@ function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {ND cell_ids[i] = i end - return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids) + return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids)#, parent.unsaved_changes) end @inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index e5086a58db6..2aa6ec24a42 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -317,7 +317,7 @@ 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 From 4b04cd20d79608f62b72e33037f09a9154755ee1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 18 Oct 2024 17:48:15 +0100 Subject: [PATCH 10/42] Added commentson next steps. --- src/meshes/p4est_mesh_view.jl | 4 ++++ src/solvers/dgsem_p4est/containers.jl | 2 ++ src/solvers/dgsem_p4est/containers_2d.jl | 2 -- 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 469f1c154e0..3bfa15cbc02 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -11,6 +11,8 @@ A view on a p4est mesh. mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: AbstractMesh{NDIMS} parent::Parent cell_ids::Vector{Int} + # SC: After some thought, we might need to create a p4est pointer to p4est data + # conatining the data from the view. # unsaved_changes::Bool end @@ -22,6 +24,8 @@ function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {ND cell_ids[i] = i end + # SC: Since we need a p4est pointer no the modified (view) p4est data, we might need a function + # like connectivity_structured that computes the connectivity. return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids)#, parent.unsaved_changes) end diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 2f534a2acf6..21a6c386ddd 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -530,6 +530,7 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView) iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh.parent) + # SC: After adding p4est ponter to the view we should change this. iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) return interfaces @@ -722,6 +723,7 @@ function count_required_surfaces(mesh::P4estMeshView) # interfaces, mortars, boundaries user_data = [0, 0, 0] + # SC: After adding p4est ponter to the view we should change this. iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) # Return counters diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index a7a806fd621..eb8d8cf85b5 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -115,8 +115,6 @@ function calc_node_coordinates!(node_coordinates, p4est_root_len = 1 << P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) - # SC: It seems that here we need to make some serious changes so it works - # with P4estMeshView without touching the p4est struct. trees = unsafe_wrap_sc(p4est_tree_t, mesh.parent.p4est.trees)[mesh.cell_ids] for tree in eachindex(trees) From a08b10fa12e1a5e9e4050647cbdb58a16a5db27b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 23 Oct 2024 19:42:10 +0100 Subject: [PATCH 11/42] Added real option of choosing elements of p4est mesh view. --- src/meshes/p4est_mesh_view.jl | 14 +++++++------- src/solvers/dgsem_p4est/containers.jl | 4 ++-- src/solvers/dgsem_p4est/containers_2d.jl | 1 + src/solvers/dgsem_p4est/dg_2d.jl | 1 + 4 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 3bfa15cbc02..31284942462 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -16,13 +16,13 @@ mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: Abs # unsaved_changes::Bool end -function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} - # SC: number of cells should be corrected. - cell_ids = Vector{Int}(undef, ncells(parent)) - # SC: do not populate this array. It needs to be given by the user. - for i in 1:ncells(parent) - cell_ids[i] = i - end +function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}, cell_ids::Vector) where {NDIMS, NDIMS_AMBIENT, RealT} +# # SC: number of cells should be corrected. +# cell_ids = Vector{Int}(undef, ncells(parent)) +# # SC: do not populate this array. It needs to be given by the user. +# for i in 1:ncells(parent) +# cell_ids[i] = i +# end # SC: Since we need a p4est pointer no the modified (view) p4est data, we might need a function # like connectivity_structured that computes the connectivity. diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 21a6c386ddd..853de3e47a7 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -530,7 +530,7 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView) iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh.parent) - # SC: After adding p4est ponter to the view we should change this. + # SC: After adding p4est pointer to the view we should change this. iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) return interfaces @@ -723,7 +723,7 @@ function count_required_surfaces(mesh::P4estMeshView) # interfaces, mortars, boundaries user_data = [0, 0, 0] - # SC: After adding p4est ponter to the view we should change this. + # SC: After adding p4est pointer to the view we should change this. iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) # Return counters diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index eb8d8cf85b5..d0220aaf6dd 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -143,6 +143,7 @@ function calc_node_coordinates!(node_coordinates, end end + @autoinfiltrate return node_coordinates end diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 24a531c3724..917f3feb996 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -64,6 +64,7 @@ function prolong2interfaces!(cache, u, index_range = eachnode(dg) @threaded for interface in eachinterface(dg, cache) + @autoinfiltrate # Copy solution data from the primary element using "delayed indexing" with # a start value and a step size to get the correct face and orientation. # Note that in the current implementation, the interface will be From 3de317a0800deba182d82f0002f720bad6995d09 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 6 Nov 2024 17:24:00 +0000 Subject: [PATCH 12/42] Little clean up. --- src/solvers/dgsem_p4est/containers.jl | 3 +-- src/solvers/dgsem_p4est/containers_2d.jl | 1 - src/solvers/dgsem_p4est/dg_2d.jl | 1 - 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 853de3e47a7..b88d7fd2886 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -530,7 +530,6 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView) iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh.parent) - # SC: After adding p4est pointer to the view we should change this. iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) return interfaces @@ -723,9 +722,9 @@ function count_required_surfaces(mesh::P4estMeshView) # interfaces, mortars, boundaries user_data = [0, 0, 0] - # SC: After adding p4est pointer to the view we should change this. iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) + @autoinfiltrate # Return counters return (interfaces = user_data[1], mortars = user_data[2], diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index d0220aaf6dd..eb8d8cf85b5 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -143,7 +143,6 @@ function calc_node_coordinates!(node_coordinates, end end - @autoinfiltrate return node_coordinates end diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 917f3feb996..24a531c3724 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -64,7 +64,6 @@ function prolong2interfaces!(cache, u, index_range = eachnode(dg) @threaded for interface in eachinterface(dg, cache) - @autoinfiltrate # Copy solution data from the primary element using "delayed indexing" with # a start value and a step size to get the correct face and orientation. # Note that in the current implementation, the interface will be From 0777a0099ed741b6fa47b64e4b7d1a11013ee2c1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 7 Nov 2024 17:33:03 +0000 Subject: [PATCH 13/42] Added first batch of corrections for P4estMeshView. --- src/meshes/p4est_mesh_view.jl | 18 ++++++++++++++++++ src/solvers/dgsem_p4est/containers.jl | 10 ++++++++-- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 31284942462..a276dc78160 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -36,4 +36,22 @@ end @inline balance!(::P4estMeshView) = nothing @inline ncells(mesh::P4estMeshView) = length(mesh.cell_ids) +function extract_interfaces!(mesh::P4estMeshView, interfaces) + u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) + u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] + u_new[:, :, :, 2:2:end] .= interfaces.u[:, :, :, mesh.cell_ids.*2] + node_indices_new = Array{eltype(interfaces.node_indices)}(undef, (2, size(mesh.cell_ids)[1]*2)) + node_indices_new[:, 1:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2 .-1)] + node_indices_new[:, 2:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2)] + neighbor_ids_new = Array{eltype(interfaces.neighbor_ids)}(undef, (2, size(mesh.cell_ids)[1]*2)) + neighbor_ids_new[:, 1:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2 .-1)] + neighbor_ids_new[:, 2:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2)] + interfaces.u = u_new + interfaces.node_indices = node_indices_new + interfaces.neighbor_ids = neighbor_ids_new + interfaces._u = vec(u_new) + interfaces._node_indices = vec(node_indices_new) + interfaces._neighbor_ids = vec(neighbor_ids_new) +end + end # @muladd diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index b88d7fd2886..914131bbefa 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -525,13 +525,17 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) return interfaces end -function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView) +function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView{2}) # 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.parent) iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) + # Extract the entry we need for this view. + extract_interfaces!(mesh, interfaces) + @autoinfiltrate + return interfaces end @@ -709,6 +713,9 @@ function count_required_surfaces(mesh::P4estMesh) iterate_p4est(mesh.p4est, user_data; iter_face_c = iter_face_c) + # Extract the entry we need for this view. + extract_interfaces!(mesh, interfaces) + # Return counters return (interfaces = user_data[1], mortars = user_data[2], @@ -724,7 +731,6 @@ function count_required_surfaces(mesh::P4estMeshView) iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) - @autoinfiltrate # Return counters return (interfaces = user_data[1], mortars = user_data[2], From 90e98bc67fc7ea470dc4b18e64aa8d29c6332315 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 14 Nov 2024 17:52:29 +0000 Subject: [PATCH 14/42] Added place holders for p4est mesh view. --- src/solvers/dgsem_p4est/containers.jl | 11 +++-------- src/solvers/dgsem_p4est/dg.jl | 26 +++++++++++++++++++++++++- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 914131bbefa..ab43f7fd858 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -532,16 +532,14 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView{2}) iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) - # Extract the entry we need for this view. - extract_interfaces!(mesh, interfaces) - @autoinfiltrate - return interfaces end # Initialization of interfaces after the function barrier function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack interfaces, interface_id, mesh = user_data + @autoinfiltrate + println("sides: ", sides_pw[1].treeid[], " ", sides_pw[2].treeid[]) user_data.interface_id += 1 # Get Tuple of local trees, one-based indexing @@ -713,9 +711,6 @@ function count_required_surfaces(mesh::P4estMesh) iterate_p4est(mesh.p4est, user_data; iter_face_c = iter_face_c) - # Extract the entry we need for this view. - extract_interfaces!(mesh, interfaces) - # Return counters return (interfaces = user_data[1], mortars = user_data[2], @@ -724,7 +719,7 @@ end function count_required_surfaces(mesh::P4estMeshView) # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face - iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh))) + iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh.parent))) # interfaces, mortars, boundaries user_data = [0, 0, 0] diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index b6259ff0207..7fd2c0ccfb9 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -8,7 +8,7 @@ # 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::Union{P4estMesh, P4estMeshView}, equations::AbstractEquations, dg::DG, ::Any, +function create_cache(mesh::P4estMesh, 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 @@ -29,6 +29,30 @@ function create_cache(mesh::Union{P4estMesh, P4estMeshView}, equations::Abstract return cache end +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) + + elements_parent = init_elements(mesh.parent, equations, dg.basis, uEltype) + interfaces_parent = init_interfaces(mesh.parent, equations, dg.basis, elements) + boundaries_parent = init_boundaries(mesh.parent, equations, dg.basis, elements) + mortars_parent = init_mortars(mesh.parent, equations, dg.basis, elements) + + # Extract data for views. + elemnts, interfaces, boundaries, mortars = extract...(elements_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 From 23e89d3f4d926fc3b1df5e6f5b532bb8feb6fc6b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 15 Nov 2024 10:26:48 +0000 Subject: [PATCH 15/42] Added p4est -> p4est mesh view extraction routine. --- src/meshes/p4est_mesh_view.jl | 23 ++++++++++++++++++++++- src/solvers/dgsem_p4est/containers.jl | 4 +--- src/solvers/dgsem_p4est/dg.jl | 17 ++++++++++++----- 3 files changed, 35 insertions(+), 9 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index a276dc78160..7a9a87d1544 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -36,7 +36,26 @@ end @inline balance!(::P4estMeshView) = nothing @inline ncells(mesh::P4estMeshView) = length(mesh.cell_ids) -function extract_interfaces!(mesh::P4estMeshView, interfaces) +# +function extract_p4est_mesh_view(elements_parent, + interfaces_parent, + boundaries_parent, + mortars_parent, + mesh) + @autoinfiltrate + 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] + interfaces = extract_interfaces(mesh, interfaces_parent) + + return elements, interfaces, boundaries_parent, mortars_parent +end + +function extract_interfaces(mesh::P4estMeshView, interfaces_parent) + interfaces = interfaces_parent u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] u_new[:, :, :, 2:2:end] .= interfaces.u[:, :, :, mesh.cell_ids.*2] @@ -52,6 +71,8 @@ function extract_interfaces!(mesh::P4estMeshView, interfaces) interfaces._u = vec(u_new) interfaces._node_indices = vec(node_indices_new) interfaces._neighbor_ids = vec(neighbor_ids_new) + + return interfaces end end # @muladd diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index ab43f7fd858..345627202a1 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -527,7 +527,7 @@ end function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView{2}) # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face - iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) + iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh.mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh.parent) iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) @@ -538,8 +538,6 @@ end # Initialization of interfaces after the function barrier function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack interfaces, interface_id, mesh = user_data - @autoinfiltrate - println("sides: ", sides_pw[1].treeid[], " ", sides_pw[2].treeid[]) user_data.interface_id += 1 # Get Tuple of local trees, one-based indexing diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 7fd2c0ccfb9..60f50ac08c4 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -29,19 +29,26 @@ 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) + balance!(mesh.parent) elements_parent = init_elements(mesh.parent, equations, dg.basis, uEltype) - interfaces_parent = init_interfaces(mesh.parent, equations, dg.basis, elements) - boundaries_parent = init_boundaries(mesh.parent, equations, dg.basis, elements) - mortars_parent = init_mortars(mesh.parent, equations, dg.basis, elements) + 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. - elemnts, interfaces, boundaries, mortars = extract...(elements_parent, ..., mesh) + elements, interfaces, boundaries, mortars = extract_p4est_mesh_view(elements_parent, + interfaces_parent, + boundaries_parent, + mortars_parent, + mesh) cache = (; elements, interfaces, boundaries, mortars) From b0edd8d3196e90964f7f2754068430abaddb6d88 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 15 Nov 2024 11:35:10 +0000 Subject: [PATCH 16/42] Added flattened arrays. --- src/meshes/p4est_mesh_view.jl | 5 +++++ src/solvers/dgsem_p4est/dg_2d.jl | 1 + 2 files changed, 6 insertions(+) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 7a9a87d1544..8b5a68a0817 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -49,6 +49,11 @@ function extract_p4est_mesh_view(elements_parent, 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 diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 24a531c3724..591a3f17fbf 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -102,6 +102,7 @@ function prolong2interfaces!(cache, u, j_secondary = j_secondary_start for i in eachnode(dg) for v in eachvariable(equations) + @autoinfiltrate interfaces.u[2, v, i, interface] = u[v, i_secondary, j_secondary, secondary_element] end From db837ded844810f706d5d546b47284c8c46db9b2 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 20 Nov 2024 18:02:16 +0000 Subject: [PATCH 17/42] Added prolong2interfaces to p4est mesh view. --- src/meshes/p4est_mesh_view.jl | 59 +++++++++++++++++++++++++++++++- src/solvers/dgsem_p4est/dg_2d.jl | 3 +- 2 files changed, 59 insertions(+), 3 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 8b5a68a0817..49b4d517b3f 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -42,7 +42,6 @@ function extract_p4est_mesh_view(elements_parent, boundaries_parent, mortars_parent, mesh) - @autoinfiltrate elements = elements_parent elements.inverse_jacobian = elements_parent.inverse_jacobian[.., mesh.cell_ids] elements.jacobian_matrix = elements_parent.jacobian_matrix[.., mesh.cell_ids] @@ -80,4 +79,62 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) return interfaces end +# We pass the `surface_integral` argument solely for dispatch +function prolong2interfaces!(cache, u, + mesh::P4estMeshView{2}, + equations, surface_integral, dg) + @unpack interfaces = cache + index_range = eachnode(dg) + + @threaded for interface in eachinterface(dg, cache) + # Copy solution data from the primary element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + @autoinfiltrate + primary_element = interfaces.neighbor_ids[1, interface] + primary_indices = interfaces.node_indices[1, interface] + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + + i_primary = i_primary_start + j_primary = j_primary_start + for i in eachnode(dg) + for v in eachvariable(equations) + interfaces.u[1, v, i, interface] = u[v, i_primary, j_primary, + primary_element] + end + i_primary += i_primary_step + j_primary += j_primary_step + end + + # Copy solution data from the secondary element using "delayed indexing" with + # a start value and a step size to get the correct face and orientation. + secondary_element = interfaces.neighbor_ids[2, interface] + secondary_indices = interfaces.node_indices[2, interface] + + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + for i in eachnode(dg) + for v in eachvariable(equations) + interfaces.u[2, v, i, interface] = u[v, i_secondary, j_secondary, + secondary_element] + end + i_secondary += i_secondary_step + j_secondary += j_secondary_step + end + end + + return nothing +end + end # @muladd diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 591a3f17fbf..611260b45f8 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -58,7 +58,7 @@ end # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack interfaces = cache index_range = eachnode(dg) @@ -102,7 +102,6 @@ function prolong2interfaces!(cache, u, j_secondary = j_secondary_start for i in eachnode(dg) for v in eachvariable(equations) - @autoinfiltrate interfaces.u[2, v, i, interface] = u[v, i_secondary, j_secondary, secondary_element] end From 9f5ae86f034eff0eaacde09f638eeae5ca45dae7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 21 Nov 2024 19:25:40 +0000 Subject: [PATCH 18/42] Minor corrections. --- src/meshes/p4est_mesh_view.jl | 2 +- src/solvers/dgsem_p4est/dg_2d.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 49b4d517b3f..cba9db45708 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -92,7 +92,6 @@ function prolong2interfaces!(cache, u, # Note that in the current implementation, the interface will be # "aligned at the primary element", i.e., the index of the primary side # will always run forwards. - @autoinfiltrate primary_element = interfaces.neighbor_ids[1, interface] primary_indices = interfaces.node_indices[1, interface] @@ -100,6 +99,7 @@ function prolong2interfaces!(cache, u, index_range) j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range) + @autoinfiltrate i_primary = i_primary_start j_primary = j_primary_start diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 611260b45f8..2c46e963a80 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -114,7 +114,7 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces From 03651bddddcafca4f1b87697e3092ad530568126 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 25 Nov 2024 19:08:50 +0000 Subject: [PATCH 19/42] Removed redundant P4estMeshView in calc_interface_flux. --- src/meshes/p4est_mesh_view.jl | 1 + src/solvers/dgsem_p4est/dg_2d.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index cba9db45708..d6cc17f3e2f 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -59,6 +59,7 @@ function extract_p4est_mesh_view(elements_parent, end function extract_interfaces(mesh::P4estMeshView, interfaces_parent) + @autoinfiltrate interfaces = interfaces_parent u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 2c46e963a80..611260b45f8 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -114,7 +114,7 @@ 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 From b36e64364cb1c39ba0bb59a852dfe74a28650c7f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 27 Nov 2024 19:31:27 +0000 Subject: [PATCH 20/42] Corrected extraction of interfaces given an array of cell ids to include in the p4est mesh view. --- src/meshes/p4est_mesh_view.jl | 64 ++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 31 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index d6cc17f3e2f..fc114108507 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -11,22 +11,11 @@ A view on a p4est mesh. mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: AbstractMesh{NDIMS} parent::Parent cell_ids::Vector{Int} - # SC: After some thought, we might need to create a p4est pointer to p4est data - # conatining the data from the view. -# unsaved_changes::Bool + unsaved_changes::Bool end function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}, cell_ids::Vector) where {NDIMS, NDIMS_AMBIENT, RealT} -# # SC: number of cells should be corrected. -# cell_ids = Vector{Int}(undef, ncells(parent)) -# # SC: do not populate this array. It needs to be given by the user. -# for i in 1:ncells(parent) -# cell_ids[i] = i -# end - - # SC: Since we need a p4est pointer no the modified (view) p4est data, we might need a function - # like connectivity_structured that computes the connectivity. - return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids)#, parent.unsaved_changes) + return P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT, typeof(parent)}(parent, cell_ids, parent.unsaved_changes) end @inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS @@ -59,23 +48,38 @@ function extract_p4est_mesh_view(elements_parent, end function extract_interfaces(mesh::P4estMeshView, interfaces_parent) - @autoinfiltrate + """ + 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 - u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) - u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] - u_new[:, :, :, 2:2:end] .= interfaces.u[:, :, :, mesh.cell_ids.*2] - node_indices_new = Array{eltype(interfaces.node_indices)}(undef, (2, size(mesh.cell_ids)[1]*2)) - node_indices_new[:, 1:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2 .-1)] - node_indices_new[:, 2:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2)] - neighbor_ids_new = Array{eltype(interfaces.neighbor_ids)}(undef, (2, size(mesh.cell_ids)[1]*2)) - neighbor_ids_new[:, 1:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2 .-1)] - neighbor_ids_new[:, 2:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2)] - interfaces.u = u_new - interfaces.node_indices = node_indices_new - interfaces.neighbor_ids = neighbor_ids_new - interfaces._u = vec(u_new) - interfaces._node_indices = vec(node_indices_new) - interfaces._neighbor_ids = vec(neighbor_ids_new) + interfaces.u = interfaces_parent.u[.., mask] + interfaces.node_indices = interfaces_parent.node_indices[.., mask] + interfaces.neighbor_ids = interfaces_parent.neighbor_ids[.., mask] + +# u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) +# u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] +# u_new[:, :, :, 2:2:end] .= interfaces.u[:, :, :, mesh.cell_ids.*2] +# node_indices_new = Array{eltype(interfaces.node_indices)}(undef, (2, size(mesh.cell_ids)[1]*2)) +# node_indices_new[:, 1:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2 .-1)] +# node_indices_new[:, 2:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2)] +# neighbor_ids_new = Array{eltype(interfaces.neighbor_ids)}(undef, (2, size(mesh.cell_ids)[1]*2)) +# neighbor_ids_new[:, 1:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2 .-1)] +# neighbor_ids_new[:, 2:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2)] +# interfaces.u = u_new +# interfaces.node_indices = node_indices_new +# interfaces.neighbor_ids = neighbor_ids_new + + # 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 @@ -100,8 +104,6 @@ function prolong2interfaces!(cache, u, index_range) j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range) - @autoinfiltrate - i_primary = i_primary_start j_primary = j_primary_start for i in eachnode(dg) From 8821ff2d8eaeb2b016e28b060be9a6e1044ccd04 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 28 Nov 2024 17:01:26 +0000 Subject: [PATCH 21/42] Write out of P4estMeshView data. Not finished yet. --- src/callbacks_step/save_solution_dg.jl | 1 + src/meshes/mesh_io.jl | 4 +-- src/meshes/p4est_mesh_view.jl | 47 +++++++++++++++++++++++++- 3 files changed, 49 insertions(+), 3 deletions(-) diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 57cdb3ef8a2..4c804a0f084 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/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 127bc420f65..1e0eea1b7c3 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -6,7 +6,7 @@ #! 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 @@ -303,7 +303,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/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index fc114108507..62fd807c86d 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -12,10 +12,13 @@ mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: Abs parent::Parent cell_ids::Vector{Int} unsaved_changes::Bool + current_filename::String end 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) + 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 @@ -140,4 +143,46 @@ function prolong2interfaces!(cache, u, return nothing 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 P4estMesh 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 From 240918f51b84ef1749e44a3ae3395f9592dcc7c1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Dec 2024 15:40:57 +0000 Subject: [PATCH 22/42] Fixed issue with non-simply connected p4est mesh views. --- src/meshes/p4est_mesh_view.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 62fd807c86d..b67b024be33 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -64,7 +64,13 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) interfaces = interfaces_parent interfaces.u = interfaces_parent.u[.., mask] interfaces.node_indices = interfaces_parent.node_indices[.., mask] - interfaces.neighbor_ids = interfaces_parent.neighbor_ids[.., 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(x->x==neighbor_ids[1, interface], mesh.cell_ids)[1] + interfaces.neighbor_ids[2, interface] = findall(x->x==neighbor_ids[2, interface], mesh.cell_ids)[1] + end # u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) # u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] From 9cae5b535c6a05470e12d488eb71b8d1418d35da Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 16 Dec 2024 16:07:43 +0000 Subject: [PATCH 23/42] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 03384151102..0f1eca3d8c3 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -30,8 +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::Union{P4estMesh{2, NDIMS_AMBIENT}, - P4estMeshView{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} From d21b6362fcfdba2994696621cf9f9a3849ba95a7 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 16 Dec 2024 16:08:32 +0000 Subject: [PATCH 24/42] Update src/meshes/mesh_io.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/meshes/mesh_io.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 1e0eea1b7c3..271e315fe28 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, P4estMeshView, 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 From 8a14d1f625362f7c76f0ad311f64f55a9cc6d915 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 16 Dec 2024 16:08:43 +0000 Subject: [PATCH 25/42] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 0f1eca3d8c3..f91c1ba7586 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -241,7 +241,8 @@ function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} - integrate_via_indices(u, mesh, equations, dg, cache; + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) return func(u_local, equations) From cc87c3adcf8801080cac6a3303134cdccfd75b24 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 16 Dec 2024 16:09:01 +0000 Subject: [PATCH 26/42] Update src/meshes/p4est_mesh_view.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/meshes/p4est_mesh_view.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index b67b024be33..b9ce5939f3c 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -8,7 +8,8 @@ A view on a p4est mesh. """ -mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: AbstractMesh{NDIMS} +mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: + AbstractMesh{NDIMS} parent::Parent cell_ids::Vector{Int} unsaved_changes::Bool From 3af246ee923b2a14d340b8a4750558cf748ed8cf Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Mon, 16 Dec 2024 16:09:15 +0000 Subject: [PATCH 27/42] Update src/meshes/p4est_mesh_view.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/meshes/p4est_mesh_view.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index b9ce5939f3c..a51b49fcf0e 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -16,7 +16,8 @@ mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: current_filename::String end -function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}, cell_ids::Vector) where {NDIMS, NDIMS_AMBIENT, RealT} +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) From fe26e98b6eef773291921480dc7f9992d83793e3 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Dec 2024 16:28:54 +0000 Subject: [PATCH 28/42] No nead for P4estMeshView having its own prolong2interfaces function. --- src/meshes/p4est_mesh_view.jl | 93 ++++++-------------------------- src/solvers/dgsem_p4est/dg_2d.jl | 2 +- 2 files changed, 18 insertions(+), 77 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index b67b024be33..ac92b43a3f1 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -1,13 +1,15 @@ +# 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, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, - NNODES} <: - AbstractMesh{NDIMS} + 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} @@ -15,6 +17,15 @@ mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: Abs 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, @@ -28,7 +39,6 @@ end @inline balance!(::P4estMeshView) = nothing @inline ncells(mesh::P4estMeshView) = length(mesh.cell_ids) -# function extract_p4est_mesh_view(elements_parent, interfaces_parent, boundaries_parent, @@ -68,23 +78,10 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) # 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(x->x==neighbor_ids[1, interface], mesh.cell_ids)[1] - interfaces.neighbor_ids[2, interface] = findall(x->x==neighbor_ids[2, interface], mesh.cell_ids)[1] + 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 -# u_new = Array{eltype(interfaces.u)}(undef, (size(interfaces.u)[1:3]..., size(mesh.cell_ids)[1]*2)) -# u_new[:, :, :, 1:2:end] .= interfaces.u[:, :, :, (mesh.cell_ids.*2 .-1)] -# u_new[:, :, :, 2:2:end] .= interfaces.u[:, :, :, mesh.cell_ids.*2] -# node_indices_new = Array{eltype(interfaces.node_indices)}(undef, (2, size(mesh.cell_ids)[1]*2)) -# node_indices_new[:, 1:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2 .-1)] -# node_indices_new[:, 2:2:end] .= interfaces.node_indices[:, (mesh.cell_ids.*2)] -# neighbor_ids_new = Array{eltype(interfaces.neighbor_ids)}(undef, (2, size(mesh.cell_ids)[1]*2)) -# neighbor_ids_new[:, 1:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2 .-1)] -# neighbor_ids_new[:, 2:2:end] .= interfaces.neighbor_ids[:, (mesh.cell_ids.*2)] -# interfaces.u = u_new -# interfaces.node_indices = node_indices_new -# interfaces.neighbor_ids = neighbor_ids_new - # Flatten the arrays. interfaces._u = vec(interfaces.u) interfaces._node_indices = vec(interfaces.node_indices) @@ -93,62 +90,6 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) return interfaces end -# We pass the `surface_integral` argument solely for dispatch -function prolong2interfaces!(cache, u, - mesh::P4estMeshView{2}, - equations, surface_integral, dg) - @unpack interfaces = cache - index_range = eachnode(dg) - - @threaded for interface in eachinterface(dg, cache) - # Copy solution data from the primary element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - # Note that in the current implementation, the interface will be - # "aligned at the primary element", i.e., the index of the primary side - # will always run forwards. - primary_element = interfaces.neighbor_ids[1, interface] - primary_indices = interfaces.node_indices[1, interface] - - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) - i_primary = i_primary_start - j_primary = j_primary_start - for i in eachnode(dg) - for v in eachvariable(equations) - interfaces.u[1, v, i, interface] = u[v, i_primary, j_primary, - primary_element] - end - i_primary += i_primary_step - j_primary += j_primary_step - end - - # Copy solution data from the secondary element using "delayed indexing" with - # a start value and a step size to get the correct face and orientation. - secondary_element = interfaces.neighbor_ids[2, interface] - secondary_indices = interfaces.node_indices[2, interface] - - i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], - index_range) - j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], - index_range) - - i_secondary = i_secondary_start - j_secondary = j_secondary_start - for i in eachnode(dg) - for v in eachvariable(equations) - interfaces.u[2, v, i, interface] = u[v, i_secondary, j_secondary, - secondary_element] - end - i_secondary += i_secondary_step - j_secondary += j_secondary_step - end - end - - return nothing -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 P4estMesh and its node coordinates are reconstructured from diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 611260b45f8..24a531c3724 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -58,7 +58,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) From 43ab4d5fd03a8086f42123afefd90c2a82ab5da0 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Dec 2024 16:33:43 +0000 Subject: [PATCH 29/42] Corrected typo in comment. --- src/meshes/p4est_mesh_view.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index ac92b43a3f1..ed25948a75e 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -92,7 +92,7 @@ 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 P4estMesh and its node coordinates are reconstructured from +# 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) From 86f984b545b45e5e7b2a3d682b3e380d230dae14 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 17 Dec 2024 12:17:31 +0000 Subject: [PATCH 30/42] Fixed typos. --- src/callbacks_step/analysis_dg2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index f91c1ba7586..0cf315be9b8 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -241,8 +241,7 @@ function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache; normalize = true) where {Func} - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, - T8codeMesh{2}}, + integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) return func(u_local, equations) From 0af718530f5fccdbb6cb242c9cc9e0754d13d47e Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 17 Dec 2024 12:33:08 +0000 Subject: [PATCH 31/42] Update src/solvers/dgsem_unstructured/dg_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/solvers/dgsem_unstructured/dg_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 2aa6ec24a42..605d1338d59 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -317,7 +317,8 @@ end # Function barrier for type stability function calc_boundary_flux!(cache, t, boundary_conditions, - mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, T8codeMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, + T8codeMesh}, equations, surface_integral, dg::DG) @unpack boundary_condition_types, boundary_indices = boundary_conditions From 8d011a68e3fad6b713c6c04a9df21c304e8e74ec Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 17 Dec 2024 12:33:28 +0000 Subject: [PATCH 32/42] Update src/solvers/dgsem_tree/dg_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/solvers/dgsem_tree/dg_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 1786c359c05..91c8ac4299d 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -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}, P4estMeshView{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 From 8bf43a9465365eda166366e1b97299d39ca39161 Mon Sep 17 00:00:00 2001 From: Simon Candelaresi <10759273+SimonCan@users.noreply.github.com> Date: Tue, 17 Dec 2024 12:40:44 +0000 Subject: [PATCH 33/42] Update src/callbacks_step/analysis_dg2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/callbacks_step/analysis_dg2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index d05d9dcee66..7d4de80b952 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -241,7 +241,8 @@ end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{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 From c51b0aeb5c7855bdae573f9ca020d940307632ce Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 17 Dec 2024 13:29:51 +0000 Subject: [PATCH 34/42] Applied autoformatter. --- src/callbacks_step/analysis_dg2d.jl | 5 ++-- src/meshes/p4est_mesh_view.jl | 26 ++++++++++++------- .../semidiscretization_hyperbolic.jl | 3 ++- src/solvers/dgsem_p4est/containers.jl | 16 +++++++----- src/solvers/dgsem_p4est/containers_2d.jl | 9 ++++--- src/solvers/dgsem_p4est/dg.jl | 14 +++++----- src/solvers/dgsem_p4est/dg_2d.jl | 15 +++++++---- src/solvers/dgsem_tree/dg_2d.jl | 6 +++-- src/solvers/dgsem_unstructured/dg_2d.jl | 6 +++-- 9 files changed, 64 insertions(+), 36 deletions(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index d05d9dcee66..ab57d785f52 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -241,9 +241,10 @@ end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{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; + integrate_via_indices(u, mesh, equations, dg, cache; normalize = normalize) do u, i, j, element, equations, dg u_local = get_node_vars(u, equations, dg, i, j, element) return func(u_local, equations) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index ed25948a75e..4c1546de06e 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -10,7 +10,8 @@ A view on a p4est mesh. """ -mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: AbstractMesh{NDIMS} +mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, Parent} <: + AbstractMesh{NDIMS} parent::Parent cell_ids::Vector{Int} unsaved_changes::Bool @@ -26,7 +27,8 @@ Create a P4estMeshView on a P4estMesh parent. - `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} +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) @@ -48,8 +50,10 @@ function extract_p4est_mesh_view(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.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) @@ -69,7 +73,7 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) 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) + (interfaces_parent.neighbor_ids[2, interface] in mesh.cell_ids) end interfaces = interfaces_parent interfaces.u = interfaces_parent.u[.., mask] @@ -78,8 +82,12 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) # 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] + 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. @@ -120,7 +128,8 @@ function save_mesh_file(mesh::P4estMeshView, output_directory, timestep, attributes(file)["ndims"] = ndims(mesh) attributes(file)["p4est_file"] = p4est_filename - file["tree_node_coordinates"] = mesh.parent.tree_node_coordinates[.., mesh.cell_ids] + 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 @@ -131,5 +140,4 @@ function save_mesh_file(mesh::P4estMeshView, output_directory, timestep, end @inline mpi_parallel(mesh::P4estMeshView) = False() - end # @muladd diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 8185d0507c9..c909196b5db 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic.jl @@ -213,7 +213,8 @@ function digest_boundary_conditions(boundary_conditions::AbstractArray, mesh, so end # No checks for these meshes yet available -function check_periodicity_mesh_boundary_conditions(mesh::Union{P4estMesh, P4estMeshView, +function check_periodicity_mesh_boundary_conditions(mesh::Union{P4estMesh, + P4estMeshView, UnstructuredMesh2D, T8codeMesh, DGMultiMesh}, diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 345627202a1..4f879678e97 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -82,8 +82,8 @@ end # Create element container and initialize element data function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT}, - P4estMeshView{NDIMS, NDIMS, RealT}, - T8codeMesh{NDIMS, RealT}}, + P4estMeshView{NDIMS, NDIMS, RealT}, + T8codeMesh{NDIMS, RealT}}, equations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} @@ -168,7 +168,8 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity) end # Create interface container and initialize interface data. -function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, basis, elements) +function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, + basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -243,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, P4estMeshView, T8codeMesh}, equations, basis, elements) +function init_boundaries(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, + basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -374,7 +376,8 @@ function Base.resize!(mortars::P4estMortarContainer, capacity) end # Create mortar container and initialize mortar data. -function init_mortars(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, basis, elements) +function init_mortars(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, + basis, elements) NDIMS = ndims(elements) uEltype = eltype(elements) @@ -528,7 +531,8 @@ end function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView{2}) # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh.mesh))) - user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh.parent) + user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, + mesh.parent) iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index eb8d8cf85b5..8fb812fd5d7 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}, P4estMeshView{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}, P4estMeshView{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. @@ -138,7 +140,8 @@ function calc_node_coordinates!(node_coordinates, multiply_dimensionwise!(view(node_coordinates, :, :, :, element), matrix1, matrix2, - view(mesh.parent.tree_node_coordinates, :, :, :, tree), + view(mesh.parent.tree_node_coordinates, :, :, :, + tree), tmp1) end end diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 60f50ac08c4..d4a1b2d3a50 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -39,16 +39,18 @@ function create_cache(mesh::P4estMeshView, equations::AbstractEquations, dg::DG, 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) + 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) + interfaces_parent, + boundaries_parent, + mortars_parent, + mesh) cache = (; elements, interfaces, boundaries, mortars) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index e1ae7a96b9c..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}, P4estMeshView{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)}, @@ -118,7 +119,8 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMeshView{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}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -382,7 +385,8 @@ end end function prolong2mortars!(cache, u, - mesh::Union{P4estMesh{2}, P4estMeshView{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 @@ -635,7 +639,8 @@ end end function calc_surface_integral!(du, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 1786c359c05..0550f0a9718 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -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}, P4estMeshView{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}, P4estMeshView{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 2aa6ec24a42..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, P4estMeshView, 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, P4estMeshView, T8codeMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, + T8codeMesh}, equations, surface_integral, dg::DG) @unpack boundary_condition_types, boundary_indices = boundary_conditions From 426471803a7448d4542be8b61ed7ea1979052023 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 18 Dec 2024 18:12:40 +0000 Subject: [PATCH 35/42] Fixed t8code issue. --- _mhd/Project.toml | 2 -- src/solvers/dgsem_tree/dg_2d.jl | 2 +- test/Project.toml.bkp | 5 ----- 3 files changed, 1 insertion(+), 8 deletions(-) delete mode 100644 _mhd/Project.toml delete mode 100644 test/Project.toml.bkp diff --git a/_mhd/Project.toml b/_mhd/Project.toml deleted file mode 100644 index cf4a19ed76e..00000000000 --- a/_mhd/Project.toml +++ /dev/null @@ -1,2 +0,0 @@ -[deps] -MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 0550f0a9718..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, - 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, diff --git a/test/Project.toml.bkp b/test/Project.toml.bkp deleted file mode 100644 index c19054c58d7..00000000000 --- a/test/Project.toml.bkp +++ /dev/null @@ -1,5 +0,0 @@ -[deps] -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" -Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95" From b56d1d7440eb57db7224a9492ae7b0f0b5e91822 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 19 Dec 2024 13:59:19 +0000 Subject: [PATCH 36/42] Added p4est mesh view into tests. --- test/test_p4est_2d.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 0e70282e77c..ce0412d7aa6 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -125,6 +125,20 @@ end end end +@trixi_testset "elixir_advection_view.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_view.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"), From ecf74d1410b458e360bbaf28023ecdbd7c216897 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 19 Dec 2024 14:00:29 +0000 Subject: [PATCH 37/42] Reduced the number of elements. --- .../elixir_advection_meshview.jl | 65 +++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_advection_meshview.jl 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() From 07e89edc8459ab0106f594b2dde3fb14dd6a70f0 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 19 Dec 2024 14:00:49 +0000 Subject: [PATCH 38/42] Added P4estMeshView to the error norms. --- src/callbacks_step/analysis_dg2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index ab57d785f52..1785a5e725e 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -140,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 From fe8b20913848655fb2b76b93f743176373f48c41 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 19 Dec 2024 14:01:58 +0000 Subject: [PATCH 39/42] Fixed typo in example run. --- test/test_p4est_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index ce0412d7aa6..9541c36c7a9 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -125,8 +125,8 @@ end end end -@trixi_testset "elixir_advection_view.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_view.jl"), +@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 From 84d22c85ded8cf997142a18d9e656a7a59acf716 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 19 Dec 2024 14:39:49 +0000 Subject: [PATCH 40/42] Applied autoformatter. --- examples/p4est_2d_dgsem/Project.toml | 9 --------- src/callbacks_step/analysis_dg2d.jl | 3 ++- 2 files changed, 2 insertions(+), 10 deletions(-) delete mode 100644 examples/p4est_2d_dgsem/Project.toml diff --git a/examples/p4est_2d_dgsem/Project.toml b/examples/p4est_2d_dgsem/Project.toml deleted file mode 100644 index 6b8f4dc58c1..00000000000 --- a/examples/p4est_2d_dgsem/Project.toml +++ /dev/null @@ -1,9 +0,0 @@ -[deps] -AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4" -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" -Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" -Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95" diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 1785a5e725e..c538d444f5a 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -140,7 +140,8 @@ end function calc_error_norms(func, u, t, analyzer, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, equations, initial_condition, dg::DGSEM, cache, cache_analysis) From b83976741a96715a79e1b4e74b9c49f09c36ce83 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 20 Dec 2024 18:11:58 +0000 Subject: [PATCH 41/42] Removed unused calc_node_coordinates for p4est mesh view. --- src/solvers/dgsem_p4est/containers_2d.jl | 62 ------------------------ 1 file changed, 62 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 8fb812fd5d7..04a2b9f9ecd 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -37,17 +37,6 @@ function calc_node_coordinates!(node_coordinates, calc_node_coordinates!(node_coordinates, mesh, basis.nodes) end -# Interpolate tree_node_coordinates to each quadrant at the nodes of the specified basis -function calc_node_coordinates!(node_coordinates, - mesh::P4estMeshView{2}, - basis::LobattoLegendreBasis) - # Hanging nodes will cause holes in the mesh if its polydeg is higher - # than the polydeg of the solver. - @assert length(basis.nodes)>=length(mesh.parent.nodes) "The solver can't have a lower polydeg than the mesh" - - calc_node_coordinates!(node_coordinates, mesh, basis.nodes) -end - # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, mesh::P4estMesh{2, NDIMS_AMBIENT}, @@ -98,57 +87,6 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end -# Interpolate tree_node_coordinates to each quadrant at the specified nodes -function calc_node_coordinates!(node_coordinates, - mesh::P4estMeshView{2, NDIMS_AMBIENT}, - nodes::AbstractVector) where {NDIMS_AMBIENT} - # 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(NDIMS_AMBIENT), static_length(nodes), - static_length(mesh.parent.nodes)) - matrix1 = StrideArray(undef, real(mesh.parent), - static_length(nodes), static_length(mesh.parent.nodes)) - matrix2 = similar(matrix1) - baryweights_in = barycentric_weights(mesh.parent.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.parent.p4est.trees)[mesh.cell_ids] - - 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.parent.nodes, nodes_out_x, - baryweights_in) - polynomial_interpolation_matrix!(matrix2, mesh.parent.nodes, nodes_out_y, - baryweights_in) - - multiply_dimensionwise!(view(node_coordinates, :, :, :, element), - matrix1, matrix2, - view(mesh.parent.tree_node_coordinates, :, :, :, - tree), - tmp1) - end - end - - return node_coordinates -end - # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{2}, faces, orientation, interface_id) From 93eab52c96570327bf312291157f1154463d21de Mon Sep 17 00:00:00 2001 From: SimonCan Date: Sat, 21 Dec 2024 17:50:27 +0000 Subject: [PATCH 42/42] Removed unused code. --- src/solvers/dgsem_p4est/containers.jl | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 4f879678e97..a070db6b701 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -528,17 +528,6 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) return interfaces end -function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMeshView{2}) - # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face - iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh.mesh))) - user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, - mesh.parent) - - iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) - - return interfaces -end - # Initialization of interfaces after the function barrier function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data) @unpack interfaces, interface_id, mesh = user_data @@ -719,21 +708,6 @@ function count_required_surfaces(mesh::P4estMesh) boundaries = user_data[3]) end -function count_required_surfaces(mesh::P4estMeshView) - # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face - iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh.parent))) - - # interfaces, mortars, boundaries - user_data = [0, 0, 0] - - iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) - - # Return counters - return (interfaces = user_data[1], - mortars = user_data[2], - boundaries = user_data[3]) -end - # Return direction of the face, which is indexed by node_indices @inline function indices2direction(indices) if indices[1] === :begin