add test

Project.toml

@@ -5,9 +5,11 @@ version = "0.1.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
 PhysicalMeshes = "97d9904f-034f-4fb7-aeaa-03a173434233"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

src/PhysicalFDM.jl

@@ -1,7 +1,8 @@
 module PhysicalFDM
 
+using LinearAlgebra
 using DocStringExtensions
-using PrecompieTools
+using PrecompileTools
 
 using SparseArrays
 using OffsetArrays
@@ -61,7 +62,7 @@ Generate differential matrix.
 - `lpoints`: Number of points at left to the target point, which differential is calculated by the fitted polynomial.
 - `rpoints`: Number of points at right to the target point. If `lpoints==rpoints`, then the differential is estimated as central finite difference. If `lpoints==0`, then it is normal forward finite difference. If `rpoints==0`, then it is backward finite difference.
 - `fitting_order`: The order of the fitted polynomial for estimating differential.
-- `boundary`: Boundary condition. Can be `Dirichlet()`(boundary value is zero), `Periodic`(assume data is periodic), `:Extrapolation`(boundary value is extrapolated according to `boundary_points` and `boundary_order`), `:None`(not deal with boundary, will return non-square matrix).
+- `boundary`: Boundary condition. Can be `:Dirichlet`(boundary value is zero), `:Periodic`(assume data is periodic), `:Extrapolation`(boundary value is extrapolated according to `boundary_points` and `boundary_order`), `:None`(not deal with boundary, will return non-square matrix).
 - `boundary_points`: Number of points for fitting polynomial to estimate differential at boundary. Normally it should not be much less than `points`, otherwise sometimes the current point may not be used to estimate the differential.
 - `boundary_order`: The order of the fitted polynomial for points determined by `boundary_points`.
 - `sparse`: If true, return sparse matrix instead of dense one.
@@ -78,7 +79,7 @@ diff_mat(k,1;points=2,lpoints=0)*x #do the normal 1st order forward differential
 diff_mat(k,1;lpoints=1,rpoints=0)*x #do the 1st order backward differential (x[n-1]-x[n]).
 ```
 """
-function diff_mat(n, order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, lpoints=div(points,2), rpoints=points-lpoints-1, fitting_order=lpoints+rpoints, boundary=Dirichlet(), boundary_points=lpoints+rpoints+1, boundary_order=boundary_points-1, sparse=false)
+function diff_mat(n, order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, lpoints=div(points,2), rpoints=points-lpoints-1, fitting_order=lpoints+rpoints, boundary=:Dirichlet, boundary_points=lpoints+rpoints+1, boundary_order=boundary_points-1, sparse=false)
     n<lpoints+rpoints+1 && throw(ArgumentError("matrix size $n must be greater than or equal to lpoints+rpoints+1 ($lpoints+$rpoints+1)"))
     x=-lpoints:rpoints
     dt=dt^order
@@ -90,23 +91,23 @@ function diff_mat(n, order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, l
         diagm1=diagm
         zeros1=zeros
     end
-    if boundary isa Dirichlet || boundary isa Vacuum # in the vacuum case, we manually compute solution on the boundaries
-        m=diagm1((x[i]=>v[i]*ones(T,n-abs(x[i])) for i=1:length(x))...)
+    if boundary == :Dirichlet || boundary == :Vacuum # in the vacuum case, we manually compute solution on the boundaries
+        m=diagm1((x[i]=>v[i]*ones(T,n-abs(x[i])) for i=eachindex(x))...)
     elseif boundary isa Periodic
         m=zeros1(T,n,n)
-        for i=1:n, j=1:length(x)
+        for i=1:n, j=eachindex(x)
             jj=mod1(x[j]+i,n)
             m[i,jj]=v[j]
         end
     #elseif boundary==:None #TODO
     #    nn=n-lpoints-rpoints
-    #    m=diagm1(nn, n, (x[i]+lpoints=>v[i]*ones(T,nn) for i=1:length(x))...)
+    #    m=diagm1(nn, n, (x[i]+lpoints=>v[i]*ones(T,nn) for i=eachindex(x))...)
     #elseif boundary==:Extrapolation #TODO
     #    m=zeros1(T,n,n)
     #    for i=1:n
     #        if i<=lpoints
     #            b=smooth_coef(T.(1-i:boundary_points-i),boundary_order,order)/dt
-    #            m[i,1:length(b)]=b
+    #            m[i,eachindex(b)]=b
     #        elseif i>n-rpoints
     #            b=smooth_coef(T.(n-i-boundary_points+1:n-i),boundary_order,order)/dt
     #            m[i,n-length(b)+1:n]=b
@@ -123,68 +124,68 @@ diff_vec(order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, lpoints=div(p
 
 
 #generate given order differential matrix for a vector which is expanded from row*col matrix
-function diff_mat2_x(row,col,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=Dirichlet(), sparse=false)
+function diff_mat2_x(row,col,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=:Dirichlet, sparse=false)
     t=diff_mat(col,order; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
     m=kron(t,I(row))
     return m
 end
-function diff_mat2_y(row,col,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=Dirichlet(), sparse=false)
+function diff_mat2_y(row,col,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=:Dirichlet, sparse=false)
     t=diff_mat(row,order; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
     m=kron(I(col),t)
     return m
 end
 
 #2D Δ(Laplacian) operator
-function delta_mat2(row,col; T=Float64, dt=one(T), points=3, boundary=Dirichlet(), sparse=false)
+function delta_mat2(row,col; T=Float64, dt=one(T), points=3, boundary=:Dirichlet, sparse=false)
     return diff_mat2_x(row,col,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)+diff_mat2_y(row,col,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
 end
 
-function delta_mat2(row,col,Δx,Δy; T=Float64, dt=one(T), points=3, boundary=Dirichlet(), sparse=false)
+function delta_mat2(row,col,Δx,Δy; T=Float64, dt=one(T), points=3, boundary=:Dirichlet, sparse=false)
     return diff_mat2_x(row,col,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)/Δx^2+diff_mat2_y(row,col,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)/Δy^2
 end
 
 #generate given order differential matrix for a vector which is expanded from row*col*page tensor
-function diff_mat3_x(row,col,page,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=Dirichlet(), sparse=false)
+function diff_mat3_x(row,col,page,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=:Dirichlet, sparse=false)
     t=diff_mat(col,order; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
     m=kron(I(page),kron(t,I(row)))
     return m
 end
-function diff_mat3_y(row,col,page,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=Dirichlet(), sparse=false)
+function diff_mat3_y(row,col,page,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=:Dirichlet, sparse=false)
     t=diff_mat(row,order; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
     m=kron(I(col*page),t) #or: m=kron(I(page),kron(I(col),t))
     return m
 end
-function diff_mat3_z(row,col,page,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=Dirichlet(), sparse=false)
+function diff_mat3_z(row,col,page,order=1; T=Float64, dt=one(T), points=2*div(order+1,2)+1, boundary=:Dirichlet, sparse=false)
     t=diff_mat(page,order; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
     m=kron(t,I(row*col)) #or: m=kron(kron(t,I(row)),I(col))
     return m
 end
 
 
 #3D Δ(Laplacian) operator
-function delta_mat3(row,col,page; T=Float64, dt=one(T), points=3, boundary=Dirichlet(), sparse=false)
+function delta_mat3(row,col,page; T=Float64, dt=one(T), points=3, boundary=:Dirichlet, sparse=false)
     return diff_mat3_x(row,col,page,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)+diff_mat3_y(row,col,page,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)+diff_mat3_z(row,col,page,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)
 end
 
-function delta_mat3(row,col,page,Δx,Δy,Δz; T=Float64, dt=one(T), points=3, boundary=Dirichlet(), sparse=false)
+function delta_mat3(row,col,page,Δx,Δy,Δz; T=Float64, dt=one(T), points=3, boundary=:Dirichlet, sparse=false)
     return diff_mat3_x(row,col,page,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)/Δx^2+diff_mat3_y(row,col,page,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)/Δy^2+diff_mat3_z(row,col,page,2; T=T,dt=dt,points=points,boundary=boundary,sparse=sparse)/Δz^2
 end
 
-function conv(kernel::AbstractArray{T,1}, d::AbstractArray{T,1}, boundary=Dirichlet(); fill = zero(T)) where T
+function conv(kernel::AbstractArray{T,1}, d::AbstractArray{T,1}, boundary=:Dirichlet; fill = zero(T)) where T
     h = div(length(kernel), 2)
     d1 = PaddedView(fill, d, (1-h:length(d)+h,))
     @tullio out[x] := d1[x-i] * kernel[i]
     return out
 end
 
-function conv(kernel::AbstractArray{T,2}, d::AbstractArray{T,2}, boundary=Dirichlet(); fill = zero(T)) where T
+function conv(kernel::AbstractArray{T,2}, d::AbstractArray{T,2}, boundary=:Dirichlet; fill = zero(T)) where T
     h=div.(size(kernel),2)
     d1=PaddedView(fill,d,(1-h[1]:size(d,1)+h[1],1-h[2]:size(d,2)+h[2]))
     @tullio out[x,y]:=d1[x-i,y-j]*kernel[i,j]
     return parent(out)
 end
 
-function conv(kernel::AbstractArray{T,3}, d::AbstractArray{T,3}, boundary=Dirichlet(); fill = zero(T)) where T
+function conv(kernel::AbstractArray{T,3}, d::AbstractArray{T,3}, boundary=:Dirichlet; fill = zero(T)) where T
     h=div.(size(kernel),2)
     d1=PaddedView(fill,d,size(d).+h.+h,h.+1)
     @tullio out[x,y,z]:=d1[x-i,y-j,z-k]*kernel[i,j,k]
@@ -304,4 +305,12 @@ end
 
 
 
+### Precompile
+@setup_workload begin
+    @compile_workload begin
+        diff_mat(3,1)
+        diff_mat(3,2)
+    end
+end
+
 end # module PhysicalFDM