Skip to content
This repository was archived by the owner on Jul 19, 2023. It is now read-only.

Multi-Dimensional BC Operator #414

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Multi-Dimensional BC Operator #414

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b266846
correcting MultiDimBC Operator
mjsheikh Jun 25, 2021
b14b160
correcting MultiDimBC operator
mjsheikh Jun 25, 2021
6541d76
ModelingToolkit fallback to previous release
mjsheikh Jun 26, 2021
f9aced4
ModelingToolkit fallback to previous release
mjsheikh Jun 26, 2021
afd25b2
support for multiplication with a constant
mjsheikh Jun 28, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ DomainSets = "0.5"
ForwardDiff = "0.10"
LazyArrays = "0.17, 0.18, 0.19, 0.20, 0.21"
LazyBandedMatrices = "0.3, 0.4, 0.5"
ModelingToolkit = "4, 5"
ModelingToolkit = "4 - 5.18"
NNlib = "0.6, 0.7"
NonlinearSolve = "0.3.7"
RuntimeGeneratedFunctions = "0.4, 0.5"
Expand Down
358 changes: 179 additions & 179 deletions src/derivative_operators/concretization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -163,190 +163,190 @@ function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
end

# Multidimensional BC operators
_concretize(Q::MultiDimDirectionalBC, M) = _concretize(Q.BCs, M)
# _concretize(Q::MultiDimDirectionalBC, M) = _concretize(Q.BCs, M)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what's the reason for commenting instead of removing?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Earlier it was an Affine operator with concretizations. Now it isn't and hence won't support this but rather than removing I thought about commenting out keeping in mind that if by any means one can come up with a way to still concretize this later.
Now when I think of it, I don't know why this should necessarily be manipulated and concretized ┐(•_•)┌


function _concretize(Q::AbstractArray{T,N}, M) where {T,N}
return (stencil.(Q, fill(M,size(Q))), affine.(Q))
end

function LinearAlgebra.Array(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
@assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
blip = zeros(Int64, N)
blip[D] = 2
s_pad = s.+ blip # extend s in the right direction
Q = _concretize(Q.BCs, s[D])
ē = unit_indices(N)
QL = zeros(T, prod(s_pad), prod(s))
Qb = zeros(T, prod(s_pad))
ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
ranges[D] = ranges[D] .+ 1

interior = CartesianIndices(Tuple(ranges))
I1 = CartesianIndex(Tuple(ones(Int64, N)))
for I in interior
i = cartesian_to_linear(I, s_pad)
j = cartesian_to_linear(I-ē[D], s)
QL[i,j] = one(T)
end
ranges[D] = 1
lower = CartesianIndices((Tuple(ranges)))
ranges[D] = s_pad[D]
upper = CartesianIndices((Tuple(ranges)))
for K in CartesianIndices(upper)
I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
il = cartesian_to_linear(lower[K], s_pad)
iu = cartesian_to_linear(upper[K], s_pad)
Qb[il] = Q[2][I][1]
Qb[iu] = Q[2][I][2]
for k in 0:s[D]-1
j = cartesian_to_linear(lower[K] + k*ē[D], s)
QL[il, j] = Q[1][I][1][k+1]
QL[iu, j] = Q[1][I][2][k+1]
end
end

return (QL, Qb)
end

"""
This is confusing, but it does work.
"""
function LinearAlgebra.Array(Q::ComposedMultiDimBC{T, B, N,M} , s::NTuple{N,G}) where {T, B, N, M, G<:Int}
for d in 1:N
@assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
end
s_pad = s.+2
Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC

QL = zeros(T, prod(s_pad), prod(s))
Qb = zeros(T, prod(s_pad))

ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
interior = CartesianIndices(Tuple(ranges))

ē = unit_indices(N) #setup unit indices in each direction
I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
for I in interior #loop over interior
i = cartesian_to_linear(I, s_pad) #find the index on the padded side
j = cartesian_to_linear(I-I1, s) #find the index on the unpadded side
QL[i,j] = one(T) #create a padded identity matrix
end
for dim in 1:N #Loop over boundaries
r_ = deepcopy(ranges)
r_[dim] = 1
lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
r_[dim] = s_pad[dim]
upper = CartesianIndices((Tuple(r_)))
for K in CartesianIndices(upper) #for every element of the boundaries
I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
iu = cartesian_to_linear(upper[K], s_pad) # ditto
Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
Qb[iu] = Q[dim][2][I][2] #ditto with upper index
for k in 1:s[dim] #loop over the direction orthogonal to the boundary
j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
QL[il, j] = Q[dim][1][I][1][k]
QL[iu, j] = Q[dim][1][I][2][k]
end
end
end

return (QL, Qb)
end

"""
See comments on the `Array` method for this type for an idea of what is going on.
"""
function SparseArrays.SparseMatrixCSC(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
@assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
blip = zeros(Int64, N)
blip[D] = 2
s_pad = s.+ blip # extend s in the right direction
Q = _concretize(Q.BCs, s[D])
ē = unit_indices(N)
QL = spzeros(T, prod(s_pad), prod(s))
Qb = spzeros(T, prod(s_pad))
ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
ranges[D] = ranges[D] .+ 1

interior = CartesianIndices(Tuple(ranges))
I1 = CartesianIndex(Tuple(ones(Int64, N)))
for I in interior
i = cartesian_to_linear(I, s_pad)
j = cartesian_to_linear(I-ē[D], s)
QL[i,j] = one(T)
end
ranges[D] = 1
lower = CartesianIndices((Tuple(ranges)))
ranges[D] = s_pad[D]
upper = CartesianIndices((Tuple(ranges)))
for K in CartesianIndices(upper)
I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
il = cartesian_to_linear(lower[K], s_pad)
iu = cartesian_to_linear(upper[K], s_pad)
Qb[il] = Q[2][I][1]
Qb[iu] = Q[2][I][2]
for k in 0:s[D]-1
j = cartesian_to_linear(lower[K] + k*ē[D], s)
QL[il, j] = Q[1][I][1][k+1]
QL[iu, j] = Q[1][I][2][k+1]
end
end

return (QL, Qb)
end


function SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC{T, B, N,M}, s::NTuple{N,G}) where {T, B, N, M, G<:Int}
for d in 1:N
@assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
end
s_pad = s.+2
Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC

QL = spzeros(T, prod(s_pad), prod(s))
Qb = spzeros(T, prod(s_pad))

ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
interior = CartesianIndices(Tuple(ranges))

ē = unit_indices(N) #setup unit indices in each direction
I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
for I in interior #loop over interior
i = cartesian_to_linear(I, s_pad) #find the index on the padded side
j = cartesian_to_linear(I-I1, s) #find the index on the unpadded side
QL[i,j] = one(T) #create a padded identity matrix
end
for dim in 1:N #Loop over boundaries
r_ = deepcopy(ranges)
r_[dim] = 1
lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
r_[dim] = s_pad[dim]
upper = CartesianIndices((Tuple(r_)))
for K in CartesianIndices(upper) #for every element of the boundaries
I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
iu = cartesian_to_linear(upper[K], s_pad) # ditto
Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
Qb[iu] = Q[dim][2][I][2] #ditto with upper index
for k in 1:s[dim] #loop over the direction orthogonal to the boundary
j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
QL[il, j] = Q[dim][1][I][1][k]
QL[iu, j] = Q[dim][1][I][2][k]
end
end
end

return (QL, Qb)
end

SparseArrays.sparse(Q::MultiDimDirectionalBC, s) = SparseMatrixCSC(Q, s)
SparseArrays.sparse(Q::ComposedMultiDimBC, s) = SparseMatrixCSC(Q, s)


function BandedMatrices.BandedMatrix(Q:: MultiDimensionalBC, M) where {T, B, D,N,K}
throw("Banded Matrix cocnretization not yet supported for MultiDimensionalBCs")
end
# function LinearAlgebra.Array(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
# @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
# blip = zeros(Int64, N)
# blip[D] = 2
# s_pad = s.+ blip # extend s in the right direction
# Q = _concretize(Q.BCs, s[D])
# ē = unit_indices(N)
# QL = zeros(T, prod(s_pad), prod(s))
# Qb = zeros(T, prod(s_pad))
# ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
# ranges[D] = ranges[D] .+ 1

# interior = CartesianIndices(Tuple(ranges))
# I1 = CartesianIndex(Tuple(ones(Int64, N)))
# for I in interior
# i = cartesian_to_linear(I, s_pad)
# j = cartesian_to_linear(I-ē[D], s)
# QL[i,j] = one(T)
# end
# ranges[D] = 1
# lower = CartesianIndices((Tuple(ranges)))
# ranges[D] = s_pad[D]
# upper = CartesianIndices((Tuple(ranges)))
# for K in CartesianIndices(upper)
# I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
# il = cartesian_to_linear(lower[K], s_pad)
# iu = cartesian_to_linear(upper[K], s_pad)
# Qb[il] = Q[2][I][1]
# Qb[iu] = Q[2][I][2]
# for k in 0:s[D]-1
# j = cartesian_to_linear(lower[K] + k*ē[D], s)
# QL[il, j] = Q[1][I][1][k+1]
# QL[iu, j] = Q[1][I][2][k+1]
# end
# end

# return (QL, Qb)
# end

# """
# This is confusing, but it does work.
# """
# function LinearAlgebra.Array(Q::ComposedMultiDimBC{T, B, N,M} , s::NTuple{N,G}) where {T, B, N, M, G<:Int}
# for d in 1:N
# @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
# end
# s_pad = s.+2
# Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC

# QL = zeros(T, prod(s_pad), prod(s))
# Qb = zeros(T, prod(s_pad))

# ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
# interior = CartesianIndices(Tuple(ranges))

# ē = unit_indices(N) #setup unit indices in each direction
# I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
# for I in interior #loop over interior
# i = cartesian_to_linear(I, s_pad) #find the index on the padded side
# j = cartesian_to_linear(I-I1, s) #find the index on the unpadded side
# QL[i,j] = one(T) #create a padded identity matrix
# end
# for dim in 1:N #Loop over boundaries
# r_ = deepcopy(ranges)
# r_[dim] = 1
# lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
# r_[dim] = s_pad[dim]
# upper = CartesianIndices((Tuple(r_)))
# for K in CartesianIndices(upper) #for every element of the boundaries
# I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
# il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
# iu = cartesian_to_linear(upper[K], s_pad) # ditto
# Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
# Qb[iu] = Q[dim][2][I][2] #ditto with upper index
# for k in 1:s[dim] #loop over the direction orthogonal to the boundary
# j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
# QL[il, j] = Q[dim][1][I][1][k]
# QL[iu, j] = Q[dim][1][I][2][k]
# end
# end
# end

# return (QL, Qb)
# end

# """
# See comments on the `Array` method for this type for an idea of what is going on.
# """
# function SparseArrays.SparseMatrixCSC(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
# @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
# blip = zeros(Int64, N)
# blip[D] = 2
# s_pad = s.+ blip # extend s in the right direction
# Q = _concretize(Q.BCs, s[D])
# ē = unit_indices(N)
# QL = spzeros(T, prod(s_pad), prod(s))
# Qb = spzeros(T, prod(s_pad))
# ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
# ranges[D] = ranges[D] .+ 1

# interior = CartesianIndices(Tuple(ranges))
# I1 = CartesianIndex(Tuple(ones(Int64, N)))
# for I in interior
# i = cartesian_to_linear(I, s_pad)
# j = cartesian_to_linear(I-ē[D], s)
# QL[i,j] = one(T)
# end
# ranges[D] = 1
# lower = CartesianIndices((Tuple(ranges)))
# ranges[D] = s_pad[D]
# upper = CartesianIndices((Tuple(ranges)))
# for K in CartesianIndices(upper)
# I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
# il = cartesian_to_linear(lower[K], s_pad)
# iu = cartesian_to_linear(upper[K], s_pad)
# Qb[il] = Q[2][I][1]
# Qb[iu] = Q[2][I][2]
# for k in 0:s[D]-1
# j = cartesian_to_linear(lower[K] + k*ē[D], s)
# QL[il, j] = Q[1][I][1][k+1]
# QL[iu, j] = Q[1][I][2][k+1]
# end
# end

# return (QL, Qb)
# end


# function SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC{T, B, N,M}, s::NTuple{N,G}) where {T, B, N, M, G<:Int}
# for d in 1:N
# @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
# end
# s_pad = s.+2
# Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC

# QL = spzeros(T, prod(s_pad), prod(s))
# Qb = spzeros(T, prod(s_pad))

# ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
# interior = CartesianIndices(Tuple(ranges))

# ē = unit_indices(N) #setup unit indices in each direction
# I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
# for I in interior #loop over interior
# i = cartesian_to_linear(I, s_pad) #find the index on the padded side
# j = cartesian_to_linear(I-I1, s) #find the index on the unpadded side
# QL[i,j] = one(T) #create a padded identity matrix
# end
# for dim in 1:N #Loop over boundaries
# r_ = deepcopy(ranges)
# r_[dim] = 1
# lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
# r_[dim] = s_pad[dim]
# upper = CartesianIndices((Tuple(r_)))
# for K in CartesianIndices(upper) #for every element of the boundaries
# I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
# il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
# iu = cartesian_to_linear(upper[K], s_pad) # ditto
# Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
# Qb[iu] = Q[dim][2][I][2] #ditto with upper index
# for k in 1:s[dim] #loop over the direction orthogonal to the boundary
# j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
# QL[il, j] = Q[dim][1][I][1][k]
# QL[iu, j] = Q[dim][1][I][2][k]
# end
# end
# end

# return (QL, Qb)
# end

# SparseArrays.sparse(Q::MultiDimDirectionalBC, s) = SparseMatrixCSC(Q, s)
# SparseArrays.sparse(Q::ComposedMultiDimBC, s) = SparseMatrixCSC(Q, s)


# function BandedMatrices.BandedMatrix(Q:: MultiDimensionalBC, M) where {T, B, D,N,K}
# throw("Banded Matrix cocnretization not yet supported for MultiDimensionalBCs")
# end

################################################################################
# Higher-Dimensional DerivativeOperator Concretizations
Expand Down
Loading