Merge branch 'master' of https://github.com/xtalax/DiffEqOperators.jl

Author: xtalax

appveyor.yml

@@ -1,13 +1,16 @@
 environment:
   matrix:
-    - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/1.0/julia-1.0-latest-win32.exe"
-    - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/1.0/julia-1.0-latest-win64.exe"
-    - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
-    - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
-matrix:
-  allow_failures:
-    - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
-    - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
+  - julia_version: 1
+
+platform:
+  - x64 # 64-bit
+
+# # Uncomment the following lines to allow failures on nightly julia
+# # (tests will run but not make your overall status red)
+# matrix:
+#   allow_failures:
+#   - julia_version: nightly
+
 branches:
   only:
     - master
@@ -20,19 +23,18 @@ notifications:
     on_build_status_changed: false
 
 install:
-  - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12"
-# Download most recent Julia Windows binary
-  - ps: (new-object net.webclient).DownloadFile(
-        $env:JULIA_URL,
-        "C:\projects\julia-binary.exe")
-# Run installer silently, output to C:\projects\julia
-  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
+  - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1"))
 
 build_script:
-# Need to convert from shallow to complete for Pkg.clone to work
-  - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
-      Pkg.clone(pwd(), \"DiffEqOperators\"); Pkg.build(\"DiffEqOperators\")"
+  - echo "%JL_BUILD_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%"
 
 test_script:
-  - C:\projects\julia\bin\julia -e "Pkg.test(\"DiffEqOperators\")"
+  - echo "%JL_TEST_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%"
+
+# # Uncomment to support code coverage upload. Should only be enabled for packages
+# # which would have coverage gaps without running on Windows
+# on_success:
+#   - echo "%JL_CODECOV_SCRIPT%"
+#   - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%"

src/DiffEqOperators.jl

@@ -20,13 +20,16 @@ include("basic_operators.jl")
 include("matrixfree_operators.jl")
 include("jacvec_operators.jl")
 
+### Boundary Operators
+include("derivative_operators/robin_bc_extended.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/boundary_operators.jl")
+include("derivative_operators/convolutions.jl")
 
 ### Composite Operators
 include("composite_operators.jl")
@@ -41,4 +44,5 @@ end
 export MatrixFreeOperator
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
+export RobinBC
 end # module

src/derivative_operators/abstract_operator_functions.jl

@@ -17,9 +17,6 @@ checkbounds(A::AbstractDerivativeOperator, k::Colon, j::Integer) =
 checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
     (0 < k ≤ size(A, 1) || throw(BoundsError(A, (k,size(A,2)))))
 
-
-# BandedMatrix{A,B,C,D}(A::DerivativeOperator{A,B,C,D}) = BandedMatrix(convert(Array, A, A.stencil_length), A.stencil_length, div(A.stencil_length,2), div(A.stencil_length,2))
-
 # ~~ getindex ~~
 @inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
     @boundscheck checkbounds(A, i, j)
@@ -113,21 +110,6 @@ end
     end
 end
 
-#=
-    This the basic mul! function which is broken up into 3 parts ie. handling the left
-    boundary, handling the interior and handling the right boundary. Finally the vector is
-    scaled appropriately according to the derivative order and the degree of discretizaiton.
-    We also update the time stamp of the DerivativeOperator inside to manage time dependent
-    boundary conditions.
-=#
-function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::Union{DerivativeOperator{T},UpwindOperator{T}}, 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)
-    rmul!(x_temp, @.(1/(A.dx^A.derivative_order)))
-    A.t[] += 1 # incrementing the internal time stamp
-end
-
 #=
     This definition of the mul! function makes it possible to apply the LinearOperator on
     a matrix and not just a vector. It basically transforms the rows one at a time.
@@ -159,58 +141,6 @@ Base.transpose(A::Union{DerivativeOperator,UpwindOperator}) = A
 Base.adjoint(A::Union{DerivativeOperator,UpwindOperator}) = A
 LinearAlgebra.issymmetric(::Union{DerivativeOperator,UpwindOperator}) = true
 
-#=
-    This function ideally should give us a matrix which when multiplied with the input vector
-    returns the derivative. But the presence of different boundary conditons complicates the
-    situation. It is not possible to incorporate the full information of the boundary conditions
-    in the matrix as they are additive in nature. So in this implementation we will just
-    return the inner stencil of the matrix transformation. The user will have to manually
-    incorporate the BCs at the ends.
-=#
-function Base.convert(::Type{Array}, A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T
-    @assert N >= A.stencil_length # stencil must be able to fit in the matrix
-    mat = zeros(T, (N, N))
-    v = zeros(T, N)
-    for i=1:N
-        v[i] = one(T)
-        #=
-            calculating the effect on a unit vector to get the matrix of transformation
-            to get the vector in the new vector space.
-        =#
-        mul!(view(mat,:,i), A, v)
-        v[i] = zero(T)
-    end
-    return mat
-end
-
-
-#=
-    This function ideally should give us a matrix which when multiplied with the input vector
-    returns the derivative. But the presence of different boundary conditons complicates the
-    situation. It is not possible to incorporate the full information of the boundary conditions
-    in the matrix as they are additive in nature. So in this implementation we will just
-    return the inner stencil of the matrix transformation. The user will have to manually
-    incorporate the BCs at the ends.
-=#
-function SparseArrays.sparse(A::AbstractDerivativeOperator{T}) where T
-    N = A.dimension
-    mat = spzeros(T, N, N)
-    v = zeros(T, N)
-    row = zeros(T, N)
-    for i=1:N
-        v[i] = one(T)
-        #=
-            calculating the effect on a unit vector to get the matrix of transformation
-            to get the vector in the new vector space.
-        =#
-        mul!(row, A, v)
-        copyto!(view(mat,:,i), row)
-        @. row = 0 * row;
-        v[i] = zero(T)
-    end
-    return mat
-end
-
 #=
     Fallback methods that use the full representation of the operator
 =#
@@ -219,3 +149,51 @@ Base.:\(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A \ convert(Array,
 Base.:\(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = convert(Array,A) \ B
 Base.:/(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A / convert(Array,B)
 Base.:/(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = convert(Array,A) / B
+
+########################################################################
+
+# Are these necessary?
+
+get_type(::AbstractDerivativeOperator{T}) where {T} = T
+
+function *(A::AbstractDerivativeOperator,x::AbstractVector)
+    #=
+        We will output a vector which is a supertype of the types of A and x
+        to ensure numerical stability
+    =#
+    get_type(A) != eltype(x) ? error("DiffEqOperator and array are not of same type!") : nothing
+    y = zeros(promote_type(eltype(A),eltype(x)), length(x))
+    LinearAlgebra.mul!(y, A::AbstractDerivativeOperator, x::AbstractVector)
+    return y
+end
+
+
+function *(A::AbstractDerivativeOperator,M::AbstractMatrix)
+    #=
+        We will output a vector which is a supertype of the types of A and x
+        to ensure numerical stability
+    =#
+    get_type(A) != eltype(M) ? error("DiffEqOperator and array are not of same type!") : nothing
+    y = zeros(promote_type(eltype(A),eltype(M)), size(M))
+    LinearAlgebra.mul!(y, A::AbstractDerivativeOperator, M::AbstractMatrix)
+    return y
+end
+
+
+function *(M::AbstractMatrix,A::AbstractDerivativeOperator)
+    #=
+        We will output a vector which is a supertype of the types of A and x
+        to ensure numerical stability
+    =#
+    get_type(A) != eltype(M) ? error("DiffEqOperator and array are not of same type!") : nothing
+    y = zeros(promote_type(eltype(A),eltype(M)), size(M))
+    LinearAlgebra.mul!(y, A::AbstractDerivativeOperator, M::AbstractMatrix)
+    return y
+end
+
+
+function *(A::AbstractDerivativeOperator,B::AbstractDerivativeOperator)
+    # TODO: it will result in an operator which calculates
+    #       the derivative of order A.dorder + B.dorder of
+    #       approximation_order = min(approx_A, approx_B)
+end