diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 12d54f239..59c9e5c1c 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -25,8 +25,6 @@ include("derivative_operators/BC_operators.jl") ### Derivative Operators include("derivative_operators/fornberg.jl") -include("derivative_operators/upwind_operator.jl") -include("derivative_operators/derivative_irreg_operator.jl") include("derivative_operators/derivative_operator.jl") include("derivative_operators/abstract_operator_functions.jl") include("derivative_operators/convolutions.jl") @@ -44,6 +42,7 @@ end export MatrixFreeOperator export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops -export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference -export RobinBC +export AbstractDerivativeOperator, DerivativeOperator, + CenteredDifference, UpwindDifference +export RobinBC, GeneralBC end # module diff --git a/src/derivative_operators/abstract_operator_functions.jl b/src/derivative_operators/abstract_operator_functions.jl index 128e87b41..2ca5e35e4 100644 --- a/src/derivative_operators/abstract_operator_functions.jl +++ b/src/derivative_operators/abstract_operator_functions.jl @@ -20,7 +20,7 @@ checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) = @inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int) @boundscheck checkbounds(A, i, j) bpc = A.boundary_point_count - N = A.dimension + N = A.len bsl = A.boundary_stencil_length slen = A.stencil_length if bpc > 0 && 1<=i<=bpc @@ -104,16 +104,16 @@ end # Base.length(A::AbstractDerivativeOperator) = A.stencil_length Base.ndims(A::AbstractDerivativeOperator) = 2 -Base.size(A::AbstractDerivativeOperator) = (A.dimension, A.dimension + 2) +Base.size(A::AbstractDerivativeOperator) = (A.len, A.len + 2) Base.size(A::AbstractDerivativeOperator,i::Integer) = size(A)[i] Base.length(A::AbstractDerivativeOperator) = reduce(*, size(A)) #= For the evenly spaced grid we have a symmetric matrix =# -Base.transpose(A::Union{DerivativeOperator,UpwindOperator}) = A -Base.adjoint(A::Union{DerivativeOperator,UpwindOperator}) = A -LinearAlgebra.issymmetric(::Union{DerivativeOperator,UpwindOperator}) = true +Base.transpose(A::DerivativeOperator) = A +Base.adjoint(A::DerivativeOperator) = A +LinearAlgebra.issymmetric(::DerivativeOperator) = true #= Fallback methods that use the full representation of the operator @@ -164,3 +164,16 @@ end function *(A::AbstractDerivativeOperator,B::AbstractDerivativeOperator) return BandedMatrix(A)*BandedMatrix(B) end + +################################################################################ + +function DiffEqBase.update_coefficients!(A::AbstractDerivativeOperator,u,p,t) + if A.coeff_func !== nothing + A.coeff_func(A.coefficients,u,p,t) + end +end + +################################################################################ + +(L::DerivativeOperator)(u,p,t) = L*u +(L::DerivativeOperator)(du,u,p,t) = mul!(du,L,u) diff --git a/src/derivative_operators/concretization.jl b/src/derivative_operators/concretization.jl index 2a6362434..badc88b32 100644 --- a/src/derivative_operators/concretization.jl +++ b/src/derivative_operators/concretization.jl @@ -1,4 +1,4 @@ -function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.dimension) where T +function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T L = zeros(T, N, N+2) bl = A.boundary_point_count stl = A.stencil_length @@ -16,7 +16,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.dimension) where return L / A.dx^A.derivative_order end -function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.dimension) where T +function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T L = spzeros(T, N, N+2) bl = A.boundary_point_count stl = A.stencil_length @@ -34,12 +34,11 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.dimensi return L / A.dx^A.derivative_order end -function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T +function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.len) where T SparseMatrixCSC(A,N) end -function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.dimension) where T - N = A.dimension +function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) where T bl = A.boundary_point_count stl = A.stencil_length bstl = A.boundary_stencil_length diff --git a/src/derivative_operators/convolutions.jl b/src/derivative_operators/convolutions.jl index 53b13e457..cccfdac82 100644 --- a/src/derivative_operators/convolutions.jl +++ b/src/derivative_operators/convolutions.jl @@ -1,5 +1,5 @@ # mul! done by convolutions -function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::Union{DerivativeOperator{T},UpwindOperator{T}}, x::AbstractVector{T}) where T<:Real +function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::DerivativeOperator, x::AbstractVector{T}) where T<:Real convolve_BC_left!(x_temp, x, A) convolve_interior!(x_temp, x, A) convolve_BC_right!(x_temp, x, A) @@ -9,102 +9,101 @@ end ################################################ # Against a standard vector, assume already padded and just apply the stencil -function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} +function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real} @assert length(x_temp)+2 == length(x) - coeffs = A.stencil_coefs + stencil = A.stencil_coefs + coeff = A.coefficients mid = div(A.stencil_length,2) for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count) xtempi = zero(T) + cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i] : stencil + cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i] : true + cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil for idx in 1:A.stencil_length - xtempi += coeffs[idx] * x[i - mid + idx] + xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx] end x_temp[i] = xtempi end end -function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} - coeffs = A.low_boundary_coefs +function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real} + stencil = A.low_boundary_coefs + coeff = A.coefficients for i in 1 : A.boundary_point_count - xtempi = coeffs[i][1]*x[1] + xtempi = stencil[i][1]*x[1] + cur_stencil = stencil[i] + cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i] : true + cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil for idx in 2:A.boundary_stencil_length - xtempi += coeffs[i][idx] * x[idx] + xtempi += cur_coeff * cur_stencil[idx] * x[idx] end x_temp[i] = xtempi end end -function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} - coeffs = A.high_boundary_coefs +function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real} + stencil = A.high_boundary_coefs + coeff = A.coefficients for i in 1 : A.boundary_point_count - xtempi = coeffs[i][end]*x[end] + xtempi = stencil[i][end]*x[end] + cur_stencil = stencil[i] + cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i] : true + cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil for idx in (A.boundary_stencil_length-1):-1:1 - xtempi += coeffs[i][end-idx] * x[end-idx] + xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx] end x_temp[end-A.boundary_point_count+i] = xtempi end end ########################################### -# NOT NEEDED ANYMORE # Against A BC-padded vector, specialize the computation to explicitly use the left, right, and middle parts -function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} - coeffs = A.stencil_coefs +function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real} + stencil = A.stencil_coefs + coeff = A.coefficients x = _x.u mid = div(A.stencil_length,2) + 1 # Just do the middle parts for i in (1+A.boundary_length) : (length(x_temp)-A.boundary_length) xtempi = zero(T) + cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_length] : stencil + cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true + cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil @inbounds for idx in 1:A.stencil_length - xtempi += coeffs[idx] * x[i - (mid-idx) + 1] + xtempi += cur_coeff * cur_stencil[idx] * x[i - (mid-idx) + 1] end x_temp[i] = xtempi end end -function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} - coeffs = A.low_boundary_coefs +function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real} + stencil = A.low_boundary_coefs + coeff = A.coefficients for i in 1 : A.boundary_length - xtempi = coeffs[i][1]*x.l + xtempi = stencil[i][1]*x.l + cur_stencil = stencil[i] + cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true + cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil @inbounds for idx in 2:A.stencil_length - xtempi += coeffs[i][idx] * x.u[idx-1] + xtempi += cur_coeff * cur_stencil[idx] * x.u[idx-1] end x_temp[i] = xtempi end end -function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} - coeffs = A.low_boundary_coefs +function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real} + stencil = A.low_boundary_coefs + coeff = A.coefficients bc_start = length(x.u) - A.stencil_length for i in 1 : A.boundary_length - xtempi = coeffs[i][end]*x.r + xtempi = stencil[i][end]*x.r + cur_stencil = stencil[i] + cur_coeff = eltype(coeff) <: AbstractVector ? coeff[i-A.boundary_length] : true + cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil @inbounds for idx in A.stencil_length:-1:2 - xtempi += coeffs[i][end-idx] * x.u[end-idx+1] + xtempi += cur_coeff * cur_stencil[end-idx] * x.u[end-idx+1] end x_temp[i] = xtempi end end - -########################################## - -#= INTERIOR CONVOLUTION =# -function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,LBC,RBC}) where {T<:Real,S<:SVector,LBC,RBC} - N = length(x) - stencil_length = length(A.up_stencil_coefs) - stencil_rem = 1 - stencil_length%2 - Threads.@threads for i in A.boundary_point_count[1]+1 : N-A.boundary_point_count[2] - xtempi = zero(T) - if A.directions[][i] == false - @inbounds for j in 1 : length(A.up_stencil_coefs) - xtempi += A.up_stencil_coefs[j] * x[i+j-1-stencil_rem] - end - else - @inbounds for j in -length(A.down_stencil_coefs)+1 : 0 - xtempi += A.down_stencil_coefs[j+stencil_length] * x[i+j+stencil_rem] - # println("i = $i, j = $j, s_idx = $(j+stencil_length), x_idx = $(i+j+stencil_rem), $(A.down_stencil_coefs[j+stencil_length]) * $(x[i+j+stencil_rem]), xtempi = $xtempi") - end - end - - x_temp[i] = xtempi - end -end diff --git a/src/derivative_operators/derivative_irreg_operator.jl b/src/derivative_operators/derivative_irreg_operator.jl deleted file mode 100644 index ef2fb1570..000000000 --- a/src/derivative_operators/derivative_irreg_operator.jl +++ /dev/null @@ -1,261 +0,0 @@ -struct FiniteDifference{T<:Real,S<:SVector,LBC,RBC} <: AbstractDerivativeOperator{T} - derivative_order :: Int - approximation_order :: Int - dx :: Vector{T} - dimension :: Int - stencil_length :: Int - stencil_coefs :: Vector{S} - boundary_point_count:: Tuple{Int,Int} - boundary_length :: Tuple{Int,Int} - low_boundary_coefs :: Ref{Vector{Vector{T}}} - high_boundary_coefs :: Ref{Vector{Vector{T}}} - boundary_condition :: Ref{Tuple{Tuple{T,T,Any},Tuple{T,T,Any}}} - t :: Ref{Int} - - function FiniteDifference{T,S,LBC,RBC}(derivative_order::Int, - approximation_order::Int, - dx::Vector{T}, - dimension::Int, BC) where - {T<:Real,S<:SVector,LBC,RBC} - # dimension = dimension - # dx = dx - stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 - bl = derivative_order + approximation_order - boundary_length = (bl,bl) - bpc = stencil_length - div(stencil_length,2) + 1 - bpc_array = [bpc,bpc] - grid_step = dx #unnecessary - - if any(x->xBC[1] : ret = BC[1] - return (one(T),zero(T),ret) - - elseif LBC == :Neumann0 - return (zero(T),one(T),zero(T)) - - elseif LBC == :periodic - return (zero(T),zero(T),zero(T)) - - else - error("Unrecognized Boundary Type!") - end -end - - -#= - This function sets us up to apply the boundary condition correctly. It returns a - 3-tuple which is basically the coefficients in the equation - a*u + b*du/dt = c(t) - The RHS ie. 'c' can be a function of time as well and therefore it is implemented - as an anonymous function. -=# -function initialize_right_boundary!(::Type{Val{:FD}},high_boundary_coefs,stencil_coefs,BC,derivative_order,grid_step::Vector{T}, - boundary_length,boundary_point_count,RBC) where T - stencil_length = length(stencil_coefs[1]) - - if RBC == :None - boundary_point_count[2] = div(stencil_length,2) - #= - Val{:FD} is type parametrization on symbols. There are two different right_None_BC! - functions, one for LinearOperator and one for UpwindOperator. So to select the correct - function without changing the name, this trick has been applied. - =# - return (zero(T),zero(T),right_None_BC!(Val{:FD},high_boundary_coefs,stencil_length,derivative_order, - grid_step,boundary_length)*BC[2]) - elseif RBC == :Neumann - error("RBC Robin not yet supported for irregular grid.") - # return (zero(T),one(T),right_Neumann_BC!(Val{:FD},high_boundary_coefs,stencil_length,derivative_order, - # grid_step,boundary_length)*BC[2]) - elseif RBC == :Robin - error("RBC Robin not yet supported for irregular grid.") - # return (BC[2][1],BC[2][2],right_Robin_BC!(Val{:FD},high_boundary_coefs,stencil_length, - # BC[2],derivative_order,grid_step, - # boundary_length,dx)*BC[2][3]) - elseif RBC == :Dirichlet0 - return (one(T),zero(T),zero(T)*BC[2]) - - elseif RBC == :Dirichlet - typeof(BC[2]) <: Real ? ret = t->BC[2] : ret = BC[2] - return (one(T),zero(T),ret) - - elseif RBC == :Neumann0 - return (zero(T),one(T),zero(T)) - - elseif RBC == :periodic - return (zero(T),zero(T),zero(T)) - - else - error("Unrecognized Boundary Type!") - end -end - - -function left_None_BC!(::Type{Val{:FD}},low_boundary_coefs,stencil_length,derivative_order, - grid_step::Vector{T},boundary_length) where T - # Fixes the problem excessive boundary points - boundary_point_count = div(stencil_length,2) - l_diff = zero(T) - mid = div(stencil_length,2) - x = [0.0; cumsum(grid_step)] - for i in 1 : boundary_point_count - # One-sided stencils require more points for same approximation order - # TODO: I don't know if this is the correct stencil length for i > 1? - push!(low_boundary_coefs, calculate_weights(derivative_order, zero(T), x[1:boundary_length].-x[i])) - end - return l_diff -end - - -function right_None_BC!(::Type{Val{:FD}},high_boundary_coefs,stencil_length,derivative_order, - grid_step::Vector{T},boundary_length) where T - boundary_point_count = div(stencil_length,2) - high_temp = zeros(T,boundary_length) - flag = derivative_order*boundary_point_count%2 - aorder = boundary_length - 1 - r_diff = zero(T) - x = [0.0; cumsum(grid_step)] - - for i in length(x) : -1 : length(x) - boundary_point_count + 1 - # One-sided stencils require more points for same approximation order - push!(high_boundary_coefs, calculate_weights(derivative_order,x[i],x[end - boundary_length + 1 : end])) - end - return r_diff -end - - -function left_Neumann_BC!(::Type{Val{:FD}},low_boundary_coefs,stencil_length,derivative_order, - grid_step::Vector{T},boundary_length) where T - boundary_point_count = stencil_length - div(stencil_length,2) + 1 - # first_order_coeffs = zeros(T,boundary_length) - # original_coeffs = zeros(T,boundary_length) - l_diff = one(T) - mid = div(stencil_length,2)+1 - x = [0.0;cumsum(grid_step)] - - first_order_coeffs = calculate_weights(1, zero(T), x[1:boundary_length]) - original_coeffs = calculate_weights(derivative_order, zero(T), x[1:boundary_length]) - l_diff = original_coeffs[end]/first_order_coeffs[end] - rmul!(first_order_coeffs, original_coeffs[end]/first_order_coeffs[end]) - # rmul!(original_coeffs, first_order_coeffs[end]/original_coeffs[end]) - @. original_coeffs = original_coeffs - first_order_coeffs - # copyto!(first_order_coeffs, first_order_coeffs[1:end-1]) - push!(low_boundary_coefs, original_coeffs[1:end-1]) - - for i in 2 : boundary_point_count - #= this means that a stencil will suffice ie. we dont't need to worry about the boundary point - being considered in the low_boundary_coefs - =# - - if i > mid - pos=i - push!(low_boundary_coefs, calculate_weights(derivative_order, x[pos], x[1:boundary_length])) - else - pos=i-1 - push!(low_boundary_coefs, append!([zero(T)],calculate_weights(derivative_order, x[pos], x[1:boundary_length]))) - end - end - - return l_diff -end diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index b6518fb67..47479c77f 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -1,55 +1,140 @@ -struct DerivativeOperator{T<:Real,S1<:SVector,S2<:SVector} <: AbstractDerivativeOperator{T} +struct DerivativeOperator{T<:Real,T2,S1,S2<:SVector,T3,F} <: AbstractDerivativeOperator{T} derivative_order :: Int approximation_order :: Int - dx :: T - dimension :: Int + dx :: T2 + len :: Int stencil_length :: Int stencil_coefs :: S1 boundary_stencil_length :: Int boundary_point_count :: Int - low_boundary_coefs :: Vector{S2} - high_boundary_coefs :: Vector{S2} - - function DerivativeOperator{T,S1,S2}(derivative_order::Int, - approximation_order::Int, dx::T, - dimension::Int) where - {T<:Real,S1<:SVector,S2<:SVector} - dimension = dimension - dx = dx - stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 - boundary_stencil_length = derivative_order + approximation_order - dummy_x = -div(stencil_length,2) : div(stencil_length,2) - boundary_x = -boundary_stencil_length+1:0 - boundary_point_count = div(stencil_length,2) - 1 # -1 due to the ghost point - # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point. - deriv_spots = (-div(stencil_length,2)+1) : -1 - boundary_deriv_spots = boundary_x[2:div(stencil_length,2)] - - stencil_coefs = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x)) - low_boundary_coefs = [convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots] - high_boundary_coefs = reverse!(copy([convert(SVector{boundary_stencil_length, T}, reverse(low_boundary_coefs[i])) for i in 1:boundary_point_count])) - - new(derivative_order, approximation_order, dx, dimension, stencil_length, - stencil_coefs, - boundary_stencil_length, - boundary_point_count, - low_boundary_coefs, - high_boundary_coefs - ) - end - DerivativeOperator{T}(dorder::Int,aorder::Int,dx::T,dim::Int) where {T<:Real} = - DerivativeOperator{T, SVector{dorder+aorder-1+(dorder+aorder)%2,T}, SVector{dorder+aorder,T}}(dorder, aorder, dx, dim) + low_boundary_coefs :: S2 + high_boundary_coefs :: S2 + coefficients :: T3 + coeff_func :: F + use_winding :: Bool end -#= - This function is used to update the boundary conditions especially if they evolve with - time. -=# -function DiffEqBase.update_coefficients!(A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector} - nothing +function CenteredDifference(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int) where {T<:Real} + + stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 + boundary_stencil_length = derivative_order + approximation_order + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + boundary_x = -boundary_stencil_length+1:0 + boundary_point_count = div(stencil_length,2) - 1 # -1 due to the ghost point + # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point. + deriv_spots = (-div(stencil_length,2)+1) : -1 + boundary_deriv_spots = boundary_x[2:div(stencil_length,2)] + + stencil_coefs = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x)) + _low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots] + low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) + high_boundary_coefs = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count])) + + DerivativeOperator{T,T,typeof(stencil_coefs), + typeof(low_boundary_coefs),Nothing, + Nothing}( + derivative_order, approximation_order, dx, len, stencil_length, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,nothing,nothing,false + ) +end + +function CenteredDifference(derivative_order::Int, + approximation_order::Int, dx::AbstractVector{T}, + len::Int) where {T<:Real} + + stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 + boundary_stencil_length = derivative_order + approximation_order + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + boundary_x = -boundary_stencil_length+1:0 + boundary_point_count = div(stencil_length,2) - 1 # -1 due to the ghost point + # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point. + deriv_spots = (-div(stencil_length,2)+1) : -1 + boundary_deriv_spots = boundary_x[2:div(stencil_length,2)] + + stencil_coefs = [convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x)) for i in 1:len(dx)] + _low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots] + low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) + high_boundary_coefs = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count])) + + DerivativeOperator{T,typeof(dx),typeof(stencil_coefs), + typeof(low_boundary_coefs),Nothing, + Nothing}( + derivative_order, approximation_order, dx, + len, stencil_length, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,nothing,nothing,false + ) +end + +function UpwindDifference(derivative_order::Int, + approximation_order::Int, dx::T, + len::Int, coeff_func) where {T<:Real} + + stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 + boundary_stencil_length = derivative_order + approximation_order + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + boundary_x = -boundary_stencil_length+1:0 + boundary_point_count = div(stencil_length,2) - 1 # -1 due to the ghost point + # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point. + deriv_spots = (-div(stencil_length,2)+1) : -1 + boundary_deriv_spots = boundary_x[2:div(stencil_length,2)] + + stencil_coefs = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x)) + _low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots] + low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) + high_boundary_coefs = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count])) + + coefficients = Vector{T}(undef,len) + + DerivativeOperator{T,T,typeof(stencil_coefs), + typeof(low_boundary_coefs),Vector{T}, + typeof(coeff_func)}( + derivative_order, approximation_order, dx, len, stencil_length, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,coefficients,coeff_func,true + ) end -################################################################################################# +function UpwindDifference(derivative_order::Int, + approximation_order::Int, dx::AbstractVector{T}, + len::Int, coeff_func) where {T<:Real} -(L::DerivativeOperator)(u,p,t) = L*u -(L::DerivativeOperator)(du,u,p,t) = mul!(du,L,u) + stencil_length = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2 + boundary_stencil_length = derivative_order + approximation_order + dummy_x = -div(stencil_length,2) : div(stencil_length,2) + boundary_x = -boundary_stencil_length+1:0 + boundary_point_count = div(stencil_length,2) - 1 # -1 due to the ghost point + # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point. + deriv_spots = (-div(stencil_length,2)+1) : -1 + boundary_deriv_spots = boundary_x[2:div(stencil_length,2)] + + stencil_coefs = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x)) + _low_boundary_coefs = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots] + low_boundary_coefs = convert(SVector{boundary_point_count},_low_boundary_coefs) + high_boundary_coefs = convert(SVector{boundary_point_count},reverse(SVector{boundary_stencil_length, T}[reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count])) + + coefficients = Vector{T}(undef,len) + + DerivativeOperator{T,typeof(dx),typeof(stencil_coefs), + typeof(low_boundary_coefs),Vector{T}, + typeof(coeff_func)}( + derivative_order, approximation_order, dx, len, stencil_length, + stencil_coefs, + boundary_stencil_length, + boundary_point_count, + low_boundary_coefs, + high_boundary_coefs,coefficients,coeff_func,true + ) +end diff --git a/src/derivative_operators/upwind_operator.jl b/src/derivative_operators/upwind_operator.jl deleted file mode 100644 index e49a11426..000000000 --- a/src/derivative_operators/upwind_operator.jl +++ /dev/null @@ -1,247 +0,0 @@ -function negate!(arr::T) where T - if size(arr,2) == 1 - rmul!(arr,-one(eltype(arr[1]))) #fix right neumann bc, eltype(Vector{T}) doesnt work. - else - for row in arr - rmul!(row,-one(eltype(arr[1]))) - end - end -end - -struct UpwindOperator{T<:Real,S<:SVector,LBC,RBC} <: AbstractDerivativeOperator{T} - derivative_order :: Int - approximation_order :: Int - dx :: T - dimension :: Int - directions :: Ref{BitArray{1}} - stencil_length :: Int - up_stencil_coefs :: S - down_stencil_coefs :: S - boundary_point_count:: Tuple{Int,Int} - boundary_length :: Tuple{Int,Int} - low_boundary_coefs :: Ref{Vector{Vector{T}}} - high_boundary_coefs :: Ref{Vector{Vector{T}}} - boundary_condition :: Ref{Tuple{Tuple{T,T,Any},Tuple{T,T,Any}}} - t :: Ref{Int} - - Base.@pure function UpwindOperator{T,S,LBC,RBC}(derivative_order::Int, approximation_order::Int, dx::T, - dimension::Int, directions::BitArray{1}, bndry_fn) where {T<:Real,S<:SVector,LBC,RBC} - dimension = dimension - dx = dx - directions = directions - stencil_length = derivative_order + approximation_order - bl = derivative_order + approximation_order - - boundary_length = (bl,bl) - grid_step = one(T) - low_boundary_coefs = Vector{T}[] - high_boundary_coefs = Vector{T}[] - - #= - We are implementing biased Upwind Operators which use a point from the other direction - also to ensure a more stable solution. - http://ac.els-cdn.com/S0378475401002889/1-s2.0-S0378475401002889-main.pdf?_tid=534e0818-8b0b-11e7-8b12-00000aab0f01&acdnat=1503826821_3bf9422abe7aa2d3613c5b644b6e258f - page 11 - =# - up_stencil_coefs = convert(SVector{stencil_length, T}, calculate_weights(derivative_order,convert(T,(stencil_length+1)%2), - grid_step .* collect(zero(T) : grid_step : stencil_length-1))) - - down_stencil_coefs = reverse(up_stencil_coefs) - - # no need to flip the sign of downwind operator, doesn't converge that way - # derivative_order%2 == 1 ? negate!(down_stencil_coefs) : nothing - - bpc_array = [length(up_stencil_coefs)-1,length(down_stencil_coefs)-1] - - left_bndry = initialize_left_boundary!(Val{:UO},low_boundary_coefs,bndry_fn,derivative_order,grid_step, - bl,bpc_array,directions,dx,LBC) - right_bndry = initialize_right_boundary!(Val{:UO},high_boundary_coefs,bndry_fn,derivative_order,grid_step, - bl,bpc_array,directions,dx,RBC) - - boundary_condition = (left_bndry, right_bndry) - boundary_point_count = (bpc_array[1],bpc_array[2]) - t = 0 - - new(derivative_order, approximation_order, dx, dimension, directions, - stencil_length, - up_stencil_coefs, - down_stencil_coefs, - boundary_point_count, - boundary_length, - low_boundary_coefs, - high_boundary_coefs, - boundary_condition, - t - ) - end - UpwindOperator{T}(dorder::Int,aorder::Int,dx::T,dim::Int,direction::BitArray{1},LBC::Symbol,RBC::Symbol;BC=(zero(T),zero(T))) where {T<:Real} = UpwindOperator{T, SVector{dorder+aorder,T}, LBC, RBC}(dorder, aorder, dx, dim, direction, BC) -end - - -(L::UpwindOperator)(t,u) = L*u -(L::UpwindOperator)(t,u,du) = mul!(du,L,u) - -function update_coefficients!(A::UpwindOperator{T,S,LBC,RBC};BC=nothing, directions=nothing) where {T<:Real,S<:SVector,RBC,LBC} - if BC != nothing - LBC == :Robin ? (length(BC[1])==3 || error("Enter the new left boundary condition as a 1-tuple")) : - (length(BC[1])==1 || error("Robin BC needs a 3-tuple for left boundary condition")) - - RBC == :Robin ? length(BC[2])==3 || error("Enter the new right boundary condition as a 1-tuple") : - length(BC[2])==1 || error("Robin BC needs a 3-tuple for right boundary condition") - - left_bndry = initialize_left_boundary!(A.low_boundary_coefs[],A.stencil_coefs,BC, - A.derivative_order,one(T),A.boundary_length,A.dx,A.directions,LBC) - - right_bndry = initialize_right_boundary!(A.high_boundary_coefs[],A.stencil_coefs,BC, - A.derivative_order,one(T),A.boundary_length,A.dx,A.directions,RBC) - - boundary_condition = (left_bndry, right_bndry) - A.boundary_condition[] = boundary_condition - end - if directions != nothing - A.directions[] = directions - end -end - - -function initialize_left_boundary!(::Type{Val{:UO}},low_boundary_coefs,BC,derivative_order,grid_step::T, - boundary_length,boundary_point_count,directions,dx,LBC) where T - approximation_order = boundary_length - derivative_order - up_stencil_length = boundary_length - - if LBC == :None - # up_stencil_length%2 == 1 ? boundary_point_count[1] = 1 : boundary_point_count[1] = 0 - return (zero(T),zero(T),left_None_BC!(Val{:UO},low_boundary_coefs,up_stencil_length,derivative_order, - grid_step,boundary_length)*BC[1]*dx) - - elseif LBC == :Neumann - return (zero(T),one(T),left_Neumann_BC!(Val{:UO},low_boundary_coefs,up_stencil_length,derivative_order, - grid_step,boundary_length)*BC[1]*dx) - - elseif LBC == :Robin - return (BC[1][1],-BC[1][2],left_Robin_BC!(Val{:UO},low_boundary_coefs,up_stencil_length, - BC[1],derivative_order,grid_step, - boundary_length,dx)*BC[1][3]*dx) - - elseif LBC == :nothing - return (zero(T),zero(T),left_nothing_BC!(Val{:UO},low_boundary_coefs,up_stencil_length,derivative_order, - grid_step,boundary_length)*BC[1]*dx) - - elseif LBC == :Dirichlet0 - return (one(T),zero(T),zero(T)*BC[1]) - - elseif LBC == :Dirichlet - return (one(T),zero(T),one(T)*BC[1]) - - elseif LBC == :Neumann0 - return (zero(T),one(T),zero(T)) - - elseif LBC == :periodic - return (zero(T),zero(T),zero(T)) - - else - error("Unrecognized Boundary Type!") - end -end - -# well it says that we have to use somewhere inside the function definition -function initialize_right_boundary!(::Type{Val{:UO}},high_boundary_coefs,BC,derivative_order,grid_step::T, - boundary_length,boundary_point_count,directions,dx,RBC) where T - approximation_order = boundary_length - derivative_order - down_stencil_length = boundary_length - - if RBC == :None - # down_stencil_length%2 == 1 ? boundary_point_count[2] = 1 : boundary_point_count[2] = 0 - return (zero(T),zero(T),right_None_BC!(Val{:UO},high_boundary_coefs,down_stencil_length,derivative_order, - grid_step,boundary_length)*BC[2]*dx) - - elseif RBC == :Neumann - return (zero(T),one(T),right_Neumann_BC!(Val{:UO},high_boundary_coefs,down_stencil_length,derivative_order, - grid_step,boundary_length)*BC[2]*dx) - - elseif RBC == :Robin - return (BC[2][1],BC[2][2],right_Robin_BC!(Val{:UO},high_boundary_coefs,down_stencil_length, - BC[2],derivative_order,grid_step, - boundary_length,dx)*BC[2][3]*dx) - - elseif RBC == :nothing - return (zero(T),zero(T),right_nothing_BC!(Val{:UO},high_boundary_coefs,down_stencil_length,derivative_order, - grid_step,boundary_length)*BC[2]*dx) - - elseif RBC == :Dirichlet0 - return (one(T),zero(T),zero(T)*BC[2]) - - elseif RBC == :Dirichlet - return (one(T),zero(T),one(T)*BC[2]) - - elseif RBC == :Neumann0 - return (zero(T),one(T),zero(T)) - - elseif RBC == :periodic - return (zero(T),zero(T),zero(T)) - - else - error("Unrecognized Boundary Type!") - end -end - - -function left_None_BC!(::Type{Val{:UO}},low_boundary_coefs,up_stencil_length,derivative_order, - grid_step::T,boundary_length) where T - # Fixes the problem excessive boundary points - boundary_point_count = up_stencil_length - l_diff = zero(T) - - for i in 1:boundary_point_count - # One-sided stencils require more points for same approximation order - # TODO: I don't know if this is the correct stencil length for i > 1? - push!(low_boundary_coefs, negate!(calculate_weights(derivative_order, (i-1)*grid_step, collect(zero(T) : grid_step : (boundary_length-1)*grid_step)))) - end - return l_diff -end - - -function right_None_BC!(::Type{Val{:UO}},high_boundary_coefs,down_stencil_length,derivative_order, - grid_step::T,boundary_length) where T - boundary_point_count = down_stencil_length - high_temp = zeros(T,boundary_length) - aorder = boundary_length - 1 - r_diff = zero(T) - - for i in 1 : boundary_point_count - push!(high_boundary_coefs, negate!(calculate_weights(derivative_order, -(i-1)*grid_step, reverse(collect(zero(T) : -grid_step : -(boundary_length-1)*grid_step))))) - end - return r_diff -end - - -function left_nothing_BC!(::Type{Val{:UO}},low_boundary_coefs,up_stencil_length,derivative_order, - grid_step::T,boundary_length) where T - # Fixes the problem excessive boundary points - boundary_point_count = up_stencil_length - l_diff = zero(T) - - push!(low_boundary_coefs, negate!(calculate_weights(derivative_order, (0)*grid_step, collect(zero(T) : grid_step : (boundary_length-1)*grid_step)))) - return l_diff -end - - -function right_nothing_BC!(::Type{Val{:UO}},high_boundary_coefs,down_stencil_length,derivative_order, - grid_step::T,boundary_length) where T - boundary_point_count = down_stencil_length - r_diff = zero(T) - - push!(high_boundary_coefs, negate!(calculate_weights(derivative_order, -(0)*grid_step, reverse(collect(zero(T) : -grid_step : -(boundary_length-1)*grid_step))))) - return r_diff -end - -#= - The Inf opnorm can be calculated easily using the stencil coeffiicents, while other opnorms - default to compute from the full matrix form. -=# -function LinearAlgebra.opnorm(A::UpwindOperator{T,S,LBC,RBC}, p::Real=2) where {T,S,LBC,RBC} - if p == Inf && LBC in [:Dirichlet0, :Neumann0, :periodic] && RBC in [:Dirichlet0, :Neumann0, :periodic] - max(sum(abs.(A.up_stencil_coefs)) / A.dx^A.derivative_order, sum(abs.(A.down_stencil_coefs)) / A.dx^A.derivative_order) - else - opnorm(convert(Array,A), p) - end -end diff --git a/test/2nd_order_check.jl b/test/2nd_order_check.jl index 38631ad62..7c018f106 100644 --- a/test/2nd_order_check.jl +++ b/test/2nd_order_check.jl @@ -5,8 +5,8 @@ using DiffEqOperators, Test, LinearAlgebra Δx = 0.025π N = Int(2*(π/Δx)) -2 x = -π:Δx:π-Δx - D1 = DerivativeOperator{Float64}(1,ord,Δx,N,:periodic,:periodic); # 1nd Derivative - D2 = DerivativeOperator{Float64}(2,ord,Δx,N,:periodic,:periodic); # 2nd Derivative + D1 = CenteredDifference(1,ord,Δx,N,:periodic,:periodic); # 1nd Derivative + D2 = CenteredDifference(2,ord,Δx,N,:periodic,:periodic); # 2nd Derivative u0 = cos.(x) du_true = -cos.(x) diff --git a/test/KdV.jl b/test/KdV.jl index 9e402bd8f..4130df58c 100644 --- a/test/KdV.jl +++ b/test/KdV.jl @@ -17,10 +17,10 @@ using DiffEqOperators, OrdinaryDiffEq const du3 = zeros(size(x)); const temp = zeros(size(x)); - # A = DerivativeOperator{Float64}(1,2,Δx,length(x),:Dirichlet0,:Dirichlet0); + # A = CenteredDifference(1,2,Δx,length(x),:Dirichlet0,:Dirichlet0); A = UpwindOperator{Float64}(1,3,Δx,length(x),true.|BitVector(undef,length(x)), :Dirichlet0,:Dirichlet0); - # C = DerivativeOperator{Float64}(3,2,Δx,length(x),:Dirichlet0,:Dirichlet0); + # C = CenteredDifference(3,2,Δx,length(x),:Dirichlet0,:Dirichlet0); C = UpwindOperator{Float64}(3,3,Δx,length(x),true.|BitVector(undef,length(x)), :Dirichlet0,:Dirichlet0); diff --git a/test/derivative_operators_interface.jl b/test/derivative_operators_interface.jl index 505da0575..2bb321e60 100644 --- a/test/derivative_operators_interface.jl +++ b/test/derivative_operators_interface.jl @@ -53,7 +53,7 @@ end # Do not modify the following test-set unless you are completely certain of your changes. @testset "Correctness of Stencils" begin N = 20 - L = DerivativeOperator{Float64}(4,4, 1.0, N) + L = CenteredDifference(4,4, 1.0, N) correct = fourth_deriv_approx_stencil(N) # Check that stencils (according to convert_by_multiplication) agree with correct @@ -64,7 +64,7 @@ end @test sparse(L) ≈ correct @test BandedMatrix(L) ≈ correct - L = DerivativeOperator{Float64}(2,4, 1.0, N) + L = CenteredDifference(2,4, 1.0, N) correct = second_deriv_fourth_approx_stencil(N) # Check that stencils (according to convert_by_multiplication) agree with correct @@ -82,7 +82,7 @@ end d_order = 2 approx_order = 2 correct = second_derivative_stencil(N) - A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N) + A = CenteredDifference(d_order,approx_order,1.0,N) @test convert_by_multiplication(Array,A,N) == correct @test Array(A) == second_derivative_stencil(N) @@ -95,7 +95,7 @@ end N = 20 d_order = 4 approx_order = 4 - A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N) + A = CenteredDifference(d_order,approx_order,1.0,N) correct = convert_by_multiplication(Array,A,N) @test Array(A) ≈ correct @@ -105,7 +105,7 @@ end N = 26 d_order = 8 approx_order = 8 - A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N) + A = CenteredDifference(d_order,approx_order,1.0,N) correct = convert_by_multiplication(Array,A,N) @test Array(A) ≈ correct @@ -119,12 +119,11 @@ end y = collect(1:1.0:N+2).^4 - 2*collect(1:1.0:N+2).^3 + collect(1:1.0:N+2).^2; y = convert(Array{BigFloat, 1}, y) - A = DerivativeOperator{BigFloat}(d_order,approx_order,one(BigFloat),N) + A = CenteredDifference(d_order,approx_order,one(BigFloat),N) correct = convert_by_multiplication(Array,A,N) @test Array(A) ≈ correct @test sparse(A) ≈ correct @test BandedMatrix(A) ≈ correct - @test A*y ≈ Array(A)*y end @@ -133,7 +132,7 @@ end d_order = 4 approx_order = 10 - A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N) + A = CenteredDifference(d_order,approx_order,1.0,N) @test A[1,1] == Array(A)[1,1] @test A[10,20] == 0 @@ -147,7 +146,7 @@ end d_order = 2 approx_order = 2 - A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N) + A = CenteredDifference(d_order,approx_order,1.0,N) M = Array(A,1000) @test A[1,1] == M[1,1] @test A[1:4,1] == M[1:4,1] @@ -167,8 +166,8 @@ end dy = yarr[2]-yarr[1] F = [x^2+y for x = xarr, y = yarr] - A = DerivativeOperator{Float64}(d_order,approx_order,dx,length(xarr)-2) - B = DerivativeOperator{Float64}(d_order,approx_order,dy,length(yarr)) + A = CenteredDifference(d_order,approx_order,dx,length(xarr)-2) + B = CenteredDifference(d_order,approx_order,dy,length(yarr)) @test A*F ≈ 2*ones(N-2,M) atol=1e-2 @@ -186,7 +185,7 @@ end # Only tests the additional functionality defined in "operator_combination.jl" N = 10 Random.seed!(0); LA = DiffEqArrayOperator(rand(N,N)) - LD = DerivativeOperator{Float64}(2,2,1.0,N) + LD = CenteredDifference(2,2,1.0,N) @test_broken begin L = 1.1*LA - 2.2*LD + 3.3*I # Builds convert(L) the brute-force way diff --git a/test/generic_operator_validation.jl b/test/generic_operator_validation.jl index 5c8106ecb..93d5a988b 100644 --- a/test/generic_operator_validation.jl +++ b/test/generic_operator_validation.jl @@ -11,7 +11,7 @@ y_ = y[2:(end-1)] @test_broken for dor in 1:6, aor in 1:8 - D1 = DerivativeOperator{Float64}(dor,aor,dx[1],length(x)) + D1 = CenteredDifference(dor,aor,dx[1],length(x)) D2 = DiffEqOperators.FiniteDifference{Float64}(dor,aor,dx,length(x)) D = (D1,D2) @test_broken convert(Array, D1) ≈ convert(Array, D2) diff --git a/test/heat_eqn.jl b/test/heat_eqn.jl index 4ba4903b4..430f84eb9 100644 --- a/test/heat_eqn.jl +++ b/test/heat_eqn.jl @@ -4,7 +4,7 @@ using OrdinaryDiffEq @testset "Parabolic Heat Equation with Dirichlet BCs" begin x = collect(-pi : 2pi/511 : pi) u0 = @. -(x - 0.5).^2 + 1/12 - A = DerivativeOperator{Float64}(2,2,2π/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end])) + A = CenteredDifference(2,2,2π/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end])) heat_eqn = ODEProblem(A, u0, (0.,10.)) soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10) @@ -20,10 +20,10 @@ end dx = 2π/(N-1) x = collect(-pi : dx : pi) u0 = @. -(x - 0.5)^2 + 1/12 - B = DerivativeOperator{Float64}(1,2,dx,N,:None,:None) + B = CenteredDifference(1,2,dx,N,:None,:None) deriv_start, deriv_end = (B*u0)[1], (B*u0)[end] - A = DerivativeOperator{Float64}(2,2,dx,N,:Neumann,:Neumann;BC=(deriv_start,deriv_end)) + A = CenteredDifference(2,2,dx,N,:Neumann,:Neumann;BC=(deriv_start,deriv_end)) heat_eqn = ODEProblem(A, u0, (0.,10.)) soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10) @@ -42,14 +42,14 @@ end dx = 2π/(N-1) x = collect(-pi : dx : pi) u0 = @. -(x - 0.5).^2 + 1/12 - B = DerivativeOperator{Float64}(1,2,dx,N,:None,:None) + B = CenteredDifference(1,2,dx,N,:None,:None) deriv_start, deriv_end = (B*u0)[1], (B*u0)[end] params = [1.0,0.5] left_RBC = params[1]*u0[1] - params[2]*deriv_start right_RBC = params[1]*u0[end] + params[2]*deriv_end - A = DerivativeOperator{Float64}(2,2,dx,N,:Robin,:Dirichlet;BC=((params[1],params[2],left_RBC),u0[end])); + A = CenteredDifference(2,2,dx,N,:Robin,:Dirichlet;BC=((params[1],params[2],left_RBC),u0[end])); heat_eqn = ODEProblem(A, u0, (0.,10.)); soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10); diff --git a/test/runtests.jl b/test/runtests.jl index 8763f06ac..942d7c4c0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,31 +11,3 @@ import Base: isapprox #@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails #@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end @time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end - -#= -using StaticArrays, DiffEqOperators -function isapprox(x::DerivativeOperator{T,S,LBC,RBC},y::FiniteDifference{T,S,LBC,RBC}; kwargs...) where {T<:Real,S<:StaticArrays.SVector,LBC,RBC} - der_order = (x,y) -> x.derivative_order == y.derivative_order - approximation_order = (x,y) -> x.approximation_order == y.approximation_order - dx = (x,y) -> !any(x.dx .!= y.dx) - dimension = (x,y) -> x.dimension == y.dimension - stencil_length = (x,y) -> x.stencil_length == y.stencil_length - stencil_coefs = (x,y) -> !any(!isapprox.(x.stencil_coefs,y.stencil_coefs; kwargs...)) - boundary_point_count= (x,y) -> x.boundary_point_count == y.boundary_point_count - boundary_length = (x,y) -> x.boundary_length == y.boundary_length - low_boundary_coefs = (x,y) -> !any(!isapprox.(x.low_boundary_coefs,y.low_boundary_coefs; kwargs...)) - high_boundary_coefs = (x,y) -> !any(!isapprox.(x.high_boundary_coefs,y.high_boundary_coefs; kwargs...)) - boundary_condition = (x,y) -> !any(!isapprox.(x.high_boundary_coefs,y.high_boundary_coefs; kwargs...)) - t = (x,y) -> x.t ≈ y.t - - #in order of fast to slow comparison - for f in (der_order,approximation_order,dimension,t,boundary_point_count,boundary_length,dx,boundary_condition,low_boundary_coefs,high_boundary_coefs,stencil_coefs) - if !f(x,y) - return false - end - end - return true -end - -isapprox(x::FiniteDifference{T,S,LBC,RBC},y::DerivativeOperator{T,S,LBC,RBC}; kwargs...) where {T<:Real,S<:StaticArrays.SVector,LBC,RBC} = isapprox(y,x; kwargs...) -=#