diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 28e0d04..58c9c80 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -18,7 +18,8 @@ jobs: fail-fast: false matrix: version: - - '1.8' + - '1.6' + - '1.10' - 'nightly' os: - ubuntu-latest @@ -34,3 +35,7 @@ jobs: - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v4 + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} diff --git a/.gitignore b/.gitignore index 99dee2c..39a195d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ .DS_Store Manifest.toml +test/tests.jl diff --git a/Project.toml b/Project.toml index 05b605a..af06939 100644 --- a/Project.toml +++ b/Project.toml @@ -1,22 +1,28 @@ name = "FMM2D" uuid = "2d63477d-9690-4b75-bcc1-c3461d43fecc" authors = ["Mikkel Paltorp Schmitt"] -version = "0.1.0" +version = "0.2.0" [deps] FMM2D_jll = "0fc7e017-fbf9-5352-9b8d-9e37f15ce1cd" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] +FMM2D_jll = "1.1" +Aqua = "0.8" SafeTestsets = "0.1.0" SpecialFunctions = "1.8.1,2" -julia = "1.6,1.7,1.8,1.9" +LinearAlgebra = "1" +Test = "1" +Random = "1" +julia = "1" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [targets] -test = ["Test"] +test = ["Test", "Aqua", "LinearAlgebra", "Random", "SafeTestsets", "SpecialFunctions"] diff --git a/README.md b/README.md index 96bf380..7493c41 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,11 @@ # FMM2D.jl [![Build Status](https://github.com/mipals/FMM2D.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/mipals/FMM2D.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Coverage](https://codecov.io/gh/mipals/FMM2D.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/mipals/FMM2D.jl) FMM2D.jl is a Julia interface for computing N-body interactions using the [Flatiron Institute's FMM2D library](https://github.com/flatironinstitute/fmm2d/). -Currently, the interface only supports Helmholtz problems. +Currently, the wrapper only wraps the Helmholtz and Laplace functionalities. ## Helmholtz FMM Let $c_j \in \mathbb{C},\ j = 1,\dots,N$ denote a collection of charge strengths, $v_j \in \mathbb{C},\ j = 1,\dots,N$ denote a collection of dipole strengths, and $\mathbf{d}_j\in\mathbb{R}^2,\ j = 1,\dots,N$ denote the corresponding dipole orientation vectors. Furthermore, $k \in \mathbb{C}$ denote the wave number. @@ -18,7 +19,7 @@ $$ where $H_0^{(1)}$ is the Hankel function of the first kind of order 0. When $\mathbf{x} = \mathbf{x}_j$ the $j$ th term is dropped from the sum. -## Example +### Example ```julia using FMM2D @@ -26,13 +27,108 @@ using FMM2D thresh = 10.0^(-5) # Tolerance zk = rand(ComplexF64) # Wavenumber -# Source-to-source, +# Source-to-source n = 200 sources = rand(2,n) charges = rand(ComplexF64,n) -vals = hfmm2d(thresh,zk,sources,charges=charges,pg=1) +pg = 3 # Evaluate potential, gradient, and Hessian at the sources +vals = hfmm2d(eps=thresh,zk=zk,sources=sources,charges=charges,pg=pg) +vals.pot +vals.grad +vals.hess + +# Source-to-target +m = 200 +targets = rand(2,m) +pgt = 3 # Evaluate potential, gradient, and Hessian at the targets +vals = hfmm2d(targets=targets,eps=thresh,zk=zk,sources=sources,charges=charges,pgt=pgt) +vals.pottarg +vals.gradtarg +vals.hesstarg +``` + +## Laplace +The Laplace problem in 2D have the following form + +$$ + u(x) = \sum_{j=1}^{N} \left[c_{j} \text{log}\left(\|x-x_{j}\|\right) - d_{j}v_{j} \cdot \nabla( \text{log}(\|x-x_{j}\|) )\right], +$$ + +In the case of complex charges and dipole strengths ($c_j, v_j \in \mathbb{C}^n$) the function call `lfmm2d` has to be used. In the case of real charges and dipole strengths ($c_j, v_j \in \mathbb{R}^n$) the function call `rfmm2d` has to be used. + + +### Example +```julia +using FMM2D + +# Simple example for the FMM2D Library +thresh = 10.0^(-5) # Tolerance + +# Source-to-source +n = 200 +sources = rand(2,n) +charges = rand(ComplexF64,n) +dipvecs = randn(2,n) +dipstr = rand(ComplexF64,n) +pg = 3 # Evaluate potential, gradient, and Hessian at the sources +vals = lfmm2d(eps=thresh,sources=sources,charges=charges,dipvecs=dipvecs,dipstr=dipstr,pg=pg) +vals.pot +vals.grad +vals.hess + +# Source-to-target +m = 100 +targets = rand(2,m) +pgt = 3 # Evaluate potential, gradient, and Hessian at the targets +vals = lfmm2d(targets=targets,eps=thresh,sources=sources,charges=charges,dipvecs=dipvecs,dipstr=dipstr,pgt=pgt) +vals.pottarg +vals.gradtarg +vals.hesstarg ``` + + + + + + + +## Stokes + +$$ + u(x) = \sum_{j=1}^NG^\text{stok}(x,x_j)c_j + d_j\cdot T^\text{stok}(x,x_j)\cdot v_j +$$ + +$$ + p(x) = \sum_{j=1}^NP^\text{stok}(x,x_j)c_j + d_j\cdot \Pi^\text{stok}(x,x_j)\cdot v_j^\top +$$ + ## Related Package -[FMMLIB2D.jl](https://github.com/ludvigak/FMMLIB2D.jl) interfaces the [FMMLIB2D](https://github.com/zgimbutas/fmmlib2d) library which the [FMM2D library improves on](https://fmm2d.readthedocs.io/en/latest/). Due to missing sufficient documentation this package currently only supports Helmholtz problems whereas FMMLIB2D.jl also includes Laplace problems. A future release of this package is expected to also support Laplace problems as well as the other supported kernels from the FMM2D library. +[FMMLIB2D.jl](https://github.com/ludvigak/FMMLIB2D.jl) interfaces the [FMMLIB2D](https://github.com/zgimbutas/fmmlib2d) library which the [FMM2D library improves on](https://fmm2d.readthedocs.io/en/latest/). diff --git a/examples/hfmm2d_example.jl b/examples/hfmm2d_example.jl index b772abd..3ddc36c 100644 --- a/examples/hfmm2d_example.jl +++ b/examples/hfmm2d_example.jl @@ -8,10 +8,10 @@ zk = rand(ComplexF64) # Wavenumber n = 200 sources = rand(2,n) charges = rand(ComplexF64,n) -vals = hfmm2d(thresh,zk,sources,charges=charges,pg=1) +vals = hfmm2d(eps=thresh,zk=zk,sources=sources,charges=charges,pg=1) # Source-to-target ntargets = 300 targets = rand(2,ntargets) -vals = hfmm2d(thresh,zk,sources;charges=charges,targets=targets,pgt=1) +vals = hfmm2d(eps=thresh,zk=zk,sources=sources,charges=charges,targets=targets,pgt=1) vals.pottarg diff --git a/src/FMM2D.jl b/src/FMM2D.jl index e8a35df..acdb1ed 100644 --- a/src/FMM2D.jl +++ b/src/FMM2D.jl @@ -9,8 +9,9 @@ All N-body codes return output in an `FMMVals` structure. See documentation of N-body codes for details. N-body interactions with the Helmholtz kernel -- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code +- [`hfmm3d`](@ref): ``O(N)`` fast mutlipole code for the Helmholtz kernel - [`h3ddir`](@ref): ``O(N^2)`` direct code +- [`lfmm3d`](@ref): ``O(N)`` fast mutlipole code for the Laplace kernel """ module FMM2D @@ -20,6 +21,9 @@ using FMM2D_jll # Exporting interface export hfmm2d +export lfmm2d, rfmm2d, cfmm2d +# export bhfmm +# export stfmm2d # Fortran input/return types Fd = Ref{Float64} @@ -29,6 +33,7 @@ Fc = Ref{ComplexF64} # Common input types TFN = Union{Array{Float64},Nothing} TCN = Union{Array{ComplexF64},Nothing} +TFCN = Union{Array{Float64},Array{ComplexF64},Nothing} # Return struct mutable struct FMMVals @@ -48,5 +53,6 @@ end include("helper_functions.jl") # Include wrappers include("helmholtz_wrappers.jl") +include("laplace_wrappers.jl") end diff --git a/src/helmholtz_wrappers.jl b/src/helmholtz_wrappers.jl index e8d041a..bad0f11 100644 --- a/src/helmholtz_wrappers.jl +++ b/src/helmholtz_wrappers.jl @@ -53,7 +53,7 @@ in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\pa Non-zero values may indicate insufficient memory available. See the documentation for the FMM2D library. If not set (`nothing`), then FMM2D library was never called. """ -function hfmm2d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float64}; +function hfmm2d(;eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float64}, charges::TCN=nothing,dipvecs::TFN=nothing,dipstr::TCN=nothing,targets::TFN=nothing, pg::Integer=0,pgt::Integer=0,nd::Integer=1) @@ -110,9 +110,9 @@ function hfmm2d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float6 if pg == 3 if nd > 1 - hess = zeros(ComplexF64,nd,4,n) + hess = zeros(ComplexF64,nd,3,n) else - hess = zeros(ComplexF64,4,n) + hess = zeros(ComplexF64,3,n) end end @@ -134,9 +134,9 @@ function hfmm2d(eps::Float64,zk::Union{Float64,ComplexF64},sources::Array{Float6 if pgt == 3 if nd > 1 - hesstarg = zeros(ComplexF64,nd,4,nt) + hesstarg = zeros(ComplexF64,nd,3,nt) else - hesstarg = zeros(ComplexF64,4,nt) + hesstarg = zeros(ComplexF64,3,nt) end end @@ -171,8 +171,7 @@ end """ ```julia - vals = h2ddir(zk,sources,targets; - charges=nothing,dipstr=nothing,dipvecs=nothing,pgt=0,nd=1,thresh=1e-16) + vals = h2ddir(zk,sources,targets, charges=nothing,dipstr=nothing,dipvecs=nothing,pgt=0,nd=1,thresh=1e-16) ``` This function computes the N-body Helmholtz interactions in two dimensions where the interaction kernel is given by ``H_0^{(1)}(kr)`` and its gradients. This is the @@ -215,7 +214,7 @@ in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\pa * `vals.gradtarg::Array{ComplexF64}` size (nd,2,nt) or (2,nt) gradient at target locations if requested * `vals.hesstarg::Array{ComplexF64}` size (nd,4,nt) or (4,nt) Hessian at target locations if requested """ -function h2ddir(zk::Union{ComplexF64,Float64},sources::Array{Float64}, targets::Array{Float64}; +function h2ddir(;zk::Union{ComplexF64,Float64},sources::Array{Float64}, targets::Array{Float64}, charges::TCN=nothing,dipvecs::TFN=nothing,dipstr::TCN=nothing, pgt::Integer=0,nd::Integer=1,thresh::Float64=1e-15) @@ -271,9 +270,9 @@ function h2ddir(zk::Union{ComplexF64,Float64},sources::Array{Float64}, targets:: if pgt == 3 if nd > 1 - hesstarg = zeros(ComplexF64,nd,4,nt) + hesstarg = zeros(ComplexF64,nd,3,nt) else - hesstarg = zeros(ComplexF64,4,nt) + hesstarg = zeros(ComplexF64,3,nt) end end diff --git a/src/laplace_wrappers.jl b/src/laplace_wrappers.jl index e69de29..aa6266e 100644 --- a/src/laplace_wrappers.jl +++ b/src/laplace_wrappers.jl @@ -0,0 +1,507 @@ +""" +```julia + vals = lfmm2d(eps,zk,sources; + charges=nothing,dipstr=nothing,dipvecs=nothing,targets=nothing,pg=0,pgt=0,nd=1) +``` +This subroutine computes the N-body Laplace interactions in two dimensions where the +interaction kernel is given by log(r) and its gradients. +```math + u(x) = \\sum_{j=1}^{N} c_{j} * log\\(\\|x-x_{j}\\|\\) - d_{j}v_{j} \\cdot \\nabla( log(\\|x-x_{j}\\|) ), +``` + +where ``c_{j}`` are the charge densities, + ``d_{j}`` are the dipole strength vectors + ``v_{j}`` are the dipole orientation vectors, + ``x_{j}`` are the source locations. +when ``x=x_{j}``, the jth term is dropped from the sum + +# Input +* `eps::Float64` precision requested +* `sources::Array{Float64}` size (2,n) source locations (``x_{j}``) + +# Optional Keyword Input +* `charges::Array{ComplexF64}` size (nd,n) +* `dipstr::Array{ComplexF64}` size (nd,n) +* `dipvecs::Array{Float64}` size (nd,2,n) or (2,n) dipole orientation vectors (``d_{j}``) +* `targets::Array{Float64}` size (2,nt) target locations (``x``) +* `pg::Integer` source eval flag. + + Do not compute any quantities at sources if `pg == 0` + + Potential (``u``) at sources evaluated if `pg == 1`. + + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2` + + Potential, gradient, and Hessian (``\\nabla \\nabla u``) at sources evaluated if `pg == 3` +* `pgt::Integer` target eval flag. + + Do not compute any quantities at targets if `pgt == 0` + + Potential at targets evaluated if `pgt == 1`. + + Potential and gradient at targets evaluated if `pgt == 2` + + Potential, gradient, and Hessian at targets evaluated if `pgt == 3` +* `nd::Integer` number of densities +Note: if all default values are used for optional input, nothing is computed. + +# Output +`vals::FMMVals` with the fields +Note that the Hessian is returned as a length 4 vector at each point with the second derivatives +in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\partial_{yx}``. +* `vals.pot::Array{ComplexF64}` size (nd,n) or (n) potential at source locations if requested +* `vals.grad::Array{ComplexF64}` size (nd,2,n) or (2,n) gradient at source locations if requested +* `vals.hess::Array{ComplexF64}` size (nd,3,n) or (3,n) Hessian at source locations if requested +* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested +* `vals.gradtarg::Array{ComplexF64}` size (nd,2,nt) or (2,nt) gradient at target locations if requested +* `vals.hesstarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) Hessian at target locations if requested +* `vals.ier` error flag as returned by FMM2D library. A value of 0 indicates a successful call. +Non-zero values may indicate insufficient memory available. See the documentation for the FMM2D library. +If not set (`nothing`), then FMM2D library was never called. +""" +function lfmm2d(;eps::Float64,sources::Array{Float64}, + charges::TCN=nothing,dipvecs::TFN=nothing,dipstr::TCN=nothing,targets::TFN=nothing, + pg::Integer=0,pgt::Integer=0,nd::Integer=1) + + + # default values + vals = FMMVals() + + zero = ComplexF64(0) + + # + pot = zero + grad = zero + hess = zero + pottarg = zero + gradtarg = zero + hesstarg = zero + + # check inputs + anyfail, n, nt, ifcharge, ifdipole = ( + scalarfmm2dinputcheck(sources,charges,dipvecs,dipstr,targets,pg,pgt,nd)) + + if anyfail + return vals + end + + if dipstr === nothing && dipvecs !== nothing + dipstr = ones(eltype(charges),n) + end + + if (ifcharge == 0); charges = zero end + if (ifdipole == 0); dipvecs = 0.0 end + if (ifdipole == 0); dipstr = zero end + if (nt == 0); targets = 0.0 end + + # allocate memory for return values + + if pg == 1 || pg == 2 || pg == 3 + if nd > 1 + pot = zeros(ComplexF64,nd,n) + else + pot = zeros(ComplexF64,n) + end + end + + if pg == 2 || pg == 3 + if nd > 1 + grad = zeros(ComplexF64,nd,2,n) + else + grad = zeros(ComplexF64,2,n) + end + end + + if pg == 3 + if nd > 1 + hess = zeros(ComplexF64,nd,3,n) + else + hess = zeros(ComplexF64,3,n) + end + end + + if pgt == 1 || pgt == 2 || pgt == 3 + if nd > 1 + pottarg = zeros(ComplexF64,nd,nt) + else + pottarg = zeros(ComplexF64,nt) + end + end + + if pgt == 2 || pgt == 3 + if nd > 1 + gradtarg = zeros(ComplexF64,nd,2,nt) + else + gradtarg = zeros(ComplexF64,2,nt) + end + end + + if pgt == 3 + if nd > 1 + hesstarg = zeros(ComplexF64,nd,3,nt) + else + hesstarg = zeros(ComplexF64,3,nt) + end + end + + + ier = Integer(0) + iper = Integer(0) + + # Calling the Fortran function + # subroutine lfmm2d(nd,eps,ns,sources,ifcharge,charge, + # ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess, + # nt,targ,ifpghtarg,pottarg,gradtarg, + # hesstarg,ier) + ccall((:lfmm2d_,libfmm2d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fc,Fi, + Fc,Fd,Fi,Fi,Fc,Fc,Fc,Fi, + Fd,Fi,Fc,Fc,Fc,Fi), + nd,eps,n,sources,ifcharge,charges,ifdipole, + dipstr,dipvecs,iper,pg,pot,grad,hess,nt,targets,pgt, + pottarg,gradtarg,hesstarg,ier) + + if (vals.ier != 0); @warn "libfmm2d had an error, see vals.ier"; end + + # Save computed values to output struct + if pg == 1 || pg == 2 || pg == 3; vals.pot = pot end + if pg == 2 || pg == 3; vals.grad = grad end + if pg == 3; vals.hess = hess end + if pgt == 1 || pgt == 2 || pgt == 3; vals.pottarg = pottarg end + if pgt == 2 || pgt == 3; vals.gradtarg = gradtarg end + if pgt == 3; vals.hesstarg = hesstarg end + + return vals + +end + +""" +```julia + vals = rfmm2d(eps,zk,sources; + charges=nothing,dipstr=nothing,dipvecs=nothing,targets=nothing,pg=0,pgt=0,nd=1) +``` +This subroutine computes the N-body Laplace interactions in two dimensions where the +interaction kernel is given by log(r) and its gradients. +```math + u(x) = \\sum_{j=1}^{N} c_{j} * log\\(\\|x-x_{j}\\|\\) - d_{j}v_{j} \\cdot \\nabla( log(\\|x-x_{j}\\|) ), +``` + +where ``c_{j}`` are the charge densities, + ``d_{j}`` are the dipole strength vectors + ``v_{j}`` are the dipole orientation vectors, + ``x_{j}`` are the source locations. +when ``x=x_{j}``, the jth term is dropped from the sum + +# Input +* `eps::Float64` precision requested +* `sources::Array{Float64}` size (2,n) source locations (``x_{j}``) + +# Optional Keyword Input +* `charges::Array{Float64}` size (nd,n) +* `dipstr::Array{Float64}` size (nd,n) +* `dipvecs::Array{Float64}` size (nd,2,n) or (2,n) dipole orientation vectors (``d_{j}``) +* `targets::Array{Float64}` size (2,nt) target locations (``x``) +* `pg::Integer` source eval flag. + + Do not compute any quantities at sources if `pg == 0` + + Potential (``u``) at sources evaluated if `pg == 1`. + + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2` + + Potential, gradient, and Hessian (``\\nabla \\nabla u``) at sources evaluated if `pg == 3` +* `pgt::Integer` target eval flag. + + Do not compute any quantities at targets if `pgt == 0` + + Potential at targets evaluated if `pgt == 1`. + + Potential and gradient at targets evaluated if `pgt == 2` + + Potential, gradient, and Hessian at targets evaluated if `pgt == 3` +* `nd::Integer` number of densities +Note: if all default values are used for optional input, nothing is computed. + +# Output +`vals::FMMVals` with the fields +Note that the Hessian is returned as a length 4 vector at each point with the second derivatives +in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\partial_{yx}``. +* `vals.pot::Array{ComplexF64}` size (nd,n) or (n) potential at source locations if requested +* `vals.grad::Array{ComplexF64}` size (nd,2,n) or (2,n) gradient at source locations if requested +* `vals.hess::Array{ComplexF64}` size (nd,3,n) or (3,n) Hessian at source locations if requested +* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested +* `vals.gradtarg::Array{ComplexF64}` size (nd,2,nt) or (2,nt) gradient at target locations if requested +* `vals.hesstarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) Hessian at target locations if requested +* `vals.ier` error flag as returned by FMM2D library. A value of 0 indicates a successful call. +Non-zero values may indicate insufficient memory available. See the documentation for the FMM2D library. +If not set (`nothing`), then FMM2D library was never called. +""" +function rfmm2d(;eps::Float64,sources::Array{Float64}, + charges::TFN=nothing,dipvecs::TFN=nothing,dipstr::TFN=nothing,targets::TFN=nothing, + pg::Integer=0,pgt::Integer=0,nd::Integer=1) + + + # default values + vals = FMMVals() + + zero = Float64(0) + + # + pot = zero + grad = zero + hess = zero + pottarg = zero + gradtarg = zero + hesstarg = zero + + # check inputs + anyfail, n, nt, ifcharge, ifdipole = ( + scalarfmm2dinputcheck(sources,charges,dipvecs,dipstr,targets,pg,pgt,nd)) + + if anyfail + return vals + end + + if dipstr === nothing && dipvecs !== nothing + dipstr = ones(n) + end + + if (ifcharge == 0); charges = zero end + if (ifdipole == 0); dipvecs = 0.0 end + if (ifdipole == 0); dipstr = zero end + if (nt == 0); targets = 0.0 end + + # allocate memory for return values + + if pg == 1 || pg == 2 || pg == 3 + if nd > 1 + pot = zeros(nd,n) + else + pot = zeros(n) + end + end + + if pg == 2 || pg == 3 + if nd > 1 + grad = zeros(nd,2,n) + else + grad = zeros(2,n) + end + end + + if pg == 3 + if nd > 1 + hess = zeros(nd,3,n) + else + hess = zeros(3,n) + end + end + + if pgt == 1 || pgt == 2 || pgt == 3 + if nd > 1 + pottarg = zeros(nd,nt) + else + pottarg = zeros(nt) + end + end + + if pgt == 2 || pgt == 3 + if nd > 1 + gradtarg = zeros(nd,2,nt) + else + gradtarg = zeros(2,nt) + end + end + + if pgt == 3 + if nd > 1 + hesstarg = zeros(nd,3,nt) + else + hesstarg = zeros(3,nt) + end + end + + + ier = Integer(0) + iper = Integer(0) + + # Calling the Fortran function + # subroutine rfmm2d(nd,eps,ns,sources,ifcharge,charge, + # ifdipole,dipstr,dipvec,iper,ifpgh,pot,grad,hess, + # nt,targ,ifpghtarg,pottarg,gradtarg, + # hesstarg,ier) + ccall((:rfmm2d_,libfmm2d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fd,Fi, + Fd,Fd,Fi,Fi,Fd,Fd,Fd,Fi, + Fd,Fi,Fd,Fd,Fd,Fi), + nd,eps,n,sources,ifcharge,charges,ifdipole, + dipstr,dipvecs,iper,pg,pot,grad,hess,nt,targets,pgt, + pottarg,gradtarg,hesstarg,ier) + + if (vals.ier != 0); @warn "libfmm2d had an error, see vals.ier"; end + + # Save computed values to output struct + if pg == 1 || pg == 2 || pg == 3; vals.pot = pot end + if pg == 2 || pg == 3; vals.grad = grad end + if pg == 3; vals.hess = hess end + if pgt == 1 || pgt == 2 || pgt == 3; vals.pottarg = pottarg end + if pgt == 2 || pgt == 3; vals.gradtarg = gradtarg end + if pgt == 3; vals.hesstarg = hesstarg end + + return vals + +end + +""" +```julia + vals = cfmm2d(eps,zk,sources; + charges=nothing,dipstr=nothing,targets=nothing,pg=0,pgt=0,nd=1) +``` +This subroutine computes the N-body Laplace interactions in two dimensions where the +interaction kernel is given by log(r) and its gradients. +```math + u(z) = \\sum_{j=1}^{N} c_{j} * log\\(\\|z - \\epsilon_j\\|\\) - \\frac{v_j}{z - \\epsilon_j}, +``` + +where ``c_{j}`` are the charge strengths, + ``v_{j}`` are the dipole strengths, + ``z_{j} = x_{j}(1) + ix_{j}(2)`` are the complex number corresponding to ``x_{j}``. +when ``z=\\epsilon_{j}``, the jth term is dropped from the sum + +# Input +* `eps::Float64` precision requested +* `sources::Array{Float64}` size (2,n) source locations (``x_{j}``) + +# Optional Keyword Input +* `charges::Array{ComplexF64}` size (nd,n) +* `dipstr::Array{ComplexF64}` size (nd,n) +* `targets::Array{Float64}` size (2,nt) target locations (``x``) +* `pg::Integer` source eval flag. + + Do not compute any quantities at sources if `pg == 0` + + Potential (``u``) at sources evaluated if `pg == 1`. + + Potential and gradient (``\\nabla u``) at sources evaluated if `pg == 2` + + Potential, gradient, and Hessian (``\\nabla \\nabla u``) at sources evaluated if `pg == 3` +* `pgt::Integer` target eval flag. + + Do not compute any quantities at targets if `pgt == 0` + + Potential at targets evaluated if `pgt == 1`. + + Potential and gradient at targets evaluated if `pgt == 2` + + Potential, gradient, and Hessian at targets evaluated if `pgt == 3` +* `nd::Integer` number of densities +Note: if all default values are used for optional input, nothing is computed. + +# Output +`vals::FMMVals` with the fields +Note that the Hessian is returned as a length 4 vector at each point with the second derivatives +in the order: ``\\partial_{xx}``, ``\\partial_{yy}``, ``\\partial_{xy}``, ``\\partial_{yx}``. +* `vals.pot::Array{ComplexF64}` size (nd,n) or (n) potential at source locations if requested +* `vals.grad::Array{ComplexF64}` size (nd,2,n) or (2,n) gradient at source locations if requested +* `vals.hess::Array{ComplexF64}` size (nd,3,n) or (3,n) Hessian at source locations if requested +* `vals.pottarg::Array{ComplexF64}` size (nd,nt) or (nt) potential at target locations if requested +* `vals.gradtarg::Array{ComplexF64}` size (nd,2,nt) or (2,nt) gradient at target locations if requested +* `vals.hesstarg::Array{ComplexF64}` size (nd,3,nt) or (3,nt) Hessian at target locations if requested +* `vals.ier` error flag as returned by FMM2D library. A value of 0 indicates a successful call. +Non-zero values may indicate insufficient memory available. See the documentation for the FMM2D library. +If not set (`nothing`), then FMM2D library was never called. +""" +function cfmm2d(;eps::Float64,sources::Array{Float64}, + charges::TCN=nothing,dipstr::TCN=nothing,targets::TFN=nothing, + pg::Integer=0,pgt::Integer=0,nd::Integer=1) + + + # default values + vals = FMMVals() + + zero = Float64(0) + + # + pot = zero + grad = zero + hess = zero + pottarg = zero + gradtarg = zero + hesstarg = zero + + # check inputs + dipvecs = nothing # No dipvectors here + anyfail, n, nt, ifcharge, ifdipole = ( + scalarfmm2dinputcheck(sources,charges,dipvecs,dipstr,targets,pg,pgt,nd)) + + # Overwriting ifdipole + if dipstr !== nothing + ifdipole = 1 + end + + # println("$(ifdipole)") + + if anyfail + return vals + end + + + if (ifcharge == 0); charges = zero end + if (ifdipole == 0); dipvecs = 0.0 end + if (ifdipole == 0); dipstr = zero end + if (nt == 0); targets = 0.0 end + + # allocate memory for return values + + if pg == 1 || pg == 2 || pg == 3 + if nd > 1 + pot = zeros(ComplexF64,nd,n) + else + pot = zeros(ComplexF64,n) + end + end + + if pg == 2 || pg == 3 + if nd > 1 + grad = zeros(ComplexF64,nd,n) + else + grad = zeros(ComplexF64,n) + end + end + + if pg == 3 + if nd > 1 + hess = zeros(ComplexF64,nd,n) + else + hess = zeros(ComplexF64,n) + end + end + + if pgt == 1 || pgt == 2 || pgt == 3 + if nd > 1 + pottarg = zeros(ComplexF64,nd,nt) + else + pottarg = zeros(ComplexF64,nt) + end + end + + if pgt == 2 || pgt == 3 + if nd > 1 + gradtarg = zeros(ComplexF64,nd,nt) + else + gradtarg = zeros(ComplexF64,nt) + end + end + + if pgt == 3 + if nd > 1 + hesstarg = zeros(ComplexF64,nd,nt) + else + hesstarg = zeros(ComplexF64,nt) + end + end + + ier = Integer(0) + iper = Integer(0) + + # Calling the Fortran function + # subroutine cfmm3d(nd,eps,ns,sources,ifcharge,charge, + # ifdipole,dipstr,iper,ifpgh,pot,grad,hess, + # nt,targ,ifpghtarg,pottarg,gradtarg, + # hesstarg,ier) + ccall((:cfmm2d_,libfmm2d),Cvoid,(Fi,Fd,Fi,Fd,Fi,Fc, + Fi,Fc,Fi,Fi,Fc,Fc,Fc, + Fi,Fd,Fi,Fc,Fc, + Fc,Fi), + nd,eps,n,sources,ifcharge,charges, + ifdipole,dipstr,iper,pg,pot,grad,hess, + nt,targets,pgt,pottarg,gradtarg, + hesstarg,ier) + + if (vals.ier != 0); @warn "libfmm2d had an error, see vals.ier"; end + + # Save computed values to output struct + if pg == 1 || pg == 2 || pg == 3; vals.pot = pot end + if pg == 2 || pg == 3; vals.grad = grad end + if pg == 3; vals.hess = hess end + if pgt == 1 || pgt == 2 || pgt == 3; vals.pottarg = pottarg end + if pgt == 2 || pgt == 3; vals.gradtarg = gradtarg end + if pgt == 3; vals.hesstarg = hesstarg end + + return vals + +end diff --git a/test/laplace.jl b/test/laplace.jl deleted file mode 100644 index e69de29..0000000 diff --git a/test/runtests.jl b/test/runtests.jl index 5c99e9e..0275c7c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,6 @@ using SafeTestsets -@safetestset "Helmholtz Wrappers" begin include("helmholtz.jl") end +@safetestset "Aqua testing " begin include("test_aqua.jl") end +@safetestset "Laplace Wrappers " begin include("test_laplace.jl") end +@safetestset "Helmholtz Wrappers" begin include("test_helmholtz.jl") end +# @safetestset "Stokes Wrappers " begin include("test_stokes.jl") end diff --git a/test/test_aqua.jl b/test/test_aqua.jl new file mode 100644 index 0000000..95e43ee --- /dev/null +++ b/test/test_aqua.jl @@ -0,0 +1,12 @@ +using Test +using Aqua +using FMM2D + +Aqua.test_undefined_exports(FMM2D) +Aqua.test_project_extras(FMM2D) +Aqua.test_unbound_args(FMM2D) +Aqua.test_ambiguities(FMM2D) +Aqua.test_deps_compat(FMM2D) +Aqua.test_stale_deps(FMM2D) +Aqua.test_piracies(FMM2D) +Aqua.test_persistent_tasks(FMM2D) diff --git a/test/test_helmholtz.jl b/test/test_helmholtz.jl new file mode 100644 index 0000000..6d9031d --- /dev/null +++ b/test/test_helmholtz.jl @@ -0,0 +1,41 @@ +using FMM2D +using SpecialFunctions +using LinearAlgebra +using Random +using Test + +Random.seed!(0) + +@testset "Helmholtz FMM 2D (complex)" begin + zk = rand(ComplexF64) + N = 1000 + sources = rand(2, N) + charges = rand(N) + 1im*rand(N) + dipstrs = rand(N) + 1im*rand(N) + dipvecs = ones(2, N)/sqrt(2) + # Direct computations (N^2) + refpot = FMM2D.h2ddir(zk=zk,sources=sources,targets=sources,dipvecs=dipvecs,dipstr=dipstrs,charges=charges,pgt=3) + + M = 100 + targets = rand(2, M) + # Direct computations (N^2) + refpottarg = FMM2D.h2ddir(zk=zk,sources=sources,targets=targets,dipvecs=dipvecs,dipstr=dipstrs,charges=charges,pgt=3) + for tol_exp=-14:-1 + tolerance = 0.49*10.0^tol_exp + U = hfmm2d(eps=tolerance,zk=zk,sources=sources, charges=charges, dipstr=dipstrs, dipvecs=dipvecs,pg=3) + pot_relerr = norm(U.pot - refpot.pottarg) / norm(refpot.pottarg) + grad_relerr = norm(U.grad - refpot.gradtarg) / norm(refpot.gradtarg) + hess_relerr = norm(U.hess - refpot.hesstarg) / norm(refpot.hesstarg) + @test pot_relerr < tolerance + @test grad_relerr < tolerance + @test hess_relerr < tolerance + + Utarg = hfmm2d(eps=tolerance,zk=zk,sources=sources, charges=charges, dipstr=dipstrs, dipvecs=dipvecs,pgt=3,targets=targets) + pot_relerrtarg = norm(Utarg.pottarg - refpottarg.pottarg) / norm(refpottarg.pottarg) + grad_relerrtarg = norm(Utarg.gradtarg - refpottarg.gradtarg) / norm(refpottarg.gradtarg) + hess_relerrtarg = norm(Utarg.hesstarg - refpottarg.hesstarg) / norm(refpottarg.hesstarg) + @test pot_relerrtarg < tolerance + @test grad_relerrtarg < tolerance + @test hess_relerrtarg < tolerance + end +end diff --git a/test/test_laplace.jl b/test/test_laplace.jl new file mode 100644 index 0000000..748a0ea --- /dev/null +++ b/test/test_laplace.jl @@ -0,0 +1,240 @@ +using FMM2D +using LinearAlgebra +using Random +using Test + +Random.seed!(0) + +function direct_self(source, charges, dipvecs, dipstrs) + if size(charges,2) == 1 + charges = transpose(charges) + end + if size(dipstrs,2) == 1 + dipstrs = transpose(dipstrs) + end + nd,n = size(charges) + pot = zeros(eltype(charges),nd,n) + N = size(source, 2) + for k=1:nd + charge = charges[k,:] + dipvec = dipvecs[2k-1:2k,:] + dipstr = dipstrs[k,:] + for i=1:N + for j=1:N + if i == j + continue + end + r1 = source[1,i] - source[1,j] + r2 = source[2,i] - source[2,j] + rsq = r1^2 + r2^2 + rdotq = r1*dipvec[1,j] + r2*dipvec[2,j] + # Divide by two since log(sqrt(r)) = log(r)/2 + pot[k,i] += charge[j]*log(rsq)/2 - dipstr[j]*rdotq/rsq + end + end + end + if nd == 1 + return transpose(pot) + else + return pot + end +end +function direct_targ(source, charges, dipvecs, dipstrs, target) + if size(charges,2) == 1 + charges = transpose(charges) + end + if size(dipstrs,2) == 1 + dipstrs = transpose(dipstrs) + end + nd = size(charges,1) + M = size(target, 2) + N = size(source, 2) + pot = zeros(eltype(charges), nd, M) + for k in 1:nd + charge = charges[k,:] + dipvec = dipvecs[2k-1:2k,:] + dipstr = dipstrs[k,:] + for i=1:M + for j=1:N + r1 = target[1,i] - source[1,j] + r2 = target[2,i] - source[2,j] + rsq = r1^2 + r2^2 + rdotq = r1*dipvec[1,j] + r2*dipvec[2,j] + # Divide by two since log(sqrt(r)) = log(r)/2 + pot[k,i] += charge[j]*log(rsq)/2 - dipstr[j]*rdotq/rsq + end + end + end + if nd == 1 + return transpose(pot) + else + return pot + end +end + +@testset "Laplace FMM 2D (complex)" begin + N = 1000 + sources = rand(2, N) + charges = rand(ComplexF64,N) + dipstrs = rand(ComplexF64,N) + dipvecs = ones(2, N)/sqrt(2) + # Direct computations (N^2) + refpot = direct_self(sources,charges,dipvecs,dipstrs) + + M = 100 + targets = rand(2, M) + # Direct computations (N^2) + refpottarg = direct_targ(sources,charges,dipvecs,dipstrs,targets) + for tol_exp=-14:-1 + tolerance = 0.49*10.0^tol_exp + U = lfmm2d(eps=tolerance,sources=sources,charges=charges, dipstr=dipstrs, dipvecs=dipvecs,pg=2) + pot_relerr = norm(U.pot - refpot) / norm(refpot) + # grad_relerr = norm(U.grad - refpot.gradtarg) / norm(refpot.gradtarg) + # hess_relerr = norm(U.hess - refpot.hesstarg) / norm(refpot.hesstarg) + @test pot_relerr < tolerance + # @test grad_relerr < tolerance + # @test hess_relerr < tolerance + + Utarg = lfmm2d(eps=tolerance,sources=sources, charges=charges, dipstr=dipstrs, dipvecs=dipvecs,targets=targets,pgt=1) + pot_relerrtarg = norm(Utarg.pottarg - refpottarg) / norm(refpottarg) + # grad_relerrtarg = norm(Utarg.gradtarg - refpottarg.gradtarg) / norm(refpottarg.gradtarg) + # hess_relerrtarg = norm(Utarg.hesstarg - refpottarg.hesstarg) / norm(refpottarg.hesstarg) + @test pot_relerrtarg < tolerance + # @test grad_relerrtarg < tolerance + # @test hess_relerrtarg < tolerance + end +end + +@testset "Laplace FMM 2D (real)" begin + N = 1000 + sources = rand(2, N) + charges = rand(N) + dipstrs = rand(N) + dipvecs = ones(2, N)/sqrt(2) + # Direct computations (N^2) + refpot = direct_self(sources,charges,dipvecs,dipstrs) + + M = 100 + targets = rand(2, M) + # Direct computations (N^2) + refpottarg = direct_targ(sources,charges,dipvecs,dipstrs,targets) + for tol_exp=-14:-1 + tolerance = 0.49*10.0^tol_exp + U = rfmm2d(eps=tolerance,sources=sources,charges=charges, dipstr=dipstrs, dipvecs=dipvecs,pg=2) + pot_relerr = norm(U.pot - refpot) / norm(refpot) + # grad_relerr = norm(U.grad - refpot.gradtarg) / norm(refpot.gradtarg) + # hess_relerr = norm(U.hess - refpot.hesstarg) / norm(refpot.hesstarg) + @test pot_relerr < tolerance + # @test grad_relerr < tolerance + # @test hess_relerr < tolerance + + Utarg = rfmm2d(eps=tolerance,sources=sources, charges=charges, dipstr=dipstrs, dipvecs=dipvecs,targets=targets,pgt=1) + pot_relerrtarg = norm(Utarg.pottarg - refpottarg) / norm(refpottarg) + # grad_relerrtarg = norm(Utarg.gradtarg - refpottarg.gradtarg) / norm(refpottarg.gradtarg) + # hess_relerrtarg = norm(Utarg.hesstarg - refpottarg.hesstarg) / norm(refpottarg.hesstarg) + @test pot_relerrtarg < tolerance + # @test grad_relerrtarg < tolerance + # @test hess_relerrtarg < tolerance + end +end + + +function direct_self_cfmm(source, charges, dipstrs) + if size(charges,2) == 1 + charges = transpose(charges) + end + if size(dipstrs,2) == 1 + dipstrs = transpose(dipstrs) + end + nd,n = size(charges) + pot = zeros(eltype(charges),nd,n) + N = size(source, 2) + for k=1:nd + charge = charges[k,:] + dipstr = dipstrs[k,:] + for i=1:N + for j=1:N + if i == j + continue + end + z = source[1,i] + im*source[2,i] + ej = source[1,j] + im*source[2,j] + r = z - ej + # r = ej - z + pot[k,i] += charge[j]*log(abs(r)) + dipstr[j]/r + # pot[k,i] += charge[j]*log((r)) + dipstr[j]/r + end + end + end + if nd == 1 + return transpose(pot) + else + return pot + end +end +function direct_targ_cfmm(source, charges, dipstrs, target) + if size(charges,2) == 1 + charges = transpose(charges) + end + if size(dipstrs,2) == 1 + dipstrs = transpose(dipstrs) + end + nd = size(charges,1) + M = size(target, 2) + N = size(source, 2) + pot = zeros(eltype(charges), nd, M) + for k in 1:nd + charge = charges[k,:] + dipstr = dipstrs[k,:] + for i=1:M + for j=1:N + z = target[1,i] + im*target[2,i] + ej = source[1,j] + im*source[2,j] + r = z - ej + # r = ej - z + pot[k,i] += charge[j]*log(abs(r)) + dipstr[j]/r + end + end + end + if nd == 1 + return transpose(pot) + else + return pot + end +end + +@testset "Cauchy FMM 2D (real)" begin + N = 1000 + sources = rand(2, N) + charges = rand(ComplexF64,N) # Only works with ones! Something is wrong + dipstrs = rand(ComplexF64,N) + charges_zeros = zeros(ComplexF64,N) + dipstrs_zeros = zeros(ComplexF64,N) + # Direct computations (N^2) + refpot_charges = direct_self_cfmm(sources,charges,dipstrs_zeros) + refpot_dipstrs = direct_self_cfmm(sources,charges_zeros,dipstrs) + + M = 100 + targets = rand(2, M) + # Direct computations (N^2) + refpottarg_charges = direct_targ_cfmm(sources,charges,dipstrs_zeros,targets) + refpottarg_dipstrs = direct_targ_cfmm(sources,charges_zeros,dipstrs,targets) + for tol_exp=-14:-1 + tolerance = 0.49*10.0^tol_exp + # FOR NOW ONLY TESTING REAL PART + U = cfmm2d(eps=tolerance,sources=sources,charges=charges,targets=targets, pg=1) + pot_relerr = norm(real.(U.pot - refpot_charges)) / norm(real.(refpot_charges)) + @test_broken pot_relerr < tolerance + U = cfmm2d(eps=tolerance,sources=sources, charges=charges_zeros, dipstr=dipstrs, pg=1) + pot_relerr = norm(real.(U.pot - refpot_dipstrs)) / norm(real.(refpot_dipstrs)) + @test pot_relerr < tolerance + + # FOR NOW ONLY TESTING REAL PART + Utarg = cfmm2d(eps=tolerance,sources=sources, charges=charges, targets=targets,pgt=1) + pot_relerrtarg = norm(real.(Utarg.pottarg - refpottarg_charges)) / norm(real.(refpottarg_charges)) + @test_broken pot_relerrtarg < tolerance + Utarg = cfmm2d(eps=tolerance,sources=sources, charges=charges_zeros, dipstr=dipstrs, targets=targets,pgt=1) + pot_relerrtarg = norm(real.(Utarg.pottarg - refpottarg_dipstrs)) / norm(real.(refpottarg_dipstrs)) + @test pot_relerrtarg < tolerance + end +end diff --git a/test/helmholtz.jl b/test/test_stokes.jl similarity index 100% rename from test/helmholtz.jl rename to test/test_stokes.jl