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 daa03a4..5b271e3 100644 --- a/src/ParticleScattering.jl +++ b/src/ParticleScattering.jl @@ -51,16 +51,22 @@ export optimize_φ #methods, optimize_rs.jl export optimize_radius #methods, incident.jl -export uinc - -# 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 +export uinc, hxinc, hyinc +#consts, incident.jl +export eta0 #temp export divideSpace, FMMtruncation, particleExpansion, FMMbuildMatrices 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! +include("optimize_phis_pwr.jl") +export optimize_pwr_φ_f, optimize_pwr_φ_g! end diff --git a/src/fmm_mvp.jl b/src/fmm_mvp.jl index 42244a9..8e74529 100644 --- a/src/fmm_mvp.jl +++ b/src/fmm_mvp.jl @@ -137,3 +137,72 @@ 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)) + 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) + if φs[ic] == 0.0 + At_mul_B!(v, scatteringMatrices[ids[ic]], view(beta, rng)) + else + #rotate without matrix - transposed + rotateMultipole!(temp, view(beta, rng), φs[ic], P) + At_mul_B!(v, scatteringMatrices[ids[ic]], temp) + 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 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_phis.jl b/src/optimize_phis.jl index 2122aa3..a6204e3 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_φ_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 - #TODO: check why this is here + 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_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.jl b/src/optimize_rs.jl index 655155b..d82da93 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) @@ -32,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) @@ -54,29 +56,24 @@ 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 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,9 +176,46 @@ 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 updateCircleScatteringDerivative!(S, dS_S, kout, kin, R, P) - #non-vectorized, reuses bessel +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::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) @@ -210,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 diff --git a/src/optimize_rs_mf.jl b/src/optimize_rs_mf.jl new file mode 100644 index 0000000..a2a9a1f --- /dev/null +++ b/src/optimize_rs_mf.jl @@ -0,0 +1,163 @@ +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) + 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) + 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 diff --git a/src/optimize_rs_pwr.jl b/src/optimize_rs_pwr.jl new file mode 100644 index 0000000..c0ff044 --- /dev/null +++ b/src/optimize_rs_pwr.jl @@ -0,0 +1,201 @@ +#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 (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 + 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 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 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..9420d1f 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)) @@ -56,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) @@ -80,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() @@ -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 @@ -257,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 @@ -271,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/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 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