From 5aaee8149d52afdba35241b2c6dfee3c5606dbaa Mon Sep 17 00:00:00 2001 From: SimonCan Date: Tue, 19 Mar 2024 16:09:40 +0000 Subject: [PATCH] Added parameters to p4est mesh struct. Added preliminary coupling methods for p4est meshes. --- src/meshes/p4est_mesh.jl | 14 +++- .../semidiscretization_coupled.jl | 71 +++++++++++++++++++ 2 files changed, 82 insertions(+), 3 deletions(-) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index abe9d9345b5..c44f3669c03 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -24,10 +24,14 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN current_filename::String unsaved_changes::Bool p4est_partition_allow_for_coarsening::Bool + coordinates_min::SVector{NDIMS, RealT} + coordinates_max::SVector{NDIMS, RealT} + trees_per_dimension::SVector{NDIMS, Int} function P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, boundary_names, current_filename, unsaved_changes, - p4est_partition_allow_for_coarsening) where {NDIMS} + p4est_partition_allow_for_coarsening, + coordinates_min, coordinates_max, trees_per_dimension) where {NDIMS} if NDIMS == 2 @assert p4est isa Ptr{p4est_t} elseif NDIMS == 3 @@ -57,7 +61,10 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN boundary_names, current_filename, unsaved_changes, - p4est_partition_allow_for_coarsening) + p4est_partition_allow_for_coarsening, + coordinates_min, + coordinates_max, + trees_per_dimension) # Destroy `p4est` structs when the mesh is garbage collected finalizer(destroy_mesh, mesh) @@ -215,7 +222,8 @@ function P4estMesh(trees_per_dimension; polydeg, return P4estMesh{NDIMS}(p4est, tree_node_coordinates, nodes, boundary_names, "", unsaved_changes, - p4est_partition_allow_for_coarsening) + p4est_partition_allow_for_coarsening, + coordinates_min, coordinates_max, trees_per_dimension) end # 2D version diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index dc21dbe9a1e..8ab3eba0ce7 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -467,6 +467,7 @@ end # In 2D function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2}, direction, mesh, equations, dg::DGSEM) + @autoinfiltrate if direction in (1, 2) cell_size = size(mesh, 2) else @@ -479,6 +480,48 @@ function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditi cell_size) end +# In 2D +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 + }, + direction, mesh::P4estMesh, equations, dg::DGSEM) + @autoinfiltrate + if direction in (1, 2) + cell_size = size(mesh, 2) + # Negative and positive y. + else + cell_size = size(mesh, 1) + end + + uEltype = eltype(boundary_condition) + boundary_condition.u_boundary = Array{uEltype, 3}(undef, nvariables(equations), + nnodes(dg), + cell_size) +end + +# In 2D for a p4est mesh. +function allocate_coupled_boundary_condition(boundary_condition::BoundaryConditionCoupled{2 + }, + direction, mesh::P4estMesh, equations, dg::DGSEM) + # Negative x. + if direction == 1 + cell_size = sum(mesh.tree_node_coordinates[1, 1, 1, :] .== minimum(mesh.tree_node_coordinates[1, 1, 1, :])) + # Positive x. + elseif direction == 2 + cell_size = sum(mesh.tree_node_coordinates[1, 1, 1, :] .== maximum(mesh.tree_node_coordinates[1, 1, 1, :])) + # Negative y. + elseif direction == 3 + cell_size = sum(mesh.tree_node_coordinates[2, 1, 1, :] .== minimum(mesh.tree_node_coordinates[2, 1, 1, :])) + # Positive y. + else + cell_size = sum(mesh.tree_node_coordinates[2, 1, 1, :] .== maximum(mesh.tree_node_coordinates[2, 1, 1, :])) + end + + uEltype = eltype(boundary_condition) + boundary_condition.u_boundary = Array{uEltype, 3}(undef, nvariables(equations), + nnodes(dg), + cell_size) +end + # Don't do anything for other BCs than BoundaryConditionCoupled function copy_to_coupled_boundary!(boundary_condition, u_ode, semi_coupled, semi) return nothing @@ -535,6 +578,26 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ cells = axes(mesh_other, 1) end + if mesh isa P4estMesh + linear_indices = LinearIndices((mesh.trees_per_dimension[1], mesh.trees_per_dimension[2])) + else + linear_indices = LinearIndices(size(mesh)) + end + + if mesh isa P4estMesh + if other_orientation == 1 + cells = mesh.trees_per_dimension[2] + else # other_orientation == 2 + cells = mesh.trees_per_dimension[1] + end + else + if other_orientation == 1 + cells = axes(mesh, 2) + else # other_orientation == 2 + cells = axes(mesh, 1) + end + end + # Copy solution data to the coupled boundary using "delayed indexing" with # a start value and a step size to get the correct face and orientation. node_index_range = eachnode(solver_other) @@ -543,10 +606,18 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh_other, 1)) j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2)) + if mesh isa P4estMesh + i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], mesh.trees_per_dimension[1]) + j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], mesh.trees_per_dimension[2]) + else + i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh, 1)) + j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh, 2)) + end i_cell = i_cell_start j_cell = j_cell_start + # @autoinfiltrate for cell in cells i_node = i_node_start j_node = j_node_start