diff --git a/docs/make.jl b/docs/make.jl index b138519d93f..c9885561004 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -151,7 +151,9 @@ makedocs( "Structured mesh" => joinpath("meshes", "structured_mesh.md"), "Unstructured mesh" => joinpath("meshes", "unstructured_quad_mesh.md"), "P4est-based mesh" => joinpath("meshes", "p4est_mesh.md"), - "DGMulti mesh" => joinpath("meshes", "dgmulti_mesh.md") + "DGMulti mesh" => joinpath("meshes", "dgmulti_mesh.md"), + "Mesh constructor comparison" => joinpath("meshes", + "mesh_constructor_comparison.md") ], "Time integration" => "time_integration.md", "Callbacks" => "callbacks.md", diff --git a/docs/src/meshes/mesh_constructor_comparison.md b/docs/src/meshes/mesh_constructor_comparison.md new file mode 100644 index 00000000000..2e8ba3c8517 --- /dev/null +++ b/docs/src/meshes/mesh_constructor_comparison.md @@ -0,0 +1,34 @@ +# Unified mesh constructor interface + +Trixi.jl provides several mesh types suited for different scenarios. +It may be useful to quickly swap between mesh types +without having to rewrite the mesh construction call. + +## Keyword-only interface + +All structured mesh types support a keyword-only constructor for rectangular domains +using `refinement_level`, yielding `2^refinement_level` cells per dimension. + +```@example mesh-swap +using Trixi + +for MeshType in (TreeMesh, StructuredMesh, P4estMesh, T8codeMesh) + mesh = MeshType(coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0), + refinement_level = 2) + display(mesh) +end +``` + +[`DGMultiMesh`](@ref) also supports the same keyword arguments, but requires a solver `dg` +as the first positional argument: + +```@example mesh-swap +dg = DGMulti(polydeg = 1, element_type = Quad()) +mesh = DGMultiMesh(dg; coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0), + refinement_level = 2) +``` + +Note: for `StructuredMesh` and `DGMultiMesh`, `refinement_level` directly sets +`cells_per_dimension = 2^refinement_level`. For `P4estMesh` and `T8codeMesh`, +a single tree per dimension is created and refined `refinement_level` times, +which also yields `2^refinement_level` leaf cells per dimension. diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index ae099d601ed..88c66645c4e 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -245,6 +245,56 @@ function P4estMesh(trees_per_dimension; polydeg, p4est_partition_allow_for_coarsening) end +""" + P4estMesh(; coordinates_min, coordinates_max, refinement_level, + polydeg = 1, RealT = Float64, periodicity = false, + unsaved_changes = true, + p4est_partition_allow_for_coarsening = true) + +Create a rectangular `P4estMesh` using keyword arguments only, for easy mesh-type swapping +with [`TreeMesh`](@ref), [`StructuredMesh`](@ref), and [`T8codeMesh`](@ref). + +A single tree per dimension is created and uniformly refined `refinement_level` times, +yielding `2^refinement_level` cells per dimension. + +# Arguments +- `coordinates_min`: vector or tuple of the coordinates of the corner in the negative direction of each dimension + to create a rectangular mesh. +- `coordinates_max`: vector or tuple of the coordinates of the corner in the positive direction of each dimension + to create a rectangular mesh. Must have the same length as `coordinates_min`. +- `refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. Default: `1`. +- `RealT::Type`: the type that should be used for coordinates. Default: `Float64`. +- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` + deciding for each dimension if the boundaries in this dimension are periodic. Default: `false`. +- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file. Default: `true`. +- `p4est_partition_allow_for_coarsening::Bool`: Must be `true` when using AMR to make mesh adaptivity + independent of domain partitioning. Should be `false` for static meshes + to permit more fine-grained partitioning. Default: `true`. +""" +function P4estMesh(; coordinates_min, + coordinates_max, + refinement_level, + polydeg = 1, + RealT = Float64, + periodicity = false, + unsaved_changes = true, + p4est_partition_allow_for_coarsening = true) + if length(coordinates_min) != length(coordinates_max) + throw(ArgumentError("coordinates_min and coordinates_max must have the same length")) + end + NDIMS = length(coordinates_min) + return P4estMesh(ntuple(_ -> 1, NDIMS); + polydeg = polydeg, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + RealT = RealT, + initial_refinement_level = refinement_level, + periodicity = periodicity, + unsaved_changes = unsaved_changes, + p4est_partition_allow_for_coarsening = p4est_partition_allow_for_coarsening) +end + # 2D version function structured_boundary_names!(boundary_names, trees_per_dimension::NTuple{2}, periodicity) diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 716015d746e..f566be8688c 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -145,6 +145,38 @@ function StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max; mapping_as_string = mapping_as_string) end +""" + StructuredMesh(; coordinates_min, coordinates_max, refinement_level, periodicity=false) + +Create a rectangular `StructuredMesh` using keyword arguments only, for easy mesh-type swapping +with [`TreeMesh`](@ref), [`P4estMesh`](@ref), and [`T8codeMesh`](@ref). + +The number of cells per dimension is `2^refinement_level`. + +Note that `StructuredMesh` does not support adaptive mesh refinement; +`refinement_level` only sets the initial uniform resolution. + +# Arguments +- `coordinates_min::NTuple{NDIMS, RealT}`: coordinate of the corner in the negative direction of each dimension. +- `coordinates_max::NTuple{NDIMS, RealT}`: coordinate of the corner in the positive direction of each dimension. + Must have the same length as `coordinates_min`. +- `refinement_level::Integer`: the number of cells in each dimension is set to `2^refinement_level`. +- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for + each dimension if the boundaries in this dimension are periodic. +""" +function StructuredMesh(; coordinates_min, + coordinates_max, + refinement_level, + periodicity = false) + if length(coordinates_min) != length(coordinates_max) + throw(ArgumentError("coordinates_min and coordinates_max must have the same length")) + end + NDIMS = length(coordinates_min) + cells_per_dimension = ntuple(_ -> 2^refinement_level, NDIMS) + return StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max; + periodicity = periodicity) +end + # Extract a string of the code that defines the mapping function function mapping2string(mapping, ndims, RealT = Float64) return string(code_string(mapping, ntuple(_ -> RealT, ndims))) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 9cfbe7016a4..5cfdf2df4b0 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -505,6 +505,46 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, mapping = mapping_) end +""" + T8codeMesh(; coordinates_min, coordinates_max, refinement_level, + polydeg = 1, RealT = Float64, periodicity = false) + +Create a rectangular `T8codeMesh` using keyword arguments only, for easy mesh-type swapping +with [`TreeMesh`](@ref), [`StructuredMesh`](@ref), and [`P4estMesh`](@ref). + +A single tree per dimension is created and uniformly refined `refinement_level` times, +yielding `2^refinement_level` cells per dimension. + +# Arguments +- `coordinates_min`: vector or tuple of the coordinates of the corner in the negative direction of each dimension + to create a rectangular mesh. Must have the same length as `coordinates_max`. +- `coordinates_max`: vector or tuple of the coordinates of the corner in the positive direction of each dimension + to create a rectangular mesh. Must have the same length as `coordinates_min`. +- `refinement_level::Integer`: refine the mesh uniformly to this level before the simulation starts. +- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh. Default: `1`. +- `RealT::Type`: the type that should be used for coordinates. Default: `Float64`. +- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` + deciding for each dimension if the boundaries in this dimension are periodic. Default: `false`. +""" +function T8codeMesh(; coordinates_min, + coordinates_max, + refinement_level, + polydeg = 1, + RealT = Float64, + periodicity = false) + if length(coordinates_min) != length(coordinates_max) + throw(ArgumentError("coordinates_min and coordinates_max must have the same length")) + end + NDIMS = length(coordinates_min) + return T8codeMesh(ntuple(_ -> 1, NDIMS); + polydeg = polydeg, + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + RealT = RealT, + initial_refinement_level = refinement_level, + periodicity = periodicity) +end + """ T8codeMesh(cmesh::Ptr{t8_cmesh}, mapping=nothing, polydeg=1, RealT=Float64, diff --git a/src/meshes/tree_mesh.jl b/src/meshes/tree_mesh.jl index 38132a1fb8f..a5a8813f9af 100644 --- a/src/meshes/tree_mesh.jl +++ b/src/meshes/tree_mesh.jl @@ -236,6 +236,50 @@ function TreeMesh(coordinates_min::Real, coordinates_max::Real; return TreeMesh((coordinates_min,), (coordinates_max,); kwargs...) end +""" + TreeMesh(; coordinates_min, coordinates_max, refinement_level, + n_cells_max = nothing, periodicity = false, + refinement_patches = (), coarsening_patches = (), RealT = Float64) + +Create a [`TreeMesh`](@ref) using keyword arguments only, for easy mesh-type swapping +with [`StructuredMesh`](@ref), [`P4estMesh`](@ref), and [`T8codeMesh`](@ref). + +# Arguments +- `coordinates_min`: coordinates of the low corner of the domain as a tuple, + e.g. `(-1.0, -1.0)` for 2D. +- `coordinates_max`: coordinates of the high corner of the domain as a tuple. + Must have the same length as `coordinates_min`. +- `refinement_level::Integer`: number of uniform refinements; + yields `2^refinement_level` cells per dimension. +- `n_cells_max`: initial capacity of the mesh data structures. If `nothing` + (default), the capacity is derived from `refinement_level`. The mesh grows + automatically beyond the initial capacity when AMR requires more cells. +- `periodicity`: either a `Bool` applied to all dimensions or an `NTuple{NDIMS, Bool}` + specifying periodicity per dimension. Default: `false`. +- `refinement_patches`: regions to additionally refine. Default: `()`. +- `coarsening_patches`: regions to coarsen. Default: `()`. +- `RealT`: floating-point type for coordinates. Default: `Float64`. +""" +function TreeMesh(; coordinates_min, + coordinates_max, + refinement_level, + n_cells_max = nothing, + periodicity = false, + refinement_patches = (), + coarsening_patches = (), + RealT = Float64) + if length(coordinates_min) != length(coordinates_max) + throw(ArgumentError("coordinates_min and coordinates_max must have the same length")) + end + return TreeMesh(coordinates_min, coordinates_max; + initial_refinement_level = refinement_level, + n_cells_max = n_cells_max, + periodicity = periodicity, + refinement_patches = refinement_patches, + coarsening_patches = coarsening_patches, + RealT = RealT) +end + function Base.show(io::IO, mesh::TreeMesh{NDIMS, TreeType}) where {NDIMS, TreeType} print(io, "TreeMesh{", NDIMS, ", ", TreeType, "} with length ", mesh.tree.length) return nothing diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index abf9eb1d9b0..fc34ff94e1a 100644 --- a/src/solvers/dgmulti/types.jl +++ b/src/solvers/dgmulti/types.jl @@ -298,6 +298,45 @@ function DGMultiMesh(dg::DGMulti{NDIMS}, cells_per_dimension; return DGMultiMesh(dg, GeometricTermsType(Cartesian(), dg), md, boundary_faces) end +""" + DGMultiMesh(dg::DGMulti{NDIMS}; coordinates_min, coordinates_max, + refinement_level, + is_on_boundary=nothing, + periodicity=ntuple(_->false, NDIMS)) + +Create a rectangular `DGMultiMesh` using keyword arguments only, for easy mesh-type swapping +with [`TreeMesh`](@ref), [`StructuredMesh`](@ref), [`P4estMesh`](@ref), and +[`T8codeMesh`](@ref). + +The number of cells per dimension is `2^refinement_level`. + +- `dg::DGMulti` contains information associated with the reference element (e.g., quadrature, + basis evaluation, differentiation, etc). +- `coordinates_min` is a vector or tuple of the coordinates of the corner in the negative direction of each dimension. + Must have the same length as `coordinates_max`. +- `coordinates_max` is a vector or tuple of the coordinates of the corner in the positive direction of each dimension. + Must have the same length as `coordinates_min`. +- `refinement_level` sets the number of cells per dimension to `2^refinement_level`. +- `is_on_boundary` specifies boundary using a `NamedTuple`. Default: `nothing`. +- `periodicity` is a tuple of booleans specifying if the domain is periodic `true`/`false` in the + (x,y,z) direction. Default: non-periodic in all dimensions. +""" +function DGMultiMesh(dg::DGMulti{NDIMS}; + coordinates_min, + coordinates_max, + refinement_level, + is_on_boundary = nothing, + periodicity = ntuple(_ -> false, NDIMS)) where {NDIMS} + if length(coordinates_min) != length(coordinates_max) + throw(ArgumentError("coordinates_min and coordinates_max must have the same length")) + end + cells_per_dimension = ntuple(_ -> 2^refinement_level, NDIMS) + return DGMultiMesh(dg, cells_per_dimension; + coordinates_min = coordinates_min, + coordinates_max = coordinates_max, + is_on_boundary = is_on_boundary, periodicity = periodicity) +end + """ DGMultiMesh(dg::DGMulti{NDIMS}, cells_per_dimension, mapping; is_on_boundary=nothing, diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index a59e319ad9b..298ce9129b1 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -1023,6 +1023,22 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@testset "Unified mesh constructor signatures (P4estMesh)" begin + # 2D: reference (trees_per_dimension) positional + mesh_ref = P4estMesh((4, 4); polydeg = 1, + coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0)) + + # 2D: using refinement_level (polydeg defaults to 1) + mesh_kw = P4estMesh(; coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0), + refinement_level = 2) + @test mesh_kw isa P4estMesh{2} + @test size(mesh_kw.tree_node_coordinates, ndims(mesh_kw) + 2) == 1 # 1 macro-tree + @test_throws ArgumentError P4estMesh(; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + refinement_level = 2) +end + @trixi_testset "elixir_euler_imex_warm_bubble.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_imex_warm_bubble.jl"), l2=[ diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index 8ed05ec4a22..c1fc4c26e1e 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -314,6 +314,20 @@ end @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13) @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13) end + +@testset "Unified mesh constructor signatures (T8codeMesh)" begin + using Trixi: T8codeMesh + # polydeg = 1 at default for T8codeMesh + mesh_ref = T8codeMesh((4, 4); + coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0)) + mesh_kw = T8codeMesh(; coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0), + refinement_level = 2) + @test mesh_kw isa T8codeMesh{2} + @test size(mesh_kw.tree_node_coordinates, ndims(mesh_kw) + 2) == 1 + @test_throws ArgumentError T8codeMesh(; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + refinement_level = 2) +end end # Clean up afterwards: delete Trixi.jl output directory diff --git a/test/test_unit.jl b/test/test_unit.jl index 6ee387102b5..8576a65ccd2 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -94,6 +94,17 @@ end initial_refinement_level = 2, n_cells_max = 10_000, periodicity = true) + + # Keyword-only constructor + mesh_ref = TreeMesh((-1.0, -1.0), (1.0, 1.0); + initial_refinement_level = 2, n_cells_max = 10_000) + mesh_kw = TreeMesh(; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0), + refinement_level = 2) + @test Trixi.ncells(mesh_kw) == Trixi.ncells(mesh_ref) + @test_throws ArgumentError TreeMesh(; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + refinement_level = 2) end @testset "helper functions" begin @@ -3564,6 +3575,58 @@ end end end +@testset "Unified mesh constructor signatures (StructuredMesh)" begin + # 1D: keyword interface (2^2 = 4 cells per dimension) + mesh_1d_ref = StructuredMesh((4,), (-1.0,), (1.0,)) + mesh_1d_kw = StructuredMesh(; coordinates_min = (-1.0,), coordinates_max = (1.0,), + refinement_level = 2) + @test mesh_1d_ref.cells_per_dimension == mesh_1d_kw.cells_per_dimension + + # 2D: keyword interface + mesh_2d_ref = StructuredMesh((4, 4), (-1.0, -1.0), (1.0, 1.0)) + mesh_2d_kw = StructuredMesh(; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0), + refinement_level = 2) + @test mesh_2d_ref.cells_per_dimension == mesh_2d_kw.cells_per_dimension + + # 3D: keyword interface + mesh_3d_ref = StructuredMesh((4, 4, 4), (-1.0, -1.0, -1.0), (1.0, 1.0, 1.0)) + mesh_3d_kw = StructuredMesh(; coordinates_min = (-1.0, -1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + refinement_level = 2) + @test mesh_3d_ref.cells_per_dimension == mesh_3d_kw.cells_per_dimension + @test_throws ArgumentError StructuredMesh(; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + refinement_level = 2) +end + +@testset "Unified mesh constructor signatures (DGMultiMesh)" begin + dg_1d = DGMulti(polydeg = 2, element_type = Line(), + approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_central), + volume_integral = VolumeIntegralFluxDifferencing(flux_central)) + + # 1D: keyword interface (2^2 = 4 elements) + mesh_1d_ref = DGMultiMesh(dg_1d, (4,)) + mesh_1d_kw = DGMultiMesh(dg_1d; coordinates_min = (-1.0,), coordinates_max = (1.0,), + refinement_level = 2) + @test mesh_1d_ref.md.num_elements == mesh_1d_kw.md.num_elements + + dg_2d = DGMulti(polydeg = 2, element_type = Quad(), + approximation_type = Polynomial(), + surface_integral = SurfaceIntegralWeakForm(flux_central), + volume_integral = VolumeIntegralFluxDifferencing(flux_central)) + + # 2D: keyword interface + mesh_2d_ref = DGMultiMesh(dg_2d, (4, 4)) + mesh_2d_kw = DGMultiMesh(dg_2d; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0), refinement_level = 2) + @test mesh_2d_ref.md.num_elements == mesh_2d_kw.md.num_elements + @test_throws ArgumentError DGMultiMesh(dg_2d; coordinates_min = (-1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), + refinement_level = 2) +end + @testset "TreeMesh without n_cells_max" begin for NDIMS in 1:3 coords_min = ntuple(_ -> -1.0, NDIMS)