From 030feff1b211b537eefe05d74b1d8f9700e45b46 Mon Sep 17 00:00:00 2001 From: Unknown <34473471+bblankrot@users.noreply.github.com> Date: Tue, 4 Sep 2018 16:52:05 +0200 Subject: [PATCH 1/9] Added transpose FMM mvp --- src/fmm_mvp.jl | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/src/fmm_mvp.jl b/src/fmm_mvp.jl index 42244a9..36b55b0 100644 --- a/src/fmm_mvp.jl +++ b/src/fmm_mvp.jl @@ -137,3 +137,71 @@ function FMM_mainMVP_pre2!(output, beta, scatteringMatrices, φs::Vector{Float64 end #inbounds return output end + +function FMM_mainMVP_transpose!(output, beta, scatteringMatrices, φs::Vector{Float64}, ids::Vector{Int}, P, mFMM, pre_agg, translated_sum) + #calculate matrix^T-vector product - devectorized with pre-preconditioning + + @inbounds begin + #here we first multiply by Xᵀ + + #Xβ storage can be moved outward, but for now this preallocation prevents most + #dynamic mem alloc by this MVP + Xβ = Array{Complex{Float64}}(length(beta)) + for ic = 1:length(ids) + rng = (ic-1)*(2*P+1) + (1:2*P+1) + v = view(Xβ, rng) + if φs[ic] == 0.0 + At_mul_B!(v, scatteringMatrices[ids[ic]], view(beta, rng)) + else + #rotate without matrix - transposed + psss.rotateMultipole!(temp, view(beta, rng), φs[ic], P) + At_mul_B!(v, scatteringMatrices[ids[ic]], temp) + psss.rotateMultipole!(v, -φs[ic], P) + end + end + + At_mul_B!(output, mFMM.Znear, Xβ) + G = length(mFMM.groups) + fill!(pre_agg,0.0) + #first preagg all + for ig2 = 1:G + for is = 1:mFMM.groups[ig2].size + indices = (mFMM.groups[ig2].point_ids[is]-1)*(2*P+1) + for ii = 1:2*P+1 + for iq = 1:mFMM.Q + pre_agg[iq,ig2] += conj(mFMM.Agg[ig2][iq,(is-1)*(2*P+1) + ii])*Xβ[indices + ii] + end + end + end + end + + for ig1 = 1:G + #translate plane waves from ig2 to ig1 + fill!(translated_sum,0.0) + for ig2 = 1:G + if isempty(mFMM.Trans[(ig2-1)*G + ig1]) #transposed + continue + else + for iQ = 1:mFMM.Q + translated_sum[iQ] -= mFMM.Trans[(ig2-1)*G + ig1][iQ]*pre_agg[iQ,ig2] #transposed + end + #minus sign because the real equation is (I-TᵀXᵀ)β=b + end + end + #disaggregate from ig1 center to ig1's scatterers + for is = 1:mFMM.groups[ig1].size + for ip = 1:2*P+1 + disagged = 0.0im + for iq = 1:mFMM.Q + disagged += mFMM.Agg[ig1][iq,(is-1)*(2*P+1) + ip]*translated_sum[iq] + end + output[(mFMM.groups[ig1].point_ids[is]-1)*(2*P+1) + ip] += disagged + end + end + end + + #add identity matrix (Iβ) + output .+= beta + end #inbounds + return output +end From a7ec09e9b19c63e7904f0507d833c605e516e238 Mon Sep 17 00:00:00 2001 From: bblankrot <34473471+bblankrot@users.noreply.github.com> Date: Wed, 5 Sep 2018 08:08:43 +0200 Subject: [PATCH 2/9] =?UTF-8?q?Created=20optimize=5F=CF=86=5Fadj?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/optimize_phis_adj.jl | 172 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 172 insertions(+) create mode 100644 src/optimize_phis_adj.jl diff --git a/src/optimize_phis_adj.jl b/src/optimize_phis_adj.jl new file mode 100644 index 0000000..18b03f7 --- /dev/null +++ b/src/optimize_phis_adj.jl @@ -0,0 +1,172 @@ +""" + optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts, + optimopts::Optim.Options, minimize = true) + +Optimize the rotation angles of a particle collection for minimization or +maximization (depending on `minimize`) of the field intensity at `points`. +`optimopts` and `method` define the optimization method, convergence criteria, +and other optimization parameters. +Returns an object of type `Optim.MultivariateOptimizationResults`. +""" +function optimize_φ_adj(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts, + optimopts::Optim.Options, method; minimize = true) + + #stuff that is done once + verify_min_distance(shapes, centers, ids, points) || error("Particles are too close.") + mFMM, scatteringMatrices, scatteringLU, buf = + prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopts) + Ns = size(centers,1) + H = optimizationHmatrix(points, centers, Ns, P, k0) + α = u2α(k0, ui, centers, P) + uinc_ = uinc(k0, points, ui) + # Allocate buffers + shared_var = OptimBuffer(Ns,P,size(points,1)) + initial_φs = copy(φs0) + last_φs = similar(initial_φs) + last_φs == initial_φs && (last_φs[1] += 1) #should never happen but has + + if minimize + df = OnceDifferentiable( + φs -> optimize_φ_f(φs, shared_var, last_φs, α, + H, points, P, uinc_, Ns, k0, centers, + scatteringMatrices, ids, mFMM, fmmopts, buf), + (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, + last_φs, α, H, points, P, uinc_, + Ns, k0, centers, scatteringMatrices, + scatteringLU, ids, mFMM, fmmopts, buf, minimize), + initial_φs) + else + df = OnceDifferentiable( + φs -> -optimize_φ_f(φs, shared_var, last_φs, α, + H, points, P, uinc_, Ns, k0, centers, + scatteringMatrices, ids, mFMM, fmmopts, buf), + (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, + last_φs, α, H, points, P, uinc_, + Ns, k0, centers, scatteringMatrices, + scatteringLU, ids, mFMM, fmmopts, buf, minimize), + initial_φs) + end + optimize(df, initial_φs, method, optimopts) +end + +function optimize_φ_adj_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, ids, mFMM, opt, buf) + if φs != last_φs + copy!(last_φs, φs) + #do whatever common calculations and save to shared_var + #construct rhs + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + if φs[ic] == 0.0 + buf.rhs[rng] = scatteringMatrices[ids[ic]]*α[rng] + else + #rotate without matrix + rotateMultipole!(view(buf.rhs,rng),view(α,rng),-φs[ic],P) + buf.rhs[rng] = scatteringMatrices[ids[ic]]*buf.rhs[rng] + rotateMultipole!(view(buf.rhs,rng),φs[ic],P) + end + end + + if opt.method == "pre" + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_pre!(output_, x_, + scatteringMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + elseif opt.method == "pre2" + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_pre2!(output_, x_, + scatteringMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + end + + fill!(shared_var.β,0.0) + shared_var.β,ch = gmres!(shared_var.β, MVP, buf.rhs, + restart = Ns*(2*P+1), tol = opt.tol, log = true, + initially_zero = true) #no restart, preconditioning + !ch.isconverged && error("FMM process did not converge") + + shared_var.f[:] = H.'*shared_var.β + shared_var.f .+= uinc_ + end +end + +function optimize_φ_adj_f(φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers,scatteringMatrices, ids, mFMM, opt, buf) + optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, + k0, centers,scatteringMatrices, ids, mFMM, opt, buf) + + func = sum(abs2, shared_var.f) +end + +function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize) + optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, + k0, centers,scatteringMatrices, ids, mFMM, opt, buf) + + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_transpose!(output_, x_, + scatteringMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + #build rhs of adjoint problem + shared_var.rhs_grad[:] = ... + #solve adjoint problem + λadj, ch = gmres(MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + tol = opt.tol, log = true) + # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + # tol = opt.tol, log = true, initially_zero = true) + if ch.isconverged == false + display("FMM process did not converge for partial derivative.") + error("..") + end + + for n = 1:Ns + #compute n-th gradient + rng = (n-1)*(2*P+1) + (1:2*P+1) + ??rotateMultipole!(tempn, shared_var.β[rng], -φs[n], P) + ??tempn[:] = scatteringLU[ids[n]]\tempn #LU decomp with pivoting + ??tempn[:] .*= -D + ??v = view(, rng) + A_mul_B!(v, scatteringMatrices[ids[n]], tempn) + rotateMultipole!(v, φs[n], P) + v[:] += D.*shared_var.β[rng] + + #prepare for next one + v[:] = 0.0im + end + + grad_stor[:] = ifelse(minimize, 2, -2)* + real(shared_var.∂β.'*(H*conj(shared_var.f))) +end + +function optimizationHmatrix(points, centers, Ns, P, k0) + points_moved = Array{Float64}(2) + H = Array{Complex{Float64}}(Ns*(2*P+1), size(points,1)) + for ic = 1:Ns, i = 1:size(points,1) + points_moved[1] = points[i,1] - centers[ic,1] + points_moved[2] = points[i,2] - centers[ic,2] + r_angle = atan2(points_moved[2], points_moved[1]) + kr = k0*hypot(points_moved[1], points_moved[2]) + + ind = (ic-1)*(2*P+1) + P + 1 + H[ind,i] = besselh(0, kr) + for l = 1:P + b = besselh(l, kr) + H[ind + l,i] = b*exp(1.0im*l*r_angle) + H[ind - l,i] = b*(-1)^l*exp(-1.0im*l*r_angle) + end + end + H +end + +function prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopt) + #setup FMM reusal + Ns = size(centers,1) + groups, boxSize = divideSpace(centers, fmmopt) + P2, Q = FMMtruncation(fmmopt.acc, boxSize, k0) + mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) + scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) + scatteringLU = [lufact(scatteringMatrices[iid]) for iid = 1:length(shapes)] + buf = FMMbuffer(Ns,P,Q,length(groups)) + return mFMM, scatteringMatrices, scatteringLU, buf +end From a35055f46927d0ba15dca91ef6579d49383bd5f2 Mon Sep 17 00:00:00 2001 From: bblankrot <34473471+bblankrot@users.noreply.github.com> Date: Wed, 5 Sep 2018 21:35:22 +0200 Subject: [PATCH 3/9] updating adjoint optimization over phis, rs --- src/fmm_mvp.jl | 5 +- src/optimize_phis.jl | 80 +++++++++++++----- src/optimize_phis_adj.jl | 172 --------------------------------------- src/optimize_rs.jl | 72 +++++++++++----- 4 files changed, 116 insertions(+), 213 deletions(-) delete mode 100644 src/optimize_phis_adj.jl diff --git a/src/fmm_mvp.jl b/src/fmm_mvp.jl index 36b55b0..8e74529 100644 --- a/src/fmm_mvp.jl +++ b/src/fmm_mvp.jl @@ -147,6 +147,7 @@ function FMM_mainMVP_transpose!(output, beta, scatteringMatrices, φs::Vector{Fl #Xβ storage can be moved outward, but for now this preallocation prevents most #dynamic mem alloc by this MVP Xβ = Array{Complex{Float64}}(length(beta)) + temp = Array{Complex{Float64}}(2*P+1) for ic = 1:length(ids) rng = (ic-1)*(2*P+1) + (1:2*P+1) v = view(Xβ, rng) @@ -154,9 +155,9 @@ function FMM_mainMVP_transpose!(output, beta, scatteringMatrices, φs::Vector{Fl At_mul_B!(v, scatteringMatrices[ids[ic]], view(beta, rng)) else #rotate without matrix - transposed - psss.rotateMultipole!(temp, view(beta, rng), φs[ic], P) + rotateMultipole!(temp, view(beta, rng), φs[ic], P) At_mul_B!(v, scatteringMatrices[ids[ic]], temp) - psss.rotateMultipole!(v, -φs[ic], P) + rotateMultipole!(v, -φs[ic], P) end end diff --git a/src/optimize_phis.jl b/src/optimize_phis.jl index 2122aa3..6a636a1 100644 --- a/src/optimize_phis.jl +++ b/src/optimize_phis.jl @@ -6,10 +6,12 @@ Optimize the rotation angles of a particle collection for minimization or maximization (depending on `minimize`) of the field intensity at `points`. `optimopts` and `method` define the optimization method, convergence criteria, and other optimization parameters. +`adjoint` dictates whether the gradient is calculated using the adjoint method +(faster) or directly. Returns an object of type `Optim.MultivariateOptimizationResults`. """ function optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts, - optimopts::Optim.Options, method; minimize = true) + optimopts::Optim.Options, method; minimize = true, adjoint = true) #stuff that is done once verify_min_distance(shapes, centers, ids, points) || error("Particles are too close.") @@ -25,27 +27,22 @@ function optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, f last_φs = similar(initial_φs) last_φs == initial_φs && (last_φs[1] += 1) #should never happen but has - if minimize - df = OnceDifferentiable( - φs -> optimize_φ_f(φs, shared_var, last_φs, α, - H, points, P, uinc_, Ns, k0, centers, - scatteringMatrices, ids, mFMM, fmmopts, buf), - (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, - last_φs, α, H, points, P, uinc_, - Ns, k0, centers, scatteringMatrices, - scatteringLU, ids, mFMM, fmmopts, buf, minimize), - initial_φs) + fopt = φs -> ifelse(minimize,1,-1)*optimize_φ_f(φs, shared_var, last_φs, α, + H, points, P, uinc_, Ns, k0, centers, + scatteringMatrices, ids, mFMM, fmmopts, + buf) + if adjoint + gopt! = (grad_stor, φs) -> optimize_φ_adj_g!(grad_stor, φs, shared_var, + last_φs, α, H, points, P, uinc_, Ns, k0, + centers, scatteringMatrices, scatteringLU, + ids, mFMM, fmmopts, buf, minimize) else - df = OnceDifferentiable( - φs -> -optimize_φ_f(φs, shared_var, last_φs, α, - H, points, P, uinc_, Ns, k0, centers, - scatteringMatrices, ids, mFMM, fmmopts, buf), - (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, - last_φs, α, H, points, P, uinc_, - Ns, k0, centers, scatteringMatrices, - scatteringLU, ids, mFMM, fmmopts, buf, minimize), - initial_φs) + gopt! = (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, + last_φs, α, H, points, P, uinc_, Ns, k0, + centers, scatteringMatrices, scatteringLU, + ids, mFMM, fmmopts, buf, minimize) end + df = OnceDifferentiable(fopt, gopt!, initial_φs) optimize(df, initial_φs, method, optimopts) end @@ -149,6 +146,49 @@ function optimize_φ_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, real(shared_var.∂β.'*(H*conj(shared_var.f))) end + +function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize) + optimize_φ_adj_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, + k0, centers,scatteringMatrices, ids, mFMM, opt, buf) + + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_transpose!(output_, x_, + scatteringMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + #build rhs of adjoint problem + shared_var.rhs_grad[:] = -(H*conj(shared_var.f)) + #solve adjoint problem + λadj, ch = gmres(MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + tol = opt.tol, log = true) + # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + # tol = opt.tol, log = true, initially_zero = true) + if ch.isconverged == false + display("FMM process did not converge for adjoint system.") + error("..") + end + + #note: this assumes that there are exactly 2P+1 non zero elements in ∂X. If + #not, v must be (2P+1)Ns × 1. + D = -1.0im*collect(-P:P) + v = Array{Complex{Float64}}(2P+1) #TODO: minimize dynamic alloc + for n = 1:Ns + #compute n-th element of gradient + rng = (n-1)*(2*P+1) + (1:2*P+1) + rotateMultipole!(v, view(shared_var.β,rng), -φs[n], P) + v[:] = scatteringLU[ids[n]]\v #LU decomp with pivoting + v[:] .*= -D + v[:] = scatteringMatrices[ids[n]]*v + rotateMultipole!(v, φs[n], P) + v[:] += D.*shared_var.β[rng] + + grad_stor[n] = ifelse(minimize, -2, 2)*real(view(λadj,rng).'*v) + #prepare for next one + v[:] = 0.0im + end +end + function optimizationHmatrix(points, centers, Ns, P, k0) points_moved = Array{Float64}(2) H = Array{Complex{Float64}}(Ns*(2*P+1), size(points,1)) diff --git a/src/optimize_phis_adj.jl b/src/optimize_phis_adj.jl deleted file mode 100644 index 18b03f7..0000000 --- a/src/optimize_phis_adj.jl +++ /dev/null @@ -1,172 +0,0 @@ -""" - optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts, - optimopts::Optim.Options, minimize = true) - -Optimize the rotation angles of a particle collection for minimization or -maximization (depending on `minimize`) of the field intensity at `points`. -`optimopts` and `method` define the optimization method, convergence criteria, -and other optimization parameters. -Returns an object of type `Optim.MultivariateOptimizationResults`. -""" -function optimize_φ_adj(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts, - optimopts::Optim.Options, method; minimize = true) - - #stuff that is done once - verify_min_distance(shapes, centers, ids, points) || error("Particles are too close.") - mFMM, scatteringMatrices, scatteringLU, buf = - prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopts) - Ns = size(centers,1) - H = optimizationHmatrix(points, centers, Ns, P, k0) - α = u2α(k0, ui, centers, P) - uinc_ = uinc(k0, points, ui) - # Allocate buffers - shared_var = OptimBuffer(Ns,P,size(points,1)) - initial_φs = copy(φs0) - last_φs = similar(initial_φs) - last_φs == initial_φs && (last_φs[1] += 1) #should never happen but has - - if minimize - df = OnceDifferentiable( - φs -> optimize_φ_f(φs, shared_var, last_φs, α, - H, points, P, uinc_, Ns, k0, centers, - scatteringMatrices, ids, mFMM, fmmopts, buf), - (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, - last_φs, α, H, points, P, uinc_, - Ns, k0, centers, scatteringMatrices, - scatteringLU, ids, mFMM, fmmopts, buf, minimize), - initial_φs) - else - df = OnceDifferentiable( - φs -> -optimize_φ_f(φs, shared_var, last_φs, α, - H, points, P, uinc_, Ns, k0, centers, - scatteringMatrices, ids, mFMM, fmmopts, buf), - (grad_stor, φs) -> optimize_φ_g!(grad_stor, φs, shared_var, - last_φs, α, H, points, P, uinc_, - Ns, k0, centers, scatteringMatrices, - scatteringLU, ids, mFMM, fmmopts, buf, minimize), - initial_φs) - end - optimize(df, initial_φs, method, optimopts) -end - -function optimize_φ_adj_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, ids, mFMM, opt, buf) - if φs != last_φs - copy!(last_φs, φs) - #do whatever common calculations and save to shared_var - #construct rhs - for ic = 1:Ns - rng = (ic-1)*(2*P+1) + (1:2*P+1) - if φs[ic] == 0.0 - buf.rhs[rng] = scatteringMatrices[ids[ic]]*α[rng] - else - #rotate without matrix - rotateMultipole!(view(buf.rhs,rng),view(α,rng),-φs[ic],P) - buf.rhs[rng] = scatteringMatrices[ids[ic]]*buf.rhs[rng] - rotateMultipole!(view(buf.rhs,rng),φs[ic],P) - end - end - - if opt.method == "pre" - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_pre!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - elseif opt.method == "pre2" - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_pre2!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - end - - fill!(shared_var.β,0.0) - shared_var.β,ch = gmres!(shared_var.β, MVP, buf.rhs, - restart = Ns*(2*P+1), tol = opt.tol, log = true, - initially_zero = true) #no restart, preconditioning - !ch.isconverged && error("FMM process did not converge") - - shared_var.f[:] = H.'*shared_var.β - shared_var.f .+= uinc_ - end -end - -function optimize_φ_adj_f(φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers,scatteringMatrices, ids, mFMM, opt, buf) - optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, - k0, centers,scatteringMatrices, ids, mFMM, opt, buf) - - func = sum(abs2, shared_var.f) -end - -function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize) - optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, - k0, centers,scatteringMatrices, ids, mFMM, opt, buf) - - MVP = LinearMap{eltype(buf.rhs)}( - (output_, x_) -> FMM_mainMVP_transpose!(output_, x_, - scatteringMatrices, φs, ids, P, mFMM, - buf.pre_agg, buf.trans), - Ns*(2*P+1), Ns*(2*P+1), ismutating = true) - - #build rhs of adjoint problem - shared_var.rhs_grad[:] = ... - #solve adjoint problem - λadj, ch = gmres(MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), - tol = opt.tol, log = true) - # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), - # tol = opt.tol, log = true, initially_zero = true) - if ch.isconverged == false - display("FMM process did not converge for partial derivative.") - error("..") - end - - for n = 1:Ns - #compute n-th gradient - rng = (n-1)*(2*P+1) + (1:2*P+1) - ??rotateMultipole!(tempn, shared_var.β[rng], -φs[n], P) - ??tempn[:] = scatteringLU[ids[n]]\tempn #LU decomp with pivoting - ??tempn[:] .*= -D - ??v = view(, rng) - A_mul_B!(v, scatteringMatrices[ids[n]], tempn) - rotateMultipole!(v, φs[n], P) - v[:] += D.*shared_var.β[rng] - - #prepare for next one - v[:] = 0.0im - end - - grad_stor[:] = ifelse(minimize, 2, -2)* - real(shared_var.∂β.'*(H*conj(shared_var.f))) -end - -function optimizationHmatrix(points, centers, Ns, P, k0) - points_moved = Array{Float64}(2) - H = Array{Complex{Float64}}(Ns*(2*P+1), size(points,1)) - for ic = 1:Ns, i = 1:size(points,1) - points_moved[1] = points[i,1] - centers[ic,1] - points_moved[2] = points[i,2] - centers[ic,2] - r_angle = atan2(points_moved[2], points_moved[1]) - kr = k0*hypot(points_moved[1], points_moved[2]) - - ind = (ic-1)*(2*P+1) + P + 1 - H[ind,i] = besselh(0, kr) - for l = 1:P - b = besselh(l, kr) - H[ind + l,i] = b*exp(1.0im*l*r_angle) - H[ind - l,i] = b*(-1)^l*exp(-1.0im*l*r_angle) - end - end - H -end - -function prepare_fmm_reusal_φs(k0, kin, P, shapes, centers, ids, fmmopt) - #setup FMM reusal - Ns = size(centers,1) - groups, boxSize = divideSpace(centers, fmmopt) - P2, Q = FMMtruncation(fmmopt.acc, boxSize, k0) - mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) - scatteringMatrices,innerExpansions = particleExpansion(k0, kin, shapes, P, ids) - scatteringLU = [lufact(scatteringMatrices[iid]) for iid = 1:length(shapes)] - buf = FMMbuffer(Ns,P,Q,length(groups)) - return mFMM, scatteringMatrices, scatteringLU, buf -end diff --git a/src/optimize_rs.jl b/src/optimize_rs.jl index 655155b..e0cc1af 100644 --- a/src/optimize_rs.jl +++ b/src/optimize_rs.jl @@ -12,12 +12,14 @@ the optimized device. `optimopts` defines the convergence criteria and other optimization parameters for both the inner and outer iterations. `method` can be either `"BFGS"` or `"LBFGS"`. See the `Optim.Fminbox` documentation for more details. +`adjoint` dictates whether the gradient is calculated using the adjoint method +(faster) or directly. Returns an object of type `Optim.MultivariateOptimizationResults`. """ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, centers, fmmopts, optimopts::Optim.Options; - minimize = true, method = "BFGS") + minimize = true, method = "BFGS", adjoint = true) Ns = size(centers,1) J = length(rs0) @@ -56,27 +58,22 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, initial_rs = copy(rs0) last_rs = similar(initial_rs) - if minimize - df = OnceDifferentiable(rs -> optimize_radius_f(rs, last_rs, shared_var, - φs, α, H, points, P, uinc_, Ns, k0, - kin, centers,scatteringMatrices, dS_S, - ids, mFMM, fmmopts, buf), - (grad_stor, rs) -> optimize_radius_g!(grad_stor, rs, last_rs, - shared_var, φs, α, H, points, P, uinc_, - Ns, k0, kin, centers, scatteringMatrices, - dS_S, ids, mFMM, fmmopts, buf, minimize), - initial_rs) + if adjoint + gopt! = (grad_stor, rs) -> optimize_radius_adj_g!(grad_stor, rs, last_rs, + shared_var, φs, α, H, points, P, uinc_, + Ns, k0, kin, centers, scatteringMatrices, + dS_S, ids, mFMM, fmmopts, buf, minimize) else - df = OnceDifferentiable(rs -> -optimize_radius_f(rs, last_rs, shared_var, + gopt! = (grad_stor, rs) -> optimize_radius_g!(grad_stor, rs, last_rs, + shared_var, φs, α, H, points, P, uinc_, + Ns, k0, kin, centers, scatteringMatrices, + dS_S, ids, mFMM, fmmopts, buf, minimize) + end + fopt = rs -> ifelse(minimize,1,-1)*optimize_radius_f(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers,scatteringMatrices, dS_S, - ids, mFMM, fmmopts, buf), - (grad_stor, rs) -> optimize_radius_g!(grad_stor, rs, last_rs, - shared_var, φs, α, H, points, P, uinc_, - Ns, k0, kin, centers, scatteringMatrices, - dS_S, ids, mFMM, fmmopts, buf, minimize), - initial_rs) - end + ids, mFMM, fmmopts, buf) + df = OnceDifferentiable(fopt, gopt!, initial_rs) if method == "LBFGS" optimize(df, r_min, r_max, initial_rs, @@ -179,6 +176,43 @@ function optimize_radius_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, poin grad_stor[:] = ifelse(minimize,2,-2)*real(shared_var.∂β.'*(H*conj(shared_var.f))) end +function optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf, minimize) + optimize_radius_common!(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) + + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_transpose!(output_, x_, + scatteringMatrices, φs, ids, P, mFMM, + buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + #build rhs of adjoint problem + shared_var.rhs_grad[:] = -(H*conj(shared_var.f)) + #solve adjoint problem + λadj, ch = gmres(MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + tol = opt.tol, log = true) + # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + # tol = opt.tol, log = true, initially_zero = true) + if ch.isconverged == false + display("FMM process did not converge for adjoint system.") + error("..") + end + + for n = 1:length(rs) + #compute n-th gradient - here we must pay the price for symmetry + #as more than one beta is affected. Overwrites rhs_grad. + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + if ids[ic] == n + shared_var.rhs_grad[rng] = dS_S[n]*shared_var.β[rng] + else + shared_var.rhs_grad[rng] = 0.0 + end + end + + grad_stor[n] = ifelse(minimize,-2,2)*real(λadj.'*shared_var.rhs_grad) + end +end + function updateCircleScatteringDerivative!(S, dS_S, kout, kin, R, P) #non-vectorized, reuses bessel From 76ce18db85cec3b0aec46cc78c79030c1b55d032 Mon Sep 17 00:00:00 2001 From: Unknown <34473471+bblankrot@users.noreply.github.com> Date: Sat, 8 Sep 2018 19:14:49 +0200 Subject: [PATCH 4/9] updated optimize_rs with adj, added gradient calc --- src/optimize_phis.jl | 2 +- src/optimize_rs.jl | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/src/optimize_phis.jl b/src/optimize_phis.jl index 6a636a1..325c852 100644 --- a/src/optimize_phis.jl +++ b/src/optimize_phis.jl @@ -148,7 +148,7 @@ end function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, P, uinc_, Ns, k0, centers, scatteringMatrices, scatteringLU, ids, mFMM, opt, buf, minimize) - optimize_φ_adj_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, + optimize_φ_common!(φs, last_φs, shared_var, α, H, points, P, uinc_, Ns, k0, centers,scatteringMatrices, ids, mFMM, opt, buf) MVP = LinearMap{eltype(buf.rhs)}( diff --git a/src/optimize_rs.jl b/src/optimize_rs.jl index e0cc1af..e81ad49 100644 --- a/src/optimize_rs.jl +++ b/src/optimize_rs.jl @@ -244,3 +244,38 @@ function updateCircleScatteringDerivative!(S, dS_S, kout, kin, R, P) pre_H = H end end + +function gradient_radius(rs, points, ids, P, ui, k0, kin, centers, fmmopts, minimize) + Ns = size(centers,1) + J = length(rs) + + assert(maximum(ids) <= J) + verify_min_distance(CircleParams.(rs), centers, ids, + points) || error("Particles are too close.") + + #setup FMM reusal + groups, boxSize = divideSpace(centers, fmmopts) + P2, Q = FMMtruncation(fmmopts.acc, boxSize, k0) + mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) + + #allocate derivative + scatteringMatrices = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] + dS_S = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] + + #stuff that is done once + H = optimizationHmatrix(points, centers, Ns, P, k0) + α = u2α(k0, ui, centers, P) + uinc_ = uinc(k0, points, ui) + φs = zeros(Float64,Ns) + + # Allocate buffers + buf = FMMbuffer(Ns,P,Q,length(groups)) + shared_var = OptimBuffer(Ns,P,size(points,1),J) + last_rs = similar(rs); all(last_rs .== rs) && (last_rs[1] += 1) + + grad_stor = Array{Float64}(J) + optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, + P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, + fmmopts, buf, minimize) + grad_stor +end From 92ea900cf06b6eb03efa7638727f56147dfcbe41 Mon Sep 17 00:00:00 2001 From: Unknown <34473471+bblankrot@users.noreply.github.com> Date: Tue, 9 Oct 2018 16:56:25 +0200 Subject: [PATCH 5/9] Adding multi-frequency optimization --- src/ParticleScattering.jl | 2 + src/optimize_rs.jl | 6 +- src/optimize_rs_mf.jl | 168 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 173 insertions(+), 3 deletions(-) create mode 100644 src/optimize_rs_mf.jl diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl index daa03a4..15bf83f 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -63,4 +63,6 @@ export plot_near_field_pgf export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices export FMMbuffer, FMMmatrices +include("optimize_rs_mf.jl") +export optimize_rs_mf end diff --git a/src/optimize_rs.jl b/src/optimize_rs.jl index e81ad49..208bf71 100644 --- a/src/optimize_rs.jl +++ b/src/optimize_rs.jl @@ -34,8 +34,8 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, else assert(J == length(r_max)) end - verify_min_distance([CircleParams(r_max[i]) for i = 1:J], centers, ids, - points) || error("Particles are too close.") + verify_min_distance(CircleParams.(r_max), centers, ids, + points) || error("Particles are too close or r_max are too large.") #setup FMM reusal groups, boxSize = divideSpace(centers, fmmopts) @@ -56,7 +56,7 @@ function optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, buf = FMMbuffer(Ns,P,Q,length(groups)) shared_var = OptimBuffer(Ns,P,size(points,1),J) initial_rs = copy(rs0) - last_rs = similar(initial_rs) + last_rs = similar(initial_rs); last_rs[1] = NaN; assert(last_rs != initial_rs) #initial_rs, last_rs must be different before first iteration if adjoint gopt! = (grad_stor, rs) -> optimize_radius_adj_g!(grad_stor, rs, last_rs, diff --git a/src/optimize_rs_mf.jl b/src/optimize_rs_mf.jl new file mode 100644 index 0000000..fd1a2d1 --- /dev/null +++ b/src/optimize_rs_mf.jl @@ -0,0 +1,168 @@ +const cₚ = 10.0 #penalty term for unwanted frequency + +function optimize_rs_mf(rs0, r_min, r_max, points, point_flags, ids, P, ui, k0, + kin, centers, fmmopts, optimopts::Optim.Options; + method = BFGS(linesearch = LineSearches.BackTracking())) + + Ns = size(centers,1) + φs = zeros(Ns) + Nk = length(k0); assert(Nk == length(kin)) + J = length(rs0) + assert(maximum(ids) <= J) + if length(r_min) == 1 + r_min = r_min*ones(Float64,J) + else + assert(J == length(r_min)) + end + if length(r_max) == 1 + r_max = r_max*ones(Float64,J) + else + assert(J == length(r_max)) + end + verify_min_distance(CircleParams.(r_max), centers, ids, + points) || error("Particles are too close or r_max are too large.") + + groups,boxSize = divideSpace(centers, fmmopts) + P2,Q = FMMtruncation(fmmopts.acc, boxSize, maximum(k0)) + mFMM = [FMMbuildMatrices(k0[ik], P, P2, Q, groups, centers, boxSize) for ik = 1:Nk] + scatteringMatrices = [[speye(Complex{Float64}, 2*P+1) for i = 1:J] for ik = 1:Nk] + dS_S = [[speye(Complex{Float64}, 2*P+1) for i = 1:J] for ik = 1:Nk] + + #calculate nd expand incident field + α = [u2α(k0[ik], ui, centers, P) for ik = 1:Nk] + uinc_ = [uinc(k0[ik], points, ui) for ik = 1:Nk] + H = [optimizationHmatrix(points, centers, Ns, P, k0[ik]) for ik = 1:Nk] + + # Allocate buffers + buf = FMMbuffer(Ns,P,Q,length(groups)) + shared_var = [OptimBuffer(Ns, P, size(points,1), J) for ik = 1:Nk] + initial_rs = copy(rs0) + last_rs = similar(initial_rs); last_rs[1] = NaN; assert(last_rs != initial_rs) #initial_rs, last_rs must be different before first iteration + df = OnceDifferentiable( + rs -> optimize_rs_mf_f(rs, last_rs, shared_var, φs, α, H, points, point_flags, + P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, + ids, mFMM, fmmopts, buf), + (grad_stor, rs) -> optimize_rs_mf_g!(grad_stor, rs, last_rs, shared_var, φs, + α, H, points, point_flags, P, uinc_, Ns, k0, kin, centers, + scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf), + initial_rs) + + optimize(df, r_min, r_max, initial_rs, Fminbox(method), optimopts) +end + +function optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, + H, points, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, + dS_S, ids, mFMM, fmmopts, buf) + if rs != last_rs + copy!(last_rs, rs) + #do whatever common calculations and save to shared_var + #construct rhs + for ik in eachindex(k0) + for id in unique(ids) + try + updateCircleScatteringDerivative!(scatteringMatrices[ik][id], + dS_S[ik][id], k0[ik], kin[ik], rs[id], P) + catch + warn("Could not calculate derivatives for id=$id,k0=$k0,kin=$kin,R=$(rs[id])") + rethrow() + end + end + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + buf.rhs[rng] = scatteringMatrices[ik][ids[ic]]*α[ik][rng] + end + + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_pre!(output_, x_, + scatteringMatrices[ik], φs, ids, P, + mFMM[ik], buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + shared_var[ik].β[:] = 0 + shared_var[ik].β,ch = gmres!(shared_var[ik].β, MVP, buf.rhs, + restart = Ns*(2*P+1), tol = fmmopts.tol, + log = true, initially_zero = true) + !ch.isconverged && error("""FMM process did not converge, normalized + residual: $(norm(MVP*shared_var[ik].β - buf.rhs)/norm(buf.rhs))""") + + shared_var[ik].f[:] = H[ik].'*shared_var[ik].β + shared_var[ik].f[:] .+= uinc_[ik] + end + end +end + +function optimize_rs_mf_f(rs, last_rs, shared_var, φs, α, H, points, point_flags, + P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, + ids, mFMM, fmmopts, buf) + optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, + Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, fmmopts, buf) + + fobj_rs_mf(shared_var, point_flags) +end + +function optimize_rs_mf_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, + point_flags, P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, + ids, mFMM, opt, buf) + optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, H, points, P, uinc_, + Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, opt, buf) + + fill!(grad_stor, 0) + for ik in eachindex(k0) + MVP = LinearMap{eltype(buf.rhs)}( + (output_, x_) -> FMM_mainMVP_transpose!(output_, x_, + scatteringMatrices[ik], φs, ids, P, + mFMM[ik], buf.pre_agg, buf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + #build rhs of adjoint problem + gobj_rs_mf!(shared_var, ik, H[ik], point_flags) + #solve adjoint problem + λadj, ch = gmres(MVP, shared_var[ik].rhs_grad, restart = Ns*(2*P+1), + tol = opt.tol, log = true) + # λadj, ch = gmres!(λadj, MVP, shared_var.rhs_grad, restart = Ns*(2*P+1), + # tol = opt.tol, log = true, initially_zero = true) + + for n = 1:length(rs) + #compute n-th gradient - here we must pay the price for symmetry + #as more than one beta is affected. Overwrites rhs_grad. + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + if ids[ic] == n + shared_var[ik].rhs_grad[rng] = dS_S[ik][n]*shared_var[ik].β[rng] + else + shared_var[ik].rhs_grad[rng] = 0.0 + end + end + grad_stor[n] += -2*real(λadj.'*shared_var[ik].rhs_grad) + end + end +end + +function fobj_rs_mf(sv, point_flags) + #emphasize field maximization + fobj = 0.0 + for i in eachindex(point_flags) + ik = point_flags[i] + #add and subtract for easy sum + sumu = sum(abs2, sv[jk].f[i] for jk in eachindex(sv)) - abs2(sv[ik].f[i]) + fobj += (1 + cₚ*sumu)/abs2(sv[ik].f[i]) + end + fobj +end + +function gobj_rs_mf!(sv, ik, H, point_flags) + #calculate -(∂f/∂β)ᵀ for adjoint method + fill!(sv[ik].rhs_grad, 0) + for i in eachindex(point_flags) + ik′ = point_flags[i] + if ik == ik′ + sumu = sum(abs2, sv[jk].f[i] for jk in eachindex(sv)) + num = 1 + cₚ*(sumu - abs2(sv[ik].f[i])) + denom = abs2(sv[ik].f[i])*sv[ik].f[i] + else + num = (-cₚ)*conj(sv[ik].f[i]) + denom = abs2(sv[ik′].f[i]) + end + sv[ik].rhs_grad .+= view(H,:,i)*(num/denom) + end +end From bd8015717e196765158a4c75171f71c5ebf326c1 Mon Sep 17 00:00:00 2001 From: Unknown <34473471+bblankrot@users.noreply.github.com> Date: Wed, 19 Dec 2018 19:26:44 +0100 Subject: [PATCH 6/9] Added power optimization, Hankel recurrence --- src/ParticleScattering.jl | 10 +- src/incident.jl | 63 ++++++++++- src/multipole.jl | 37 +++++- src/optimize_rs.jl | 72 ++++++------ src/optimize_rs_mf.jl | 9 +- src/optimize_rs_pwr.jl | 229 ++++++++++++++++++++++++++++++++++++++ src/poynting.jl | 112 +++++++++++++++++++ src/shapes.jl | 4 +- src/visualization.jl | 35 +++--- 9 files changed, 505 insertions(+), 66 deletions(-) create mode 100644 src/optimize_rs_pwr.jl create mode 100644 src/poynting.jl diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl index 15bf83f..b0ccafe 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -51,7 +51,9 @@ export optimize_φ #methods, optimize_rs.jl export optimize_radius #methods, incident.jl -export uinc +export uinc, hxinc, hyinc +#consts, incident.jl +export eta0 # For advanced plotting with pgfplots import DataFrames, CSV, PGFPlotsX; const pgf = PGFPlotsX @@ -65,4 +67,10 @@ export FMMbuffer, FMMmatrices include("optimize_rs_mf.jl") export optimize_rs_mf + +include("poynting.jl") +export calc_power + +include("optimize_rs_pwr.jl") +export PowerBuffer, OptimProblemBuffer, optMatrixPwr, optimize_pwr_rs_f, optimize_pwr_rs_g! end diff --git a/src/incident.jl b/src/incident.jl index 4f703fb..a415973 100644 --- a/src/incident.jl +++ b/src/incident.jl @@ -1,3 +1,5 @@ +const eta0 = 4π*299792458e-7 + #plane wave function u2α(k, u::PlaneWave, centers::Array{T,2}, P) where T <: Real #build incoming coefficients for plane wave incident field @@ -57,7 +59,7 @@ function uinc(k, points::Array{T,1}, u::LineSource) where T <: Real 0.25im*besselh(0, 1, k*r) end -#current source +#current source TODO: simple avg->trap rule function u2α(k, u::CurrentSource, centers::Array{T,2}, P) where T <: Real α = zeros(Complex{Float64}, size(centers,1)*(2P + 1)) c = 0.25im*u.len/length(u.σ) @@ -123,3 +125,62 @@ function err_α(k, ui::Einc, P, shape, center) u_α = sum(α[p + P + 1]*besselj(p, k*shape.R)*exp.(1im*p*θ) for p = -P:P) norm(u_α - u_exact)/norm(u_exact) end + + +function hxinc(k, p, u::PlaneWave) + (sin(u.θi)/eta0)*exp.(1.0im*k*(cos(u.θi)*p[:,1] + sin(u.θi)*p[:,2])) +end + +function hyinc(k, p, u::PlaneWave) + (-cos(u.θi)/eta0)*exp.(1.0im*k*(cos(u.θi)*p[:,1] + sin(u.θi)*p[:,2])) +end + +function hxinc(k, p, u::LineSource) + R = sqrt.((view(p,:,1) - u.x).^2 + (view(p,:,2) - u.y).^2) + if any(R .== 0) + warn("Hinc: encountered singularity in incident field, returned NaN") + R[R.==0] = NaN + end + h = (0.25/eta0./R).*besselh.(1,k*R).*(u.y - view(p,:,2)) +end + +function hyinc(k, p, u::LineSource) + R = sqrt.((view(p,:,1) - u.x).^2 + (view(p,:,2) - u.y).^2) + if any(R .== 0) + warn("Hinc: encountered singularity in incident field, returned NaN") + R[R.==0] = NaN + end + h = (0.25/eta0./R).*besselh.(1,k*R).*(view(p,:,1) - u.x) +end + +#these untested +function hxinc(k, points, u::CurrentSource) + res = Array{Complex{Float64}}(size(points, 1)) + r = Array{Float64}(size(u.p, 1)) + y = Array{Float64}(size(u.p, 1)) + for i = 1:size(points, 1) + y .= points[i,2] - u.p[:,2] + r .= hypot.(points[i,1] - u.p[:,1], y) + if any(r .== 0) #point is on source + warn("Hinc: encountered singularity in incident field, returned NaN") + res[i] = NaN + end + res[i] = (-0.25*u.len/length(u.σ)/eta0)*sum(besselh.(1, 1, k*r).*u.σ.*y./r) + end + res +end +function hyinc(k, points, u::CurrentSource) + res = Array{Complex{Float64}}(size(points, 1)) + r = Array{Float64}(size(u.p, 1)) + x = Array{Float64}(size(u.p, 1)) + for i = 1:size(points, 1) + x .= points[i,1] - u.p[:,1] + r .= hypot.(x, points[i,2] - u.p[:,2]) + if any(r .== 0) #point is on source + warn("Hinc: encountered singularity in incident field, returned NaN") + res[i] = NaN + end + res[i] = (0.25*u.len/length(u.σ)/eta0)*sum(besselh.(1, 1, k*r).*u.σ.*x./r) + end + res +end diff --git a/src/multipole.jl b/src/multipole.jl index 1272303..627204f 100644 --- a/src/multipole.jl +++ b/src/multipole.jl @@ -96,13 +96,17 @@ function M2Lmatrix!(T, k, P, d) end end -function scattered_field_multipole(k0, beta, centers::Array{Float64,2}, points::Array{Float64,2}) +function scattered_field_multipole(k0, beta, centers::Array{Float64,2}, points::Array{Float64,2}; recurrence = false) Ns = size(centers,1) P = div(div(length(beta),Ns)-1,2) len = size(points,1) Esc = zeros(Complex{Float64}, len) - scattered_field_multipole!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len) + if recurrence + scattered_field_multipole_recurrence!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len) + else + scattered_field_multipole!(Esc, k0, beta, P, centers, 1:Ns, points, 1:len) + end return Esc end @@ -115,11 +119,10 @@ function scattered_field_multipole!(Esc::Array{Complex{Float64},1}, k0, beta, P, rs_moved = hypot(points_moved1, points_moved2) ts_moved = atan2(points_moved2, points_moved1) - try - Esc[ip] += beta[ind]*besselh(0, 1, k0*rs_moved) - catch + if rs_moved == 0 error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip") end + Esc[ip] += beta[ind]*besselh(0, 1, k0*rs_moved) for p = 1:P Esc[ip] += besselh(p, 1, k0*rs_moved)*(beta[p + ind]*exp(1im*p*ts_moved) + (-1)^p*beta[-p + ind]*exp(-1im*p*ts_moved)) end @@ -127,6 +130,30 @@ function scattered_field_multipole!(Esc::Array{Complex{Float64},1}, k0, beta, P, end end +function scattered_field_multipole_recurrence!(Esc::Array{Complex{Float64},1}, k0, beta, P, centers::Array{Float64,2}, ind_centers, points::Array{Float64,2}, ind_points) + for ic in ind_centers + ind = (ic-1)*(2*P+1) + P + 1 + for ip in ind_points + points_moved1 = points[ip,1] - centers[ic,1] + points_moved2 = points[ip,2] - centers[ic,2] + + rs_moved = hypot(points_moved1, points_moved2) + ts_moved = atan2(points_moved2, points_moved1) + if rs_moved == 0 + error("rs_moved == 0, center=$(centers[ic,:]), point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip") + end + Hₚ₋₂ = besselh(-1, 1, k0*rs_moved) + Hₚ₋₁ = besselh(0, 1, k0*rs_moved) + Esc[ip] += beta[ind]*Hₚ₋₁ + for p = 1:P + Hₚ = (2*(p-1)/(k0*rs_moved))*Hₚ₋₁ - Hₚ₋₂ + Esc[ip] += Hₚ*(beta[p + ind]*exp(1im*p*ts_moved) + (-1)^p*beta[-p + ind]*exp(-1im*p*ts_moved)) + Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ + end + end + end +end + function circleScatteringMatrix(kout, kin, R, P; gamma = false) #non-vectorized, reuses bessel S = Array{Complex{Float64}}(2*P+1) diff --git a/src/optimize_rs.jl b/src/optimize_rs.jl index 208bf71..d82da93 100644 --- a/src/optimize_rs.jl +++ b/src/optimize_rs.jl @@ -213,9 +213,9 @@ function optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, end end -function updateCircleScatteringDerivative!(S, dS_S, kout, kin, R, P) +function updateCircleScatteringDerivative!(S, dS_S, kout, kin, R::Real, P) #non-vectorized, reuses bessel - + R > 0 || throw(DomainError(R, "`R` must be positive.")) pre_J0 = besselj(-1,kout*R) pre_J1 = besselj(-1,kin*R) pre_H = besselh(-1,kout*R) @@ -245,37 +245,37 @@ function updateCircleScatteringDerivative!(S, dS_S, kout, kin, R, P) end end -function gradient_radius(rs, points, ids, P, ui, k0, kin, centers, fmmopts, minimize) - Ns = size(centers,1) - J = length(rs) - - assert(maximum(ids) <= J) - verify_min_distance(CircleParams.(rs), centers, ids, - points) || error("Particles are too close.") - - #setup FMM reusal - groups, boxSize = divideSpace(centers, fmmopts) - P2, Q = FMMtruncation(fmmopts.acc, boxSize, k0) - mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) - - #allocate derivative - scatteringMatrices = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] - dS_S = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] - - #stuff that is done once - H = optimizationHmatrix(points, centers, Ns, P, k0) - α = u2α(k0, ui, centers, P) - uinc_ = uinc(k0, points, ui) - φs = zeros(Float64,Ns) - - # Allocate buffers - buf = FMMbuffer(Ns,P,Q,length(groups)) - shared_var = OptimBuffer(Ns,P,size(points,1),J) - last_rs = similar(rs); all(last_rs .== rs) && (last_rs[1] += 1) - - grad_stor = Array{Float64}(J) - optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, - P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, - fmmopts, buf, minimize) - grad_stor -end +# function gradient_radius(rs, points, ids, P, ui, k0, kin, centers, fmmopts, minimize) +# Ns = size(centers,1) +# J = length(rs) +# +# assert(maximum(ids) <= J) +# verify_min_distance(CircleParams.(rs), centers, ids, +# points) || error("Particles are too close.") +# +# #setup FMM reusal +# groups, boxSize = divideSpace(centers, fmmopts) +# P2, Q = FMMtruncation(fmmopts.acc, boxSize, k0) +# mFMM = FMMbuildMatrices(k0, P, P2, Q, groups, centers, boxSize, tri = true) +# +# #allocate derivative +# scatteringMatrices = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] +# dS_S = [speye(Complex{Float64}, 2*P+1) for ic = 1:J] +# +# #stuff that is done once +# H = optimizationHmatrix(points, centers, Ns, P, k0) +# α = u2α(k0, ui, centers, P) +# uinc_ = uinc(k0, points, ui) +# φs = zeros(Float64,Ns) +# +# # Allocate buffers +# buf = FMMbuffer(Ns,P,Q,length(groups)) +# shared_var = OptimBuffer(Ns,P,size(points,1),J) +# last_rs = similar(rs); all(last_rs .== rs) && (last_rs[1] += 1) +# +# grad_stor = Array{Float64}(J) +# optimize_radius_adj_g!(grad_stor, rs, last_rs, shared_var, φs, α, H, points, +# P, uinc_, Ns, k0, kin, centers, scatteringMatrices, dS_S, ids, mFMM, +# fmmopts, buf, minimize) +# grad_stor +# end diff --git a/src/optimize_rs_mf.jl b/src/optimize_rs_mf.jl index fd1a2d1..a2a9a1f 100644 --- a/src/optimize_rs_mf.jl +++ b/src/optimize_rs_mf.jl @@ -59,13 +59,8 @@ function optimize_rs_mf_common!(rs, last_rs, shared_var, φs, α, #construct rhs for ik in eachindex(k0) for id in unique(ids) - try - updateCircleScatteringDerivative!(scatteringMatrices[ik][id], - dS_S[ik][id], k0[ik], kin[ik], rs[id], P) - catch - warn("Could not calculate derivatives for id=$id,k0=$k0,kin=$kin,R=$(rs[id])") - rethrow() - end + updateCircleScatteringDerivative!(scatteringMatrices[ik][id], + dS_S[ik][id], k0[ik], kin[ik], rs[id], P) end for ic = 1:Ns rng = (ic-1)*(2*P+1) + (1:2*P+1) diff --git a/src/optimize_rs_pwr.jl b/src/optimize_rs_pwr.jl new file mode 100644 index 0000000..c1f01ab --- /dev/null +++ b/src/optimize_rs_pwr.jl @@ -0,0 +1,229 @@ +#for now just the main functions, no wrapper +#remove all unnecessary allocations +mutable struct PowerBuffer #one for each *power calculation* (can be > 1 for each frequency/incident wave) + iₚ::Integer #which problem (β) is this?? + Ez::Vector{Complex{Float64}} + Hx::Vector{Complex{Float64}} + Hy::Vector{Complex{Float64}} + Ez_inc::Vector{Complex{Float64}} + Hx_inc::Vector{Complex{Float64}} + Hy_inc::Vector{Complex{Float64}} + pow::Float64 #power + HEz::Vector{Vector{Complex{Float64}}} + HHx::Vector{Vector{Complex{Float64}}} + HHy::Vector{Vector{Complex{Float64}}} + nhat::Array{Float64,2} + S::Array{Complex{Float64},2} + ∂pow::Vector{Complex{Float64}} # ∂power/∂βᵀ + l::Float64 #arc length + PowerBuffer(Ns::Integer, P::Integer, nhat, iₚ) = + new(iₚ, + Array{Complex{Float64}}(size(nhat, 1)), + Array{Complex{Float64}}(size(nhat, 1)), + Array{Complex{Float64}}(size(nhat, 1)), + Array{Complex{Float64}}(size(nhat, 1)), + Array{Complex{Float64}}(size(nhat, 1)), + Array{Complex{Float64}}(size(nhat, 1)), + 0.0, + [Array{Complex{Float64}}(Ns*(2P+1)) for i=1:size(nhat, 1)], + [Array{Complex{Float64}}(Ns*(2P+1)) for i=1:size(nhat, 1)], + [Array{Complex{Float64}}(Ns*(2P+1)) for i=1:size(nhat, 1)], + nhat, + Array{Complex{Float64},2}(size(nhat)), + Array{Complex{Float64}}(Ns*(2P+1)), + 0.0) +end + +#TODO: take k0,ui from opb +function optMatrixPwr(points, centers, Ns, P, k0, ui, iₚ, nhat, l) + sv = PowerBuffer(Ns, P, nhat, iₚ) + cf = 1/(1im*k0*eta0) + pt = Array{Float64}(2) + for ip = 1:size(points,1) + for ic = 1:Ns + pt[1] = points[ip,1] - centers[ic,1] + pt[2] = points[ip,2] - centers[ic,2] + θ = atan2(pt[2], pt[1]) + R = hypot(pt[1], pt[2]) + + ind = (ic-1)*(2*P+1) + P + 1 + Hₚ₋₂ = besselh(-1, 1, k0*R) + Hₚ₋₁ = besselh(0, 1, k0*R) + + sv.HEz[ip][ind] = Hₚ₋₁ + ∂H∂R = k0*Hₚ₋₂ + sv.HHx[ip][ind] = cf*∂H∂R*pt[2]/R + sv.HHy[ip][ind] = -cf*∂H∂R*pt[1]/R + for p = 1:P + Hₚ = (2*(p-1)/(k0*R))*Hₚ₋₁ - Hₚ₋₂ + ∂H∂R = k0*(Hₚ₋₁ - (p/k0/R)*Hₚ) + sv.HEz[ip][ind - p] = (-1)^p*Hₚ*exp(-1im*p*θ) + sv.HHy[ip][ind - p] = -cf*exp(-1im*p*θ)*(-1)^p*(∂H∂R*pt[1] - Hₚ*1im*p*(-pt[2])/R)/R + sv.HHx[ip][ind - p] = cf*exp(-1im*p*θ)*(-1)^p*(∂H∂R*pt[2] - Hₚ*1im*p*(pt[1])/R)/R + + sv.HEz[ip][ind + p] = Hₚ*exp(1im*p*θ) + sv.HHy[ip][ind + p] = -cf*exp(1im*p*θ)*(∂H∂R*pt[1] + Hₚ*1im*p*(-pt[2])/R)/R + sv.HHx[ip][ind + p] = cf*exp(1im*p*θ)*(∂H∂R*pt[2] + Hₚ*1im*p*(pt[1])/R)/R + Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ + end + end + end + sv.l = l + sv.Ez_inc = uinc(k0, points, ui) + sv.Hx_inc = hxinc(k0, points, ui) + sv.Hy_inc = hyinc(k0, points, ui) + sv +end + +mutable struct OptimProblemBuffer#one for each problem, e.g. 3 if same incident wave but 3 different frequencies, or 2 if same frequency but different incident wave + k0::Float64 + kin::Float64 + α::Vector{Complex{Float64}} + ui::Einc + β::Vector{Complex{Float64}} + rhs_grad::Vector{Complex{Float64}} + λadj::Vector{Complex{Float64}} +end + +function OptimProblemBuffer(k0::Float64, kin::Float64, centers::Array{Float64,2}, ui::Einc, P::Integer) + Ns = size(centers, 1) + α = ParticleScattering.u2α(k0, ui, centers, P) + OptimProblemBuffer(k0, kin, α, ui, Array{Complex{Float64}}(Ns*(2*P+1)), + Array{Complex{Float64}}(Ns*(2*P+1)), Array{Complex{Float64}}(Ns*(2*P+1))) +end + +function optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM) + if rs != last_rs + copy!(last_rs, rs) + #do whatever common calculations and save to shared_var + for i in eachindex(opb) + #TODO: if multiple problems have same k0, these calcs are duplicated + for id in unique(ids) + ParticleScattering.updateCircleScatteringDerivative!(scatteringMatrices[i][id], dS_S[i][id], opb[i].k0, opb[i].kin, rs[id], P) + end + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*opb[i].α[rng] + end + + MVP = LinearMap{eltype(fmmbuf.rhs)}( + (output_, x_) -> ParticleScattering.FMM_mainMVP_pre!(output_, x_, + scatteringMatrices[i], φs, ids, P, mFMM[i], + fmmbuf.pre_agg, fmmbuf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + opb[i].β[:] = 0 + opb[i].β,ch = gmres!(opb[i].β, MVP, fmmbuf.rhs, + restart = Ns*(2*P+1), tol = fmmopts.tol, + log = true, initially_zero = true) #no restart + !ch.isconverged && error("""FMM process did not converge, normalized + residual: $(norm(MVP*opb[i].β - fmmbuf.rhs)/norm(fmmbuf.rhs))""") + end + #calculate power for each power calculation arc + calc_multi_pwr!(power_buffer, opb) + end +end + +function optimize_pwr_rs_f(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM, fobj_pwr) + optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM) + + fobj_pwr(power_buffer) +end + +function optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM, gobj_pwr!) + optimize_pwr_rs_common!(rs, last_rs, opb, power_buffer, ids, scatteringMatrices, dS_S, P, Ns, fmmbuf, fmmopts, φs, mFMM) + + fill!(grad_stor, 0) #gradient is the sum of both adjoint gradients + #compute all ∂P/∂βᵀ + dPdβ_pwr!.(power_buffer) + #build rhs of all adjoint problems (⁠-∂f/∂βᵀ) + gobj_pwr!(power_buffer, opb) + + #for each problem, solve adjoint problem and add to total gradient + for i in eachindex(opb) + MVP = LinearMap{eltype(fmmbuf.rhs)}( + (output_, x_) -> ParticleScattering.FMM_mainMVP_transpose!(output_, x_, + scatteringMatrices[i], φs, ids, P, mFMM[i], + fmmbuf.pre_agg, fmmbuf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + #solve adjoint problem + opb[i].λadj[:] = 0 + opb[i].λadj, ch = gmres!(opb[i].λadj, MVP, opb[i].rhs_grad, + restart = Ns*(2*P+1), tol = fmmopts.tol, log = true, + initially_zero = true) #no restart + for n = 1:length(rs) + #compute n-th gradient - here we must pay the price for symmetry + #as more than one β is affected. Overwrites rhs_grad. TODO: + #rhs_grad is still usually sparse - utilize this to reduce complexity + #here O(N^2) -> O(N) + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + if ids[ic] == n + opb[i].rhs_grad[rng] = dS_S[i][n]*opb[i].β[rng] + else + opb[i].rhs_grad[rng] = 0.0 + end + end + grad_stor[n] += -2*real(opb[i].λadj.'*opb[i].rhs_grad) + end + end +end + +function calc_multi_pwr!(power_buffer, opb) + for sv in power_buffer + len = length(sv.HEz) + Sn = Array{Complex{Float64}}(len) # S⋅n + βi = opb[sv.iₚ].β + for ip = 1:len + sv.Ez[ip] = sv.HEz[ip].'*βi + sv.Ez_inc[ip] + sv.Hx[ip] = sv.HHx[ip].'*βi + sv.Hx_inc[ip] + sv.Hy[ip] = sv.HHy[ip].'*βi + sv.Hy_inc[ip] + Sn[ip] = -0.5*real(sv.Ez[ip]*conj(sv.Hy[ip]))*sv.nhat[ip,1] + Sn[ip] += 0.5*real(sv.Ez[ip]*conj(sv.Hx[ip]))*sv.nhat[ip,2] + end + sv.pow = (sum(Sn) - 0.5*Sn[1] - 0.5*Sn[end])*sv.l/(len-1) #trapezoidal rule, correct up to arc length for simple integrals + end +end + +# function calc_multi_pwr_bad_quad!(power_buffer, opb) +# for sv in power_buffer +# len = length(sv.HEz) +# Sn = Array{Complex{Float64}}(len) # S⋅n +# βi = opb[sv.iₚ].β +# for ip = 1:len +# sv.Ez[ip] = sv.HEz[ip].'*βi + sv.Ez_inc[ip] +# sv.Hx[ip] = sv.HHx[ip].'*βi + sv.Hx_inc[ip] +# sv.Hy[ip] = sv.HHy[ip].'*βi + sv.Hy_inc[ip] +# Sn[ip] = -0.5*real(sv.Ez[ip]*conj(sv.Hy[ip]))*sv.nhat[ip,1] +# Sn[ip] += 0.5*real(sv.Ez[ip]*conj(sv.Hx[ip]))*sv.nhat[ip,2] +# end +# sv.pow = sum(Sn)*sv.l/len +# end +# end + +function dPdβ_pwr!(sv::PowerBuffer) + len = length(sv.Ez) #num of points + fill!(sv.∂pow, 0.0) + #this is a sum of real and imaginary parts, hence the additional 1/2 factor + for ip = 1:len + cf = ((ip == 1 || ip == len) ? 0.5 : 1.0)*sv.l/(len-1) #trapezoidal rule constant + sv.∂pow .+= (-0.25*cf*sv.nhat[ip,1])*(conj(sv.Hy[ip])*sv.HEz[ip] + + conj(sv.Ez[ip])*sv.HHy[ip]) + sv.∂pow .+= (0.25*cf*sv.nhat[ip,2])*(conj(sv.Hx[ip])*sv.HEz[ip] + + conj(sv.Ez[ip])*sv.HHx[ip]) + end +end +# +# function dPdβ_pwr_bad_quad!(sv::PowerBuffer) +# len = length(sv.Ez) #num of points +# fill!(sv.∂pow, 0.0) +# #this is a sum of real and imaginary parts, hence the additional 1/2 factor +# for ip = 1:len +# cf = 1.0*sv.l/len +# sv.∂pow .+= (-0.25*cf*sv.nhat[ip,1])*(conj(sv.Hy[ip])*sv.HEz[ip] + +# conj(sv.Ez[ip])*sv.HHy[ip]) +# sv.∂pow .+= (0.25*cf*sv.nhat[ip,2])*(conj(sv.Hx[ip])*sv.HEz[ip] + +# conj(sv.Ez[ip])*sv.HHx[ip]) +# end +# end diff --git a/src/poynting.jl b/src/poynting.jl new file mode 100644 index 0000000..f2d9475 --- /dev/null +++ b/src/poynting.jl @@ -0,0 +1,112 @@ +function poynting_vector(k0, beta, centers, points, Ez_inc, Hx_inc, Hy_inc) + #assumes points lie outside all scattering disks + Ns = size(centers,1) + P = div(div(length(beta),Ns)-1,2) + len = size(points,1) + pyntg = Array{Float64}(len, 2) + Ez = copy(Ez_inc) + for ip in 1:len + ∂Ez∂x = 0.0im + ∂Ez∂y = 0.0im + for ic in 1:Ns + ind = (ic-1)*(2*P+1) + P + 1 + pt1 = points[ip,1] - centers[ic,1] + pt2 = points[ip,2] - centers[ic,2] + + R = hypot(pt1, pt2) + θ = atan2(pt2, pt1) + R == 0 && error("R == 0, center=$(centers[ic,:]), + point=$(points[ip,:]), k0 = $k0, ic=$ic, ip=$ip") + + Hₚ₋₂ = besselh(-1, 1, k0*R) + Hₚ₋₁ = besselh(0, 1, k0*R) + + Ez[ip] += beta[ind]*Hₚ₋₁ + ∂H∂R = k0*Hₚ₋₂ + ∂Ez∂x += beta[ind]*∂H∂R*pt1/R + ∂Ez∂y += beta[ind]*∂H∂R*pt2/R + + for p = 1:P + Hₚ = (2*(p-1)/(k0*R))*Hₚ₋₁ - Hₚ₋₂ + ∂H∂R = k0*(Hₚ₋₁ - p/(k0*R)*Hₚ) + + Ez[ip] += beta[ind - p]*(-1)^p*Hₚ*exp(-1im*p*θ) + ∂Ez∂x += beta[ind - p]*(-1)^p*exp(-1im*p*θ)*(∂H∂R*pt1/R + Hₚ*1im*(-p)*(-pt2)/R^2) + ∂Ez∂y += beta[ind - p]*(-1)^p*exp(-1im*p*θ)*(∂H∂R*pt2/R + Hₚ*1im*(-p)*(pt1)/R^2) + + Ez[ip] += beta[ind + p]*Hₚ*exp(1im*p*θ) + ∂Ez∂x += beta[ind + p]*exp(1im*p*θ)*(∂H∂R*pt1/R + Hₚ*1im*p*(-pt2)/R^2) + ∂Ez∂y += beta[ind + p]*exp(1im*p*θ)*(∂H∂R*pt2/R + Hₚ*1im*p*(pt1)/R^2) + Hₚ₋₂ = Hₚ₋₁; Hₚ₋₁ = Hₚ + end + end + Hx = Hx_inc[ip] + (1/1im/k0/eta0)*∂Ez∂y + Hy = Hy_inc[ip] + (-1/1im/k0/eta0)*∂Ez∂x + pyntg[ip,1] = (-0.5)*real(Ez[ip]*conj(Hy)) + pyntg[ip,2] = 0.5*real(Ez[ip]*conj(Hx)) + end + pyntg +end + +function poynting_vector(k0, points, Ez_inc, Hx_inc, Hy_inc) + #assumes points lie outside all scattering disks + len = size(points,1) + pyntg = Array{Float64}(len, 2) + for ip in 1:len + pyntg[ip,1] = (-0.5)*real(Ez_inc[ip]*conj(Hy_inc[ip])) + pyntg[ip,2] = 0.5*real(Ez_inc[ip]*conj(Hx_inc[ip])) + end + pyntg +end + +function calc_power(k0::Array, kin, P, sp, points, nhat, ui) + #this assumes points are equidistant. correct result only after multiplying by arc length + power = Array{Float64}(length(k0)) + for i in eachindex(k0) + Ez_inc = uinc(k0[i], points, ui) + Hx_inc = hxinc(k0[i], points, ui) + Hy_inc = hyinc(k0[i], points, ui) + if sp == nothing + pyntg = poynting_vector(k0[i], points, Ez_inc, Hx_inc, Hy_inc) + else + beta = solve_particle_scattering(k0[i], kin[i], P, sp, ui; + get_inner = false, verbose = false) + pyntg = poynting_vector(k0[i], beta, sp.centers, points, Ez_inc, Hx_inc, Hy_inc) + end + power[i] = power_quadrature(nhat, pyntg) + end + power +end + +function calc_power(k0::Number, kin, P, sp, points, nhat, ui) + #this assumes points are equidistant. correct result only after multiplying by arc length + Ez_inc = uinc(k0, points, ui) + Hx_inc = hxinc(k0, points, ui) + Hy_inc = hyinc(k0, points, ui) + if sp == nothing + pyntg = poynting_vector(k0, points, Ez_inc, Hx_inc, Hy_inc) + else + beta = solve_particle_scattering(k0, kin, P, sp, ui; + get_inner = false, verbose = false) + pyntg = poynting_vector(k0, beta, sp.centers, points, Ez_inc, Hx_inc, Hy_inc) + end + power_quadrature(nhat, pyntg) +end + +function power_quadrature(n, S) + #trapezoidal rule + len = size(S,1) + (dot(n, S) - 0.5*(n[1,1]*S[1,1] + n[1,2]*S[1,2] + + n[len,1]*S[len,1] + n[len,2]*S[len,2]))/(len - 1) +end + +function rect_border(b, Nx, Ny) + #clockwise from top + points1 = [linspace(b[1], b[2], Nx) b[4]*ones(Nx)] + points2 = [b[2]*ones(Ny) linspace(b[4], b[3], Ny)] + points3 = [linspace(b[2], b[1], Nx) b[3]*ones(Nx)] + points4 = [b[1]*ones(Ny) linspace(b[3], b[4], Ny)] + + n = [[zeros(Nx) ones(Nx)], [ones(Ny) zeros(Ny)], [zeros(Nx) -ones(Nx)], [-ones(Ny) zeros(Ny)]] + [points1, points2, points3, points4], n +end diff --git a/src/shapes.jl b/src/shapes.jl index fbd5915..cf4ab13 100644 --- a/src/shapes.jl +++ b/src/shapes.jl @@ -77,8 +77,8 @@ end """ hex_grid(a::Integer, rows::Integer, d; minus1 = false) -Return `centers`, an `(M,2)` array containing the points on a hexagonal lattice -with horizontal rows, with `a` points in each row and `rows` rows, distanced `d`. +Return `centers`, an `(M,2)` array containing points on a hexagonal lattice +with horizontal rows, with `a` points distanced `d` in each row and `rows` rows. If `minus1` is true, the last point in every odd row is omitted. """ function hex_grid(a::Integer, rows::Integer, d; minus1 = false) diff --git a/src/visualization.jl b/src/visualization.jl index 05477ac..5a2fd2a 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -1,6 +1,6 @@ """ plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc; - opt::FMMoptions = FMMoptions(), use_multipole = true, + opt::FMMoptions = FMMoptions(), method = "multipole", x_points = 201, y_points = 201, border = find_border(sp), normalize = 1.0) @@ -10,9 +10,10 @@ TM field `ui` scattering from the ScatteringProblem `sp`, using matplotlib's bounding box or calculate automatically. Uses the FMM options given by `opt` (FMM is disabled by default); -`use_multipole` dictates whether electric field is calculated using the -multipole/cylindrical harmonics (true) or falls back on potential densities -(false). Either way, the multiple-scattering system is solved in the cylindrical +`method = "multipole"` dictates whether electric field is calculated using the +multipole/cylindrical harmonics, uses a faster but less accurate Hankel recurrence +formula (`"recurrence"`), or falls back on potential densities (`"density"`). +Either way, the multiple-scattering system is solved in the cylindrical harmonics space. Normalizes all distances and sizes in plot (but not output) by `normalize`. @@ -22,7 +23,7 @@ Returns the calculated field in two formats: contains the field at `(mean(xgrid[i, j:j+1]), mean(ygrid[i:i+1, j]))`. """ function plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc; - opt::FMMoptions = FMMoptions(), use_multipole = true, + opt::FMMoptions = FMMoptions(), method = "multipole", x_points = 201, y_points = 201, border = find_border(sp), normalize = 1.0) @@ -40,7 +41,7 @@ function plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc; vec(ygrid[1:y_points, 1:x_points]) + dy) Ez = calc_near_field(k0, kin, P, sp, points, ui, - use_multipole = use_multipole, opt = opt) + method = method, opt = opt) zgrid = reshape(Ez, y_points, x_points) figure() pcolormesh(xgrid/normalize, ygrid/normalize, abs.(zgrid)) @@ -127,21 +128,22 @@ draw_shapes(sp; ax = gca(), normalize = 1.0) = draw_shapes(sp.shapes, """ calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc; - opt::FMMoptions = FMMoptions(), use_multipole = true, + opt::FMMoptions = FMMoptions(), method = "multipole", verbose = true) Calculates the total electric field as a result of a plane wave with incident field `ui` scattering from the ScatteringProblem `sp`, at `points`. -Uses the FMM options given by `opt` (default behavious is disabled FMM); -`use_multipole` dictates whether electric field is calculated using the -multipole/cylindrical harmonics (true) or falls back on potential densities -(false). Either way, the multiple-scattering system is solved in the cylindrical +Uses the FMM options given by `opt` (default behaviour is disabled FMM); +`method = "multipole"` dictates whether electric field is calculated using the +multipole/cylindrical harmonics, uses a faster but less accurate Hankel recurrence +formula (`"recurrence"`), or falls back on potential densities (`"density"`). +Either way, the multiple-scattering system is solved in the cylindrical harmonics space, and the field by a particular scatterer inside its own scattering discs is calculated by potential densities, as the cylindrical harmonics approximation is not valid there. """ function calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc; - opt::FMMoptions = FMMoptions(), use_multipole = true, + opt::FMMoptions = FMMoptions(), method = "multipole", verbose = true) shapes = sp.shapes; ids = sp.ids; centers = sp.centers; φs = sp.φs @@ -199,9 +201,12 @@ function calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc; u[rng_out] += uinc(k0, points[rng_out,:], ui) for ic2 = 1:size(sp) ic == ic2 && continue - if use_multipole + if method == "multipole" scattered_field_multipole!(u, k0, beta, P, centers, ic2, points, find(rng_out)) + elseif method == "recurrence" + scattered_field_multipole_recurrence!(u, k0, beta, P, centers, ic2, points, + find(rng_out)) else if φs[ic2] == 0.0 ft_rot2 = shapes[ids[ic2]].ft .+ centers[ic2,:]' @@ -224,8 +229,10 @@ function calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc; rng = (tags .== 0) #incident field u[rng] = uinc(k0, points[rng,:], ui) - if use_multipole + if method == "multipole" scattered_field_multipole!(u, k0, beta, P, centers, 1:size(sp), points, find(rng)) + elseif method == "recurrence" + scattered_field_multipole_recurrence!(u, k0, beta, P, centers, 1:size(sp), points, find(rng)) else for ic = 1:size(centers,1) if typeof(shapes[ids[ic]]) == ShapeParams From 5a02b620f2e0b1f67fe6cc228ade7176b70ae8c2 Mon Sep 17 00:00:00 2001 From: Unknown <34473471+bblankrot@users.noreply.github.com> Date: Fri, 15 Mar 2019 12:04:23 +0100 Subject: [PATCH 7/9] initial power optimization commit --- src/PS_types.jl | 2 +- src/ParticleScattering.jl | 2 + src/optimize_phis.jl | 2 +- src/optimize_phis_pwr.jl | 90 +++++++++++++++++++++++++++++++++++++++ src/optimize_rs_pwr.jl | 32 +------------- 5 files changed, 96 insertions(+), 32 deletions(-) create mode 100644 src/optimize_phis_pwr.jl diff --git a/src/PS_types.jl b/src/PS_types.jl index 942a55f..7706652 100644 --- a/src/PS_types.jl +++ b/src/PS_types.jl @@ -191,7 +191,7 @@ end """ LineSource(x, y) -Constructor for the `LineSource` type, where `(x0,y0)` is the coordinate of the +Constructor for the `LineSource` type, where `(x,y)` is the coordinate of the current filament. """ struct LineSource <: Einc diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl index b0ccafe..81a5326 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -73,4 +73,6 @@ export calc_power include("optimize_rs_pwr.jl") export PowerBuffer, OptimProblemBuffer, optMatrixPwr, optimize_pwr_rs_f, optimize_pwr_rs_g! +include("optimize_phis_pwr.jl") +export optimize_pwr_φ_f, optimize_pwr_φ_g! end diff --git a/src/optimize_phis.jl b/src/optimize_phis.jl index 325c852..a6204e3 100644 --- a/src/optimize_phis.jl +++ b/src/optimize_phis.jl @@ -184,7 +184,7 @@ function optimize_φ_adj_g!(grad_stor, φs, shared_var, last_φs, α, H, points, v[:] += D.*shared_var.β[rng] grad_stor[n] = ifelse(minimize, -2, 2)*real(view(λadj,rng).'*v) - #prepare for next one + #prepare for next one - #TODO: check why this is here v[:] = 0.0im end end diff --git a/src/optimize_phis_pwr.jl b/src/optimize_phis_pwr.jl new file mode 100644 index 0000000..9e94707 --- /dev/null +++ b/src/optimize_phis_pwr.jl @@ -0,0 +1,90 @@ +#for now just the main functions, no wrapper +function optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM) + if φs != last_φs + copy!(last_φs, φs) + #do whatever common calculations and save to shared_var + for i in eachindex(opb) + for ic = 1:Ns + rng = (ic-1)*(2*P+1) + (1:2*P+1) + if φs[ic] == 0.0 + fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*opb[i].α[rng] + else + #rotate without matrix + rotateMultipole!(view(fmmbuf.rhs,rng), view(opb[i].α,rng), -φs[ic], P) + fmmbuf.rhs[rng] = scatteringMatrices[i][ids[ic]]*fmmbuf.rhs[rng] + rotateMultipole!(view(fmmbuf.rhs,rng), φs[ic], P) + end + end + + MVP = LinearMap{eltype(fmmbuf.rhs)}( + (output_, x_) -> ParticleScattering.FMM_mainMVP_pre!(output_, x_, + scatteringMatrices[i], φs, ids, P, mFMM[i], + fmmbuf.pre_agg, fmmbuf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + opb[i].β[:] = 0 + opb[i].β,ch = gmres!(opb[i].β, MVP, fmmbuf.rhs, + restart = Ns*(2*P+1), tol = fmmopts.tol, + log = true, initially_zero = true) #no restart + !ch.isconverged && error("""FMM process did not converge, normalized + residual: $(norm(MVP*opb[i].β - fmmbuf.rhs)/norm(fmmbuf.rhs))""") + end + #calculate power for each power calculation arc + calc_multi_pwr!(power_buffer, opb) + end +end + +function optimize_pwr_φ_f(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM, fobj_pwr) + optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM) + + fobj_pwr(power_buffer) +end + +function optimize_pwr_φ_g!(grad_stor, φs, last_φs, opb, power_buffer, ids, scatteringMatrices, scatteringLU, P, Ns, fmmbuf, fmmopts, mFMM, gobj_pwr!) + optimize_pwr_φ_common!(φs, last_φs, opb, power_buffer, ids, scatteringMatrices, P, Ns, fmmbuf, fmmopts, mFMM) + + fill!(grad_stor, 0) #gradient is the sum of all adjoint gradients + #compute all ∂P/∂βᵀ + dPdβ_pwr!.(power_buffer) + #build rhs of all adjoint problems (⁠-∂f/∂βᵀ) + gobj_pwr!(power_buffer, opb) + + #for each problem, solve adjoint problem and add to total gradient + for i in eachindex(opb) + MVP = LinearMap{eltype(fmmbuf.rhs)}( + (output_, x_) -> ParticleScattering.FMM_mainMVP_transpose!(output_, x_, + scatteringMatrices[i], φs, ids, P, mFMM[i], + fmmbuf.pre_agg, fmmbuf.trans), + Ns*(2*P+1), Ns*(2*P+1), ismutating = true) + + #solve adjoint problem + opb[i].λadj[:] = 0 + opb[i].λadj, ch = gmres!(opb[i].λadj, MVP, opb[i].rhs_grad, + restart = Ns*(2*P+1), tol = fmmopts.tol, + log=true, initially_zero = true) #no restart + !ch.isconverged && error("""FMM process did not converge for adjoint, + normalized residual: + $(norm(MVP*opb[i].λadj - opb[i].rhs_grad)/norm(opb[i].rhs_grad))""") + D = -1.0im*collect(-P:P) + v = Array{Complex{Float64}}(2*P+1) #TODO: minimize dynamic alloc + for n = 1:Ns + rng = (n-1)*(2*P+1) + (1:2*P+1) + # rotateMultipole!(v, view(opb[i].β,rng), -φs[n], P) + # v[:] = scatteringLU[i][ids[n]]\v #LU decomp with pivoting + # v[:] .*= -D + # v[:] = scatteringMatrices[i][ids[n]]*v + # rotateMultipole!(v, φs[n], P) + # v[:] += D.*opb[i].β[rng] + # grad_stor[n] += (-2)*real(view(opb[i].λadj,rng).'*v) + # #prepare for next one - #TODO: check why this is here + # v[:] = 0.0im + + #dumb way + Φ = spdiagm(exp.(-1.0im*(-P:P)*φs[n])) + temppp = Φ*(scatteringLU[i][ids[n]]\(conj(Φ)*opb[i].β[rng])) + temppp2 = spdiagm(D)*(Φ*(scatteringMatrices[i][ids[n]]*(conj(Φ)*temppp))) - + Φ*(scatteringMatrices[i][ids[n]]*(conj(Φ)*(spdiagm(D)*temppp))) + grad_stor[n] += (-2)*real(opb[i].λadj[rng].'*temppp2) + end + end +end diff --git a/src/optimize_rs_pwr.jl b/src/optimize_rs_pwr.jl index c1f01ab..c0ff044 100644 --- a/src/optimize_rs_pwr.jl +++ b/src/optimize_rs_pwr.jl @@ -154,7 +154,8 @@ function optimize_pwr_rs_g!(grad_stor, rs, last_rs, opb, power_buffer, ids, scat initially_zero = true) #no restart for n = 1:length(rs) #compute n-th gradient - here we must pay the price for symmetry - #as more than one β is affected. Overwrites rhs_grad. TODO: + #as more than one β is affected. Overwrites rhs_grad (ok because + #there is seperate one for each i). TODO: #rhs_grad is still usually sparse - utilize this to reduce complexity #here O(N^2) -> O(N) for ic = 1:Ns @@ -186,22 +187,6 @@ function calc_multi_pwr!(power_buffer, opb) end end -# function calc_multi_pwr_bad_quad!(power_buffer, opb) -# for sv in power_buffer -# len = length(sv.HEz) -# Sn = Array{Complex{Float64}}(len) # S⋅n -# βi = opb[sv.iₚ].β -# for ip = 1:len -# sv.Ez[ip] = sv.HEz[ip].'*βi + sv.Ez_inc[ip] -# sv.Hx[ip] = sv.HHx[ip].'*βi + sv.Hx_inc[ip] -# sv.Hy[ip] = sv.HHy[ip].'*βi + sv.Hy_inc[ip] -# Sn[ip] = -0.5*real(sv.Ez[ip]*conj(sv.Hy[ip]))*sv.nhat[ip,1] -# Sn[ip] += 0.5*real(sv.Ez[ip]*conj(sv.Hx[ip]))*sv.nhat[ip,2] -# end -# sv.pow = sum(Sn)*sv.l/len -# end -# end - function dPdβ_pwr!(sv::PowerBuffer) len = length(sv.Ez) #num of points fill!(sv.∂pow, 0.0) @@ -214,16 +199,3 @@ function dPdβ_pwr!(sv::PowerBuffer) conj(sv.Ez[ip])*sv.HHx[ip]) end end -# -# function dPdβ_pwr_bad_quad!(sv::PowerBuffer) -# len = length(sv.Ez) #num of points -# fill!(sv.∂pow, 0.0) -# #this is a sum of real and imaginary parts, hence the additional 1/2 factor -# for ip = 1:len -# cf = 1.0*sv.l/len -# sv.∂pow .+= (-0.25*cf*sv.nhat[ip,1])*(conj(sv.Hy[ip])*sv.HEz[ip] + -# conj(sv.Ez[ip])*sv.HHy[ip]) -# sv.∂pow .+= (0.25*cf*sv.nhat[ip,2])*(conj(sv.Hx[ip])*sv.HEz[ip] + -# conj(sv.Ez[ip])*sv.HHx[ip]) -# end -# end From 58bd7f2e6428be3a7ff0b9e8d9b93aeb3ec3664f Mon Sep 17 00:00:00 2001 From: Boaz Blankrot <34473471+bblankrot@users.noreply.github.com> Date: Thu, 23 May 2019 14:15:33 +0200 Subject: [PATCH 8/9] Removed PGFPlotsX support Built-in plotting now only though PyPlot --- src/ParticleScattering.jl | 6 --- src/visualization_pgf.jl | 91 --------------------------------------- 2 files changed, 97 deletions(-) delete mode 100644 src/visualization_pgf.jl diff --git a/src/ParticleScattering.jl b/src/ParticleScattering.jl index 81a5326..5b271e3 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -55,12 +55,6 @@ export uinc, hxinc, hyinc #consts, incident.jl export eta0 -# For advanced plotting with pgfplots -import DataFrames, CSV, PGFPlotsX; const pgf = PGFPlotsX -include("visualization_pgf.jl") -#methods, visualization_pgf.jl -export plot_near_field_pgf - #temp export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices export FMMbuffer, FMMmatrices diff --git a/src/visualization_pgf.jl b/src/visualization_pgf.jl deleted file mode 100644 index f8302bc..0000000 --- a/src/visualization_pgf.jl +++ /dev/null @@ -1,91 +0,0 @@ -""" - plot_near_field_pgf(filename, k0, kin, P, sp::ScatteringProblem, ui::Einc; - opt::FMMoptions = FMMoptions(), use_multipole = true, - x_points = 201, y_points = 201, border = find_border(sp), - downsample = 1, include_preamble = false, normalize = 1.0) - -Plots the total electric field as a result of a plane wave with incident TM -field `ui` scattering from the ScatteringProblem `sp`, using pgfplots's -`surf`. Can accept number of sampling points in each direction, and either -a given `border` or calculate it automatically. The plots of the shapes (but not -the field) can be downsampled by setting an integer `downsample`, since pgfplots -slows down dramatically when drawing many shapes with many nodes. - -Uses the FMM options given by `opt` (FMM is disabled by default); -`use_multipole` dictates whether electric field is calculated using the -multipole/cylindrical harmonics (true) or falls back on potential densities -(false). Either way, the multiple-scattering system is solved in the cylindrical -harmonics space. Normalizes all distances and sizes in plot by `normalize`. - -Saves the generated pgfplots file to `filename`, with just a surrounding `tikzpicture` -environment if `include_preamble=false`, and a compilable tandalone document -otherwise. -""" -function plot_near_field_pgf(filename, k0, kin, P, sp::ScatteringProblem, ui::Einc; - opt::FMMoptions = FMMoptions(), use_multipole = true, - x_points = 201, y_points = 201, border = find_border(sp), - downsample = 1, include_preamble = false, normalize = 1.0) - #important difference: here function is measured on the grid, unlike pcolormesh. - shapes = sp.shapes; ids = sp.ids; centers = sp.centers; φs = sp.φs - x_min, x_max, y_min, y_max = border - - x = linspace(x_min, x_max, x_points) - y = linspace(y_min, y_max, y_points) - - xgrid = repmat(x',y_points,1) - ygrid = repmat(y,1,x_points) - points = [xgrid[:] ygrid[:]] - - Ez = calc_near_field(k0, kin, P, sp, points, ui, - use_multipole=use_multipole, opt = opt) - - aspect = (y_max - y_min)/(x_max - x_min) - - dt = DataFrames.DataFrame(x = points[:,1]/normalize, y = points[:,2]/normalize, z = abs.(Ez)) - CSV.write(filename * ".dat", dt, delim = '\t') - pgf.@pgf begin - ax = pgf.Axis({ xmin = x_min/normalize, - xmax = x_max/normalize, - ymin = y_min/normalize, - ymax = y_max/normalize, - xlabel = "\$x/$normalize\$", - ylabel = "\$y/$normalize\$", - scale_only_axis, - width = "\\linewidth", - height = "$aspect*\\linewidth", - view = "{0}{90}", - "mesh/ordering" = "y varies", - colorbar, - point_meta_max = maximum(abs.(Ez))}) - #new release of pgfplotsx isn't printing table correctly (too many newlines for pgfplots) - push!(ax, "\\addplot3[surf, no markers, shader={interp}, - mesh/rows={$(y_points)}] table{$(basename(filename)).dat};") - - end - draw_shapes_pgf(shapes, ids, centers, φs, ax, 1.1*maximum(abs.(Ez)), downsample, normalize = normalize) - pgf.save(filename, ax ,include_preamble = include_preamble) -end - -function draw_shapes_pgf(shapes, ids, centers, φs, ax, floating, downsample; normalize = 1.0) - #draw in 3d so it is "above" surf plot - for ic = 1:size(centers,1) - if typeof(shapes[ids[ic]]) == ShapeParams - if φs[ic] == 0.0 - ft_rot = shapes[ids[ic]].ft[1:downsample:end,:] .+ centers[ic,:]' - else - Rot = [cos(φs[ic]) sin(φs[ic]);-sin(φs[ic]) cos(φs[ic])] - ft_rot = shapes[ids[ic]].ft[1:downsample:end,:]*Rot .+ centers[ic,:]' - end - co = pgf.Coordinates([ft_rot[:,1];ft_rot[1,1]]/normalize, - [ft_rot[:,2];ft_rot[1,2]]/normalize, - floating*ones(size(ft_rot,1) + 1)) - pgf.@pgf push!(ax, pgf.Plot3({black, no_markers, thick}, co)) - else - x = centers[ic,1]/normalize - y = centers[ic,2]/normalize - R = shapes[ids[ic]].R/normalize - push!(ax, "\\addplot3 [black, thick, domain=0:2*pi,samples=100] - ({$x+$R*cos(deg(x))},{$y+$R*sin(deg(x))},$floating);") - end - end -end From b8670caf10d5b5eb8bc224706df33a2004ad5e48 Mon Sep 17 00:00:00 2001 From: Boaz Blankrot <34473471+bblankrot@users.noreply.github.com> Date: Thu, 23 May 2019 14:15:33 +0200 Subject: [PATCH 9/9] added recurrence option to calc_far_field and test --- src/visualization.jl | 16 ++++++++++------ test/fmm_test.jl | 4 ++-- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/visualization.jl b/src/visualization.jl index 5a2fd2a..9420d1f 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -57,14 +57,14 @@ end """ plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave; - opt::FMMoptions = FMMoptions(), use_multipole = true, + opt::FMMoptions = FMMoptions(), method = "multipole", plot_points = 200) Plots and returns the echo width (radar cross section in two dimensions) for a -given scattering problem. `opt`, `use_multipole` are as in `plot_near_field`. +given scattering problem. `opt`, `method` are as in `plot_near_field`. """ function plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave; - opt::FMMoptions = FMMoptions(), use_multipole = true, + opt::FMMoptions = FMMoptions(), method = "multipole", plot_points = 200) Rmax = maximum(s.R for s in sp.shapes) @@ -81,7 +81,7 @@ function plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave; points = [x_far y_far] Ez = calc_far_field(k0, kin, P, points, sp, pw, - use_multipole = use_multipole, opt = opt) + method = method, opt = opt) Ez[:] = (k0*Rfar)*abs2.(Ez) #plot echo width figure() @@ -264,7 +264,7 @@ function calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc; end function calc_far_field(k0, kin, P, points, sp::ScatteringProblem, pw::PlaneWave; - opt::FMMoptions = FMMoptions(), use_multipole = true) + opt::FMMoptions = FMMoptions(), method = "multipole") #calc only scattered field + assumes all points are outside shapes shapes = sp.shapes; centers = sp.centers; ids = sp.ids; φs = sp.φs if opt.FMM @@ -278,9 +278,13 @@ function calc_far_field(k0, kin, P, points, sp::ScatteringProblem, pw::PlaneWave beta, sigma_mu = solve_particle_scattering(k0, kin, P, sp, pw) end Ez = zeros(Complex{Float64}, size(points,1)) - if use_multipole + + if method == "multipole" scattered_field_multipole!(Ez, k0, beta, P, centers, 1:size(sp), points, 1:size(points,1)) + elseif method == "recurrence" + scattered_field_multipole_recurrence!(Ez, k0, beta, P, centers, + 1:size(sp), points, 1:size(points,1)) else for ic = 1:size(sp) if typeof(shapes[ids[ic]]) == ShapeParams diff --git a/test/fmm_test.jl b/test/fmm_test.jl index fb06213..4771c9d 100644 --- a/test/fmm_test.jl +++ b/test/fmm_test.jl @@ -43,12 +43,12 @@ u1 = calc_near_field(k0, kin, P, sp, points, PlaneWave(θ_i); opt = fmm_options, verbose = false) u2 = calc_near_field(k0, kin, P, sp, points, PlaneWave(θ_i); opt = fmm_options, - verbose = true, use_multipole = false) + verbose = true, method = "density") @test norm(u1 - u2)/norm(u1) < 1e-6 u3 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0); opt = fmm_options, verbose = false) u4 = calc_near_field(k0, kin, P, sp, points, LineSource(-0.8λ0,-λ0); opt = fmm_options, - verbose = true, use_multipole = false) + verbose = true, method = "recurrence") @test norm(u3 - u4)/norm(u3) < 1e-6 end