From bb074193b9f8b1f66be706ef2a5db8cad7fa6b6b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 26 Jun 2023 16:31:21 +0100 Subject: [PATCH 01/43] Added coupling converters. --- src/Trixi.jl | 2 ++ .../coupling_converters.jl | 13 ++++++++++ .../coupling_converters_2d.jl | 25 +++++++++++++++++++ 3 files changed, 40 insertions(+) create mode 100644 src/coupling_converters/coupling_converters.jl create mode 100644 src/coupling_converters/coupling_converters_2d.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 66878f4b45..85b536352b 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -189,6 +189,8 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled +export coupling_converter_heaviside_2d + export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic, diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl new file mode 100644 index 0000000000..a82b849f28 --- /dev/null +++ b/src/coupling_converters/coupling_converters.jl @@ -0,0 +1,13 @@ +# 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 + +#################################################################################################### +# Include files with actual implementations for different systems of equations. + +include("coupling_converters_2d.jl") + +end # @muladd diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl new file mode 100644 index 0000000000..f18058632c --- /dev/null +++ b/src/coupling_converters/coupling_converters_2d.jl @@ -0,0 +1,25 @@ +# 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 + +@doc raw""" + Coupling converter function for a system of two LinearScalarAdvectionEquation2D. + +The coupling is given as a Heaviside step. +```math +c(x) = {c_0, for x \ge x_0 \times s + 0, for x < x_0 \times s} +``` +Here, `s` is the sign of the step function, x_0 the save_position +of the step and c_0 the amplitude. +""" +function coupling_converter_heaviside_2d(x, x_0, c_0, s, + equations_left::LinearScalarAdvectionEquation2D, + equation_right::LinearScalarAdvectionEquation2D) + return c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 +end + +end # @muladd From 95892bf8348689fa19e303f0d37c71cb563c3a05 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 29 Jun 2023 14:23:06 +0100 Subject: [PATCH 02/43] Added generic converter_function for structured 2d meshes. --- src/Trixi.jl | 1 + .../coupling_converters_2d.jl | 20 +++++++++++++++---- .../semidiscretization_coupled.jl | 14 ++++++++----- 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 85b536352b..0759ba129d 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -123,6 +123,7 @@ include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") include("time_integration/time_integration.jl") +include("coupling_converters/coupling_converters.jl") # `trixi_include` and special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl index f18058632c..e360ffc4ea 100644 --- a/src/coupling_converters/coupling_converters_2d.jl +++ b/src/coupling_converters/coupling_converters_2d.jl @@ -5,6 +5,7 @@ @muladd begin #! format: noindent + @doc raw""" Coupling converter function for a system of two LinearScalarAdvectionEquation2D. @@ -16,10 +17,21 @@ c(x) = {c_0, for x \ge x_0 \times s Here, `s` is the sign of the step function, x_0 the save_position of the step and c_0 the amplitude. """ -function coupling_converter_heaviside_2d(x, x_0, c_0, s, - equations_left::LinearScalarAdvectionEquation2D, - equation_right::LinearScalarAdvectionEquation2D) - return c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 +function coupling_converter_heaviside_2d(x_0, c_0, s) + return (x, u) -> c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 * u +end + + +@doc raw""" + Coupling converter function for a system of two LinearScalarAdvectionEquation2D. + +The coupling is given as a linear function. +```math +c(x) = x * u(x) +``` +""" +function coupling_converter_linear_2d() + return (x, u) -> x[2]*u end end # @muladd diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index b7adff7842..babf4527a5 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -371,8 +371,9 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_semi_index :: Int other_orientation :: Int indices :: Indices + coupling_converter :: Function - function BoundaryConditionCoupled(other_semi_index, indices, uEltype) + function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) NDIMS = length(indices) u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1)) @@ -385,7 +386,7 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic end new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index, - other_orientation, indices) + other_orientation, indices, coupling_converter) end end @@ -495,9 +496,12 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for i in eachnode(solver) for v in 1:size(u, 1) - boundary_condition.u_boundary[v, i, cell] = u[v, i_node, j_node, - linear_indices[i_cell, - j_cell]] + x = cache.elements.node_coordinates[:, i_node, j_node, linear_indices[i_cell, j_cell]] + converted_u = boundary_condition.coupling_converter(x, u[:, i_node, j_node, linear_indices[i_cell, j_cell]]) + boundary_condition.u_boundary[v, i, cell] = converted_u[v] + # boundary_condition.u_boundary[v, i, cell] = u[v, i_node, j_node, + # linear_indices[i_cell, + # j_cell]] end i_node += i_node_step j_node += j_node_step From 6fdc7147966276b8030a2d01a9364c6852df25d7 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 28 Nov 2023 14:20:54 +0000 Subject: [PATCH 03/43] Added p4est-view. --- src/meshes/p4est_mesh_view.jl | 126 ++++++++++++++++++++++++++++++++++ 1 file changed, 126 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 0000000000..2146a604f9 --- /dev/null +++ b/src/meshes/p4est_mesh_view.jl @@ -0,0 +1,126 @@ +@muladd begin +#! format: noindent + +mutable struct P4estMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} + parent::P4estMesh{NDIMS, RealT} + mapping::Any # Not relevant for performance + index_min::NTuple{NDIMS, Int} + index_max::NTuple{NDIMS, Int} +end + +function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; + index_min = ntuple(_ -> 1, Val(NDIMS)), + index_max = size(parent)) where {NDIMS, RealT} + @assert index_min <= index_max + @assert all(index_min .> 0) + @assert index_max <= size(parent) + + return P4estMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, + index_max) +end + +# Check if mesh is periodic +function isperiodic(mesh::P4estMeshView) + @unpack parent = mesh + return isperiodic(parent) && size(parent) == size(mesh) +end + +function isperiodic(mesh::P4estMeshView, dimension) + @unpack parent, index_min, index_max = mesh + return (isperiodic(parent, dimension) && + index_min[dimension] == 1 && + index_max[dimension] == size(parent, dimension)) +end + +@inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS +@inline Base.real(::P4estMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT +function Base.size(mesh::P4estMeshView) + @unpack index_min, index_max = mesh + return index_max .- index_min .+ 1 +end +function Base.size(mesh::P4estMeshView, i) + @unpack index_min, index_max = mesh + return index_max[i] - index_min[i] + 1 +end +Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) +Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) + +# function calc_node_coordinates!(node_coordinates, element, +# cell_x, cell_y, mapping, +# mesh::P4estMeshView{2}, +# # basis::LobattoLegendreBasis) +# basis) +# @unpack nodes = basis +# @unpack parent, index_min, index_max = mesh + +# # Get cell length in reference mesh +# dx = 2 / size(parent, 1) +# dy = 2 / size(parent, 2) + +# # Calculate index offsets with respect to parent +# parent_offset_x = index_min[1] - 1 +# parent_offset_y = index_min[2] - 1 + +# # Calculate node coordinates of reference mesh +# cell_x_offset = -1 + (cell_x - 1 + parent_offset_x) * dx + dx / 2 +# cell_y_offset = -1 + (cell_y - 1 + parent_offset_y) * dy + dy / 2 +# p4est_tree_t +# for j in eachnode(basis), i in eachnode(basis) +# # node_coordinates are the mapped reference node_coordinates +# node_coordinates[:, i, j, element] .= mapping(cell_x_offset + dx / 2 * nodes[i], +# cell_y_offset + dy / 2 * nodes[j]) +# end +# end + +# Interpolate tree_node_coordinates to each quadrant at the specified nodes +function calc_node_coordinates!(node_coordinates, + mesh::P4estMeshView{2}, + nodes::AbstractVector) + @unpack parent, index_min, index_max = mesh + + # 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(parent), + StaticInt(2), static_length(nodes), static_length(parent.nodes)) + matrix1 = StrideArray(undef, real(parent), + static_length(nodes), static_length(parent.nodes)) + matrix2 = similar(matrix1) + baryweights_in = barycentric_weights(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, parent.p4est.trees) + @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, parent.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, parent.nodes, nodes_out_y, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, element), + matrix1, matrix2, + view(parent.tree_node_coordinates, :, :, :, tree), + tmp1) + end + end + + return node_coordinates +end +end # @muladd From 1a38c0c8ba57b8f71c97a237cf53c604c1aa3d1b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 28 Nov 2023 16:40:00 +0000 Subject: [PATCH 04/43] Added p4estMeshView to various routines. More to come, especially bug fixing. --- examples/p4est_2d_dgsem/Project.toml | 1 + src/Trixi.jl | 3 +- src/callbacks_step/analysis_dg2d.jl | 16 ++++----- src/meshes/meshes.jl | 1 + src/meshes/p4est_mesh_view.jl | 46 +++++++++++++++++++----- src/solvers/dg.jl | 2 +- src/solvers/dgsem_p4est/containers.jl | 5 +-- src/solvers/dgsem_p4est/containers_2d.jl | 6 ++-- src/solvers/dgsem_p4est/dg.jl | 2 +- src/solvers/dgsem_p4est/dg_2d.jl | 10 +++--- src/solvers/dgsem_tree/dg_2d.jl | 19 +++++----- 11 files changed, 72 insertions(+), 39 deletions(-) diff --git a/examples/p4est_2d_dgsem/Project.toml b/examples/p4est_2d_dgsem/Project.toml index 094d841150..6b8f4dc58c 100644 --- a/examples/p4est_2d_dgsem/Project.toml +++ b/examples/p4est_2d_dgsem/Project.toml @@ -1,5 +1,6 @@ [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" diff --git a/src/Trixi.jl b/src/Trixi.jl index ad81fb33fd..1b96e361a5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -220,7 +220,8 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, export lake_at_rest_error export ncomponents, eachcomponent -export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, T8codeMesh +export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, P4estMeshView, + T8codeMesh export DG, DGSEM, LobattoLegendreBasis, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index a9e0cf87b0..e3609b9211 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -31,7 +31,7 @@ end function create_cache_analysis(analyzer, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -108,7 +108,7 @@ end function calc_error_norms(func, u, t, analyzer, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, equations, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @unpack node_coordinates, inverse_jacobian = cache.elements @@ -176,7 +176,7 @@ end function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, equations, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis @@ -204,7 +204,7 @@ end function integrate(func::Func, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + 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 @@ -214,7 +214,7 @@ function integrate(func::Func, u, end function integrate(func::Func, u, - mesh::Union{TreeMesh{2}, P4estMesh{2}}, + mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}}, equations, equations_parabolic, dg::DGSEM, cache, cache_parabolic; normalize = true) where {Func} @@ -233,7 +233,7 @@ end function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, @@ -278,7 +278,7 @@ end function analyze(::Val{:l2_divb}, du, u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) @unpack contravariant_vectors = cache.elements integrate_via_indices(u, mesh, equations, dg, cache, cache, @@ -346,7 +346,7 @@ end function analyze(::Val{:linf_divb}, du, u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) @unpack derivative_matrix, weights = dg.basis @unpack contravariant_vectors = cache.elements diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index ed2158b169..d5150ee202 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -12,6 +12,7 @@ include("unstructured_mesh.jl") include("face_interpolant.jl") include("transfinite_mappings_3d.jl") include("p4est_mesh.jl") +include("p4est_mesh_view.jl") include("t8code_mesh.jl") include("mesh_io.jl") include("dgmulti_meshes.jl") diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 2146a604f9..daed597410 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -3,20 +3,19 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::P4estMesh{NDIMS, RealT} - mapping::Any # Not relevant for performance - index_min::NTuple{NDIMS, Int} - index_max::NTuple{NDIMS, Int} + nodes::SVector + index_min::Int + index_max::Int end function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; - index_min = ntuple(_ -> 1, Val(NDIMS)), - index_max = size(parent)) where {NDIMS, RealT} + index_min = 1::Int, + index_max = sizeof(unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees))::Int) where {NDIMS, RealT} @assert index_min <= index_max @assert all(index_min .> 0) - @assert index_max <= size(parent) + @assert index_max <= sizeof(unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees)) - return P4estMeshView{NDIMS, RealT}(parent, parent.mapping, index_min, - index_max) + return P4estMeshView{NDIMS, RealT}(parent, parent.nodes, index_min, index_max) end # Check if mesh is periodic @@ -45,6 +44,36 @@ end Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) +function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL) + p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) + # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes + # See https://github.com/cburstedde/p4est/issues/112 + p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) +end + +@inline ncells(mesh::P4estMeshView) = Int(mesh.parent.p4est.local_num_quadrants[]) + +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 + +function init_interfaces!(interfaces, mesh::P4estMeshView) + init_surfaces!(interfaces, nothing, nothing, mesh.parent) + + return interfaces +end + # function calc_node_coordinates!(node_coordinates, element, # cell_x, cell_y, mapping, # mesh::P4estMeshView{2}, @@ -93,7 +122,6 @@ function calc_node_coordinates!(node_coordinates, p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) trees = unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees) - @autoinfiltrate for tree in eachindex(trees) offset = trees[tree].quadrants_offset diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 9e5ebc7f9b..5fb8774db7 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -436,7 +436,7 @@ function get_node_variables!(node_variables, mesh, equations, dg::DG, cache) end const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, - T8codeMesh} + 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 5fe68e0671..084c587877 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -81,7 +81,8 @@ function Base.resize!(elements::P4estElementContainer, capacity) end # Create element container and initialize element data -function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, T8codeMesh{NDIMS, RealT}}, +function init_elements(mesh::Union{P4estMesh{NDIMS, RealT}, P4estMeshView{NDIMS, RealT}, + T8codeMesh{NDIMS, RealT}}, equations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} @@ -166,7 +167,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) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 236d7d24c0..fe196562cb 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 @@ -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. @@ -37,7 +37,7 @@ end # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, - mesh::P4estMesh{2}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, nodes::AbstractVector) # We use `StrideArray`s here since these buffers are used in performance-critical # places and the additional information passed to the compiler makes them faster diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index ec50627d3e..55ac41c8f3 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 diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 36624f2ce8..f868306866 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)}, @@ -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}, P4estMesh{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}, P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -206,7 +206,7 @@ end # Inlined version of the interface flux computation for equations with conservative and nonconservative terms @inline function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 7ecf4c0003..a2f34737a7 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -37,14 +37,14 @@ end # The methods below are specialized on the volume integral type # and called from the basic `create_cache` method at the top. function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DG, uEltype) NamedTuple() end function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, equations, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] element_ids_dgfv = Int[] @@ -70,7 +70,7 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMe end function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, equations, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype) A3dp1_x = Array{uEltype, 3} @@ -92,7 +92,7 @@ end # The methods below are specialized on the mortar type # and called from the basic `create_cache` method at the top. function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal performance using different types MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, @@ -110,7 +110,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, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -181,7 +181,7 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -235,7 +235,7 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) @@ -332,7 +332,7 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DGSEM, cache) @@ -391,7 +391,8 @@ end @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2} + UnstructuredMesh2D, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2} }, nonconservative_terms, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) From 65a00ba33539efaa3e779a275ca2db8b3cdaf228 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 30 Nov 2023 12:54:04 +0000 Subject: [PATCH 05/43] Added P4estMeshView in further functions. --- src/callbacks_step/stepsize_dg2d.jl | 4 ++-- src/solvers/dgsem_p4est/containers.jl | 4 ++-- src/solvers/dgsem_p4est/dg_2d.jl | 12 ++++++------ src/solvers/dgsem_structured/dg_2d.jl | 2 +- src/solvers/dgsem_tree/dg_2d.jl | 2 +- src/solvers/dgsem_tree/indicators_2d.jl | 2 +- src/solvers/dgsem_unstructured/dg_2d.jl | 8 ++++---- 7 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 673c3ba6aa..7655463006 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -77,7 +77,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -113,7 +113,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, constant_speed::True, equations, dg::DG, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 084c587877..c1903a4be7 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -242,7 +242,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) @@ -373,7 +373,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) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index f868306866..a9f735a64b 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}, P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -182,7 +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}, P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -247,7 +247,7 @@ end end function prolong2boundaries!(cache, u, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack boundaries = cache index_range = eachnode(dg) @@ -276,7 +276,7 @@ function prolong2boundaries!(cache, u, end function calc_boundary_flux!(cache, t, boundary_condition::BC, boundary_indexing, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) where {BC} @unpack boundaries = cache @unpack surface_flux_values = cache.elements @@ -312,7 +312,7 @@ end # inlined version of the boundary flux calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -343,7 +343,7 @@ end # inlined version of the boundary flux with nonconservative terms calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, i_index, j_index, diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 25a0eea096..18d9c62d82 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -59,7 +59,7 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 @inline function weak_form_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, + P4estMesh{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] diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index a2f34737a7..28f7b6de72 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -181,7 +181,7 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - P4estMeshView{2}, T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 8333bb515d..1b8823d662 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -208,7 +208,7 @@ end end # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha -function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, alpha, +function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, alpha, alpha_tmp, dg, cache) # Copy alpha values such that smoothing is indpedenent of the element access order diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index b12a96c4c3..41eca25c04 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -307,14 +307,14 @@ end # TODO: Taal dimension agnostic function calc_boundary_flux!(cache, t, boundary_condition::BoundaryConditionPeriodic, - mesh::Union{UnstructuredMesh2D, P4estMesh, T8codeMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, T8codeMesh}, equations, surface_integral, dg::DG) @assert isempty(eachboundary(dg, cache)) end # Function barrier for type stability function calc_boundary_flux!(cache, t, boundary_conditions, - mesh::Union{UnstructuredMesh2D, P4estMesh, T8codeMesh}, + mesh::Union{UnstructuredMesh2D, P4estMesh, P4estMeshView, T8codeMesh}, equations, surface_integral, dg::DG) @unpack boundary_condition_types, boundary_indices = boundary_conditions @@ -328,7 +328,7 @@ end function calc_boundary_flux_by_type!(cache, t, BCs::NTuple{N, Any}, BC_indices::NTuple{N, Vector{Int}}, mesh::Union{UnstructuredMesh2D, P4estMesh, - T8codeMesh}, + P4estMeshView, T8codeMesh}, equations, surface_integral, dg::DG) where {N} # Extract the boundary condition type and index vector boundary_condition = first(BCs) @@ -352,7 +352,7 @@ end # terminate the type-stable iteration over tuples function calc_boundary_flux_by_type!(cache, t, BCs::Tuple{}, BC_indices::Tuple{}, mesh::Union{UnstructuredMesh2D, P4estMesh, - T8codeMesh}, + P4estMeshView, T8codeMesh}, equations, surface_integral, dg::DG) nothing end From 40c4e0aff036596cfe049d3bc4b1323f0c3b5f76 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 27 Jun 2024 16:34:22 +0100 Subject: [PATCH 06/43] Removed redundant toml files. --- examples/structured_2d_dgsem/Project.toml | 5 ----- examples/structured_3d_dgsem/Project.toml | 5 ----- 2 files changed, 10 deletions(-) delete mode 100644 examples/structured_2d_dgsem/Project.toml delete mode 100644 examples/structured_3d_dgsem/Project.toml diff --git a/examples/structured_2d_dgsem/Project.toml b/examples/structured_2d_dgsem/Project.toml deleted file mode 100644 index c19054c58d..0000000000 --- a/examples/structured_2d_dgsem/Project.toml +++ /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" diff --git a/examples/structured_3d_dgsem/Project.toml b/examples/structured_3d_dgsem/Project.toml deleted file mode 100644 index c19054c58d..0000000000 --- a/examples/structured_3d_dgsem/Project.toml +++ /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 8f5216edf0149523851045118430122cfb36d60b Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 27 Jun 2024 16:41:41 +0100 Subject: [PATCH 07/43] Removed redundant files. --- examples/tree_1d_dgsem/Project.toml | 5 ----- examples/tree_2d_dgsem/Project.toml | 8 -------- examples/tree_3d_dgsem/Project.toml | 5 ----- examples/unstructured_2d_dgsem/Project.toml | 5 ----- src/equations/.ideal_glm_mhd_2d.jl.swp | Bin 24576 -> 0 bytes test/.runtests.jl.swp | Bin 16384 -> 0 bytes tmp/Project.toml | 7 ------- utils/Project.toml | 6 ------ 8 files changed, 36 deletions(-) delete mode 100644 examples/tree_1d_dgsem/Project.toml delete mode 100644 examples/tree_2d_dgsem/Project.toml delete mode 100644 examples/tree_3d_dgsem/Project.toml delete mode 100644 examples/unstructured_2d_dgsem/Project.toml delete mode 100644 src/equations/.ideal_glm_mhd_2d.jl.swp delete mode 100644 test/.runtests.jl.swp delete mode 100644 tmp/Project.toml delete mode 100644 utils/Project.toml diff --git a/examples/tree_1d_dgsem/Project.toml b/examples/tree_1d_dgsem/Project.toml deleted file mode 100644 index c19054c58d..0000000000 --- a/examples/tree_1d_dgsem/Project.toml +++ /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" diff --git a/examples/tree_2d_dgsem/Project.toml b/examples/tree_2d_dgsem/Project.toml deleted file mode 100644 index 094d841150..0000000000 --- a/examples/tree_2d_dgsem/Project.toml +++ /dev/null @@ -1,8 +0,0 @@ -[deps] -AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4" -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/examples/tree_3d_dgsem/Project.toml b/examples/tree_3d_dgsem/Project.toml deleted file mode 100644 index c19054c58d..0000000000 --- a/examples/tree_3d_dgsem/Project.toml +++ /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" diff --git a/examples/unstructured_2d_dgsem/Project.toml b/examples/unstructured_2d_dgsem/Project.toml deleted file mode 100644 index c19054c58d..0000000000 --- a/examples/unstructured_2d_dgsem/Project.toml +++ /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" diff --git a/src/equations/.ideal_glm_mhd_2d.jl.swp b/src/equations/.ideal_glm_mhd_2d.jl.swp deleted file mode 100644 index 74c1493caf0289d7fcaa5516f68827b57982b4c2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 24576 zcmeI3Ym6jUb;m1qAjTLRL8K@k!DUAZ-MiD>(>=R8vtyiiW_P{oWgo*njOAtBQ`22N zUD@iYc2#xH46DVF6Dto}*hrD0L{T0>j`AW7au6bk9YugbLS#q4Cj^plL?V%w2*HVo zjPpD9-m32DnVpTvDj!0%^zW*=&vVZ`_trW0oLZf^_nuQ~zS?m3{jlSF<`Dfmd$B%kRDCyVvbq0wah4lm1Ddv1gH&|;$Z7->vkNpk5`cpwIiCZ=A!GW97MpEV32ky^KWuy8}%6}I<@Cr>aB0NaJ^bZe!+o) z0|f^P4ip?HI8bn);6TBFf&>4LIFR;lcOGIqH`y_@?DYpGuCLnP{WiTlk^fzr|2~`k z>_q-6Hvha$zuk6__47KL{tkO>P2~TeO`o>chbQv?%eH@y&Hva${&#HteK!B4iTsys z{sT7uv}K4b=KQD4&&K};6Z!vP^CMgT^*5WFto}FA-?sV7w*K=I`TuV7-*5BZ zX&I8WZ;Ly(+Wg@}zD+JJ1qTWa6dWixP;j8&K*5260|f^P4ip?HIPiVS0T)C57RLTZ zVo2iuZ&^P4@$HV&0e6E}-|RTw1YZST0e=hr5PS@H-~gzDm*3<#e+fPe-Ubb+jgIq+;0$;#I00@0pT~%O2wVXR;Ja^dobP}?0}p{K;Dg`}FbDpX z0sISi3j98J04#%_1Qqaha4UG_hbRkP1YZJ=fu9FI2Yv?3f(H1;>mBFs!AHPZa2t4u z1By?BUjdJT9#{r@!8^d~z&~MgzW}}feir;RcnDkp(_kOC4SWS-{b6td+ztMXvxEbCV?WVs+y1vg*5i z?iw?V8C5#xwW{i1b++tu24PELO@t|_Pp}E3wz-SKrnKJFtv99h(#Et(r&S)DvaoQE zM5i1JdZ)VWqglvG<4D=rvu6+Ezv#yI`NIWzoSCoI8#DX&Un+Iew4W^0YHqyY535nU zR&!UA+NgxEby4->Xy~p6n~F(_qyDB!QXV~(3|3oS5U8XZwJxg`Y#GFUxTf4bHQZKL zb)r~x{k5)Q9vExV5ByXTd+NwBr@~~d!YbblM?9KV$6{}dYA32Wvr59Fva_fr+3VVw z5H>q{hHo&_94~C^c4LpBg)qYJHJFg3bnafSl}2$X?ncvUSO*O-JFWWDYKg}s9+%`X z@yla#WWTMx=WnBlto2RrL+N8lC)k3fGu_-!`l;l(BT{rOqb0kZ5pLSZx|Fil z=+bmg)s2g+XG2{+LcIy>9UHaoN@r7JrrJ25N1;=vR6FXqEZN?GHMWWfvfk53U1Sfa z6S|d@p8I(DlG?BdLXsp*tLL7ktphc~Glw?O@FdTV7MzdTr}-G2BXxn`A^+mEcR zu`ZQ#nN+W}48}@+55!HbgDujJP&8!)=qqBcuZX@}z3r<(&*Kg&V5WLNXDx|%R_cWO zDSk^LkF!?bbc-^N+0rpxN%^d{PSGtOIJyO$BCW76vIL$)h^#82<_vLR%Cok#bC6k=^d&YVn zx^c)`U;0Wn`(D&3$s9FgD~IaWK!6b>>C*^Vey4ktH1B^ykWjTNL0#VYs8OL5G%~|8 zN>gtaHX`WfyVQ-mG3XRjnS6ZeN%H`;Veqrd7)yRFfS>leIAg ztz?f8cL!QpDDO@E3Ss`%V?X6B8bs@;#Y10|S$2p-bwr2lv1)yEHnT?<{Z}(C>tG9o zYD&ejnTRqBc`d(#=G^t%)KnZRySU2s=eE+F*m8w8$f2uFZ(vo-G#<24wNuH)kt>E( z?^#fCb5*@(L94o3JI}>Z6#exG(%hEV06j!kMq72xqyGD%IwM~_*^v`pv|IRH5gk8Q zu=NHB`%Kw2xTLLwY|!S_yECO`3d+u1?WhHp)VK-e+VdRHxHd^N{c<#2_ zOwpPjI;x`>WFnsVrLLz(Zs6lrLG?^ zz$r9mRsu=`GamKR3as!pWMcg^*p&QpK3a>cP}cX|rxp_}CZzF_B1P;Y^hZ(LT`9)| zbF-z(k6ZDJ%3uA7zDmsJ9+ZnQ>bnjS7C&KMw>+CLywiP z&5DdUq7VHxN>oXRC|4b^|1V>^KLy18$K2MxFJkw97W^LgZSW*G3*G|$g0j!pI&uzB z{0a^f94I(YaG>Bo!GVGU1qTWa6dWixP;j8&zzyaACr@fhSsOk-LAA~!eeNg@Pw{_* z$;E@0aPR25**tA7PuriT&0i|3nkk6)^DaLO{7~#?>nCi*QIa&fp0?&Ubtk(WYfh=d zZqVY8fny)B>veUGX-suh`oPd@T^U||@zF>0-O|;vIrh*5WXF!;0HpqkQ2*}H(PDPQ zqX%{95}x0hYK%_RB)j5dh2-kt3tii6W`0&ne170k*7pBr@LfNFk6Y~jX**$0W9L5% zeDGHAZS4Ff!DHYsxE1^h%@rIdI8bn) z;K2Va2eOw}XAaj2K4;$LrD45E49s3aKIB>5^51o=oT7VgW~UGoy-;@OGBTpDt&XB^0+A!DL&U!*~1#c;%io^6f~Id3nME+a(} z4`k|d)w9a1>-M;2?OCCLQ+8G)fkH5uo4k}}?Wl|`k>^`+$_to|AYjxHG>F8FifGz2 zjae^dg2sxBFA6h!abU`n_RXZ?UNs}?%K-OkZ*36jJc+g%<(UW*T}En{#-@-)T6po{ zUsnCEed418IINXQ&wo677EPK*$&oT;`A=_j5VXPvt~@D zZASataC$0NJIZuiA@NBq&*m=%Yu;+?iaSjr?H(KuXb~XJ;!#YOJ}@z>!tw9!T$4nA=!X3N5?2Do#h ze);Za$dggTv4q!|i%7f#+3AanR)UMn#boIk#;7;OVim*c2~_^{C?;xoSI z@R@Uq$IqN10<}Irial*dzK%St6QDISGgmt_H(!~p%*@VIW)B_QKUbM;?$<`7SgS_| z0TOsL3vN{xSV3#*+<=%B*{Cb4o0VlZP7`%T`;wzrSaj4ztsX`kuwbQzF!H^a(GuKy zR2!jc*@R+A?&;jGeacI_QJZjxD0V6mrP_}D%j&Fg#7j(|1f~*NyWBmZz2|;1tX%{RmabqHZy)`f!QN^SYF)gQG8iOFYn*6_~21ILC^4A z+y4i!#s3P3{hxh*|4D5A$G{~p58eR2j?Mo#I1AngK84NyI9LEb2EK&d|2N?C;FrJ? z;BoLMSOPbLo4}uQ7T|*}kh6g*_^PGh8Lms948F)&z$btQPJ)}kUvmy{1+0T<@B(K5 z&wws?Gx!we0Zs4#cpLZ@XB$s~N5D^je(5pB_%c81 zFgER^OE>!ts>A0`9Rk-oeD>7Y4-SBXU=GX!7|wKy!_x%q%f0z=YAk1H9A$cR7p3A1 zHHV%%9tNA@6UwXAUOA4-D(`c$_RW=3;{1SkAEfbwN-G`nG$&7n44&rY$s8r~bcm;aQ)h=I5m(nDfQ4 zb5p52g}kSPy;}j+E7T`wp3^jK2#e{9U5sWJ1w#(eWNGkW!O)Ked{(FraHK3Bo`hpx zN5Jgu9son5BfZLx>EUD1OlIsn!7MXj**rSBtP~eN#faHEI~BC&VbrFP3)(a{+uPjk haj=Zc1ZUL!q;Ema`St1IJ}a~N+=whzGBM2i{BKoKLFE7d diff --git a/test/.runtests.jl.swp b/test/.runtests.jl.swp deleted file mode 100644 index 5e1034243b2009478d035b481dfcd9e870587aea..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeI3Ux*|}9mgx@pIr3f6|f3e1*`&A0jt0}QUM(wBcF%o#|rPCE%tZp+P|y#?G*d# zyUrIs?anG-6|f3e1*`&A0jq#jz$#!BunJfOtO8bn|3L+qN60Do^)0jD$NT@$`v31I z33(Iz349$~2PeTBA0Xs=AOh#WS#TeC^=?AG3w&?^{PP4MuYezcAA)P(e(?JH2^oS9 zfuFvQkc;2}@aMY-`3dNP{~RaekKhULDNqL`@CopjI|=zE_y$-9C&9npOUS>#Hjp3! zHE;sF^&Ue00Nw!K0SVXw0(=tO0sirBLf!_y0AB-N0bd4}!Nf)~LHU;t9E3QmE$z;W<8w0;qM8yMZc3HELj8;?9<_GS~28>AjD zIp#x&wtmEXM$vA1_+r9YCZ@o7EkaI(PdBfweRi$Ad1+&lwp!F-LEzBmpWl7Bun@CE zHyn~#xOsK)L#EP|!`GG@FqJiLv9cpHZBp6%mMNX2u1*u~%{0#lI;;3p!D-hHEZ8`_ zV3X?R>ZKFT;e6YuyK?w^FjkC-hZm}_dGrcb7pkyv^a>xkwF;ZJP12QxD!dJn&hKPd z(iAP6l5}r6X9bUt$l3Kb8)fj_SE+NR^e9ujx|IfkwcXGwEicnbMfVfVJk>gr4`=I} zm#<#>^u-OFar@l8QR&Jw!s#Jb<`4!${4?MQ>+x13BTmnhSXaeN_eUdUW{>ARa=Ejd z4|%4&=go*pF%EEIewuBGx_*`RwT{)fTCFE^KkZgs8P*;*PF`UN22~!^p1iV#Y^=p; z5Y$%AHb{1DcKMo5b)VBXkuFz?@+6UoqB2E2IGla!#x>oRY2f82vvW+*Z5{;WDMJ?} zWzc{PKDH;)VQfAube2t!jUS^8 z$N3I-OUu-g6mF5x3A10K*RpU8gvX(Hr&5`%szWs*F_6NmVXjDBZUR=3ekJ6gOopQf zdubw~9=))Z#Tb8%E5XEkx-?2r`rw0fW2mY^ZHvg0+sfoTj_p^1#c{xKer(COTWqJL z>!yhkQOVh_6F6nl&dHVI>Ra@t7Z%sR8#1cJgs??1BcVNegElUWqCjmg)bSyyV*A8sxQrBaSA7xG9~8ugXcv-L*3S)ZzH8{LwQ znpCDrqVrJ5k>HA8Xo4I`^f%MFk#=;ZIE8CZh$!8u($C0*B82doYbF9U-TT(i0M8Q| zEh$r>X&~hdMFVjoZ)awKJ|^_#v}@8m;^^eqXPsI4w8)JjDU32J7xT_P)r`%)z(rY3 zsZ2e$&m!a(@~w|58C5ZoN~Snzn&QJp+^*2>5N+VJ&jP&5ys5bWh2d8 zL9XSwIBjm9A@HWL%Oa*~h|{z`eGX?C*tZ!N6LG^@d<&e1$ixcKDV#xEIsbncs}~%8 zF<9lp_LAvB0cQw3R8bQ5_-d#$4^fdJ6thpBlSqnW zTc~WE(AhO>n1w>8OmauhkLD<69G&clL%X~uB0^jlp_*NjbBnqG3*(X*+b&#rx?N1# zW!fy$mF1ahre=Z#uF6+6w8wSn()u%H>P*^zfy)_9d{`1ex5&%845K|MJ<3AdGC)Ni zKa;m4Q2=Ya8P2iVOjU=^?mSOu&CRspMkRlq9nj#psx zz-+Xq*k)fi#4P4XTk*he_YyV`8mo(_*X}XZ_J$D)g=_x1?BQPi+=A`7kF4RHFO$$b zt7xm3y8>&d=inpOg)^^G7F_rqWY#=Jp8sh3APdv4Zf>(wNteCZmPSHP&Mzp<*rE9vn5Ar22Ev`mk-h#*%1wvtpmIleM=`E! rZqj?b5RV4L9Op-+4t9N1YTEI=4NUQvW&DbA+TemQnz2( Date: Fri, 28 Jun 2024 12:23:25 +0100 Subject: [PATCH 08/43] Removed elixir on p4est coupling. This should be part of a coupling branch. --- .../elixir_advection_coupled.jl | 142 ------------------ 1 file changed, 142 deletions(-) delete mode 100644 examples/p4est_2d_dgsem/elixir_advection_coupled.jl diff --git a/examples/p4est_2d_dgsem/elixir_advection_coupled.jl b/examples/p4est_2d_dgsem/elixir_advection_coupled.jl deleted file mode 100644 index 5a62d2af2e..0000000000 --- a/examples/p4est_2d_dgsem/elixir_advection_coupled.jl +++ /dev/null @@ -1,142 +0,0 @@ -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Coupled semidiscretization of two linear advection systems using converter functions such that -# the upper half of the domain is coupled periodically, while the lower half is not coupled -# and any incoming wave is completely absorbed. -# -# In this elixir, we have a square domain that is divided into a left half and a right half. On each -# half of the domain, a completely independent SemidiscretizationHyperbolic is created for the -# linear advection equations. The two systems are coupled in the x-direction and have periodic -# boundaries in the y-direction. For a high-level overview, see also the figure below: -# -# (-1, 1) ( 1, 1) -# ┌────────────────────┬────────────────────┐ -# │ ↑ periodic ↑ │ ↑ periodic ↑ │ -# │ │ │ -# │ │ │ -# │ ========= │ ========= │ -# │ system #1 │ system #2 │ -# │ ========= │ ========= │ -# │ │ │ -# │ │ │ -# │ │ │ -# │ │ │ -# │ coupled -->│<-- coupled │ -# │ │ │ -# │<-- coupled │ coupled -->│ -# │ │ │ -# │ │ │ -# │ ↓ periodic ↓ │ ↓ periodic ↓ │ -# └────────────────────┴────────────────────┘ -# (-1, -1) ( 1, -1) - -advection_velocity = (0.2, -0.7) -equations = LinearScalarAdvectionEquation2D(advection_velocity) - -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) - -# First mesh is the left half of a [-1,1]^2 square -coordinates_min1 = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) -coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y)) - -# Define identical resolution as a variable such that it is easier to change from `trixi_include` -trees_per_dimension = (8, 16) - -trees_per_dimension1 = trees_per_dimension - -# Create P4estMesh with 8 x 16 trees and 16 x 32 elements -mesh1 = P4estMesh(trees_per_dimension1, polydeg = 3, - coordinates_min=coordinates_min1, coordinates_max=coordinates_max1, - initial_refinement_level = 1) - - -# The user can define their own coupling functions. -coupling_function1 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * u - -boundary_conditions_x_neg1 = BoundaryConditionCoupled(2, (:end, :i_forward), Float64, - coupling_function1) -boundary_conditions_x_pos1 = BoundaryConditionCoupled(2, (:begin, :i_forward), Float64, - coupling_function1) - -# A semidiscretization collects data structures and functions for the spatial discretization -semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, - solver, - boundary_conditions = ( - # Connect left boundary with right boundary of right mesh - x_neg = boundary_conditions_x_neg1, - # Connect right boundary with left boundary of right mesh - x_pos = boundary_conditions_x_pos1, - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic)) - -# Second mesh is the right half of a [-1,1]^2 square -coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y)) -coordinates_max2 = (1.0, 1.0) # maximum coordinates (max(x), max(y)) - -trees_per_dimension2 = trees_per_dimension - -# Create P4estMesh with 8 x 16 trees and 16 x 32 elements -mesh2 = P4estMesh(trees_per_dimension2, polydeg = 3, - coordinates_min=coordinates_min2, coordinates_max=coordinates_max2, - initial_refinement_level = 1) - - -coupling_function2 = (x, u) -> (sign(x[2] - 0.0)*0.1 + 1.0)/1.1 * u - -boundary_conditions_x_neg2 = BoundaryConditionCoupled(1, (:end, :i_forward), Float64, - coupling_function2) -boundary_conditions_x_pos2 = BoundaryConditionCoupled(1, (:begin, :i_forward), Float64, - coupling_function2) - -semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, - solver, - boundary_conditions = ( - # Connect left boundary with right boundary of left mesh - x_neg = boundary_conditions_x_neg2, - # Connect right boundary with left boundary of left mesh - x_pos = boundary_conditions_x_pos2, - y_neg = boundary_condition_periodic, - y_pos = boundary_condition_periodic)) - -# Create a semidiscretization that bundles semi1 and semi2 -semi = SemidiscretizationCoupled(semi1, semi2) - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span from 0.0 to 2.0 -ode = semidiscretize(semi, (0.0, 20.0)); - -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers -summary_callback = SummaryCallback() - -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback1 = AnalysisCallback(semi1, interval = 100) -analysis_callback2 = AnalysisCallback(semi2, interval = 100) -analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) - -# The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval = 1, - solution_variables = cons2prim) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 1.6) - -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, - stepsize_callback) - -############################################################################### -# run the simulation - -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), - dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback - save_everystep = false, callback = callbacks); - -# Print the timer summary -summary_callback() From e13094698c1dcf6a602981400889ea04f16f4d0d Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 28 Jun 2024 12:26:07 +0100 Subject: [PATCH 09/43] Removed default coupling converter functions. --- src/Trixi.jl | 3 -- .../coupling_converters.jl | 13 ------- .../coupling_converters_2d.jl | 37 ------------------- 3 files changed, 53 deletions(-) delete mode 100644 src/coupling_converters/coupling_converters.jl delete mode 100644 src/coupling_converters/coupling_converters_2d.jl diff --git a/src/Trixi.jl b/src/Trixi.jl index 4b9e056205..741fcae49e 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -131,7 +131,6 @@ include("semidiscretization/semidiscretization_hyperbolic.jl") include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("semidiscretization/semidiscretization_coupled.jl") -include("time_integration/time_integration.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") @@ -209,8 +208,6 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled -export coupling_converter_heaviside_2d - export initial_condition_convergence_test, source_terms_convergence_test export source_terms_harmonic export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic, diff --git a/src/coupling_converters/coupling_converters.jl b/src/coupling_converters/coupling_converters.jl deleted file mode 100644 index a82b849f28..0000000000 --- a/src/coupling_converters/coupling_converters.jl +++ /dev/null @@ -1,13 +0,0 @@ -# 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 - -#################################################################################################### -# Include files with actual implementations for different systems of equations. - -include("coupling_converters_2d.jl") - -end # @muladd diff --git a/src/coupling_converters/coupling_converters_2d.jl b/src/coupling_converters/coupling_converters_2d.jl deleted file mode 100644 index e360ffc4ea..0000000000 --- a/src/coupling_converters/coupling_converters_2d.jl +++ /dev/null @@ -1,37 +0,0 @@ -# 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 - - -@doc raw""" - Coupling converter function for a system of two LinearScalarAdvectionEquation2D. - -The coupling is given as a Heaviside step. -```math -c(x) = {c_0, for x \ge x_0 \times s - 0, for x < x_0 \times s} -``` -Here, `s` is the sign of the step function, x_0 the save_position -of the step and c_0 the amplitude. -""" -function coupling_converter_heaviside_2d(x_0, c_0, s) - return (x, u) -> c_0 * (s*sign(x[2] - x_0) + 1.0)/2.0 * u -end - - -@doc raw""" - Coupling converter function for a system of two LinearScalarAdvectionEquation2D. - -The coupling is given as a linear function. -```math -c(x) = x * u(x) -``` -""" -function coupling_converter_linear_2d() - return (x, u) -> x[2]*u -end - -end # @muladd From 6b9783b8f3554b214608c893e89790f06db265ce Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 28 Jun 2024 12:43:05 +0100 Subject: [PATCH 10/43] Fixed some merge errors. --- src/Trixi.jl | 6 +--- src/callbacks_step/analysis_dg2d.jl | 35 +++---------------- src/callbacks_step/stepsize_dg2d.jl | 12 ++----- .../semidiscretization_coupled.jl | 17 --------- 4 files changed, 8 insertions(+), 62 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 741fcae49e..98c7542b12 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -225,12 +225,8 @@ export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, export lake_at_rest_error export ncomponents, eachcomponent -<<<<<<< HEAD -export TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, P4estMeshView, -======= export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh, ->>>>>>> main - T8codeMesh + P4estMeshView, T8codeMesh export DG, DGSEM, LobattoLegendreBasis, diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 434b045bed..34d1ab8582 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -30,14 +30,9 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2}, end function create_cache_analysis(analyzer, -<<<<<<< HEAD - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, -======= mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, ->>>>>>> main + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -113,14 +108,9 @@ function calc_error_norms(func, u, t, analyzer, end function calc_error_norms(func, u, t, analyzer, -<<<<<<< HEAD - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, -======= mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, ->>>>>>> main initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @unpack node_coordinates, inverse_jacobian = cache.elements @@ -187,15 +177,10 @@ function integrate_via_indices(func::Func, u, end function integrate_via_indices(func::Func, u, -<<<<<<< HEAD - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, -======= mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, ->>>>>>> main dg::DGSEM, cache, args...; normalize = true) where {Func} @unpack weights = dg.basis @@ -222,13 +207,8 @@ function integrate_via_indices(func::Func, u, end function integrate(func::Func, u, -<<<<<<< HEAD - mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, -======= mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, ->>>>>>> main + 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 @@ -256,13 +236,8 @@ function integrate(func::Func, u, end function analyze(::typeof(entropy_timederivative), du, u, t, -<<<<<<< HEAD - mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, -======= mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, ->>>>>>> main + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DG, cache) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 3901eda3b1..4263451cb8 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -77,11 +77,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, -<<<<<<< HEAD - P4estMeshView{2}, T8codeMesh{2}}, -======= - T8codeMesh{2}, StructuredMeshView{2}}, ->>>>>>> main + T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}} constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -117,11 +113,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, -<<<<<<< HEAD - P4estMeshView{2}, T8codeMesh{2}}, -======= - T8codeMesh{2}, StructuredMeshView{2}}, ->>>>>>> main + T8codeMesh{2}, StructuredMeshView{2} P4estMeshView{2}}, constant_speed::True, equations, dg::DG, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index d71370a31e..30d5da0742 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -438,7 +438,6 @@ mutable struct BoundaryConditionCoupled{NDIMS, uEltype <: Real, Indices, CouplingConverter} # NDIMST2M1 == NDIMS * 2 - 1 # Buffer for boundary values: [variable, nodes_i, nodes_j, cell_i, cell_j] -<<<<<<< HEAD u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 other_semi_index :: Int other_orientation :: Int @@ -446,15 +445,6 @@ mutable struct BoundaryConditionCoupled{NDIMS, coupling_converter :: Function function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) -======= - u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 - other_orientation :: Int - indices :: Indices - coupling_converter :: CouplingConverter - - function BoundaryConditionCoupled(other_semi_index, indices, uEltype, - coupling_converter) ->>>>>>> main NDIMS = length(indices) u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1)) @@ -466,15 +456,8 @@ mutable struct BoundaryConditionCoupled{NDIMS, other_orientation = 3 end -<<<<<<< HEAD new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index, other_orientation, indices, coupling_converter) -======= - new{NDIMS, other_semi_index, NDIMS * 2 - 1, uEltype, typeof(indices), - typeof(coupling_converter)}(u_boundary, - other_orientation, - indices, coupling_converter) ->>>>>>> main end end From b8b2573eb3b7846789ca816618eb10b944fc4c03 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 1 Jul 2024 10:33:17 +0100 Subject: [PATCH 11/43] Fixed merging issues. --- src/Trixi.jl | 2 +- .../semidiscretization_coupled.jl | 29 +++++++------------ src/solvers/dg.jl | 7 +---- src/solvers/dgsem_structured/dg_2d.jl | 7 +---- src/solvers/dgsem_tree/dg_2d.jl | 14 ++------- 5 files changed, 15 insertions(+), 44 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 98c7542b12..0578e045ad 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -131,10 +131,10 @@ include("semidiscretization/semidiscretization_hyperbolic.jl") include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl") include("semidiscretization/semidiscretization_euler_acoustics.jl") include("semidiscretization/semidiscretization_coupled.jl") +include("time_integration/time_integration.jl") include("callbacks_step/callbacks_step.jl") include("callbacks_stage/callbacks_stage.jl") include("semidiscretization/semidiscretization_euler_gravity.jl") -include("time_integration/time_integration.jl") # Special elixirs such as `convergence_test` include("auxiliary/special_elixirs.jl") diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 30d5da0742..6b009cfad2 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -438,13 +438,13 @@ mutable struct BoundaryConditionCoupled{NDIMS, uEltype <: Real, Indices, CouplingConverter} # NDIMST2M1 == NDIMS * 2 - 1 # Buffer for boundary values: [variable, nodes_i, nodes_j, cell_i, cell_j] - u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 - other_semi_index :: Int - other_orientation :: Int - indices :: Indices - coupling_converter :: Function + u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 + other_orientation :: Int + indices :: Indices + coupling_converter :: CouplingConverter - function BoundaryConditionCoupled(other_semi_index, indices, uEltype, coupling_converter) + function BoundaryConditionCoupled(other_semi_index, indices, uEltype, + coupling_converter) NDIMS = length(indices) u_boundary = Array{uEltype, NDIMS * 2 - 1}(undef, ntuple(_ -> 0, NDIMS * 2 - 1)) @@ -456,8 +456,10 @@ mutable struct BoundaryConditionCoupled{NDIMS, other_orientation = 3 end - new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices)}(u_boundary, other_semi_index, - other_orientation, indices, coupling_converter) + new{NDIMS, other_semi_index, NDIMS * 2 - 1, uEltype, typeof(indices), + typeof(coupling_converter)}(u_boundary, + other_orientation, + indices, coupling_converter) end end @@ -606,16 +608,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ j_node = j_node_start element_id = linear_indices[i_cell, j_cell] -<<<<<<< HEAD - for i in eachnode(solver) - for v in 1:size(u, 1) - x = cache.elements.node_coordinates[:, i_node, j_node, linear_indices[i_cell, j_cell]] - converted_u = boundary_condition.coupling_converter(x, u[:, i_node, j_node, linear_indices[i_cell, j_cell]]) - boundary_condition.u_boundary[v, i, cell] = converted_u[v] - # boundary_condition.u_boundary[v, i, cell] = u[v, i_node, j_node, - # linear_indices[i_cell, - # j_cell]] -======= for element_id in eachnode(solver_other) x_other = get_node_coords(node_coordinates_other, equations_other, solver_other, @@ -628,7 +620,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ for i in eachindex(u_node_converted) u_boundary[i, element_id, cell] = u_node_converted[i] ->>>>>>> main end i_node += i_node_step diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 4fefec32bc..d9b34a1536 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -436,14 +436,9 @@ function get_node_variables!(node_variables, mesh, equations, dg::DG, cache) get_node_variables!(node_variables, mesh, equations, dg.volume_integral, dg, cache) end -<<<<<<< HEAD -const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh, - P4estMeshView, T8codeMesh} -======= const MeshesDGSEM = Union{TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, - P4estMesh, T8codeMesh} ->>>>>>> main + P4estMesh, P4estMeshView, T8codeMesh} @inline function ndofs(mesh::MeshesDGSEM, dg::DG, cache) nelements(cache.elements) * nnodes(dg)^ndims(mesh) diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 5498da8e8f..36bed30868 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -58,14 +58,9 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 =# @inline function weak_form_kernel!(du, u, element, -<<<<<<< HEAD - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, -======= mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, + UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, ->>>>>>> main nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) # true * [some floating point value] == [exactly the same floating point value] diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 216ca40b34..2e0a69eb6b 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -180,13 +180,8 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, -<<<<<<< HEAD - UnstructuredMesh2D, P4estMesh{2}, - P4estMeshView{2}, T8codeMesh{2}}, -======= StructuredMeshView{2}, UnstructuredMesh2D, - P4estMesh{2}, T8codeMesh{2}}, ->>>>>>> main + P4estMesh{2}, P4estMeshView{2},T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -397,12 +392,7 @@ end @inline function fv_kernel!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, -<<<<<<< HEAD - P4estMeshView{2}, T8codeMesh{2} - }, -======= - T8codeMesh{2}}, ->>>>>>> main + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_flux_fv, dg::DGSEM, cache, element, alpha = true) @unpack fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded = cache From bedfd555c154ac11fc15e0b2633cc87b79ece22f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 2 Jul 2024 21:53:40 +0100 Subject: [PATCH 12/43] Removed redundant p4est-view code. Added a few comments, mostly todos. --- src/meshes/p4est_mesh_view.jl | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index daed597410..83dd53c1b0 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -3,21 +3,17 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} parent::P4estMesh{NDIMS, RealT} + # We need to store information about which cells are part of this mesh view. nodes::SVector - index_min::Int - index_max::Int end -function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; - index_min = 1::Int, - index_max = sizeof(unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees))::Int) where {NDIMS, RealT} - @assert index_min <= index_max - @assert all(index_min .> 0) - @assert index_max <= sizeof(unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees)) +function P4estMeshView(parent::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} - return P4estMeshView{NDIMS, RealT}(parent, parent.nodes, index_min, index_max) + return P4estMeshView{NDIMS, RealT}(parent, parent.nodes) end +# TODO: Check if this is still needed. +# At the end we will have every cell boundary with a boundary condition. # Check if mesh is periodic function isperiodic(mesh::P4estMeshView) @unpack parent = mesh @@ -34,12 +30,13 @@ end @inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::P4estMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT function Base.size(mesh::P4estMeshView) - @unpack index_min, index_max = mesh - return index_max .- index_min .+ 1 + # TODO: Implement size function. This used to be with checking index_min and inde_max. + return 0 end function Base.size(mesh::P4estMeshView, i) - @unpack index_min, index_max = mesh - return index_max[i] - index_min[i] + 1 +# @unpack index_min, index_max = mesh + # TODO: Implement size function. This used to be with checking index_min and inde_max. + return 0 end Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) From 7b7900ff2dec9a131bd496b652a6f9e54bafd025 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 4 Jul 2024 16:35:47 +0100 Subject: [PATCH 13/43] Added P4estMeshView. Runs without syntax errors. No save_solution yet. --- src/callbacks_step/save_solution_dg.jl | 2 +- src/callbacks_step/stepsize_dg2d.jl | 4 +- src/meshes/p4est_mesh_view.jl | 213 +++++++++++------- .../semidiscretization_hyperbolic.jl | 2 +- src/solvers/dgsem_p4est/containers.jl | 4 +- src/solvers/dgsem_p4est/dg_2d.jl | 6 +- src/solvers/dgsem_structured/dg_2d.jl | 2 +- 7 files changed, 141 insertions(+), 92 deletions(-) diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index deae8f7c93..1e3800a62a 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -9,7 +9,7 @@ function save_solution_file(u, time, dt, timestep, mesh::Union{SerialTreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, SerialP4estMesh, - SerialT8codeMesh}, + P4estMeshView, SerialT8codeMesh}, equations, dg::DG, cache, solution_callback, element_variables = Dict{Symbol, Any}(), diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 4263451cb8..3399f718cb 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -77,7 +77,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}} + T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection @@ -113,7 +113,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}, StructuredMeshView{2} P4estMeshView{2}}, + T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}}, constant_speed::True, equations, dg::DG, cache) @unpack contravariant_vectors, inverse_jacobian = cache.elements diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 83dd53c1b0..cb44a35f3f 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -1,15 +1,35 @@ @muladd begin #! format: noindent -mutable struct P4estMeshView{NDIMS, RealT <: Real} <: AbstractMesh{NDIMS} +mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: + AbstractMesh{NDIMS} + # Attributes from the original P4est mesh + p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t} + is_parallel :: IsParallel + ghost :: Ghost # Either PointerWrapper{p4est_ghost_t} or PointerWrapper{p8est_ghost_t} + # Coordinates at the nodes specified by the tensor product of `nodes` (NDIMS times). + # This specifies the geometry interpolation for each tree. + tree_node_coordinates::Array{RealT, NDIMSP2} # [dimension, i, j, k, tree] + nodes::SVector{NNODES, RealT} + boundary_names::Array{Symbol, 2} # [face direction, tree] + current_filename::String + unsaved_changes::Bool + p4est_partition_allow_for_coarsening::Bool + + # Attributes partaining the views. parent::P4estMesh{NDIMS, RealT} - # We need to store information about which cells are part of this mesh view. - nodes::SVector +# # We need to store information about which cells are part of this mesh view. +# nodes::SVector end function P4estMeshView(parent::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} + ghost = ghost_new_p4est(parent.p4est) + ghost_pw = PointerWrapper(ghost) - return P4estMeshView{NDIMS, RealT}(parent, parent.nodes) + return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), + typeof(parent.p4est), typeof(parent.ghost), NDIMS + 2, length(parent.nodes)}(parent.p4est, parent.is_parallel, parent.ghost, parent.tree_node_coordinates, + parent.nodes,parent.boundary_names,parent.current_filename, parent.unsaved_changes, + parent.p4est_partition_allow_for_coarsening, parent) end # TODO: Check if this is still needed. @@ -21,31 +41,51 @@ function isperiodic(mesh::P4estMeshView) end function isperiodic(mesh::P4estMeshView, dimension) - @unpack parent, index_min, index_max = mesh - return (isperiodic(parent, dimension) && - index_min[dimension] == 1 && - index_max[dimension] == size(parent, dimension)) + @unpack parent = mesh + + return (isperiodic(parent, dimension)) end @inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::P4estMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT +@inline function ntrees(mesh::P4estMeshView) + return mesh.p4est.trees.elem_count[] +end +@inline ncellsglobal(mesh::P4estMeshView) = Int(mesh.p4est.global_num_quadrants[]) function Base.size(mesh::P4estMeshView) - # TODO: Implement size function. This used to be with checking index_min and inde_max. return 0 end function Base.size(mesh::P4estMeshView, i) -# @unpack index_min, index_max = mesh # TODO: Implement size function. This used to be with checking index_min and inde_max. return 0 end Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) +function Base.show(io::IO, mesh::P4estMesh) + print(io, "P4estMeshView{", ndims(mesh), ", ", real(mesh), "}") +end + +function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMeshView) + if get(io, :compact, false) + show(io, mesh) + else + setup = [ + "#trees" => ntrees(mesh), + "current #cells" => ncellsglobal(mesh), + "polydeg" => length(mesh.nodes) - 1, + ] + summary_box(io, + "P4estMeshView{" * string(ndims(mesh)) * ", " * string(real(mesh)) * + "}", setup) + end +end + function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL) - p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) + p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn) # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes # See https://github.com/cburstedde/p4est/issues/112 - p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) + p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn) end @inline ncells(mesh::P4estMeshView) = Int(mesh.parent.p4est.local_num_quadrants[]) @@ -71,81 +111,90 @@ function init_interfaces!(interfaces, mesh::P4estMeshView) return interfaces end -# function calc_node_coordinates!(node_coordinates, element, -# cell_x, cell_y, mapping, -# mesh::P4estMeshView{2}, -# # basis::LobattoLegendreBasis) -# basis) -# @unpack nodes = basis -# @unpack parent, index_min, index_max = mesh - -# # Get cell length in reference mesh -# dx = 2 / size(parent, 1) -# dy = 2 / size(parent, 2) - -# # Calculate index offsets with respect to parent -# parent_offset_x = index_min[1] - 1 -# parent_offset_y = index_min[2] - 1 - -# # Calculate node coordinates of reference mesh -# cell_x_offset = -1 + (cell_x - 1 + parent_offset_x) * dx + dx / 2 -# cell_y_offset = -1 + (cell_y - 1 + parent_offset_y) * dy + dy / 2 -# p4est_tree_t -# for j in eachnode(basis), i in eachnode(basis) -# # node_coordinates are the mapped reference node_coordinates -# node_coordinates[:, i, j, element] .= mapping(cell_x_offset + dx / 2 * nodes[i], -# cell_y_offset + dy / 2 * nodes[j]) -# end -# end - # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, mesh::P4estMeshView{2}, nodes::AbstractVector) - @unpack parent, index_min, index_max = mesh - - # 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(parent), - StaticInt(2), static_length(nodes), static_length(parent.nodes)) - matrix1 = StrideArray(undef, real(parent), - static_length(nodes), static_length(parent.nodes)) - matrix2 = similar(matrix1) - baryweights_in = barycentric_weights(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, parent.p4est.trees) - - for tree in eachindex(trees) - offset = trees[tree].quadrants_offset - quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants) - - for i in eachindex(quadrants) - element = offset + i - quad = quadrants[i] - - quad_length = p4est_quadrant_len(quad.level) / p4est_root_len - - nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ - quad.x / p4est_root_len) .- 1 - nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ - quad.y / p4est_root_len) .- 1 - polynomial_interpolation_matrix!(matrix1, parent.nodes, nodes_out_x, - baryweights_in) - polynomial_interpolation_matrix!(matrix2, parent.nodes, nodes_out_y, - baryweights_in) - - multiply_dimensionwise!(view(node_coordinates, :, :, :, element), - matrix1, matrix2, - view(parent.tree_node_coordinates, :, :, :, tree), - tmp1) - end - end + @unpack parent = mesh + +# # 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(parent), +# StaticInt(2), static_length(nodes), static_length(parent.nodes)) +# matrix1 = StrideArray(undef, real(parent), +# static_length(nodes), static_length(parent.nodes)) +# matrix2 = similar(matrix1) +# baryweights_in = barycentric_weights(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, parent.p4est.trees) +# +# for tree in eachindex(trees) +# offset = trees[tree].quadrants_offset +# quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants) +# +# for i in eachindex(quadrants) +# element = offset + i +# quad = quadrants[i] +# +# quad_length = p4est_quadrant_len(quad.level) / p4est_root_len +# +# nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ +# quad.x / p4est_root_len) .- 1 +# nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ +# quad.y / p4est_root_len) .- 1 +# polynomial_interpolation_matrix!(matrix1, parent.nodes, nodes_out_x, +# baryweights_in) +# polynomial_interpolation_matrix!(matrix2, parent.nodes, nodes_out_y, +# baryweights_in) +# +# multiply_dimensionwise!(view(node_coordinates, :, :, :, element), +# matrix1, matrix2, +# view(parent.tree_node_coordinates, :, :, :, tree), +# tmp1) +# end +# end return node_coordinates end + +function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; + system = "") + # Create output directory (if it does not exist) + mkpath(output_directory) + + # Determine file name based on existence of meaningful time step + if timestep > 0 + filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) + p4est_filename = @sprintf("p4est_data_%s_%06d", system, 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.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.tree_node_coordinates + file["nodes"] = Vector(mesh.nodes) # the mesh uses `SVector`s for the nodes + # to increase the runtime performance + # but HDF5 can only handle plain arrays + file["boundary_names"] = mesh.boundary_names .|> String + end + + return filename +end end # @muladd diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 7c82a132a0..609109875e 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/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index e0af3be794..b61da518cc 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -271,7 +271,7 @@ function init_boundaries(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa return boundaries end -function init_boundaries!(boundaries, mesh::P4estMesh) +function init_boundaries!(boundaries, mesh::Union{P4estMesh, P4estMeshView}) init_surfaces!(nothing, nothing, boundaries, mesh) return boundaries @@ -514,7 +514,7 @@ function init_surfaces_iter_face_inner(info, user_data) return nothing end -function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh) +function init_surfaces!(interfaces, mortars, boundaries, mesh::Union{P4estMesh, P4estMeshView}) # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index a9f735a64b..6980799f4e 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -385,7 +385,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 @@ -452,7 +452,7 @@ function prolong2mortars!(cache, u, end function calc_mortar_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -626,7 +626,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_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 36bed30868..6b462c6b21 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -620,7 +620,7 @@ 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 From 11190d657b68b78ad1bcf9fe73ac09de7013e000 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 8 Jul 2024 18:25:32 +0100 Subject: [PATCH 14/43] Added P4estMeshView to a few more function as union with P4estMesh. This version now runs. --- src/meshes/p4est_mesh_view.jl | 107 ++++++++------------------ src/solvers/dgsem_p4est/containers.jl | 4 +- 2 files changed, 34 insertions(+), 77 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index cb44a35f3f..76b5603083 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -42,7 +42,6 @@ end function isperiodic(mesh::P4estMeshView, dimension) @unpack parent = mesh - return (isperiodic(parent, dimension)) end @@ -52,13 +51,13 @@ end return mesh.p4est.trees.elem_count[] end @inline ncellsglobal(mesh::P4estMeshView) = Int(mesh.p4est.global_num_quadrants[]) -function Base.size(mesh::P4estMeshView) - return 0 -end -function Base.size(mesh::P4estMeshView, i) - # TODO: Implement size function. This used to be with checking index_min and inde_max. - return 0 -end +# function Base.size(mesh::P4estMeshView) +# return 0 +# end +# function Base.size(mesh::P4estMeshView, i) +# # TODO: Implement size function. This used to be with checking index_min and inde_max. +# return 0 +# end Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) @@ -88,79 +87,37 @@ function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL) p4est_balance(mesh.p4est, P4EST_CONNECT_FACE, init_fn) end -@inline ncells(mesh::P4estMeshView) = Int(mesh.parent.p4est.local_num_quadrants[]) - -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 - -function init_interfaces!(interfaces, mesh::P4estMeshView) - init_surfaces!(interfaces, nothing, nothing, mesh.parent) - - return interfaces -end - -# Interpolate tree_node_coordinates to each quadrant at the specified nodes -function calc_node_coordinates!(node_coordinates, - mesh::P4estMeshView{2}, - nodes::AbstractVector) - @unpack parent = mesh +@inline ncells(mesh::P4estMeshView) = Int(mesh.p4est.local_num_quadrants[]) -# # 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(parent), -# StaticInt(2), static_length(nodes), static_length(parent.nodes)) -# matrix1 = StrideArray(undef, real(parent), -# static_length(nodes), static_length(parent.nodes)) -# matrix2 = similar(matrix1) -# baryweights_in = barycentric_weights(parent.nodes) -# -# # Macros from `p4est` -# p4est_root_len = 1 << P4EST_MAXLEVEL -# p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) +# 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))) # -# trees = unsafe_wrap_sc(p4est_tree_t, parent.p4est.trees) +# # interfaces, mortars, boundaries +# user_data = [0, 0, 0] # -# for tree in eachindex(trees) -# offset = trees[tree].quadrants_offset -# quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants) +# iterate_p4est(mesh.parent.p4est, user_data; iter_face_c = iter_face_c) # -# for i in eachindex(quadrants) -# element = offset + i -# quad = quadrants[i] +# # Return counters +# return (interfaces = user_data[1], +# mortars = user_data[2], +# boundaries = user_data[3]) +# end + +# function init_interfaces!(interfaces, mesh::P4estMeshView) +# init_surfaces!(interfaces, nothing, nothing, mesh.parent) # -# quad_length = p4est_quadrant_len(quad.level) / p4est_root_len +# return interfaces +# end + +# # Interpolate tree_node_coordinates to each quadrant at the specified nodes +# function calc_node_coordinates!(node_coordinates, +# mesh::P4estMeshView{2}, +# nodes::AbstractVector) +# @unpack parent = mesh # -# 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, parent.nodes, nodes_out_x, -# baryweights_in) -# polynomial_interpolation_matrix!(matrix2, parent.nodes, nodes_out_y, -# baryweights_in) -# -# multiply_dimensionwise!(view(node_coordinates, :, :, :, element), -# matrix1, matrix2, -# view(parent.tree_node_coordinates, :, :, :, tree), -# tmp1) -# end -# end - - return node_coordinates -end +# return node_coordinates +# end function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; system = "") diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index b61da518cc..72b0741bd7 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -197,7 +197,7 @@ function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa 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 @@ -689,7 +689,7 @@ function cfunction(::typeof(count_surfaces_iter_face), ::Val{3}) (Ptr{p8est_iter_face_info_t}, Ptr{Cvoid})) end -function count_required_surfaces(mesh::P4estMesh) +function count_required_surfaces(mesh::Union{P4estMesh, P4estMeshView}) # Let `p4est` iterate over all interfaces and call count_surfaces_iter_face iter_face_c = cfunction(count_surfaces_iter_face, Val(ndims(mesh))) From ba965eb9a350d635bf19e478f57e8dec30461520 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 9 Jul 2024 17:40:54 +0100 Subject: [PATCH 15/43] Removed commented out code. --- src/meshes/p4est_mesh_view.jl | 39 ----------------------------------- 1 file changed, 39 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 76b5603083..1ba28c3e3c 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -18,8 +18,6 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2 # Attributes partaining the views. parent::P4estMesh{NDIMS, RealT} -# # We need to store information about which cells are part of this mesh view. -# nodes::SVector end function P4estMeshView(parent::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} @@ -51,13 +49,6 @@ end return mesh.p4est.trees.elem_count[] end @inline ncellsglobal(mesh::P4estMeshView) = Int(mesh.p4est.global_num_quadrants[]) -# function Base.size(mesh::P4estMeshView) -# return 0 -# end -# function Base.size(mesh::P4estMeshView, i) -# # TODO: Implement size function. This used to be with checking index_min and inde_max. -# return 0 -# end Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) @@ -89,36 +80,6 @@ end @inline ncells(mesh::P4estMeshView) = Int(mesh.p4est.local_num_quadrants[]) -# 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 - -# function init_interfaces!(interfaces, mesh::P4estMeshView) -# init_surfaces!(interfaces, nothing, nothing, mesh.parent) -# -# return interfaces -# end - -# # Interpolate tree_node_coordinates to each quadrant at the specified nodes -# function calc_node_coordinates!(node_coordinates, -# mesh::P4estMeshView{2}, -# nodes::AbstractVector) -# @unpack parent = mesh -# -# return node_coordinates -# end - function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; system = "") # Create output directory (if it does not exist) From 7583e09f24f1c897b63ebe2411819b4670635a25 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 9 Jul 2024 17:43:26 +0100 Subject: [PATCH 16/43] Spelling error. --- 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 1ba28c3e3c..f5393a6958 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -16,7 +16,7 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2 unsaved_changes::Bool p4est_partition_allow_for_coarsening::Bool - # Attributes partaining the views. + # Attributes pertaining the views. parent::P4estMesh{NDIMS, RealT} end From cf0d84ccd4560485d97d20fd6572e091269fd4e6 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 9 Jul 2024 18:22:40 +0100 Subject: [PATCH 17/43] Added test eleixir for 0th version of p4estview. --- .../elixir_advection_p4estview.jl | 64 +++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_advection_p4estview.jl diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl new file mode 100644 index 0000000000..29205ed404 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -0,0 +1,64 @@ +# The same setup as tree_2d_dgsem/elixir_advection_basic.jl +# to verify the StructuredMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) + +trees_per_dimension = (8, 8) + +# Create P4estMesh with 8 x 8 trees and 16 x 16 elements +parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 1) +mesh = P4estMeshView(parent_mesh) + +# 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, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() From 05d2697c196c07523e9346d9c5df61db962e5daf Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 9 Jul 2024 18:23:17 +0100 Subject: [PATCH 18/43] Added p4estview test in 2d. --- 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 f4924bdcf7..085456974e 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -125,6 +125,20 @@ end end end +@trixi_testset "elixir_advection_p4estview.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4estview.jl"), + l2=[8.311947672995606e-6], + linf=[6.627000277448225e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"), From 15ccd82e84d64ddd84eeb3460cabac36cfae4893 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 9 Jul 2024 18:32:15 +0100 Subject: [PATCH 19/43] Applied autoformatter. --- .../elixir_advection_p4estview.jl | 3 +- src/callbacks_step/analysis_dg2d.jl | 21 ++++++++----- src/callbacks_step/stepsize_dg2d.jl | 2 +- src/meshes/p4est_mesh_view.jl | 23 +++++++++----- .../semidiscretization_hyperbolic.jl | 3 +- src/solvers/dgsem_p4est/containers.jl | 14 +++++---- src/solvers/dgsem_p4est/containers_2d.jl | 8 +++-- src/solvers/dgsem_p4est/dg.jl | 3 +- src/solvers/dgsem_p4est/dg_2d.jl | 30 ++++++++++++------- src/solvers/dgsem_structured/dg_2d.jl | 6 ++-- src/solvers/dgsem_tree/dg_2d.jl | 18 ++++++----- src/solvers/dgsem_tree/indicators_2d.jl | 3 +- src/solvers/dgsem_unstructured/dg_2d.jl | 6 ++-- test/Project.toml | 6 ++++ 14 files changed, 96 insertions(+), 50 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl index 29205ed404..f142fa4dc4 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -20,7 +20,8 @@ trees_per_dimension = (8, 8) # Create P4estMesh with 8 x 8 trees and 16 x 16 elements parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3, - coordinates_min = coordinates_min, coordinates_max = coordinates_max, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, initial_refinement_level = 1) mesh = P4estMeshView(parent_mesh) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 34d1ab8582..83532332ae 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -32,7 +32,8 @@ end function create_cache_analysis(analyzer, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, equations, dg::DG, cache, RealT, uEltype) @@ -109,7 +110,8 @@ end function calc_error_norms(func, u, t, analyzer, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + UnstructuredMesh2D, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, equations, initial_condition, dg::DGSEM, cache, cache_analysis) @unpack vandermonde, weights = analyzer @@ -178,7 +180,8 @@ end function integrate_via_indices(func::Func, u, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @@ -208,7 +211,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 @@ -218,7 +222,7 @@ function integrate(func::Func, u, end function integrate(func::Func, u, - mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}}, + mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}}, equations, equations_parabolic, dg::DGSEM, cache, cache_parabolic; normalize = true) where {Func} @@ -237,7 +241,8 @@ end function analyze(::typeof(entropy_timederivative), du, u, t, 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) # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ integrate_via_indices(u, mesh, equations, dg, cache, @@ -282,7 +287,7 @@ end function analyze(::Val{:l2_divb}, du, u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - P4estMeshView{2}, T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) @unpack contravariant_vectors = cache.elements integrate_via_indices(u, mesh, equations, dg, cache, cache, @@ -350,7 +355,7 @@ end function analyze(::Val{:linf_divb}, du, u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - P4estMeshView{2}, T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) @unpack derivative_matrix, weights = dg.basis @unpack contravariant_vectors = cache.elements diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 3399f718cb..3dcad46729 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -77,7 +77,7 @@ end function max_dt(u, t, mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}}, + T8codeMesh{2}, StructuredMeshView{2}, P4estMeshView{2}}, constant_speed::False, equations, dg::DG, cache) # to avoid a division by zero if the speed vanishes everywhere, # e.g. for steady-state linear advection diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index f5393a6958..51714eae58 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -1,8 +1,9 @@ @muladd begin #! format: noindent -mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: - AbstractMesh{NDIMS} +mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, + NNODES} <: + AbstractMesh{NDIMS} # Attributes from the original P4est mesh p4est :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t} is_parallel :: IsParallel @@ -24,10 +25,17 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} ghost = ghost_new_p4est(parent.p4est) ghost_pw = PointerWrapper(ghost) - return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), - typeof(parent.p4est), typeof(parent.ghost), NDIMS + 2, length(parent.nodes)}(parent.p4est, parent.is_parallel, parent.ghost, parent.tree_node_coordinates, - parent.nodes,parent.boundary_names,parent.current_filename, parent.unsaved_changes, - parent.p4est_partition_allow_for_coarsening, parent) + return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), + typeof(parent.is_parallel), + typeof(parent.p4est), typeof(parent.ghost), NDIMS + 2, + length(parent.nodes)}(parent.p4est, parent.is_parallel, + parent.ghost, + parent.tree_node_coordinates, + parent.nodes, parent.boundary_names, + parent.current_filename, + parent.unsaved_changes, + parent.p4est_partition_allow_for_coarsening, + parent) end # TODO: Check if this is still needed. @@ -87,7 +95,8 @@ function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; # Determine file name based on existence of meaningful time step if timestep > 0 - filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) + filename = joinpath(output_directory, + @sprintf("mesh_%s_%06d.h5", system, timestep)) p4est_filename = @sprintf("p4est_data_%s_%06d", system, timestep) else filename = joinpath(output_directory, "mesh.h5") diff --git a/src/semidiscretization/semidiscretization_hyperbolic.jl b/src/semidiscretization/semidiscretization_hyperbolic.jl index 609109875e..50c5aef114 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 72b0741bd7..293a4907d3 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, RealT}, P4estMeshView{NDIMS, RealT}, - T8codeMesh{NDIMS, RealT}}, + T8codeMesh{NDIMS, RealT}}, equations, basis, ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} @@ -167,7 +167,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) @@ -242,7 +243,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) @@ -373,7 +375,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) @@ -514,7 +517,8 @@ function init_surfaces_iter_face_inner(info, user_data) return nothing end -function init_surfaces!(interfaces, mortars, boundaries, mesh::Union{P4estMesh, P4estMeshView}) +function init_surfaces!(interfaces, mortars, boundaries, + mesh::Union{P4estMesh, P4estMeshView}) # Let `p4est` iterate over all interfaces and call init_surfaces_iter_face iter_face_c = cfunction(init_surfaces_iter_face, Val(ndims(mesh))) user_data = InitSurfacesIterFaceUserData(interfaces, mortars, boundaries, mesh) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index fe196562cb..5fe6ac0fb0 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. @@ -37,7 +39,7 @@ end # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, - mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, nodes::AbstractVector) # We use `StrideArray`s here since these buffers are used in performance-critical # places and the additional information passed to the compiler makes them faster diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index 48265c17da..4999c5f059 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -8,7 +8,8 @@ # This method is called when a SemidiscretizationHyperbolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. -function create_cache(mesh::Union{P4estMesh, P4estMeshView}, equations::AbstractEquations, dg::DG, ::Any, +function create_cache(mesh::Union{P4estMesh, P4estMeshView}, + equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype <: Real} # Make sure to balance the `p4est` before creating any containers # in case someone has tampered with the `p4est` after creating the mesh diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 6980799f4e..8c77c02461 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)}, @@ -58,7 +59,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}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) @unpack interfaces = cache index_range = eachnode(dg) @@ -114,7 +115,8 @@ function prolong2interfaces!(cache, u, end function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms, equations, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.interfaces @@ -182,7 +184,8 @@ end # Inlined version of the interface flux computation for conservation laws @inline function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -206,7 +209,8 @@ end # Inlined version of the interface flux computation for equations with conservative and nonconservative terms @inline function calc_interface_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMesh{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMesh{2}, + T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, interface_index, normal_direction, @@ -276,7 +280,7 @@ function prolong2boundaries!(cache, u, end function calc_boundary_flux!(cache, t, boundary_condition::BC, boundary_indexing, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral, dg::DG) where {BC} @unpack boundaries = cache @unpack surface_flux_values = cache.elements @@ -312,7 +316,8 @@ end # inlined version of the boundary flux calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -343,7 +348,8 @@ end # inlined version of the boundary flux with nonconservative terms calculation along a physical interface @inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, i_index, j_index, @@ -385,7 +391,8 @@ end end function prolong2mortars!(cache, u, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DGSEM) @unpack neighbor_ids, node_indices = cache.mortars @@ -452,7 +459,7 @@ function prolong2mortars!(cache, u, end function calc_mortar_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, mortar_l2::LobattoLegendreMortarL2, surface_integral, dg::DG, cache) @@ -626,7 +633,8 @@ end end function calc_surface_integral!(du, u, - mesh::Union{P4estMesh{2}, 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_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 03a7094a42..e834c1ad86 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -59,7 +59,8 @@ See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-17 @inline function weak_form_kernel!(du, u, element, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2}, + UnstructuredMesh2D, P4estMesh{2}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, dg::DGSEM, cache, alpha = true) @@ -625,7 +626,8 @@ end function apply_jacobian!(du, mesh::Union{StructuredMesh{2}, StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{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 7743c80de2..92dc6e758e 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -45,7 +45,8 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMesh end function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype) element_ids_dg = Int[] element_ids_dgfv = Int[] @@ -71,7 +72,8 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMe end function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG, uEltype) A3dp1_x = Array{uEltype, 3} @@ -93,7 +95,7 @@ end # The methods below are specialized on the mortar type # and called from the basic `create_cache` method at the top. function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal performance using different types MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, @@ -111,7 +113,8 @@ end # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer? function rhs!(du, u, t, - mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, + mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, + equations, initial_condition, boundary_conditions, source_terms::Source, dg::DG, cache) where {Source} # Reset du @@ -182,7 +185,8 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, - P4estMesh{2}, P4estMeshView{2},T8codeMesh{2}}, + P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralWeakForm, dg::DGSEM, cache) @@ -237,7 +241,7 @@ function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, - P4estMeshView{2}, T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralFluxDifferencing, dg::DGSEM, cache) @@ -334,7 +338,7 @@ end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}, - P4estMeshView{2}, T8codeMesh{2}}, + P4estMeshView{2}, T8codeMesh{2}}, nonconservative_terms, equations, volume_integral::VolumeIntegralShockCapturingHG, dg::DGSEM, cache) diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 52de8408a0..1b540378a8 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -98,7 +98,8 @@ end end # Diffuse alpha values by setting each alpha to at least 50% of neighboring elements' alpha -function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, alpha, +function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, P4estMeshView{2}, + T8codeMesh{2}}, alpha, alpha_tmp, dg, cache) # Copy alpha values such that smoothing is indpedenent of the element access order diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 04b7aad2df..7bdd5066cf 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -311,14 +311,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 diff --git a/test/Project.toml b/test/Project.toml index c8ae33a40a..35eaa83dd4 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,7 @@ [deps] +AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Convex = "f65535da-76fb-5f13-bab9-19810c17039a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" @@ -8,13 +10,17 @@ ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" +Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95" [compat] Aqua = "0.8" From 836fceeedfbd15ad7c42de68a08085ad6a89f323 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 11 Jul 2024 18:55:34 +0100 Subject: [PATCH 20/43] Moved calculation of node coordinates for P4estMeshView. --- src/meshes/p4est_mesh_view.jl | 50 ++++++++++++++++++++++++ src/solvers/dgsem_p4est/containers_2d.jl | 2 +- 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 51714eae58..6ec4d84b4b 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -19,12 +19,14 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2 # Attributes pertaining the views. parent::P4estMesh{NDIMS, RealT} + mesh_cells::Array{Int32} end function P4estMeshView(parent::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} ghost = ghost_new_p4est(parent.p4est) ghost_pw = PointerWrapper(ghost) + @autoinfiltrate return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), typeof(parent.p4est), typeof(parent.ghost), NDIMS + 2, @@ -88,6 +90,54 @@ end @inline ncells(mesh::P4estMeshView) = Int(mesh.p4est.local_num_quadrants[]) +function calc_node_coordinates!(node_coordinates, + mesh::P4estMeshView{2}, + nodes::AbstractVector) + # We use `StrideArray`s here since these buffers are used in performance-critical + # places and the additional information passed to the compiler makes them faster + # than native `Array`s. + tmp1 = StrideArray(undef, real(mesh), + StaticInt(2), static_length(nodes), static_length(mesh.nodes)) + matrix1 = StrideArray(undef, real(mesh), + static_length(nodes), static_length(mesh.nodes)) + matrix2 = similar(matrix1) + baryweights_in = barycentric_weights(mesh.nodes) + + # Macros from `p4est` + p4est_root_len = 1 << P4EST_MAXLEVEL + p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l) + + trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees) + + for tree in eachindex(trees) + offset = trees[tree].quadrants_offset + quadrants = unsafe_wrap_sc(p4est_quadrant_t, trees[tree].quadrants) + + for i in eachindex(quadrants) + element = offset + i + quad = quadrants[i] + + quad_length = p4est_quadrant_len(quad.level) / p4est_root_len + + nodes_out_x = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.x / p4est_root_len) .- 1 + nodes_out_y = 2 * (quad_length * 1 / 2 * (nodes .+ 1) .+ + quad.y / p4est_root_len) .- 1 + polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, + baryweights_in) + polynomial_interpolation_matrix!(matrix2, mesh.nodes, nodes_out_y, + baryweights_in) + + multiply_dimensionwise!(view(node_coordinates, :, :, :, element), + matrix1, matrix2, + view(mesh.tree_node_coordinates, :, :, :, tree), + tmp1) + end + end + + return node_coordinates +end + function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; system = "") # Create output directory (if it does not exist) diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 5fe6ac0fb0..8e842a9e8d 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -39,7 +39,7 @@ end # Interpolate tree_node_coordinates to each quadrant at the specified nodes function calc_node_coordinates!(node_coordinates, - mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + mesh::P4estMesh{2}, nodes::AbstractVector) # We use `StrideArray`s here since these buffers are used in performance-critical # places and the additional information passed to the compiler makes them faster From 7e6c25f8e36be2a6bde9729dc1beedd778bc82c3 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 15 Jul 2024 18:36:06 +0100 Subject: [PATCH 21/43] Added features to p4estMeshView. Currently still buggy. --- src/meshes/p4est_mesh.jl | 1 + src/meshes/p4est_mesh_view.jl | 47 ++++++++++++++++--- .../sort_boundary_conditions.jl | 1 + 3 files changed, 42 insertions(+), 7 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 2046220aec..5fc3215b2d 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -88,6 +88,7 @@ end @inline Base.ndims(::P4estMesh{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT +@inline Base.size(mesh::P4estMesh) = size(mesh.tree_node_coordinates)[end] @inline function ntrees(mesh::P4estMesh) return mesh.p4est.trees.elem_count[] diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 6ec4d84b4b..524260590f 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -19,25 +19,51 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2 # Attributes pertaining the views. parent::P4estMesh{NDIMS, RealT} - mesh_cells::Array{Int32} + indices_min::NTuple{NDIMS, Int} + indices_max::NTuple{NDIMS, Int} +# view_cells::Array{Bool} end -function P4estMeshView(parent::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} +# function P4estMeshView(parent::P4estMesh{NDIMS, RealT}, view_cells::Array{Bool}) where {NDIMS, RealT} +function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; + indices_min = ntuple(_ -> 1, Val(NDIMS)), + indices_max = size(parent), + coordinates_min = nothing, coordinates_max = nothing) where {NDIMS, RealT} + @assert indices_min <= indices_max + @assert all(indices_min .> 0) +# @assert indices_max <= size(parent) + ghost = ghost_new_p4est(parent.p4est) ghost_pw = PointerWrapper(ghost) - @autoinfiltrate + # Extract mapping + mapping = coordinates2mapping(coordinates_min, coordinates_max) + + trees_per_dimension = indices_max .- indices_min + tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, + ntuple(_ -> length(parent.nodes), + NDIMS)..., + prod(trees_per_dimension)) + calc_tree_node_coordinates!(tree_node_coordinates, parent.nodes, mapping, + trees_per_dimension) + + connectivity = connectivity_structured(trees_per_dimension..., (false, false)) + + # TODO: The initial refinment level of 1 should not be hard-coded. + p4est = new_p4est(connectivity, 1) + p4est_pw = PointerWrapper(p4est) + return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), - typeof(parent.p4est), typeof(parent.ghost), NDIMS + 2, - length(parent.nodes)}(parent.p4est, parent.is_parallel, + typeof(p4est_pw), typeof(parent.ghost), NDIMS + 2, + length(parent.nodes)}(PointerWrapper(p4est), parent.is_parallel, parent.ghost, - parent.tree_node_coordinates, + tree_node_coordinates, parent.nodes, parent.boundary_names, parent.current_filename, parent.unsaved_changes, parent.p4est_partition_allow_for_coarsening, - parent) + parent, indices_min, indices_max) end # TODO: Check if this is still needed. @@ -138,6 +164,13 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL) + p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) + # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes + # See https://github.com/cburstedde/p4est/issues/112 + p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) +end + function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; system = "") # Create output directory (if it does not exist) diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index 2c2c6876d7..de88617756 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -71,6 +71,7 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N MPI.Gatherv!(send_buffer, nothing, mpi_root(), mpi_comm()) end else + @autoinfiltrate for key in keys(boundary_dictionary) if !(key in unique_names) error("Key $(repr(key)) is not a valid boundary name. " * From 8825868d7c44e385d77b49b3e90b777057c46ab6 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 16 Jul 2024 17:28:44 +0100 Subject: [PATCH 22/43] Added correct boundary names for p4estMeshViews. --- src/meshes/p4est_mesh_view.jl | 11 ++++++++--- .../dgsem_unstructured/sort_boundary_conditions.jl | 1 - 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 524260590f..9042996a37 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -47,19 +47,24 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; calc_tree_node_coordinates!(tree_node_coordinates, parent.nodes, mapping, trees_per_dimension) - connectivity = connectivity_structured(trees_per_dimension..., (false, false)) + connectivity = connectivity_structured(trees_per_dimension..., (false, true)) # TODO: The initial refinment level of 1 should not be hard-coded. p4est = new_p4est(connectivity, 1) p4est_pw = PointerWrapper(p4est) + # Non-periodic boundaries + boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) + + structured_boundary_names!(boundary_names, trees_per_dimension, (false, true)) + return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), typeof(p4est_pw), typeof(parent.ghost), NDIMS + 2, length(parent.nodes)}(PointerWrapper(p4est), parent.is_parallel, parent.ghost, tree_node_coordinates, - parent.nodes, parent.boundary_names, + parent.nodes, boundary_names, parent.current_filename, parent.unsaved_changes, parent.p4est_partition_allow_for_coarsening, @@ -88,7 +93,7 @@ end Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) -function Base.show(io::IO, mesh::P4estMesh) +function Base.show(io::IO, mesh::P4estMeshView) print(io, "P4estMeshView{", ndims(mesh), ", ", real(mesh), "}") end diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index de88617756..2c2c6876d7 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -71,7 +71,6 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N MPI.Gatherv!(send_buffer, nothing, mpi_root(), mpi_comm()) end else - @autoinfiltrate for key in keys(boundary_dictionary) if !(key in unique_names) error("Key $(repr(key)) is not a valid boundary name. " * From 897437044560cb4a1700242e992330a33015844a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 17 Jul 2024 16:33:26 +0100 Subject: [PATCH 23/43] Reverted Project.toml from teststo an older version. --- test/Project.toml | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 35eaa83dd4..c8ae33a40a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,5 @@ [deps] -AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Convex = "f65535da-76fb-5f13-bab9-19810c17039a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" @@ -10,17 +8,13 @@ ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" -Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95" [compat] Aqua = "0.8" From b11698aa4c9f908bf1604185000e1d10ea4a8557 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 17 Jul 2024 16:58:27 +0100 Subject: [PATCH 24/43] Cleaning up of p4estMeshViews. Updated example and test results. --- .../elixir_advection_p4estview.jl | 11 +++++--- src/meshes/p4est_mesh_view.jl | 26 ++++++++++++------- test/test_p4est_2d.jl | 4 +-- 3 files changed, 27 insertions(+), 14 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl index f142fa4dc4..b4efb5e4ee 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -22,12 +22,17 @@ trees_per_dimension = (8, 8) parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3, coordinates_min = coordinates_min, coordinates_max = coordinates_max, - initial_refinement_level = 1) -mesh = P4estMeshView(parent_mesh) + initial_refinement_level = 1, periodicity = (false, true)) +mesh = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (4, 8), + coordinates_min = coordinates_min, coordinates_max = (0.0, 1.0), + periodicity = (false, true)) + +boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_convergence_test), + :x_pos => BoundaryConditionDirichlet(initial_condition_convergence_test)) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, - solver) + solver, boundary_conditions=boundary_conditions) ############################################################################### # ODE solvers, callbacks etc. diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 9042996a37..a308839889 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -21,25 +21,33 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2 parent::P4estMesh{NDIMS, RealT} indices_min::NTuple{NDIMS, Int} indices_max::NTuple{NDIMS, Int} -# view_cells::Array{Bool} end # function P4estMeshView(parent::P4estMesh{NDIMS, RealT}, view_cells::Array{Bool}) where {NDIMS, RealT} function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; indices_min = ntuple(_ -> 1, Val(NDIMS)), indices_max = size(parent), - coordinates_min = nothing, coordinates_max = nothing) where {NDIMS, RealT} + coordinates_min = nothing, coordinates_max = nothing, + periodicity = (true, true)) where {NDIMS, RealT} + trees_per_dimension = indices_max .- indices_min + @assert indices_min <= indices_max @assert all(indices_min .> 0) -# @assert indices_max <= size(parent) + @assert prod(trees_per_dimension) <= size(parent) ghost = ghost_new_p4est(parent.p4est) ghost_pw = PointerWrapper(ghost) # Extract mapping + if isnothing(coordinates_min) + coordinates_min = minimum(parent.tree_node_coordinates) + end + if isnothing(coordinates_max) + coordinates_max = maximum(parent.tree_node_coordinates) + end + mapping = coordinates2mapping(coordinates_min, coordinates_max) - trees_per_dimension = indices_max .- indices_min tree_node_coordinates = Array{RealT, NDIMS + 2}(undef, NDIMS, ntuple(_ -> length(parent.nodes), NDIMS)..., @@ -47,16 +55,16 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; calc_tree_node_coordinates!(tree_node_coordinates, parent.nodes, mapping, trees_per_dimension) - connectivity = connectivity_structured(trees_per_dimension..., (false, true)) + connectivity = connectivity_structured(trees_per_dimension..., periodicity) - # TODO: The initial refinment level of 1 should not be hard-coded. + # TODO: The initial refinement level of 1 should not be hard-coded. p4est = new_p4est(connectivity, 1) p4est_pw = PointerWrapper(p4est) # Non-periodic boundaries boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) - structured_boundary_names!(boundary_names, trees_per_dimension, (false, true)) + structured_boundary_names!(boundary_names, trees_per_dimension, periodicity) return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), @@ -184,8 +192,8 @@ function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; # Determine file name based on existence of meaningful time step if timestep > 0 filename = joinpath(output_directory, - @sprintf("mesh_%s_%06d.h5", system, timestep)) - p4est_filename = @sprintf("p4est_data_%s_%06d", system, timestep) + @sprintf("mesh_%s_%09d.h5", system, timestep)) + p4est_filename = @sprintf("p4est_data_%s_%09d", system, timestep) else filename = joinpath(output_directory, "mesh.h5") p4est_filename = "p4est_data" diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 085456974e..fbd7b065ac 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -127,8 +127,8 @@ end @trixi_testset "elixir_advection_p4estview.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4estview.jl"), - l2=[8.311947672995606e-6], - linf=[6.627000277448225e-5]) + l2=[1.842677076881425e-5], + linf=[0.00013873258442609337]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 8f435533812f48931f07689a2fff0763eda9ac4f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 17 Jul 2024 17:00:52 +0100 Subject: [PATCH 25/43] Applied autoformatter. --- examples/p4est_2d_dgsem/elixir_advection_p4estview.jl | 2 +- src/meshes/p4est_mesh_view.jl | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl index b4efb5e4ee..7187634b76 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -32,7 +32,7 @@ boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_conditio # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, - solver, boundary_conditions=boundary_conditions) + solver, boundary_conditions = boundary_conditions) ############################################################################### # ODE solvers, callbacks etc. diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index a308839889..e55bbde9fa 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -69,7 +69,8 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), typeof(p4est_pw), typeof(parent.ghost), NDIMS + 2, - length(parent.nodes)}(PointerWrapper(p4est), parent.is_parallel, + length(parent.nodes)}(PointerWrapper(p4est), + parent.is_parallel, parent.ghost, tree_node_coordinates, parent.nodes, boundary_names, From b8b992f79377db577a94172e7f1c6aa0eb605c37 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 18 Jul 2024 13:51:25 +0100 Subject: [PATCH 26/43] Changed mesh anem so that trixi2vtk can read it. --- 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 e55bbde9fa..4c208ac0ce 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -208,7 +208,7 @@ function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; # Open file (clobber existing content) h5open(filename, "w") do file # Add context information as attributes - attributes(file)["mesh_type"] = get_name(mesh) + attributes(file)["mesh_type"] = get_name(mesh.parent) attributes(file)["ndims"] = ndims(mesh) attributes(file)["p4est_file"] = p4est_filename From 09c8ee2e32f8ddd349fd148ef1849a80a76a2c51 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 22 Jul 2024 15:36:58 +0100 Subject: [PATCH 27/43] Removed redundant isperiodic checking function. --- src/meshes/p4est_mesh_view.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 4c208ac0ce..25c0ba7c68 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -80,19 +80,6 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; parent, indices_min, indices_max) end -# TODO: Check if this is still needed. -# At the end we will have every cell boundary with a boundary condition. -# Check if mesh is periodic -function isperiodic(mesh::P4estMeshView) - @unpack parent = mesh - return isperiodic(parent) && size(parent) == size(mesh) -end - -function isperiodic(mesh::P4estMeshView, dimension) - @unpack parent = mesh - return (isperiodic(parent, dimension)) -end - @inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS @inline Base.real(::P4estMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT @inline function ntrees(mesh::P4estMeshView) From a43e82d9409d0e7054a3963de7274872752a81b5 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 22 Jul 2024 15:38:04 +0100 Subject: [PATCH 28/43] Removed double occurence of function balance!(mesh::P4estMeshView{2}). --- src/meshes/p4est_mesh_view.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 25c0ba7c68..ecbafa7cbb 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -165,13 +165,6 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end -function balance!(mesh::P4estMeshView{2}, init_fn = C_NULL) - p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) - # Due to a bug in `p4est`, the forest needs to be rebalanced twice sometimes - # See https://github.com/cburstedde/p4est/issues/112 - p4est_balance(mesh.parent.p4est, P4EST_CONNECT_FACE, init_fn) -end - function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; system = "") # Create output directory (if it does not exist) From e087369a2ececa97188be48a16a06b51718d8d50 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Wed, 24 Jul 2024 20:04:11 +0100 Subject: [PATCH 29/43] Added optin for p4est write out in coupled simulations. --- src/meshes/p4est_mesh_view.jl | 2 +- src/semidiscretization/semidiscretization_coupled.jl | 9 ++++++++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index ecbafa7cbb..b94ec82e03 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -165,7 +165,7 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end -function save_mesh_file(mesh::P4estMeshView, output_directory, timestep = 0; +function save_mesh_file(mesh::P4estMeshView, output_directory; timestep = 0, system = "") # Create output directory (if it does not exist) mkpath(output_directory) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 6b009cfad2..5c895864cc 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -510,10 +510,17 @@ end function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization) n_boundaries = 2 * ndims(semi) mesh, equations, solver, _ = mesh_equations_solver_cache(semi) + boundary_dictionary_names = [:x_neg, :x_pos, :y_neg, :y_pos] for direction in 1:n_boundaries + boundary_condition = nothing + if typeof(semi.boundary_conditions) <: Trixi.UnstructuredSortedBoundaryTypes + if boundary_dictionary_names[direction] in keys(semi.boundary_conditions.boundary_dictionary) + boundary_condition = semi.boundary_conditions.boundary_dictionary[boundary_dictionary_names[direction]] + end + else boundary_condition = semi.boundary_conditions[direction] - + end allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, solver) From 455d53ecb632fcb38da64507fbdc457851ff78c1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 25 Jul 2024 11:50:22 +0100 Subject: [PATCH 30/43] Added system to the mesh file neames for p4est mesh views. --- src/meshes/p4est_mesh_view.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index b94ec82e03..6d6ccd9fd0 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -176,8 +176,8 @@ function save_mesh_file(mesh::P4estMeshView, output_directory; timestep = 0, @sprintf("mesh_%s_%09d.h5", system, timestep)) p4est_filename = @sprintf("p4est_data_%s_%09d", system, timestep) else - filename = joinpath(output_directory, "mesh.h5") - p4est_filename = "p4est_data" + filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) + p4est_filename = @sprintf("p4est_data_%s", system) end p4est_file = joinpath(output_directory, p4est_filename) From f058b570232ee45dab91d0256a4685faa0c2ec72 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 26 Jul 2024 15:12:06 +0100 Subject: [PATCH 31/43] Applied autoformatter. --- src/semidiscretization/semidiscretization_coupled.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 5c895864cc..cf1f80dcee 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -514,13 +514,13 @@ function allocate_coupled_boundary_conditions(semi::AbstractSemidiscretization) for direction in 1:n_boundaries boundary_condition = nothing - if typeof(semi.boundary_conditions) <: Trixi.UnstructuredSortedBoundaryTypes - if boundary_dictionary_names[direction] in keys(semi.boundary_conditions.boundary_dictionary) - boundary_condition = semi.boundary_conditions.boundary_dictionary[boundary_dictionary_names[direction]] + if typeof(semi.boundary_conditions) <: Trixi.UnstructuredSortedBoundaryTypes + if boundary_dictionary_names[direction] in keys(semi.boundary_conditions.boundary_dictionary) + boundary_condition = semi.boundary_conditions.boundary_dictionary[boundary_dictionary_names[direction]] + end + else + boundary_condition = semi.boundary_conditions[direction] end - else - boundary_condition = semi.boundary_conditions[direction] - end allocate_coupled_boundary_condition(boundary_condition, direction, mesh, equations, solver) From 1cdfd88e3ed6a94eabefad3c4c4054673894e65f Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 26 Jul 2024 15:57:08 +0100 Subject: [PATCH 32/43] Added indexed size function for p4est mesh views. --- src/meshes/p4est_mesh_view.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 6d6ccd9fd0..22144b93f5 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -88,6 +88,11 @@ end @inline ncellsglobal(mesh::P4estMeshView) = Int(mesh.p4est.global_num_quadrants[]) Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) +Base.size(mesh::P4estMeshView) = size(mesh.tree_node_coordinates)[end] +function Base.size(mesh::P4estMeshView, i) + @unpack indices_min, indices_max = mesh + return indices_max[i] - indices_min[i] + 1 +end function Base.show(io::IO, mesh::P4estMeshView) print(io, "P4estMeshView{", ndims(mesh), ", ", real(mesh), "}") From 4ee4a738367ba8e3f32b24b99dd9077ac47c9c7a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Jul 2024 19:06:20 +0100 Subject: [PATCH 33/43] Changed the p4est mesh view example to two mesh views that sit side by side without coupling. --- .../elixir_advection_p4estview.jl | 29 ++++++++++++------- test/test_p4est_2d.jl | 4 +-- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl index 7187634b76..37c491f965 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -18,21 +18,26 @@ coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) trees_per_dimension = (8, 8) -# Create P4estMesh with 8 x 8 trees and 16 x 16 elements -parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3, - coordinates_min = coordinates_min, - coordinates_max = coordinates_max, - initial_refinement_level = 1, periodicity = (false, true)) -mesh = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (4, 8), - coordinates_min = coordinates_min, coordinates_max = (0.0, 1.0), - periodicity = (false, true)) +# Create P4estMesh with 8 x 8 trees and 16 x 16 elements and the mesh views +parent_mesh = P4estMesh(trees_per_dimension, polydeg=3, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + initial_refinement_level=1, periodicity=(false, true)) +mesh1 = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (4, 8), + coordinates_min = coordinates_min, coordinates_max = (0.0, 1.0), + periodicity = (false, true)) +mesh2 = P4estMeshView(parent_mesh; indices_min = (5, 1), indices_max = (8, 8), + coordinates_min = (0.0, -1.0), coordinates_max = coordinates_max, + periodicity = (false, true)) boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_convergence_test), :x_pos => BoundaryConditionDirichlet(initial_condition_convergence_test)) # A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, - solver, boundary_conditions = boundary_conditions) +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, boundary_conditions=boundary_conditions) +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver, boundary_conditions=boundary_conditions) + +# Create a semidiscretization that bundles semi1 and semi2 +semi = SemidiscretizationCoupled(semi1, semi2) ############################################################################### # ODE solvers, callbacks etc. @@ -45,7 +50,9 @@ ode = semidiscretize(semi, (0.0, 1.0)); summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback = AnalysisCallback(semi, interval = 100) +analysis_callback1 = AnalysisCallback(semi1, interval=100) +analysis_callback2 = AnalysisCallback(semi2, interval=100) +analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) # The SaveSolutionCallback allows to save the solution to a file in regular intervals save_solution = SaveSolutionCallback(interval = 100, diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index fbd7b065ac..9340de7fa6 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -127,8 +127,8 @@ end @trixi_testset "elixir_advection_p4estview.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4estview.jl"), - l2=[1.842677076881425e-5], - linf=[0.00013873258442609337]) + l2 = [1.842677076881425e-5, 1.842677076883415e-5], + linf = [0.00013873258442609337, 0.00013873258442687053]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 7ab4485a29c70c997393e1b0bfa23a4c479d676a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 29 Jul 2024 20:24:00 +0100 Subject: [PATCH 34/43] Applied autoformatter. --- .../elixir_advection_p4estview.jl | 17 ++++++++++------- test/test_p4est_2d.jl | 4 ++-- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl index 37c491f965..5e942f0ccd 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_p4estview.jl @@ -19,9 +19,10 @@ coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) trees_per_dimension = (8, 8) # Create P4estMesh with 8 x 8 trees and 16 x 16 elements and the mesh views -parent_mesh = P4estMesh(trees_per_dimension, polydeg=3, - coordinates_min=coordinates_min, coordinates_max=coordinates_max, - initial_refinement_level=1, periodicity=(false, true)) +parent_mesh = P4estMesh(trees_per_dimension, polydeg = 3, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + initial_refinement_level = 1, periodicity = (false, true)) mesh1 = P4estMeshView(parent_mesh; indices_min = (1, 1), indices_max = (4, 8), coordinates_min = coordinates_min, coordinates_max = (0.0, 1.0), periodicity = (false, true)) @@ -33,8 +34,10 @@ boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_conditio :x_pos => BoundaryConditionDirichlet(initial_condition_convergence_test)) # A semidiscretization collects data structures and functions for the spatial discretization -semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, solver, boundary_conditions=boundary_conditions) -semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, solver, boundary_conditions=boundary_conditions) +semi1 = SemidiscretizationHyperbolic(mesh1, equations, initial_condition_convergence_test, + solver, boundary_conditions = boundary_conditions) +semi2 = SemidiscretizationHyperbolic(mesh2, equations, initial_condition_convergence_test, + solver, boundary_conditions = boundary_conditions) # Create a semidiscretization that bundles semi1 and semi2 semi = SemidiscretizationCoupled(semi1, semi2) @@ -50,8 +53,8 @@ ode = semidiscretize(semi, (0.0, 1.0)); summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback1 = AnalysisCallback(semi1, interval=100) -analysis_callback2 = AnalysisCallback(semi2, interval=100) +analysis_callback1 = AnalysisCallback(semi1, interval = 100) +analysis_callback2 = AnalysisCallback(semi2, interval = 100) analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2) # The SaveSolutionCallback allows to save the solution to a file in regular intervals diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 9340de7fa6..5fe108f84f 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -127,8 +127,8 @@ end @trixi_testset "elixir_advection_p4estview.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4estview.jl"), - l2 = [1.842677076881425e-5, 1.842677076883415e-5], - linf = [0.00013873258442609337, 0.00013873258442687053]) + l2=[1.842677076881425e-5, 1.842677076883415e-5], + linf=[0.00013873258442609337, 0.00013873258442687053]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 7cca87a2f6c2a08d223cbe67573cc5e8fbeb3b62 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 1 Aug 2024 16:38:22 +0100 Subject: [PATCH 35/43] Corrected trees_per_dimension in p4est mesh views. --- 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 22144b93f5..caf3cc183a 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -29,7 +29,7 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; indices_max = size(parent), coordinates_min = nothing, coordinates_max = nothing, periodicity = (true, true)) where {NDIMS, RealT} - trees_per_dimension = indices_max .- indices_min + trees_per_dimension = indices_max .- indices_min .+ (1, 1) @assert indices_min <= indices_max @assert all(indices_min .> 0) From c57364d73c7c11f1fcf1c1f2b3cff9c48191b088 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 1 Aug 2024 16:40:37 +0100 Subject: [PATCH 36/43] Updated errors in p4est view test. --- 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 6a1a542091..fd2b6d1765 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -127,8 +127,8 @@ end @trixi_testset "elixir_advection_p4estview.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_p4estview.jl"), - l2=[1.842677076881425e-5, 1.842677076883415e-5], - linf=[0.00013873258442609337, 0.00013873258442687053]) + l2=[8.286821350926049e-6, 8.286821351032878e-6], + linf=[6.674634841197236e-5, 6.674634841874472e-5]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 2724abe044eee113b487951a40ecfc18ebca3eb1 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 15 Aug 2024 13:27:23 +0100 Subject: [PATCH 37/43] Commented constructor of P4estMeshView. --- src/meshes/p4est_mesh_view.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index caf3cc183a..3f0775d33a 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -1,6 +1,13 @@ @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, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: AbstractMesh{NDIMS} From 55b3ee745eaf46f7e64a53d5e0e128246a65d89c Mon Sep 17 00:00:00 2001 From: SimonCan Date: Fri, 23 Aug 2024 16:18:46 +0100 Subject: [PATCH 38/43] Removed unused code. --- src/meshes/p4est_mesh_view.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 3f0775d33a..e266e2116d 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -93,17 +93,6 @@ end return mesh.p4est.trees.elem_count[] end @inline ncellsglobal(mesh::P4estMeshView) = Int(mesh.p4est.global_num_quadrants[]) -Base.axes(mesh::P4estMeshView) = map(Base.OneTo, size(mesh)) -Base.axes(mesh::P4estMeshView, i) = Base.OneTo(size(mesh, i)) -Base.size(mesh::P4estMeshView) = size(mesh.tree_node_coordinates)[end] -function Base.size(mesh::P4estMeshView, i) - @unpack indices_min, indices_max = mesh - return indices_max[i] - indices_min[i] + 1 -end - -function Base.show(io::IO, mesh::P4estMeshView) - print(io, "P4estMeshView{", ndims(mesh), ", ", real(mesh), "}") -end function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMeshView) if get(io, :compact, false) From ad63556750b2493ad7e5e62ce9173fb5167506e0 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Mon, 16 Sep 2024 16:24:35 +0100 Subject: [PATCH 39/43] Applied autoformatter. --- src/solvers/dgsem_p4est/containers.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 23cdf863ff..452f913560 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -81,7 +81,8 @@ function Base.resize!(elements::P4estElementContainer, capacity) end # Create element container and initialize element data -function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT}, P4estMeshView{NDIMS, RealT}, +function init_elements(mesh::Union{P4estMesh{NDIMS, NDIMS, RealT}, + P4estMeshView{NDIMS, RealT}, T8codeMesh{NDIMS, RealT}}, equations, basis, From 89370737ab68579091f2c172d5fb54c104c140a5 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 17 Sep 2024 18:05:05 +0100 Subject: [PATCH 40/43] Changes to reflect NDIMS_AMBIENT in p4est. --- src/meshes/p4est_mesh_view.jl | 11 ++++++----- src/solvers/dgsem_p4est/containers.jl | 2 +- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index e266e2116d..cfb67eb9e2 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -8,7 +8,7 @@ A view on a p4est mesh. """ -mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, +mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <: AbstractMesh{NDIMS} # Attributes from the original P4est mesh @@ -25,13 +25,13 @@ mutable struct P4estMeshView{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2 p4est_partition_allow_for_coarsening::Bool # Attributes pertaining the views. - parent::P4estMesh{NDIMS, RealT} + parent::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT} indices_min::NTuple{NDIMS, Int} indices_max::NTuple{NDIMS, Int} end # function P4estMeshView(parent::P4estMesh{NDIMS, RealT}, view_cells::Array{Bool}) where {NDIMS, RealT} -function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; +function P4estMeshView(parent::P4estMesh{NDIMS, NDIMS, RealT}; indices_min = ntuple(_ -> 1, Val(NDIMS)), indices_max = size(parent), coordinates_min = nothing, coordinates_max = nothing, @@ -73,7 +73,7 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; structured_boundary_names!(boundary_names, trees_per_dimension, periodicity) - return P4estMeshView{NDIMS, eltype(parent.tree_node_coordinates), + return P4estMeshView{NDIMS, NDIMS, eltype(parent.tree_node_coordinates), typeof(parent.is_parallel), typeof(p4est_pw), typeof(parent.ghost), NDIMS + 2, length(parent.nodes)}(PointerWrapper(p4est), @@ -88,7 +88,8 @@ function P4estMeshView(parent::P4estMesh{NDIMS, RealT}; end @inline Base.ndims(::P4estMeshView{NDIMS}) where {NDIMS} = NDIMS -@inline Base.real(::P4estMeshView{NDIMS, RealT}) where {NDIMS, RealT} = RealT +@inline Base.real(::P4estMeshView{NDIMS, NDIMS, RealT}) where {NDIMS, RealT} = RealT +@inline ndims_ambient(::P4estMesh{NDIMS, NDIMS}) where {NDIMS} = NDIMS @inline function ntrees(mesh::P4estMeshView) return mesh.p4est.trees.elem_count[] end diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 452f913560..397ac54a38 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, From 234f11f40ac7cbc144ea217dfa29676daca1683a Mon Sep 17 00:00:00 2001 From: SimonCan Date: Thu, 19 Sep 2024 16:07:42 +0100 Subject: [PATCH 41/43] Applied autoformatter. --- 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 cfb67eb9e2..ec30529ab5 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, IsParallel, P, Ghost, NDIMSP2, +mutable struct P4estMeshView{NDIMS, NDIMS_AMBIENT, RealT <: Real, IsParallel, P, Ghost, + NDIMSP2, NNODES} <: AbstractMesh{NDIMS} # Attributes from the original P4est mesh From 81d836e9df1bb7b855c5eeec4fd6f14115894236 Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 8 Oct 2024 11:29:47 +0100 Subject: [PATCH 42/43] Applied autoformatter. --- src/meshes/p4est_mesh_view.jl | 2 +- src/solvers/dgsem_tree/dg_2d.jl | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index ec30529ab5..f0480d81ab 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -103,7 +103,7 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMeshView) setup = [ "#trees" => ntrees(mesh), "current #cells" => ncellsglobal(mesh), - "polydeg" => length(mesh.nodes) - 1, + "polydeg" => length(mesh.nodes) - 1 ] summary_box(io, "P4estMeshView{" * string(ndims(mesh)) * ", " * string(real(mesh)) * diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index c49ffb9c52..176aa2fa00 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -113,7 +113,8 @@ end # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer? function rhs!(du, u, t, - mesh::Union{TreeMesh{2}, P4estMesh{2}, 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 820cf1390ac38dd38985b1c0f17ff6d3fc0980db Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 8 Oct 2024 13:14:05 +0100 Subject: [PATCH 43/43] Removed the option of writing mesh files for different time steps. This is currently not used. --- src/meshes/p4est_mesh_view.jl | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index f0480d81ab..e6110f2036 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -173,15 +173,8 @@ function save_mesh_file(mesh::P4estMeshView, output_directory; timestep = 0, # 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_%s_%09d.h5", system, timestep)) - p4est_filename = @sprintf("p4est_data_%s_%09d", system, timestep) - else - filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) - p4est_filename = @sprintf("p4est_data_%s", system) - end + filename = joinpath(output_directory, @sprintf("mesh_%s.h5", system)) + p4est_filename = @sprintf("p4est_data_%s", system) p4est_file = joinpath(output_directory, p4est_filename)