From 0c98fefe915b27c9ad6ab5b586acd85b6670304d Mon Sep 17 00:00:00 2001 From: Tomas Lycken Date: Fri, 26 Dec 2014 06:08:54 +0100 Subject: [PATCH] Infrastructure for padding the equation system The extra equations from the boundary conditions will not always allow indexing as we did before, by e.g. doing cm_d = cp_d at the lower boundary for some BC:s. In order to generalize better, we instead pad the equation system with one line for each extra equation to allow more general BC specs. This had the advantage of significantly simplifying some of the code (and, even better, some of the concepts) required to implement new functionality, so even though some of the underlying infrastructure is now more complicated than before, new contributors will hopefully have a lower threshold to start adding features. --- src/Interpolations.jl | 30 +++++-- src/quadratic.jl | 177 +++++++++++------------------------------- 2 files changed, 70 insertions(+), 137 deletions(-) diff --git a/src/Interpolations.jl b/src/Interpolations.jl index 33242f58..3d741474 100644 --- a/src/Interpolations.jl +++ b/src/Interpolations.jl @@ -2,7 +2,7 @@ module Interpolations using Base.Cartesian -import Base: size, eltype, getindex +import Base: size, eltype, getindex, ndims export Interpolation, @@ -50,8 +50,10 @@ end # However, all prefilters copy the array, so do that here as well prefilter{T,N,IT<:InterpolationType}(A::AbstractArray{T,N}, ::IT) = copy(A) -size(itp::Interpolation, d::Integer) = size(itp.coefs, d) -size(itp::Interpolation) = size(itp.coefs) +size{T,N,IT<:InterpolationType}(itp::Interpolation{T,N,IT}, d::Integer) = + size(itp.coefs, d) - 2*padding(IT()) +size(itp::Interpolation) = tuple(collect([size(itp,i) for i in 1:ndims(itp)])) +ndims(itp::Interpolation) = ndims(itp.coefs) eltype(itp::Interpolation) = eltype(itp.coefs) offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) : @@ -68,6 +70,14 @@ include("constant.jl") include("linear.jl") include("quadratic.jl") +# If nothing else is specified, don't pad at all +padding(::InterpolationType) = 0 +padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(IT()) +function similar_with_padding(A, it::InterpolationType) + pad = padding(it) + coefs = Array(eltype(A), [size(A,i)+2pad for i in 1:ndims(A)]...) + coefs, pad +end # This creates getindex methods for all supported combinations for IT in ( @@ -130,22 +140,28 @@ for IT in ( Quadratic{LinearBC,OnCell} ) @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT) - ret = similar(A) + ret, pad = similar_with_padding(A,it) szs = collect(size(A)) strds = collect(strides(A)) + strdsR = collect(strides(ret)) for dim in 1:N n = szs[dim] szs[dim] = 1 - M = prefiltering_system_matrix(eltype(A), n, it) + M, b = prefiltering_system(eltype(A), n+2pad, it) @nloops N i d->1:szs[d] begin cc = @ntuple N i strt = 1 + sum([(cc[i]-1)*strds[i] for i in 1:length(cc)]) + strtR = 1 + sum([(cc[i]-1)*strdsR[i] for i in 1:length(cc)]) rng = range(strt, strds[dim], n) - ret[rng] = M \ vec(A[rng]) - end + rngR = range(strtR, strdsR[dim], size(ret,dim)) + + b[1+pad:end-pad] += A[rng] + ret[rngR] = M \ b + b[1+pad:end-pad], A[rng] + end szs[dim] = n end ret diff --git a/src/quadratic.jl b/src/quadratic.jl index fd2b80c8..48a7c845 100644 --- a/src/quadratic.jl +++ b/src/quadratic.jl @@ -12,106 +12,12 @@ end function indices(q::Quadratic, N) quote + pad = padding($q) @nexprs $N d->begin ixp_d = ix_d + 1 ixm_d = ix_d - 1 - # fx_d is the distance from the middle point to the interpolation point - fx_d = x_d - convert(typeof(x_d), ix_d) - end - end -end - -function bc_gen(::Quadratic{Flat,OnGrid}, N) - quote - # After extrapolation has been handled, 1 <= x_d <= size(itp,d) - # The index is simply the closest integer. - @nexprs $N d->(ix_d = iround(x_d)) - # end - end -end - -function bc_gen(::Quadratic{Flat,OnCell}, N) - quote - # After extrapolation has been handled, 0.5 <= x_d <= size(itp,d)+.5 - @nexprs $N d->begin - # The index is the closest integer... - ix_d = iround(x_d) - - #...except in the case where x_d is actually at the upper edge; - # since size(itp,d)+.5 is rounded toward size(itp,d)+1, - # it needs special treatment*. - if x_d == size(itp,d)+.5 - ix_d -= 1 - end - end - end -end -function bc_gen(::Quadratic{LinearBC,OnCell}, N) - quote - @nexprs $N d->begin - if x_d < 1 - # extrapolate towards -∞ - fx_d = x_d - convert(typeof(x_d), 1) - k = itp[2] - itp[1] - - return itp[1] + k * fx_d - end - if x_d > size(itp, d) - # extrapolate towards ∞ - s_d = size(itp, d) - fx_d = x_d - convert(typeof(x_d), s_d) - k = itp[s_d] - itp[s_d - 1] - - return itp[s_d] + k * fx_d - end - - ix_d = iround(x_d) - end - end -end - -function indices(::Quadratic{Flat,OnGrid}, N) - quote - @nexprs $N d->begin - ixp_d = ix_d + 1 - ixm_d = ix_d - 1 - - fx_d = x_d - convert(typeof(x_d), ix_d) - - if ix_d == size(itp,d) - ixp_d = ixm_d - end - if ix_d == 1 - ixm_d = ixp_d - end - end - end -end -function indices(::Quadratic{Flat,OnCell}, N) - quote - @nexprs $N d->begin - ixp_d = ix_d + 1 - ixm_d = ix_d - 1 - - fx_d = x_d - convert(typeof(x_d), ix_d) - - if ix_d == size(itp,d) - ixp_d = ix_d - end - if ix_d == 1 - ixm_d = ix_d - end - end - end -end -function indices(::Quadratic{LinearBC,OnCell}, N) - quote - @nexprs $N d->begin - ixp_d = ix_d + 1 - ixm_d = ix_d - 1 - - fx_d = x_d - convert(typeof(x_d), ix_d) + fx_d = x_d - convert(typeof(x_d), ix_d - pad) end end end @@ -126,7 +32,7 @@ function coefficients(::Quadratic, N) end end -# This assumes integral values ixm_d, ix_d, and ixp_d (typically ixm_d = ix_d-1, ixp_d = ix_d+1, except at boundaries), +# This assumes integral values ixm_d, ix_d, and ixp_d, # coefficients cm_d, c_d, and cp_d, and an array itp.coefs function index_gen(degree::QuadraticDegree, N::Integer, offsets...) if length(offsets) < N @@ -140,47 +46,58 @@ function index_gen(degree::QuadraticDegree, N::Integer, offsets...) end end -function unmodified_system_matrix{T}(::Type{T}, n::Int, ::Quadratic) +# Quadratic interpolation has 1 extra coefficient at each boundary +# Therefore, make the coefficient array 1 larger in each direction, +# in each dimension. +padding(::Quadratic) = 1 + +function inner_system_diags{T}(::Type{T}, n::Int, ::Quadratic) du = fill(convert(T,1/8), n-1) d = fill(convert(T,3/4),n) dl = copy(du) + d[1] = d[end] = du[1] = dl[end] = convert(T,0) (dl,d,du) end -function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{ExtendInner}) - dl,d,du = unmodified_system_matrix(T,n,q) - d[1] = d[end] = 9/8 - du[1] = dl[end] = -1/4 - MT = lufact!(Tridiagonal(dl, d, du)) - U = zeros(T,n,2) - V = zeros(T,2,n) - C = zeros(T,2,2) - - C[1,1] = C[2,2] = 1/8 - U[1,1] = U[n,2] = 1. - V[1,3] = V[2,n-2] = 1. - - Woodbury(MT, U, C, V) -end - -function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell}) - dl,d,du = unmodified_system_matrix(T,n,q) - d[1] += 1/8 - d[end] += 1/8 - lufact!(Tridiagonal(dl, d, du)) +function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnCell}) + dl,d,du = inner_system_diags(T,n,q) + d[1] = d[end] = -1 + du[1] = dl[end] = 1 + lufact!(Tridiagonal(dl, d, du)), zeros(T, n) end -function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnGrid}) - dl,d,du = unmodified_system_matrix(T,n,q) - du[1] += 1/8 - dl[end] += 1/8 - lufact!(Tridiagonal(dl, d, du)) +function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{Flat,OnGrid}) + dl,d,du = inner_system_diags(T,n,q) + d[1] = d[end] = convert(T, -1) + du[1] = dl[end] = convert(T, 0) + + rowspec = zeros(T,n,2) + # first row last row + rowspec[1,1] = rowspec[n,2] = convert(T, 1) + colspec = zeros(T,2,n) + # third col third-to-last col + colspec[1,3] = colspec[2,n-2] = convert(T, 1) + valspec = zeros(T,2,2) + # val for [1,3], val for [n,n-2] + valspec[1,1] = valspec[2,2] = convert(T, 1) + + Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n) end -function prefiltering_system_matrix{T}(::Type{T}, n::Int, q::Quadratic{LinearBC,OnCell}) - dl,d,du = unmodified_system_matrix(T,n,q) - - d[1] = d[end] = 1 - du[1] = dl[end] = 0 - lufact!(Tridiagonal(dl, d, du)) +function prefiltering_system{T}(::Type{T}, n::Int, q::Quadratic{LinearBC}) + dl,d,du = inner_system_diags(T,n,q) + d[1] = d[end] = convert(T, 1) + du[1] = dl[end] = convert(T, -2) + + rowspec = zeros(T,n,2) + # first row last row + rowspec[1,1] = rowspec[n,2] = convert(T, 1) + colspec = zeros(T,2,n) + # third col third-to-last col + colspec[1,3] = colspec[2,n-2] = convert(T, 1) + valspec = zeros(T,2,2) + # val for [1,3], val for [n,n-2] + valspec[1,1] = valspec[2,2] = convert(T, 1) + + Woodbury(lufact!(Tridiagonal(dl, d, du)), rowspec, valspec, colspec), zeros(T, n) end