From aa418a39879dba16934090d774bdc42250775406 Mon Sep 17 00:00:00 2001 From: Lucas C Wilcox Date: Mon, 22 Apr 2024 14:00:18 -0700 Subject: [PATCH] Move cell size out of the type domain This will allow code reuse between cells of different orders, which will be particularly useful when precompiling code. --- examples/advection/semdg_advection_2d.jl | 2 +- examples/grids/grid.jl | 2 +- examples/grids/hohqmeshimport.jl | 2 +- examples/grids/sphereshellgrid.jl | 2 +- src/cells.jl | 21 ++--- src/eye.jl | 11 ++- src/gausscells.jl | 84 ++++++++++---------- src/gridarrays.jl | 8 +- src/lobattocells.jl | 98 ++++++++++++------------ test/cells.jl | 25 +++--- test/gridarrays.jl | 3 +- test/gridnumbering.jl | 12 +-- test/grids.jl | 54 ++++++------- test/kron.jl | 28 +++---- test/mpitest_n2_commgridarrays.jl | 2 +- test/mpitest_n2_gridnumbering.jl | 4 +- test/mpitest_n3_gridarrays.jl | 2 +- test/mpitest_n3_gridnumbering.jl | 4 +- 18 files changed, 178 insertions(+), 186 deletions(-) diff --git a/examples/advection/semdg_advection_2d.jl b/examples/advection/semdg_advection_2d.jl index a09e1aa..d49df42 100644 --- a/examples/advection/semdg_advection_2d.jl +++ b/examples/advection/semdg_advection_2d.jl @@ -279,7 +279,7 @@ function run( comm = MPI.COMM_WORLD, ) rank = MPI.Comm_rank(comm) - cell = LobattoCell{Tuple{(N .+ 1)...},FT,AT}() + cell = LobattoCell{FT,AT}((N .+ 1)...) coordinates = ntuple(_ -> range(FT(0), stop = FT(2π), length = K + 1), 2) periodicity = (true, true) gm = GridManager(cell, brick(coordinates, periodicity); comm = comm, min_level = L) diff --git a/examples/grids/grid.jl b/examples/grids/grid.jl index f14b08a..561d308 100644 --- a/examples/grids/grid.jl +++ b/examples/grids/grid.jl @@ -35,7 +35,7 @@ K = (2, 3, 4) coordinates = ntuple(d -> range(start = -1.0, stop = 1.0, length = K[d] + 1), length(K)) gm = GridManager( - LobattoCell{Tuple{N...},Float64,AT}(), + LobattoCell{Float64,AT}(N...), Raven.brick(coordinates); comm = comm, min_level = 2, diff --git a/examples/grids/hohqmeshimport.jl b/examples/grids/hohqmeshimport.jl index 43083fb..fdb4050 100644 --- a/examples/grids/hohqmeshimport.jl +++ b/examples/grids/hohqmeshimport.jl @@ -60,7 +60,7 @@ coarse_grid = coarsegrid("out/IceCreamCone.inp") # coarse_grid = coarsegrid("examples/grids/Pond/Pond.inp") N = (4, 4) -gm = GridManager(LobattoCell{Tuple{N...},Float64,AT}(), coarse_grid, min_level = 1) +gm = GridManager(LobattoCell{Float64,AT}(N...), coarse_grid, min_level = 1) grid = generate(gm) diff --git a/examples/grids/sphereshellgrid.jl b/examples/grids/sphereshellgrid.jl index 7fc134d..3f96a3f 100644 --- a/examples/grids/sphereshellgrid.jl +++ b/examples/grids/sphereshellgrid.jl @@ -35,7 +35,7 @@ R = 1 coarse_grid = Raven.cubeshell2dgrid(R) -gm = GridManager(LobattoCell{Tuple{N...},Float64,Array}(), coarse_grid, min_level = 2) +gm = GridManager(LobattoCell{Float64,Array}(N...), coarse_grid, min_level = 2) indicator = rand((Raven.AdaptNone, Raven.AdaptRefine), length(gm)) adapt!(gm, indicator) diff --git a/src/cells.jl b/src/cells.jl index 9d1a8fd..38a5142 100644 --- a/src/cells.jl +++ b/src/cells.jl @@ -1,18 +1,13 @@ -abstract type AbstractCell{S<:Tuple,T,A<:AbstractArray,N} end +abstract type AbstractCell{T,A<:AbstractArray,N} end -floattype(::Type{<:AbstractCell{S,T}}) where {S,T} = T -arraytype(::Type{<:AbstractCell{S,T,A}}) where {S,T,A} = A -Base.ndims(::Type{<:AbstractCell{S,T,A,N}}) where {S,T,A,N} = N -Base.size(::Type{<:AbstractCell{S,T,A}}) where {S,T,A} = size_to_tuple(S) -Base.size(::Type{<:AbstractCell{S,T,A}}, i::Integer) where {S,T,A} = size_to_tuple(S)[i] -Base.length(::Type{<:AbstractCell{S,T,A}}) where {S,T,A} = tuple_prod(S) -Base.strides(::Type{<:AbstractCell{S,T,A}}) where {S,T,A} = - Base.size_to_strides(1, size_to_tuple(S)...) +floattype(::Type{<:AbstractCell{T}}) where {T} = T +arraytype(::Type{<:AbstractCell{T,A}}) where {T,A} = A +Base.ndims(::Type{<:AbstractCell{T,A,N}}) where {T,A,N} = N floattype(cell::AbstractCell) = floattype(typeof(cell)) arraytype(cell::AbstractCell) = arraytype(typeof(cell)) Base.ndims(cell::AbstractCell) = Base.ndims(typeof(cell)) -Base.size(cell::AbstractCell) = Base.size(typeof(cell)) -Base.size(cell::AbstractCell, i::Integer) = Base.size(typeof(cell), i) -Base.length(cell::AbstractCell) = Base.length(typeof(cell)) -Base.strides(cell::AbstractCell) = Base.strides(typeof(cell)) + +Base.size(cell::AbstractCell, i::Integer) = size(cell)[i] +Base.length(cell::AbstractCell) = prod(size(cell)) +Base.strides(cell::AbstractCell) = Base.size_to_strides(1, size(cell)...) diff --git a/src/eye.jl b/src/eye.jl index 0ab5cb8..e39e510 100644 --- a/src/eye.jl +++ b/src/eye.jl @@ -1,7 +1,10 @@ -struct Eye{T,N} <: AbstractArray{T,2} end +struct Eye{T} <: AbstractArray{T,2} + N::Int +end -Base.size(::Eye{T,N}) where {T,N} = (N, N) -Base.IndexStyle(::Eye{T,N}) where {T,N} = IndexCartesian() -function Base.getindex(::Eye{T,N}, I::Vararg{Int,2}) where {T,N} +Base.size(eye::Eye{T}) where {T} = (eye.N, eye.N) +Base.IndexStyle(::Eye) = IndexCartesian() +@inline function Base.getindex(eye::Eye{T}, I::Vararg{Int,2}) where {T} + @boundscheck checkbounds(eye, I...) return (I[1] == I[2]) ? one(T) : zero(T) end diff --git a/src/gausscells.jl b/src/gausscells.jl index bd241cc..7977df8 100644 --- a/src/gausscells.jl +++ b/src/gausscells.jl @@ -46,7 +46,8 @@ function gaussoperators_1d(::Type{T}, M) where {T} ) end -struct GaussCell{S,T,A,N,O,P,D,WD,SWD,M,FM,E,H,TG,TL,TB} <: AbstractCell{S,T,A,N} +struct GaussCell{T,A,N,S,O,P,D,WD,SWD,M,FM,E,H,TG,TL,TB} <: AbstractCell{T,A,N} + size::S points_1d::O weights_1d::O points::P @@ -62,18 +63,18 @@ struct GaussCell{S,T,A,N,O,P,D,WD,SWD,M,FM,E,H,TG,TL,TB} <: AbstractCell{S,T,A,N toboundary::TB end -function Base.show(io::IO, ::GaussCell{S,T,A}) where {S,T,A} +function Base.show(io::IO, ::GaussCell{T,A,N}) where {T,A,N} print(io, "GaussCell{") - Base.show(io, S) - print(io, ", ") Base.show(io, T) print(io, ", ") Base.show(io, A) + print(io, ", ") + Base.show(io, N) print(io, "}") end -function GaussCell{Tuple{S1},T,A}() where {S1,T,A} - o = adapt(A, gaussoperators_1d(T, S1)) +function GaussCell{T,A}(m) where {T,A} + o = adapt(A, gaussoperators_1d(T, m)) points_1d = (o.points,) weights_1d = (o.weights,) @@ -90,6 +91,7 @@ function GaussCell{Tuple{S1},T,A}() where {S1,T,A} toboundary = Kron((o.toboundary,)) args = ( + (m,), points_1d, weights_1d, points, @@ -104,24 +106,23 @@ function GaussCell{Tuple{S1},T,A}() where {S1,T,A} tolobatto, toboundary, ) - GaussCell{Tuple{S1},T,A,1,typeof.(args[2:end])...}(args...) + GaussCell{T,A,1,typeof(args[1]),typeof.(args[3:end])...}(args...) end -function GaussCell{Tuple{S1,S2},T,A}() where {S1,S2,T,A} - o = adapt(A, (gaussoperators_1d(T, S1), gaussoperators_1d(T, S2))) +function GaussCell{T,A}(m1, m2) where {T,A} + o = adapt(A, (gaussoperators_1d(T, m1), gaussoperators_1d(T, m2))) - points_1d = (reshape(o[1].points, (S1, 1)), reshape(o[2].points, (1, S2))) - weights_1d = (reshape(o[1].weights, (S1, 1)), reshape(o[2].weights, (1, S2))) + points_1d = (reshape(o[1].points, (m1, 1)), reshape(o[2].points, (1, m2))) + weights_1d = (reshape(o[1].weights, (m1, 1)), reshape(o[2].weights, (1, m2))) points = vec(SVector.(points_1d...)) - derivatives = - (Kron((Eye{T,S2}(), o[1].derivative)), Kron((o[2].derivative, Eye{T,S1}()))) + derivatives = (Kron((Eye{T}(m2), o[1].derivative)), Kron((o[2].derivative, Eye{T}(m1)))) weightedderivatives = ( - Kron((Eye{T,S2}(), o[1].weightedderivative)), - Kron((o[2].weightedderivative, Eye{T,S1}())), + Kron((Eye{T}(m2), o[1].weightedderivative)), + Kron((o[2].weightedderivative, Eye{T}(m1))), ) skewweightedderivatives = ( - Kron((Eye{T,S2}(), o[1].skewweightedderivative)), - Kron((o[2].skewweightedderivative, Eye{T,S1}())), + Kron((Eye{T}(m2), o[1].skewweightedderivative)), + Kron((o[2].skewweightedderivative, Eye{T}(m1))), ) mass = Diagonal(vec(.*(weights_1d...))) ω1, ω2 = weights_1d @@ -134,6 +135,7 @@ function GaussCell{Tuple{S1,S2},T,A}() where {S1,S2,T,A} toboundary = Kron((o[2].toboundary, o[1].toboundary)) args = ( + (m1, m2), points_1d, weights_1d, points, @@ -148,40 +150,40 @@ function GaussCell{Tuple{S1,S2},T,A}() where {S1,S2,T,A} tolobatto, toboundary, ) - GaussCell{Tuple{S1,S2},T,A,2,typeof.(args[2:end])...}(args...) + GaussCell{T,A,2,typeof(args[1]),typeof.(args[3:end])...}(args...) end -function GaussCell{Tuple{S1,S2,S3},T,A}() where {S1,S2,S3,T,A} +function GaussCell{T,A}(m1, m2, m3) where {T,A} o = adapt( A, - (gaussoperators_1d(T, S1), gaussoperators_1d(T, S2), gaussoperators_1d(T, S3)), + (gaussoperators_1d(T, m1), gaussoperators_1d(T, m2), gaussoperators_1d(T, m3)), ) points_1d = ( - reshape(o[1].points, (S1, 1, 1)), - reshape(o[2].points, (1, S2, 1)), - reshape(o[3].points, (1, 1, S3)), + reshape(o[1].points, (m1, 1, 1)), + reshape(o[2].points, (1, m2, 1)), + reshape(o[3].points, (1, 1, m3)), ) weights_1d = ( - reshape(o[1].weights, (S1, 1, 1)), - reshape(o[2].weights, (1, S2, 1)), - reshape(o[3].weights, (1, 1, S3)), + reshape(o[1].weights, (m1, 1, 1)), + reshape(o[2].weights, (1, m2, 1)), + reshape(o[3].weights, (1, 1, m3)), ) points = vec(SVector.(points_1d...)) derivatives = ( - Kron((Eye{T,S3}(), Eye{T,S2}(), o[1].derivative)), - Kron((Eye{T,S3}(), o[2].derivative, Eye{T,S1}())), - Kron((o[3].derivative, Eye{T,S2}(), Eye{T,S1}())), + Kron((Eye{T}(m3), Eye{T}(m2), o[1].derivative)), + Kron((Eye{T}(m3), o[2].derivative, Eye{T}(m1))), + Kron((o[3].derivative, Eye{T}(m2), Eye{T}(m1))), ) weightedderivatives = ( - Kron((Eye{T,S3}(), Eye{T,S2}(), o[1].weightedderivative)), - Kron((Eye{T,S3}(), o[2].weightedderivative, Eye{T,S1}())), - Kron((o[3].weightedderivative, Eye{T,S2}(), Eye{T,S1}())), + Kron((Eye{T}(m3), Eye{T}(m2), o[1].weightedderivative)), + Kron((Eye{T}(m3), o[2].weightedderivative, Eye{T}(m1))), + Kron((o[3].weightedderivative, Eye{T}(m2), Eye{T}(m1))), ) skewweightedderivatives = ( - Kron((Eye{T,S3}(), Eye{T,S2}(), o[1].skewweightedderivative)), - Kron((Eye{T,S3}(), o[2].skewweightedderivative, Eye{T,S1}())), - Kron((o[3].skewweightedderivative, Eye{T,S2}(), Eye{T,S1}())), + Kron((Eye{T}(m3), Eye{T}(m2), o[1].skewweightedderivative)), + Kron((Eye{T}(m3), o[2].skewweightedderivative, Eye{T}(m1))), + Kron((o[3].skewweightedderivative, Eye{T}(m2), Eye{T}(m1))), ) mass = Diagonal(vec(.*(weights_1d...))) ω1, ω2, ω3 = weights_1d @@ -200,6 +202,7 @@ function GaussCell{Tuple{S1,S2,S3},T,A}() where {S1,S2,S3,T,A} toboundary = Kron((o[3].toboundary, o[2].toboundary, o[1].toboundary)) args = ( + (m1, m2, m3), points_1d, weights_1d, points, @@ -214,24 +217,25 @@ function GaussCell{Tuple{S1,S2,S3},T,A}() where {S1,S2,S3,T,A} tolobatto, toboundary, ) - GaussCell{Tuple{S1,S2,S3},T,A,3,typeof.(args[2:end])...}(args...) + GaussCell{T,A,3,typeof(args[1]),typeof.(args[3:end])...}(args...) end -GaussCell{S,T}() where {S,T} = GaussCell{S,T,Array}() -GaussCell{S}() where {S} = GaussCell{S,Float64}() +GaussCell{T}(args...) where {T} = GaussCell{T,Array}(args...) +GaussCell(args...) = GaussCell{Float64}(args...) -function Adapt.adapt_structure(to, cell::GaussCell{S,T,A,N}) where {S,T,A,N} +function Adapt.adapt_structure(to, cell::GaussCell{T,A,N}) where {T,A,N} names = fieldnames(GaussCell) args = ntuple(j -> adapt(to, getfield(cell, names[j])), length(names)) B = arraytype(to) - GaussCell{S,T,B,N,typeof.(args[2:end])...}(args...) + GaussCell{T,B,N,typeof(args[1]),typeof.(args[3:end])...}(args...) end const GaussLine{T,A} = GaussCell{Tuple{B},T,A} where {B} const GaussQuad{T,A} = GaussCell{Tuple{B,C},T,A} where {B,C} const GaussHex{T,A} = GaussCell{Tuple{B,C,D},T,A} where {B,C,D} +Base.size(cell::GaussCell) = cell.size points_1d(cell::GaussCell) = cell.points_1d weights_1d(cell::GaussCell) = cell.weights_1d points(cell::GaussCell) = cell.points diff --git a/src/gridarrays.jl b/src/gridarrays.jl index 4398e11..d1d10a9 100644 --- a/src/gridarrays.jl +++ b/src/gridarrays.jl @@ -89,7 +89,7 @@ end Create an array containing elements of type `T` for each point in the grid (including the ghost cells). The dimensions of the array is -`(size(celltype(grid))..., length(grid))` as the ghost cells are hidden by +`(size(referencecell(grid))..., length(grid))` as the ghost cells are hidden by default. The type `T` is assumed to be able to be interpreted into an `NTuple{M,L}`. @@ -141,9 +141,9 @@ cells. The one after is associated with the number of cells. """ function GridArray{T}(::UndefInitializer, grid::Grid) where {T} A = arraytype(grid) - dims = (size(celltype(grid))..., Int(numcells(grid, Val(false)))) - dimswithghosts = (size(celltype(grid))..., Int(numcells(grid, Val(true)))) - F = ndims(celltype(grid)) + 1 + dims = (size(referencecell(grid))..., Int(numcells(grid, Val(false)))) + dimswithghosts = (size(referencecell(grid))..., Int(numcells(grid, Val(true)))) + F = ndims(referencecell(grid)) + 1 return GridArray{T}(undef, A, dims, dimswithghosts, comm(grid), false, F) end diff --git a/src/lobattocells.jl b/src/lobattocells.jl index d79d92a..34118d0 100644 --- a/src/lobattocells.jl +++ b/src/lobattocells.jl @@ -46,7 +46,8 @@ function lobattooperators_1d(::Type{T}, M) where {T} ) end -struct LobattoCell{S,T,A,N,O,P,D,WD,SWD,M,FM,E,H,TG,TL,TB} <: AbstractCell{S,T,A,N} +struct LobattoCell{T,A,N,S,O,P,D,WD,SWD,M,FM,E,H,TG,TL,TB} <: AbstractCell{T,A,N} + size::S points_1d::O weights_1d::O points::P @@ -62,18 +63,18 @@ struct LobattoCell{S,T,A,N,O,P,D,WD,SWD,M,FM,E,H,TG,TL,TB} <: AbstractCell{S,T,A toboundary::TB end -function Base.show(io::IO, ::LobattoCell{S,T,A}) where {S,T,A} +function Base.show(io::IO, ::LobattoCell{T,A,N}) where {T,A,N} print(io, "LobattoCell{") - Base.show(io, S) - print(io, ", ") Base.show(io, T) print(io, ", ") Base.show(io, A) + print(io, ", ") + Base.show(io, N) print(io, "}") end -function LobattoCell{Tuple{S1},T,A}() where {S1,T,A} - o = adapt(A, lobattooperators_1d(T, S1)) +function LobattoCell{T,A}(m) where {T,A} + o = adapt(A, lobattooperators_1d(T, m)) points_1d = (o.points,) weights_1d = (o.weights,) @@ -90,6 +91,7 @@ function LobattoCell{Tuple{S1},T,A}() where {S1,T,A} toboundary = Kron((o.toboundary,)) args = ( + (m,), points_1d, weights_1d, points, @@ -104,24 +106,23 @@ function LobattoCell{Tuple{S1},T,A}() where {S1,T,A} tolobatto, toboundary, ) - LobattoCell{Tuple{S1},T,A,1,typeof.(args[2:end])...}(args...) + LobattoCell{T,A,1,typeof(args[1]),typeof.(args[3:end])...}(args...) end -function LobattoCell{Tuple{S1,S2},T,A}() where {S1,S2,T,A} - o = adapt(A, (lobattooperators_1d(T, S1), lobattooperators_1d(T, S2))) +function LobattoCell{T,A}(m1, m2) where {T,A} + o = adapt(A, (lobattooperators_1d(T, m1), lobattooperators_1d(T, m2))) - points_1d = (reshape(o[1].points, (S1, 1)), reshape(o[2].points, (1, S2))) - weights_1d = (reshape(o[1].weights, (S1, 1)), reshape(o[2].weights, (1, S2))) + points_1d = (reshape(o[1].points, (m1, 1)), reshape(o[2].points, (1, m2))) + weights_1d = (reshape(o[1].weights, (m1, 1)), reshape(o[2].weights, (1, m2))) points = vec(SVector.(points_1d...)) - derivatives = - (Kron((Eye{T,S2}(), o[1].derivative)), Kron((o[2].derivative, Eye{T,S1}()))) + derivatives = (Kron((Eye{T}(m2), o[1].derivative)), Kron((o[2].derivative, Eye{T}(m1)))) weightedderivatives = ( - Kron((Eye{T,S2}(), o[1].weightedderivative)), - Kron((o[2].weightedderivative, Eye{T,S1}())), + Kron((Eye{T}(m2), o[1].weightedderivative)), + Kron((o[2].weightedderivative, Eye{T}(m1))), ) skewweightedderivatives = ( - Kron((Eye{T,S2}(), o[1].skewweightedderivative)), - Kron((o[2].skewweightedderivative, Eye{T,S1}())), + Kron((Eye{T}(m2), o[1].skewweightedderivative)), + Kron((o[2].skewweightedderivative, Eye{T}(m1))), ) mass = Diagonal(vec(.*(weights_1d...))) ω1, ω2 = weights_1d @@ -135,6 +136,7 @@ function LobattoCell{Tuple{S1,S2},T,A}() where {S1,S2,T,A} args = ( + (m1, m2), points_1d, weights_1d, points, @@ -149,44 +151,44 @@ function LobattoCell{Tuple{S1,S2},T,A}() where {S1,S2,T,A} tolobatto, toboundary, ) - LobattoCell{Tuple{S1,S2},T,A,2,typeof.(args[2:end])...}(args...) + LobattoCell{T,A,2,typeof(args[1]),typeof.(args[3:end])...}(args...) end -function LobattoCell{Tuple{S1,S2,S3},T,A}() where {S1,S2,S3,T,A} +function LobattoCell{T,A}(m1, m2, m3) where {T,A} o = adapt( A, ( - lobattooperators_1d(T, S1), - lobattooperators_1d(T, S2), - lobattooperators_1d(T, S3), + lobattooperators_1d(T, m1), + lobattooperators_1d(T, m2), + lobattooperators_1d(T, m3), ), ) points_1d = ( - reshape(o[1].points, (S1, 1, 1)), - reshape(o[2].points, (1, S2, 1)), - reshape(o[3].points, (1, 1, S3)), + reshape(o[1].points, (m1, 1, 1)), + reshape(o[2].points, (1, m2, 1)), + reshape(o[3].points, (1, 1, m3)), ) weights_1d = ( - reshape(o[1].weights, (S1, 1, 1)), - reshape(o[2].weights, (1, S2, 1)), - reshape(o[3].weights, (1, 1, S3)), + reshape(o[1].weights, (m1, 1, 1)), + reshape(o[2].weights, (1, m2, 1)), + reshape(o[3].weights, (1, 1, m3)), ) points = vec(SVector.(points_1d...)) derivatives = ( - Kron((Eye{T,S3}(), Eye{T,S2}(), o[1].derivative)), - Kron((Eye{T,S3}(), o[2].derivative, Eye{T,S1}())), - Kron((o[3].derivative, Eye{T,S2}(), Eye{T,S1}())), + Kron((Eye{T}(m3), Eye{T}(m2), o[1].derivative)), + Kron((Eye{T}(m3), o[2].derivative, Eye{T}(m1))), + Kron((o[3].derivative, Eye{T}(m2), Eye{T}(m1))), ) weightedderivatives = ( - Kron((Eye{T,S3}(), Eye{T,S2}(), o[1].weightedderivative)), - Kron((Eye{T,S3}(), o[2].weightedderivative, Eye{T,S1}())), - Kron((o[3].weightedderivative, Eye{T,S2}(), Eye{T,S1}())), + Kron((Eye{T}(m3), Eye{T}(m2), o[1].weightedderivative)), + Kron((Eye{T}(m3), o[2].weightedderivative, Eye{T}(m1))), + Kron((o[3].weightedderivative, Eye{T}(m2), Eye{T}(m1))), ) skewweightedderivatives = ( - Kron((Eye{T,S3}(), Eye{T,S2}(), o[1].skewweightedderivative)), - Kron((Eye{T,S3}(), o[2].skewweightedderivative, Eye{T,S1}())), - Kron((o[3].skewweightedderivative, Eye{T,S2}(), Eye{T,S1}())), + Kron((Eye{T}(m3), Eye{T}(m2), o[1].skewweightedderivative)), + Kron((Eye{T}(m3), o[2].skewweightedderivative, Eye{T}(m1))), + Kron((o[3].skewweightedderivative, Eye{T}(m2), Eye{T}(m1))), ) mass = Diagonal(vec(.*(weights_1d...))) ω1, ω2, ω3 = weights_1d @@ -205,6 +207,7 @@ function LobattoCell{Tuple{S1,S2,S3},T,A}() where {S1,S2,S3,T,A} toboundary = Kron((o[3].toboundary, o[2].toboundary, o[1].toboundary)) args = ( + (m1, m2, m3), points_1d, weights_1d, points, @@ -219,24 +222,25 @@ function LobattoCell{Tuple{S1,S2,S3},T,A}() where {S1,S2,S3,T,A} tolobatto, toboundary, ) - LobattoCell{Tuple{S1,S2,S3},T,A,3,typeof.(args[2:end])...}(args...) + LobattoCell{T,A,3,typeof(args[1]),typeof.(args[3:end])...}(args...) end -LobattoCell{S,T}() where {S,T} = LobattoCell{S,T,Array}() -LobattoCell{S}() where {S} = LobattoCell{S,Float64}() +LobattoCell{T}(args...) where {T} = LobattoCell{T,Array}(args...) +LobattoCell(args...) = LobattoCell{Float64}(args...) -function Adapt.adapt_structure(to, cell::LobattoCell{S,T,A,N}) where {S,T,A,N} +function Adapt.adapt_structure(to, cell::LobattoCell{T,A,N}) where {T,A,N} names = fieldnames(LobattoCell) args = ntuple(j -> adapt(to, getfield(cell, names[j])), length(names)) B = arraytype(to) - LobattoCell{S,T,B,N,typeof.(args[2:end])...}(args...) + LobattoCell{T,B,N,typeof(args[1]),typeof.(args[3:end])...}(args...) end -const LobattoLine{T,A} = LobattoCell{Tuple{B},T,A} where {B} -const LobattoQuad{T,A} = LobattoCell{Tuple{B,C},T,A} where {B,C} -const LobattoHex{T,A} = LobattoCell{Tuple{B,C,D},T,A} where {B,C,D} +const LobattoLine{T,A} = LobattoCell{T,A,1} +const LobattoQuad{T,A} = LobattoCell{T,A,2} +const LobattoHex{T,A} = LobattoCell{T,A,3} +Base.size(cell::LobattoCell) = cell.size points_1d(cell::LobattoCell) = cell.points_1d weights_1d(cell::LobattoCell) = cell.weights_1d points(cell::LobattoCell) = cell.points @@ -424,7 +428,7 @@ function materializepoints( FT = floattype(referencecell) AT = arraytype(referencecell) N = (interpolation_degree + 1, interpolation_degree + 1) - interp_r = vec.(points_1d(LobattoCell{Tuple{N...},FT,AT}())) + interp_r = vec.(points_1d(LobattoCell{FT,AT}(N...))) Q = max(512 ÷ prod(length.(r)), 1) @@ -1428,7 +1432,7 @@ function materializepoints( FT = floattype(referencecell) AT = arraytype(referencecell) N = Tuple((interpolation_degree + 1) * ones(Int64, 3)) - interp_r = vec.(points_1d(LobattoCell{Tuple{N...},FT,AT}())) + interp_r = vec.(points_1d(LobattoCell{FT,AT}(N...))) IntType = typeof(length(r)) num_local = IntType(P4estTypes.lengthoflocalquadrants(forest)) diff --git a/test/cells.jl b/test/cells.jl index fb1a6fa..b75fd6a 100644 --- a/test/cells.jl +++ b/test/cells.jl @@ -1,12 +1,10 @@ using Raven.StaticArrays: size_to_tuple function cells_testsuite(AT, FT) - cell = LobattoCell{Tuple{3,3},FT,AT}() + cell = LobattoCell{FT,AT}(3, 3) @test floattype(typeof(cell)) == FT @test arraytype(typeof(cell)) <: AT @test ndims(typeof(cell)) == 2 - @test size(typeof(cell)) == (3, 3) - @test length(typeof(cell)) == 9 @test floattype(cell) == FT @test arraytype(cell) <: AT @test ndims(cell) == 2 @@ -43,15 +41,15 @@ function cells_testsuite(AT, FT) @test Array(h1d[2][1]) == I1 @test Array(h1d[2][2]) == I2 end - @test adapt(Array, cell) isa LobattoCell{S,FT,Array} where {S} + @test adapt(Array, cell) isa LobattoCell{FT,Array} - S = Tuple{3,4,2} - cell = LobattoCell{S,FT,AT}() + S = (3, 4, 2) + cell = LobattoCell{FT,AT}(S...) @test floattype(cell) == FT @test arraytype(cell) <: AT @test Base.ndims(cell) == 3 - @test size(cell) == size_to_tuple(S) - @test length(cell) == prod(size_to_tuple(S)) + @test size(cell) == S + @test length(cell) == prod(S) @test sum(mass(cell)) .≈ 8 @test mass(cell) isa Diagonal @test sum(facemass(cell)) .≈ 24 @@ -59,12 +57,9 @@ function cells_testsuite(AT, FT) @test faceoffsets(cell) == (0, 8, 16, 22, 28, 40, 52) @test strides(cell) == (1, 3, 12) D = derivatives(cell) - @test Array(D[1] * points(cell)) ≈ - fill(SVector(one(FT), zero(FT), zero(FT)), prod(size_to_tuple(S))) - @test Array(D[2] * points(cell)) ≈ - fill(SVector(zero(FT), one(FT), zero(FT)), prod(size_to_tuple(S))) - @test Array(D[3] * points(cell)) ≈ - fill(SVector(zero(FT), zero(FT), one(FT)), prod(size_to_tuple(S))) + @test Array(D[1] * points(cell)) ≈ fill(SVector(one(FT), zero(FT), zero(FT)), prod(S)) + @test Array(D[2] * points(cell)) ≈ fill(SVector(zero(FT), one(FT), zero(FT)), prod(S)) + @test Array(D[3] * points(cell)) ≈ fill(SVector(zero(FT), zero(FT), one(FT)), prod(S)) D1d = derivatives_1d(cell) @test length(D1d) == ndims(cell) @test all( @@ -86,7 +81,7 @@ function cells_testsuite(AT, FT) end end - cell = LobattoCell{Tuple{5},FT,AT}() + cell = LobattoCell{FT,AT}(5) @test floattype(cell) == FT @test arraytype(cell) <: AT @test Base.ndims(cell) == 1 diff --git a/test/gridarrays.jl b/test/gridarrays.jl index ca05381..6f20a3b 100644 --- a/test/gridarrays.jl +++ b/test/gridarrays.jl @@ -3,8 +3,7 @@ function gridarrays_testsuite(AT, FT) N = (3, 2) K = (2, 3) L = 1 - gm = - GridManager(LobattoCell{Tuple{N...},FT,AT}(), Raven.brick(FT, K); min_level = L) + gm = GridManager(LobattoCell{FT,AT}(N...), Raven.brick(FT, K); min_level = L) grid = generate(gm) x = points(grid) diff --git a/test/gridnumbering.jl b/test/gridnumbering.jl index ddf9d73..3013829 100644 --- a/test/gridnumbering.jl +++ b/test/gridnumbering.jl @@ -115,7 +115,7 @@ @test dtoc_degree_3_local == dtoc_degree_3_global @test dtoc_degree_3_local == P4estTypes.unsafe_element_nodes(nodes) .+ 0x1 - cell_degree_3 = LobattoCell{Tuple{4,4},Float64,Array}() + cell_degree_3 = LobattoCell{Float64,Array}(4, 4) dtoc_degree_3 = Raven.materializedtoc(cell_degree_3, dtoc_degree_3_local, dtoc_degree_3_global) @test isisomorphic(dtoc_degree_3, dtoc_degree_3_global) @@ -237,7 +237,7 @@ @test dtoc_degree_3_local == dtoc_degree_3_global @test dtoc_degree_3_local == P4estTypes.unsafe_element_nodes(nodes) .+ 0x1 - cell_degree_3 = LobattoCell{Tuple{4,4,4},Float64,Array}() + cell_degree_3 = LobattoCell{Float64,Array}(4, 4, 4) dtoc_degree_3 = Raven.materializedtoc(cell_degree_3, dtoc_degree_3_local, dtoc_degree_3_global) @test isisomorphic(dtoc_degree_3, dtoc_degree_3_global) @@ -280,7 +280,7 @@ AT = Array N = 4 - cell = LobattoCell{Tuple{N,N},FT,AT}() + cell = LobattoCell{FT,AT}(N, N) # 4--------5--------6 # | | | @@ -355,7 +355,7 @@ end let - cell = LobattoCell{Tuple{2,2,2}}() + cell = LobattoCell(2, 2, 2) cg = brick((1, 1, 1), (true, true, true)) gm = GridManager(cell, cg) grid = generate(gm) @@ -365,7 +365,7 @@ end let - cell = LobattoCell{Tuple{2,2}}() + cell = LobattoCell(2, 2) cg = brick((1, 1), (true, true)) gm = GridManager(cell, cg) grid = generate(gm) @@ -378,7 +378,7 @@ AT = Array N = 5 - cell = LobattoCell{Tuple{N,N,N}}() + cell = LobattoCell(N, N, N) # 10-------11-------12 # /| /| /| diff --git a/test/grids.jl b/test/grids.jl index d70fcf6..abf0f40 100644 --- a/test/grids.jl +++ b/test/grids.jl @@ -26,11 +26,7 @@ function grids_testsuite(AT, FT) coordinates = ntuple(d -> range(-one(FT), stop = one(FT), length = K[d] + 1), length(K)) - gm = GridManager( - LobattoCell{Tuple{N...},FT,AT}(), - Raven.brick(coordinates); - min_level = 2, - ) + gm = GridManager(LobattoCell{FT,AT}(N...), Raven.brick(coordinates); min_level = 2) indicator = rand(rng, (Raven.AdaptNone, Raven.AdaptRefine), length(gm)) @@ -64,7 +60,7 @@ function grids_testsuite(AT, FT) coarse_grid = Raven.cubeshellgrid(R, r) - gm = GridManager(LobattoCell{Tuple{N...},FT,AT}(), coarse_grid, min_level = 2) + gm = GridManager(LobattoCell{FT,AT}(N...), coarse_grid, min_level = 2) indicator = rand(rng, (Raven.AdaptNone, Raven.AdaptRefine), length(gm)) adapt!(gm, indicator) @@ -93,7 +89,7 @@ function grids_testsuite(AT, FT) coarse_grid = Raven.cubeshell2dgrid(R) - gm = GridManager(LobattoCell{Tuple{N...},FT,AT}(), coarse_grid, min_level = 2) + gm = GridManager(LobattoCell{FT,AT}(N...), coarse_grid, min_level = 2) indicator = rand(rng, (Raven.AdaptNone, Raven.AdaptRefine), length(gm)) adapt!(gm, indicator) @@ -121,7 +117,7 @@ function grids_testsuite(AT, FT) K = (2, 1) coordinates = ntuple(d -> range(-one(FT), stop = one(FT), length = K[d] + 1), length(K)) - cell = LobattoCell{Tuple{N...},FT,AT}() + cell = LobattoCell{FT,AT}(N...) gm = GridManager(cell, brick(coordinates, (true, true))) grid = generate(gm) @test all(boundarycodes(grid) .== 0) @@ -133,11 +129,7 @@ function grids_testsuite(AT, FT) coordinates = ntuple(d -> range(-one(FT), stop = one(FT), length = K[d] + 1), length(K)) - gm = GridManager( - LobattoCell{Tuple{N...},FT,AT}(), - Raven.brick(coordinates); - min_level = 1, - ) + gm = GridManager(LobattoCell{FT,AT}(N...), Raven.brick(coordinates); min_level = 1) indicator = rand(rng, (Raven.AdaptNone, Raven.AdaptRefine), length(gm)) @@ -164,7 +156,7 @@ function grids_testsuite(AT, FT) let - cell = LobattoCell{Tuple{4,4},FT,AT}() + cell = LobattoCell{FT,AT}(4, 4) vertices = [ SVector{2,FT}(0, 0), # 1 @@ -197,7 +189,7 @@ function grids_testsuite(AT, FT) end let - cell = LobattoCell{Tuple{3,3,3},FT,AT}() + cell = LobattoCell{FT,AT}(3, 3, 3) vertices = [ SVector{3,FT}(0, 0, 0), # 1 @@ -261,7 +253,7 @@ function grids_testsuite(AT, FT) level = 2 gm = GridManager( - LobattoCell{Tuple{L,M},FT,AT}(), + LobattoCell{FT,AT}(L, M), Raven.brick(coordinates); min_level = level, ) @@ -334,7 +326,7 @@ function grids_testsuite(AT, FT) # Polynomial orders need to be the same to match up faces. # We currently do not support different polynomial orders for # joining faces. - LobattoCell{Tuple{L,L},FT,AT}(), + LobattoCell{FT,AT}(L, L), Raven.cubeshell2dgrid(R); min_level = level, ) @@ -371,7 +363,7 @@ function grids_testsuite(AT, FT) level = 0 gm = GridManager( - LobattoCell{Tuple{L,M},FT,AT}(), + LobattoCell{FT,AT}(L, M), Raven.brick(coordinates); min_level = level, ) @@ -424,7 +416,7 @@ function grids_testsuite(AT, FT) level = 0 gm = GridManager( - LobattoCell{Tuple{L,M,N},FT,AT}(), + LobattoCell{FT,AT}(L, M, N), Raven.brick(coordinates); min_level = level, ) @@ -540,7 +532,7 @@ function grids_testsuite(AT, FT) level = 0 gm = GridManager( - LobattoCell{Tuple{L,M,N},FT,AT}(), + LobattoCell{FT,AT}(L, M, N), Raven.brick(coordinates); min_level = level, ) @@ -562,7 +554,7 @@ function grids_testsuite(AT, FT) end @testset "2D uniform brick grid" begin - cell = LobattoCell{Tuple{4,5},FT,AT}() + cell = LobattoCell{FT,AT}(4, 5) xrange = range(-FT(1000), stop = FT(1000), length = 21) yrange = range(-FT(2000), stop = FT(2000), length = 11) grid = generate(GridManager(cell, Raven.brick((xrange, yrange)))) @@ -607,7 +599,7 @@ function grids_testsuite(AT, FT) end @testset "3D uniform brick grid" begin - cell = LobattoCell{Tuple{4,5,6},FT,AT}() + cell = LobattoCell{FT,AT}(4, 5, 6) xrange = range(-FT(1000), stop = FT(1000), length = 21) yrange = range(-FT(2000), stop = FT(2000), length = 11) zrange = range(-FT(3000), stop = FT(3000), length = 6) @@ -732,7 +724,7 @@ function grids_testsuite(AT, FT) end gm = GridManager( - LobattoCell{Tuple{Nq...},FT,AT}(), + LobattoCell{FT,AT}(Nq...), Raven.brick(brickrange, ntuple(_ -> true, dim)), ) grid = generate(min_node_dist_warpfun, gm) @@ -763,18 +755,18 @@ function grids_testsuite(AT, FT) cg = coarsegrid(vertices, cells) coarse_grid = coarsegrid("curvedboxmesh2d.inp") - gm = GridManager(LobattoCell{Tuple{N...},FT,AT}(), cg, min_level = 1) - gmcurved = GridManager(LobattoCell{Tuple{N...},FT,AT}(), coarse_grid, min_level = 1) + gm = GridManager(LobattoCell{FT,AT}(N...), cg, min_level = 1) + gmcurved = GridManager(LobattoCell{FT,AT}(N...), coarse_grid, min_level = 1) grid = generate(gm) gridcurved = generate(gmcurved) @test coarse_grid.vertices ≈ cg.vertices cg1 = coarsegrid("flatGingerbreadMan.inp") - gm1 = GridManager(LobattoCell{Tuple{N...},FT,AT}(), cg1, min_level = 1) + gm1 = GridManager(LobattoCell{FT,AT}(N...), cg1, min_level = 1) grid = generate(gm1) cg2 = coarsegrid("GingerbreadMan.inp") - gm2 = GridManager(LobattoCell{Tuple{N...},FT,AT}(), cg2, min_level = 1) + gm2 = GridManager(LobattoCell{FT,AT}(N...), cg2, min_level = 1) grid2 = generate(gm2) @test cg1.vertices ≈ cg2.vertices @@ -801,19 +793,19 @@ function grids_testsuite(AT, FT) cg = coarsegrid(vertices, cells) coarse_grid = coarsegrid("curvedboxmesh3d.inp") - gm = GridManager(LobattoCell{Tuple{N...},FT,AT}(), cg, min_level = 1) - gmcurved = GridManager(LobattoCell{Tuple{N...},FT,AT}(), coarse_grid, min_level = 1) + gm = GridManager(LobattoCell{FT,AT}(N...), cg, min_level = 1) + gmcurved = GridManager(LobattoCell{FT,AT}(N...), coarse_grid, min_level = 1) grid = generate(gm) gridcurved = generate(gmcurved) @test coarse_grid.vertices ≈ cg.vertices cg1 = coarsegrid("flatHalfCircle3DRot.inp") - gm1 = GridManager(LobattoCell{Tuple{N...},FT,AT}(), cg1, min_level = 1) + gm1 = GridManager(LobattoCell{FT,AT}(N...), cg1, min_level = 1) grid = generate(gm1) cg2 = coarsegrid("HalfCircle3DRot.inp") - gm2 = GridManager(LobattoCell{Tuple{N...},FT,AT}(), cg2, min_level = 1) + gm2 = GridManager(LobattoCell{FT,AT}(N...), cg2, min_level = 1) grid2 = generate(gm2) @test cg1.vertices ≈ cg2.vertices diff --git a/test/kron.jl b/test/kron.jl index 0da633e..1db3ce6 100644 --- a/test/kron.jl +++ b/test/kron.jl @@ -6,12 +6,12 @@ function kron_testsuite(AT, FT) for args in ( (a,), - (a, Raven.Eye{FT,5}()), - (Raven.Eye{FT,2}(), b), + (a, Raven.Eye{FT}(5)), + (Raven.Eye{FT}(2), b), (a, b), - (Raven.Eye{FT,3}(), Raven.Eye{FT,2}(), c), - (Raven.Eye{FT,2}(), b, Raven.Eye{FT,7}()), - (a, Raven.Eye{FT,4}(), Raven.Eye{FT,7}()), + (Raven.Eye{FT}(3), Raven.Eye{FT}(2), c), + (Raven.Eye{FT}(2), b, Raven.Eye{FT}(7)), + (a, Raven.Eye{FT}(4), Raven.Eye{FT}(7)), (a, b, c), ) K = adapt(AT, collect(Raven.Kron(adapt(Array, args)))) @@ -44,7 +44,7 @@ function kron2dgridarray_testsuite(AT, FT) A = adapt(AT, rand(rng, FT, 5, 2)) B = adapt(AT, rand(rng, FT, 2, 3)) - cell = LobattoCell{Tuple{2,3},FT,AT}() + cell = LobattoCell{FT,AT}(2, 3) gm = GridManager(cell, Raven.brick(2, 1); min_level = 1) grid = generate(gm) @@ -52,7 +52,7 @@ function kron2dgridarray_testsuite(AT, FT) v = GridArray{EntryType}(undef, grid) v .= adapt(AT, rand(rng, EntryType, size(v))) - for args in ((Raven.Eye{FT,3}(), A), (B, Raven.Eye{FT,2}()), (B, A)) + for args in ((Raven.Eye{FT}(3), A), (B, Raven.Eye{FT}(2)), (B, A)) K = adapt(AT, collect(Raven.Kron(adapt(Array, args)))) Kv1 = adapt(Array, Raven.Kron(args) * v) dimsIn = (prod(size(v)[1:end-1]), size(v)[end]) @@ -68,7 +68,7 @@ function kron3dgridarray_testsuite(AT, FT) B = adapt(AT, rand(rng, FT, 2, 3)) C = adapt(AT, rand(rng, FT, 3, 5)) - cell = LobattoCell{Tuple{2,3,5},FT,AT}() + cell = LobattoCell{FT,AT}(2, 3, 5) gm = GridManager(cell, Raven.brick(2, 1, 1); min_level = 1) grid = generate(gm) @@ -77,12 +77,12 @@ function kron3dgridarray_testsuite(AT, FT) v .= adapt(AT, rand(rng, EntryType, size(v))) for args in ( - (Raven.Eye{FT,5}(), Raven.Eye{FT,3}(), A), - (Raven.Eye{FT,5}(), B, Raven.Eye{FT,2}()), - (C, Raven.Eye{FT,3}(), Raven.Eye{FT,2}()), - (Raven.Eye{FT,5}(), B, A), - (C, Raven.Eye{FT,3}(), A), - (C, B, Raven.Eye{FT,2}()), + (Raven.Eye{FT}(5), Raven.Eye{FT}(3), A), + (Raven.Eye{FT}(5), B, Raven.Eye{FT}(2)), + (C, Raven.Eye{FT}(3), Raven.Eye{FT}(2)), + (Raven.Eye{FT}(5), B, A), + (C, Raven.Eye{FT}(3), A), + (C, B, Raven.Eye{FT}(2)), (C, B, A), ) K = adapt(AT, collect(Raven.Kron(adapt(Array, args)))) diff --git a/test/mpitest_n2_commgridarrays.jl b/test/mpitest_n2_commgridarrays.jl index 7c2e6fd..cef4dae 100644 --- a/test/mpitest_n2_commgridarrays.jl +++ b/test/mpitest_n2_commgridarrays.jl @@ -12,7 +12,7 @@ function test(::Type{FT}, ::Type{AT}) where {FT,AT} N = (3, 2) K = (2, 1) min_level = 1 - cell = LobattoCell{Tuple{N...},Float64,AT}() + cell = LobattoCell{Float64,AT}(N...) gm = GridManager(cell, Raven.brick(K...); min_level) grid = generate(gm) diff --git a/test/mpitest_n2_gridnumbering.jl b/test/mpitest_n2_gridnumbering.jl index 9316048..40654a4 100644 --- a/test/mpitest_n2_gridnumbering.jl +++ b/test/mpitest_n2_gridnumbering.jl @@ -14,7 +14,7 @@ let AT = Array N = 4 - cell = LobattoCell{Tuple{N,N},FT,AT}() + cell = LobattoCell{FT,AT}(N, N) # 4--------5--------6 # | | | @@ -130,7 +130,7 @@ let AT = Array N = 3 - cell = LobattoCell{Tuple{N,N,N}}() + cell = LobattoCell(N, N, N) # 10-------11-------12 # /| /| /| diff --git a/test/mpitest_n3_gridarrays.jl b/test/mpitest_n3_gridarrays.jl index 74f5cf8..3657c27 100644 --- a/test/mpitest_n3_gridarrays.jl +++ b/test/mpitest_n3_gridarrays.jl @@ -19,7 +19,7 @@ end function test(N, K, ::Type{FT}, ::Type{AT}) where {FT,AT} @testset "GridArray ($N, $AT, $FT)" begin minlvl = 1 - cell = LobattoCell{Tuple{N...},Float64,AT}() + cell = LobattoCell{Float64,AT}(N...) gm = GridManager(cell, Raven.brick(K...); min_level = minlvl) grid = generate(gm) diff --git a/test/mpitest_n3_gridnumbering.jl b/test/mpitest_n3_gridnumbering.jl index 1b64d44..9504506 100644 --- a/test/mpitest_n3_gridnumbering.jl +++ b/test/mpitest_n3_gridnumbering.jl @@ -159,7 +159,7 @@ let ) end - cell_degree_3 = LobattoCell{Tuple{4,4},Float64,Array}() + cell_degree_3 = LobattoCell{Float64,Array}(4, 4) dtoc_degree_3 = Raven.materializedtoc(cell_degree_3, dtoc_degree_3_local, dtoc_degree_3_global) @@ -362,7 +362,7 @@ let ) end - cell_degree_3 = LobattoCell{Tuple{4,4,4},Float64,Array}() + cell_degree_3 = LobattoCell{Float64,Array}(4, 4, 4) dtoc_degree_3 = Raven.materializedtoc(cell_degree_3, dtoc_degree_3_local, dtoc_degree_3_global) if rank == 0