From 4a7d81d05b6407091b53bb6bb251c0dc154a84a1 Mon Sep 17 00:00:00 2001 From: Xuanzhao Gao Date: Sun, 7 Jan 2024 17:03:14 +0800 Subject: [PATCH] revise --- Manifest.toml | 12 +++---- README.md | 10 +++--- src/ChebParticleMesh.jl | 11 +++--- src/chebpoly.jl | 28 +++++++++++++-- src/grid.jl | 19 ++++++++-- src/index.jl | 24 +++++++++++-- src/interpolate.jl | 79 +++++++++++++++++++++++++++++++++++++++++ src/interpolation.jl | 0 src/types.jl | 5 +++ test/chebpoly.jl | 2 +- test/index.jl | 35 +++++++++++++++++- test/interpolate.jl | 5 +++ test/runtests.jl | 2 ++ 13 files changed, 208 insertions(+), 24 deletions(-) create mode 100644 src/interpolate.jl delete mode 100644 src/interpolation.jl create mode 100644 test/interpolate.jl diff --git a/Manifest.toml b/Manifest.toml index 94e3768..7d3e44d 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -144,9 +144,9 @@ version = "1.15.0" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d" +git-tree-sha1 = "ac67408d9ddf207de5cfa9a97e114352430f01ed" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.18.15" +version = "0.18.16" [[deps.Dates]] deps = ["Printf"] @@ -164,9 +164,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[deps.Distributions]] deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"] -git-tree-sha1 = "9242eec9b7e2e14f9952e8ea1c7e31a50501d587" +git-tree-sha1 = "a4532d110ce91bd744b99280193a317310960c46" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.25.104" +version = "0.25.106" [deps.Distributions.extensions] DistributionsChainRulesCoreExt = "ChainRulesCore" @@ -378,9 +378,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.LoopVectorization]] deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "fc712dd664097440f19a91a704299cca02134ca0" +git-tree-sha1 = "0f5648fbae0d015e3abe5867bca2b362f67a5894" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.168" +version = "0.12.166" [deps.LoopVectorization.extensions] ForwardDiffExt = ["ChainRulesCore", "ForwardDiff"] diff --git a/README.md b/README.md index 311f016..0ec252b 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # ChebParticleMesh -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://Xuanzhao Gao.github.io/ChebParticleMesh.jl/stable/) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://Xuanzhao Gao.github.io/ChebParticleMesh.jl/dev/) -[![Build Status](https://github.com/Xuanzhao Gao/ChebParticleMesh.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/Xuanzhao Gao/ChebParticleMesh.jl/actions/workflows/CI.yml?query=branch%3Amain) -[![Build Status](https://travis-ci.com/Xuanzhao Gao/ChebParticleMesh.jl.svg?branch=main)](https://travis-ci.com/Xuanzhao Gao/ChebParticleMesh.jl) -[![Coverage](https://codecov.io/gh/Xuanzhao Gao/ChebParticleMesh.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/Xuanzhao Gao/ChebParticleMesh.jl) + + +[![Build Status](https://github.com/HPMolSim/ChebParticleMesh.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/HPMolSim/ChebParticleMesh.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Build Status](https://travis-ci.com/HPMolSim/ChebParticleMesh.jl.svg?branch=main)](https://travis-ci.com/HPMolSim/ChebParticleMesh.jl) +[![Coverage](https://codecov.io/gh/HPMolSim/ChebParticleMesh.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/HPMolSim/ChebParticleMesh.jl) diff --git a/src/ChebParticleMesh.jl b/src/ChebParticleMesh.jl index 6262a31..7da8e84 100644 --- a/src/ChebParticleMesh.jl +++ b/src/ChebParticleMesh.jl @@ -1,10 +1,13 @@ module ChebParticleMesh -using LinearAlgebra, ExTinyMD, FastTransforms, SpecialFunctions, LoopVectorization +using LinearAlgebra, ExTinyMD, FastTransforms, SpecialFunctions, LoopVectorization, FFTW export horner -export ChebCoef, pwcheb_eval, pwcheb_eval! -export GridInfo, PadIndex, ImageIndex, id_image2pad, id_pad2image, grid_revise_image!, grid_revise_pad!, grid_image2pad!, grid_pad2image! +export ChebCoef, pwcheb_eval, pwcheb_eval!, f_eval, f_eval! +export GridInfo, PadIndex, ImageIndex +export id_image2pad, id_pad2image, grid_revise_image!, grid_revise_pad!, grid_image2pad!, grid_pad2image! +export nearest_grid_id, pad_grid_id, image_grid_id, image_grid_pos, pad_grid_pos +export interpolate_single!, interpolate! include("types.jl") include("utils.jl") @@ -13,7 +16,7 @@ include("chebpoly.jl") include("grid.jl") include("index.jl") -include("interpolation.jl") +include("interpolate.jl") include("scaling.jl") include("gather.jl") diff --git a/src/chebpoly.jl b/src/chebpoly.jl index 919e7a3..be69815 100644 --- a/src/chebpoly.jl +++ b/src/chebpoly.jl @@ -10,7 +10,7 @@ function ChebCoef(f::Function, h::T, w::Int, ν::Int) where{T<:Union{Float32, Fl """ function ChebCoef(f::Function, h::T, w::Int, ν::Int) where{T<:Union{Float32, Float64}} coef = pwcheb_coef(f, h, w, ν) - return ChebCoef{T}(h, w, ν, coef) + return ChebCoef{T}(f, h, w, ν, coef) end function pwcheb_coef(f::Function, h::T, window_cut::Int, ν::Int) where{T<:Union{Float32, Float64}} @@ -50,7 +50,7 @@ function pwcheb_eval!(x0::T, cheb_approx::Array{T, 1}, chebcoef::ChebCoef{T}) wh h = chebcoef.h C = chebcoef.coef - @assert abs(x0) ≤ h / 2 + @assert (abs(x0) ≤ h / 2) | (abs(x0) ≈ h / 2) @assert size(cheb_approx, 1) == size(C, 1) - 1 if x0 ≥ zero(T) @@ -73,5 +73,29 @@ function pwcheb_eval(x0::T, chebcoef::ChebCoef{T}) where{T<:Union{Float32, Float cheb_approx = zeros(T, size(chebcoef.coef, 1) - 1) cheb_approx = pwcheb_eval!(x0, cheb_approx, chebcoef) + return cheb_approx +end + +function f_eval!(x0::T, cheb_approx::Array{T, 1}, chebcoef::ChebCoef{T}) where{T<:Union{Float32, Float64}} + + h = chebcoef.h + C = chebcoef.coef + w = chebcoef.w + @assert (abs(x0) ≤ h / 2) | (abs(x0) ≈ h / 2) + @assert size(cheb_approx, 1) == size(C, 1) - 1 + + for i in - w : w + x = - x0 + i * h + cheb_approx[i + w + 1] = chebcoef.f(x) + end + + return cheb_approx +end + +function f_eval(x0::T, chebcoef::ChebCoef{T}) where{T<:Union{Float32, Float64}} + + cheb_approx = zeros(T, 2 * chebcoef.w + 1) + cheb_approx = f_eval!(x0, cheb_approx, chebcoef) + return cheb_approx end \ No newline at end of file diff --git a/src/grid.jl b/src/grid.jl index 9fd59fd..f691d8d 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -21,13 +21,28 @@ function GridInfo(N_real::NTuple{3, Int}, w::NTuple{3, Int}, periodicity::NTuple trans_info = TransInfo(N_real, periodicity, image, pad) - return GridInfo{T}(N_real, w, periodicity, image, pad, L, h, N_image, N_pad, trans_info) + k = Vector{Array{T, 1}}() + + for i in 1:3 + ki = [2π * j / (L[i] + h[i] * 2 * pad[i]) for j in 0:(N_pad[i] - 1)] + + for j in 1:N_pad[i] + if j - 1 > ceil(N_pad[i] / 2) + ki[j] -= 2π * N_pad[i] / (L[i] + h[i] * 2 * pad[i]) + end + end + + push!(k, ki) + end + + return GridInfo{T}(N_real, w, periodicity, image, pad, L, h, N_image, N_pad, trans_info, k) end function GridBox(grid_info::GridInfo{T}) where{T<:Union{Float32, Float64}} pad_grid = zeros(T, grid_info.N_pad...) image_grid = zeros(T, grid_info.N_image...) - return GridBox{T}(pad_grid, image_grid) + cheb_value = [zeros(T, 2 * grid_info.w[i] + 1) for i in 1:3] + return GridBox{T}(pad_grid, image_grid, cheb_value) end """ diff --git a/src/index.jl b/src/index.jl index b6b9257..64a1611 100644 --- a/src/index.jl +++ b/src/index.jl @@ -21,17 +21,35 @@ function pad_grid_id(x::Point{3, T}, gridinfo::GridInfo{T}) where{T} return PadIndex(idx, idy, idz) end +function pad_grid_pos(id::PadIndex, gridinfo::GridInfo{T}) where{T} + + x = (id.idx - gridinfo.pad[1] - T(0.5)) * gridinfo.h[1] + y = (id.idy - gridinfo.pad[2] - T(0.5)) * gridinfo.h[2] + z = (id.idz - gridinfo.pad[3] - T(0.5)) * gridinfo.h[3] + + return Point(x, y, z) +end + function image_grid_id(x::Point{3, T}, gridinfo::GridInfo{T}) where{T} idx, idy, idz = nearest_grid_id(x, gridinfo.h) - idx = idx + gridinfo.w[1] + 1 - idy = idy + gridinfo.w[2] + 1 - idz = idz + gridinfo.w[3] + 1 + idx = idx + gridinfo.image[1] + idy = idy + gridinfo.image[2] + idz = idz + gridinfo.image[3] return ImageIndex(idx, idy, idz) end +function image_grid_pos(id::ImageIndex, gridinfo::GridInfo{T}) where{T} + + x = (id.idx - gridinfo.image[1] - T(0.5)) * gridinfo.h[1] + y = (id.idy - gridinfo.image[2] - T(0.5)) * gridinfo.h[2] + z = (id.idz - gridinfo.image[3] - T(0.5)) * gridinfo.h[3] + + return Point(x, y, z) +end + function id_image2pad_single(id_image_i::Int, N_i::T, image_i::Int, pad_i::Int, periodicity::Bool) where{T} id_pad_i = 0 diff --git a/src/interpolate.jl b/src/interpolate.jl new file mode 100644 index 0000000..4945786 --- /dev/null +++ b/src/interpolate.jl @@ -0,0 +1,79 @@ +function interpolate_single!(q::T, pos::Point{3, T}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T} + + cheb_value = gridbox.cheb_value + near_id_image = image_grid_id(pos, gridinfo) + near_pos_image = image_grid_pos(near_id_image, gridinfo) + for i in 1:3 + dx = pos[i] - near_pos_image[i] + pwcheb_eval!(dx, cheb_value[i], chebcoefs[i]) + end + + @turbo for i in - gridinfo.w[1] : gridinfo.w[1] + for j in - gridinfo.w[2] : gridinfo.w[2] + for k in - gridinfo.w[3] : gridinfo.w[3] + gridbox.image_grid[ + near_id_image.idx + i, + near_id_image.idy + j, + near_id_image.idz + k] += + q * + cheb_value[1][i + gridinfo.w[1] + 1] * + cheb_value[2][j + gridinfo.w[2] + 1] * + cheb_value[3][k + gridinfo.w[3] + 1] + end + end + end + + return nothing +end + +function interpolate!(qs::Vector{T}, poses::Vector{Point{3, T}}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T} + + @assert length(qs) == length(poses) + grid_revise_image!(gridbox) + + for i in 1:length(qs) + interpolate_single!(qs[i], poses[i], gridinfo, gridbox, chebcoefs) + end + + return nothing +end + +function interpolate_single_direct!(q::T, pos::Point{3, T}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T} + + cheb_value = gridbox.cheb_value + near_id_image = image_grid_id(pos, gridinfo) + near_pos_image = image_grid_pos(near_id_image, gridinfo) + for i in 1:3 + dx = pos[i] - near_pos_image[i] + f_eval!(dx, cheb_value[i], chebcoefs[i]) + end + + @turbo for i in - gridinfo.w[1] : gridinfo.w[1] + for j in - gridinfo.w[2] : gridinfo.w[2] + for k in - gridinfo.w[3] : gridinfo.w[3] + gridbox.image_grid[ + near_id_image.idx + i, + near_id_image.idy + j, + near_id_image.idz + k] += + q * + cheb_value[1][i + gridinfo.w[1] + 1] * + cheb_value[2][j + gridinfo.w[2] + 1] * + cheb_value[3][k + gridinfo.w[3] + 1] + end + end + end + + return nothing +end + +function interpolate_direct!(qs::Vector{T}, poses::Vector{Point{3, T}}, gridinfo::GridInfo{T}, gridbox::GridBox{T}, chebcoefs::NTuple{3, ChebCoef{T}}) where{T} + + @assert length(qs) == length(poses) + grid_revise_image!(gridbox) + + for i in 1:length(qs) + interpolate_single_direct!(qs[i], poses[i], gridinfo, gridbox, chebcoefs) + end + + return nothing +end \ No newline at end of file diff --git a/src/interpolation.jl b/src/interpolation.jl deleted file mode 100644 index e69de29..0000000 diff --git a/src/types.jl b/src/types.jl index 7dd90b2..3ea46e1 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,4 +1,5 @@ struct ChebCoef{T} + f::Function h::T w::Int ν::Int @@ -18,11 +19,15 @@ struct GridInfo{T} N_pad::NTuple{3, Int} trans_info::Vector{Tuple{NTuple{3, Int}, NTuple{3, Int}, NTuple{3, Int}}} + + k::Vector{Array{T, 1}} end mutable struct GridBox{T} pad_grid::Array{T, 3} image_grid::Array{T, 3} + + cheb_value::Vector{Array{T, 1}} end abstract type AbstractIndex end diff --git a/test/chebpoly.jl b/test/chebpoly.jl index 8790ca5..010da2e 100644 --- a/test/chebpoly.jl +++ b/test/chebpoly.jl @@ -6,6 +6,6 @@ c = ChebCoef(f, h, w, ν) for x0 in [- h / 2 : 0.01 : h / 2...] p = [i * h for i in -w:w] - @test pwcheb_eval(x0, c) ≈ f.(p .- x0) + @test pwcheb_eval(x0, c) ≈ f_eval(x0, c) ≈ f.(p .- x0) end end \ No newline at end of file diff --git a/test/index.jl b/test/index.jl index 5dec588..fed9318 100644 --- a/test/index.jl +++ b/test/index.jl @@ -1,4 +1,37 @@ -@testset "grid transform: image to pad" begin +@testset "pos/id transform" begin + N_real = (20, 30, 40) + w = (12, 23, 3) + extra_pad = (24, 15, 6) + L = (100.0, 100.0, 100.0) + + for periodicity_x in [true, false] + for periodicity_y in [true, false] + for periodicity_z in [true, false] + periodicity = (periodicity_x, periodicity_y, periodicity_z) + @testset "peridoicity is $(periodicity)" begin + gridinfo = GridInfo(N_real, w, periodicity, extra_pad, L) + for _ = 1:10 + pos = Point(L[1] * rand(), L[2] * rand(), L[3] * rand()) + + id_image = image_grid_id(pos, gridinfo) + pos_trans_image = image_grid_pos(id_image, gridinfo) + for i in 1:3 + @test abs(pos[i] - pos_trans_image[i]) ≤ gridinfo.h[i] / 2 + end + + id_pad = pad_grid_id(pos, gridinfo) + pos_trans_pad = pad_grid_pos(id_pad, gridinfo) + for i in 1:3 + @test abs(pos[i] - pos_trans_pad[i]) ≤ gridinfo.h[i] / 2 + end + end + end + end + end + end +end + +@testset "grid transform" begin N_real = (5, 6, 7) w = (1, 3, 4) extra_pad = (4, 5, 6) diff --git a/test/interpolate.jl b/test/interpolate.jl new file mode 100644 index 0000000..b50910c --- /dev/null +++ b/test/interpolate.jl @@ -0,0 +1,5 @@ +@testset "Chebyshev interpolation vs direct interpolation" begin + f = x -> exp(-x^2 * 10.0) + + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index ada096d..fbc2098 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,8 +1,10 @@ using ChebParticleMesh using Test +using ExTinyMD @testset "ChebParticleMesh.jl" begin include("utils.jl") include("chebpoly.jl") include("index.jl") + include("interpolate.jl") end