add Readme, CI

.github/CI.yml

@@ -0,0 +1,35 @@
+name: CI
+
+on: [push, pull_request, workflow_dispatch]
+
+# 64-bit Julia only
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    runs-on: ${{ matrix.os }}
+    if: "!contains(github.event.head_commit.message, '[skip ci]') && !contains(github.event.head_commit.message, 'CompatHelper')"
+    strategy:
+      matrix:
+        version:
+          - '1.9.4'
+        os:
+          - ubuntu-latest
+          - macOS-latest
+          - windows-latest
+        arch:
+          - x64
+    steps:
+      - run: echo ACTIONS_RUNNER_DEBUG true
+      - uses: actions/[email protected]
+      - name: "Set up Julia"
+        uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      #- uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-runtest@latest
+        #env: 
+        #  ACTIONS_RUNNER_DEBUG: true
+      - uses: julia-actions/julia-uploadcodecov@latest
+        env:
+          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

.github/TagBot.yml

@@ -0,0 +1,15 @@
+name: TagBot
+on:
+  issue_comment:
+    types:
+      - created
+  workflow_dispatch:
+jobs:
+  TagBot:
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          ssh: ${{ secrets.DOCUMENTER_KEY }}

Project.toml

@@ -8,8 +8,18 @@ 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"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
+
+[compat]
+DocStringExtensions = ">=0.8"
+OffsetArrays = ">=1"
+PaddedViews = ">=0.5"
+PrecompileTools = ">=1"
+SparseArrays = ">=1"
+StaticArrays = ">=1"
+Tullio = ">=0.3"
+julia = ">=1.6"

README.md

@@ -1 +1,138 @@
 # PhysicalFDM.jl
+
+Finite Differencing Method
+
+WARNING: *This package is under development!!!*
+
+## Installation
+
+```julia
+]add PhysicalFDM
+```
+
+or
+
+```julia
+using Pkg; Pkg.add("PhysicalFDM")
+```
+
+or
+
+```julia
+using Pkg; Pkg.add("https://github.com/JuliaAstroSim/PhysicalFDM.jl")
+```
+
+To test the Package:
+```julia
+]test PhysicalFDM
+```
+
+## User Guide
+
+This package is extracted from [AstroNbodySim.jl](https://github.com/JuliaAstroSim/AstroNbodySim.jl). You may find more advanced examples there.
+
+### Differencing
+
+Central differencing scheme is defined as
+$$
+\left(\frac{\partial u}{\partial x}\right)_{i, j}=\frac{1}{2 \Delta x}\left(u_{i+1, j}-u_{i-1, j}\right)+O(\Delta x)
+$$
+Suppose we have an 1D data, and $\delta x = 5$:
+```julia
+julia> d1 = [1,2,1]
+3-element Vector{Int64}:
+ 1
+ 2
+ 1
+
+julia> grad_central(5, d1)
+3-element Vector{Float64}:
+  0.2
+  0.0
+ -0.2
+```
+
+2D example, where $(\nabla_x(d2), \nabla_y(d2))$ is returned as `Tuple`:
+```julia
+julia> d2 = [1 1 1; 1 2 1; 1 1 1]
+3×3 Matrix{Int64}:
+ 1  1  1
+ 1  2  1
+ 1  1  1
+
+julia> grad_central(5, 5, d2)
+([0.1 0.2 0.1; 0.0 0.0 0.0; -0.1 -0.2 -0.1], [0.1 0.0 -0.1; 0.2 0.0 -0.2; 0.1 0.0 -0.1])
+```
+
+3D example
+```julia
+julia> d3 = ones(3,3,3);
+       d3[2,2,2] = 2;
+
+julia> grad_central(1, 1, 1, d3)
+([0.5 0.5 0.5; 0.0 0.0 0.0; -0.5 -0.5 -0.5;;; 0.5 1.0 0.5; 0.0 0.0 0.0; -0.5 -1.0 -0.5;;; 0.5 0.5 0.5; 0.0 0.0 0.0; -0.5 -0.5 -0.5], [0.5 0.0 -0.5; 0.5 0.0 -0.5; 0.5 0.0 -0.5;;; 0.5 0.0 -0.5; 1.0 0.0 -1.0; 0.5 0.0 -0.5;;; 0.5 0.0 -0.5; 0.5 0.0 -0.5; 0.5 0.0 -0.5], [0.5 0.5 0.5; 0.5 1.0 0.5; 0.5 0.5 0.5;;; 0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0;;; -0.5 -0.5 -0.5; -0.5 -1.0 -0.5; -0.5 -0.5 -0.5])
+```
+
+### FDM solver
+
+#### Poisson equation
+$$
+\delta u = f
+$$
+can be discretized to
+$$
+\frac{u_{i+1}-2 u_{i}+u_{i-1}}{\Delta x^{2}}=f_{i}, i = 1, 2, \cdots, N
+$$
+In Dirichlet boundary conditions,
+$$
+\frac{1}{\Delta x^2} \begin{pmatrix}
+        -2 &      1 &      0 & \cdots & \cdots &      0 \\
+        1 &     -2 &      1 &      0 & \cdots & \vdots \\
+        0 &      1 &     -2 &      1 & \cdots & \vdots \\
+    \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\
+    \vdots & \ddots &      1 &     -2 &      1 &      0 \\
+    \vdots & \ddots &      0 &      1 &     -2 &      1 \\
+            0 & \cdots & \cdots &      0 &      1 &     -2
+\end{pmatrix}
+\cdot \begin{pmatrix}
+    u_1 \\ u_2 \\ \vdots \\ \vdots \\ \vdots \\ u_{N-1} \\ u_N
+\end{pmatrix}
+= \begin{pmatrix}
+    f_1 \\ f_2 \\ \vdots \\ \vdots \\ \vdots \\ f_{N-1} \\ f_N
+\end{pmatrix}
+$$
+
+The Poisson problem is converted to solving a matrix equation
+$$
+\mathbf{A} \cdot \mathbf{u} = \mathbf{f}
+$$
+
+`PhysicalFDM.jl` is here to provide user-friendly tools to generate `\mathbf{A}` matrix.
+
+For example, suppose we have a 1D mesh with 3 points, and for the Poisson problem we have a 2nd-order differencing operator:
+```julia
+julia> A = diff_mat(3,2)
+3×3 Matrix{Float64}:
+ -2.0   1.0   0.0
+  1.0  -2.0   1.0
+  0.0   1.0  -2.0
+```
+
+Here's a MWE:
+```julia
+using PhysicalFDM
+f = [0, 1, 0]
+A = diff_mat(3,2)
+u = f \ A
+```
+
+## Package ecosystem
+
+- Basic data structure: [PhysicalParticles.jl](https://github.com/JuliaAstroSim/PhysicalParticles.jl)
+- File I/O: [AstroIO.jl](https://github.com/JuliaAstroSim/AstroIO.jl)
+- Initial Condition: [AstroIC.jl](https://github.com/JuliaAstroSim/AstroIC.jl)
+- Parallelism: [ParallelOperations.jl](https://github.com/JuliaAstroSim/ParallelOperations.jl)
+- Trees: [PhysicalTrees.jl](https://github.com/JuliaAstroSim/PhysicalTrees.jl)
+- Meshes: [PhysicalMeshes.jl](https://github.com/JuliaAstroSim/PhysicalMeshes.jl)
+- Plotting: [AstroPlot.jl](https://github.com/JuliaAstroSim/AstroPlot.jl)
+- Simulation: [ISLENT](https://github.com/JuliaAstroSim/ISLENT)

src/PhysicalFDM.jl

@@ -5,14 +5,16 @@ using DocStringExtensions
 using PrecompileTools
 
 using SparseArrays
+using StaticArrays
 using OffsetArrays
 using PaddedViews
 using Tullio
 
 export diff_mat
 export diff_mat2_x, diff_mat2_y, diff_mat2
 export diff_mat3_x, diff_mat3_y, diff_mat3_z, diff_mat3
-export diff_central_x, diff_central_y, diff_central_z, diff_central
+export diff_central_x, diff_central_y, diff_central_z
+export grad_central
 export laplace_conv
 
 """
@@ -24,7 +26,7 @@ Suppose the data is `[u,v,w]`, we want to smooth `v` by fitting a line (1 order
 
 Using `Rational` type or `Sym` type of `SymPy` can get the exact coefficients. For example, `smooth_coef(Sym[0,1,2],2,1)` get `Sym[-3/2, 2, -1/2]`, which means the first order differential of `[u,v,w]` at `u` is `-1.5u+2v-0.5w`. Using this way can generate all data in <https://en.wikipedia.org/wiki/Finite_difference_coefficient>
 
-Author: Qian, Long. 2021-09 (龙潜，github: longqian95)
+Author: Qian, Long. 2021-09 (github: longqian95)
 
 # Examples
 
@@ -67,7 +69,7 @@ Generate differential matrix.
 - `boundary_order`: The order of the fitted polynomial for points determined by `boundary_points`.
 - `sparse`: If true, return sparse matrix instead of dense one.
 
-Author: Qian, Long. 2021-09 (龙潜，github: longqian95)
+Author: Qian, Long. 2021-09 (github: longqian95)
 
 # Examples
 
@@ -224,7 +226,7 @@ function diff_central_x(Δx, u::AbstractArray{T,2}, pad = zero(T)) where T
 end
 
 function diff_central_y(Δy, u::AbstractArray{T,2}, pad = zero(T)) where T
-    LenX, LenY, LenZ = size(u)
+    LenX, LenY = size(u)
     kernel = diff_central_op(Δy)
     h = div(length(kernel), 2)
     d1 = PaddedView(pad, u, (1:LenX,     1-h:LenY+h))
@@ -259,17 +261,17 @@ function diff_central_z(Δz, u::AbstractArray{T,3}, pad = zero(T)) where T
     return parent(out)
 end
 
-function diff_central(Δx, u::AbstractArray{T,1}, pad = zero(T)) where T
-    diff_central(Δx, u, pad)
+function grad_central(Δx, u::AbstractArray{T,1}, pad = zero(T)) where T
+    diff_central_x(Δx, u, pad)
 end
 
-function diff_central(Δx, Δy, u::AbstractArray{T,2}, pad = zero(T)) where T
+function grad_central(Δx, Δy, u::AbstractArray{T,2}, pad = zero(T)) where T
     dx = diff_central_x(Δx, u, pad)
     dy = diff_central_y(Δy, u, pad)
     return dx, dy
 end
 
-function diff_central(Δx, Δy, Δz, u::AbstractArray{T,3}, pad = zero(T)) where T
+function grad_central(Δx, Δy, Δz, u::AbstractArray{T,3}, pad = zero(T)) where T
     dx = diff_central_x(Δx, u, pad)
     dy = diff_central_y(Δy, u, pad)
     dz = diff_central_z(Δz, u, pad)
@@ -298,12 +300,7 @@ function laplace_conv_op(Δx, Δy, Δz)
      0    0    0], dims = 3))
 end
 
-function laplace_conv(a, Δ...)
-    kernel = laplace_conv_op(Δ...)
-    return imfilter(a, kernel, Fill(0))
-end
-
-
+#TODO: a better way to conv
 
 ### Precompile
 @setup_workload begin

test/runtests.jl

@@ -232,10 +232,18 @@ using PhysicalFDM
     end
 end
 
-@testset "Divergence" begin
-
+d1 = [1,2,1]
+d2 = [1 1 1; 1 2 1; 1 1 1]
+d3 = ones(3,3,3)
+d3[2,2,2] = 2
+
+@testset "Gradience" begin
+    # zero padding by default
+    @test grad_central(1, d1) == [1, 0, -1]
+    @test grad_central(1, 1, d2) == ([0.5 1.0 0.5; 0.0 0.0 0.0; -0.5 -1.0 -0.5], [0.5 0.0 -0.5; 1.0 0.0 -1.0; 0.5 0.0 -0.5])
+    @test grad_central(1, 1, 1, d3) == ([0.5 0.5 0.5; 0.0 0.0 0.0; -0.5 -0.5 -0.5;;; 0.5 1.0 0.5; 0.0 0.0 0.0; -0.5 -1.0 -0.5;;; 0.5 0.5 0.5; 0.0 0.0 0.0; -0.5 -0.5 -0.5], [0.5 0.0 -0.5; 0.5 0.0 -0.5; 0.5 0.0 -0.5;;; 0.5 0.0 -0.5; 1.0 0.0 -1.0; 0.5 0.0 -0.5;;; 0.5 0.0 -0.5; 0.5 0.0 -0.5; 0.5 0.0 -0.5], [0.5 0.5 0.5; 0.5 1.0 0.5; 0.5 0.5 0.5;;; 0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0;;; -0.5 -0.5 -0.5; -0.5 -1.0 -0.5; -0.5 -0.5 -0.5])
 end
 
-@testset "Laplacian" begin
-
+@testset "Laplace" begin
+    @test laplace(d1, 1)
 end