diff --git a/.github/CI.yml b/.github/CI.yml new file mode 100644 index 0000000..6e5b240 --- /dev/null +++ b/.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/checkout@v1.0.0 + - 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 }} \ No newline at end of file diff --git a/.github/TagBot.yml b/.github/TagBot.yml new file mode 100644 index 0000000..623860f --- /dev/null +++ b/.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 }} \ No newline at end of file diff --git a/Project.toml b/Project.toml index c4ab2fe..b700194 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index f341a54..444c147 100644 --- a/README.md +++ b/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) \ No newline at end of file diff --git a/src/PhysicalFDM.jl b/src/PhysicalFDM.jl index 09598ab..c859c54 100644 --- a/src/PhysicalFDM.jl +++ b/src/PhysicalFDM.jl @@ -5,6 +5,7 @@ using DocStringExtensions using PrecompileTools using SparseArrays +using StaticArrays using OffsetArrays using PaddedViews using Tullio @@ -12,7 +13,8 @@ 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 -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 diff --git a/test/runtests.jl b/test/runtests.jl index bc1a120..61237d1 100644 --- a/test/runtests.jl +++ b/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 \ No newline at end of file