rederive uniform source local error

scripts/choose_p.jl

@@ -30,7 +30,7 @@ function generate_multipole(center, ρs, θs, ϕs, qs, expansion_order)
     system = build_system(center, ρs, θs, ϕs, qs)
 
     # build branch
-    branch = Branch(1:1, 0, 1:0, 0, 1, center, 0.0, 0.0, SVector{6}(0.0,0,0,0,0,0), SVector{3}(0.0,0,0), expansion_order)
+    branch = Branch(1:length(ρs), 0, 1:0, 0, 1, center, 0.0, 0.0, SVector{6}(0.0,0,0,0,0,0), SVector{3}(0.0,0,0), expansion_order)
 
     # multipole coefficients
     body_to_multipole!(branch, system, branch.harmonics, Val(expansion_order))
@@ -120,19 +120,43 @@ function solve_p_new(ε, ρ, am, al)
     return 51
 end
 
+function solve_p_new_multipole(ε, ρ, am, al)
+    target = ε * (ρ-am-al)
+    f1 = am / (ρ-al)
+    t1 = 3 * am
+    for p in 0:50
+        this_target = target * (2*p+8)# * sqrt(4*p+6)
+        t1 < this_target && (return p)
+        t1 *= f1
+    end
+    return 51
+end
+
 function solve_p_old(ε, ρ, am)
     c = ρ/am - 1.0
     p_old = Int(ceil(-log(c,ε)))
     return p_old
 end
 
-function choose_p(epsilon::Number, dx=0.0, pmax=51)
+function solve_p_old_multipole(ε, ρ, am, al)
+    target = ε * (ρ - am - al)
+    f2 = am / (ρ-al)
+    t2 = am
+    for p in 0:50
+        t2 < target && (return p)
+        t2 *= f2
+    end
+    return 51
+end
+
+function choose_p(epsilon::Number, dx=0.0, pmax=51; multipole_only=false,
+        ρs = [1.0],
+        θs = [π/2],
+        ϕs = [0.0],
+        qs = [1.0],
+    )
     # define multipole (source)
     multipole_center = SVector{3}(0.0,0,0)
-    ρs = [1.0]
-    θs = [π/2]
-    ϕs = [0.0]
-    qs = [1.0]
     multipole, system = generate_multipole(multipole_center, ρs, θs, ϕs, qs, pmax)
 
     # define target
@@ -144,12 +168,20 @@ function choose_p(epsilon::Number, dx=0.0, pmax=51)
 
     # determine expansion order
     ρ = norm(local_center - multipole.center)
+    al = 1.0 - dx
     am = ρs[1]
-    p_old = solve_p_old(epsilon, ρ, am)
+    if multipole_only
+        p_old = solve_p_old_multipole(epsilon, ρ, am, al)
+    else
+        p_old = solve_p_old(epsilon, ρ, am)
+    end
 
     # use new error function
-    al = 1.0 - dx
-    p_new = solve_p_new(epsilon, ρ, am, al)
+    if multipole_only
+        p_new = solve_p_new_multipole(epsilon, ρ, am, al)
+    else
+        p_new = solve_p_new(epsilon, ρ, am, al)
+    end
 
     # potential
     u_multipole = evaluate_multipole(xt, multipole, p_old)
@@ -163,15 +195,20 @@ function choose_p(epsilon::Number, dx=0.0, pmax=51)
     return abs((u_multipole - u_direct)/u_direct), abs((u_local - u_direct)/u_direct), p_old, abs((u_multipole_new - u_direct)/u_direct), abs((u_local_new - u_direct)/u_direct), p_new
 end
 
-function choose_p(εs::Vector, dx=0.0, pmax=50)
+function choose_p(εs::Vector, dx=0.0, pmax=51; multipole_only=false,
+        ρs = [1.0],
+        θs = [π/2],
+        ϕs = [0.0],
+        qs = [1.0],
+    )
     multipole_error = zeros(length(εs))
     local_error = zeros(length(εs))
     new_multipole_error = zeros(length(εs))
     new_local_error = zeros(length(εs))
     old_ps = zeros(Int,length(εs))
     new_ps = zeros(Int,length(εs))
     for (i,ε) in enumerate(εs)
-        m, l, p_old, m_new, l_new, p_new = choose_p(ε, dx, pmax)
+        m, l, p_old, m_new, l_new, p_new = choose_p(ε, dx, pmax; multipole_only, ρs, θs, ϕs, qs)
         multipole_error[i] = m
         local_error[i] = l
         old_ps[i] = p_old
@@ -182,6 +219,48 @@ function choose_p(εs::Vector, dx=0.0, pmax=50)
     return multipole_error, local_error, old_ps, new_multipole_error, new_local_error, new_ps
 end
 
+function uniform_sphere(nx, A)
+    ρs = Float64[]
+    θs = Float64[]
+    ϕs = Float64[]
+    for x in range(-1, stop=1.0, length= nx)
+        for y in range(-1, stop=1.0, length=nx)
+            for z in range(-1, stop=1.0, length=nx)
+                ρ = x*x+y*y+z*z
+                if ρ <= 1.0
+                    r, θ, ϕ = FastMultipole.cartesian_to_spherical(SVector{3}(x,y,z))
+                    push!(ρs, r)
+                    push!(θs, θ)
+                    push!(ϕs, ϕ)
+                end
+            end
+        end
+    end
+    qs = zeros(length(ρs))
+    qs .= A / length(qs)
+    return ρs, θs, ϕs, qs
+end
+
+function uniform_box(nx, A)
+    ρs = Float64[]
+    θs = Float64[]
+    ϕs = Float64[]
+    for x in range(-1, stop=1.0, length= nx)
+        for y in range(-1, stop=1.0, length=nx)
+            for z in range(-1, stop=1.0, length=nx)
+                r, θ, ϕ = FastMultipole.cartesian_to_spherical(SVector{3}(x,y,z))
+                push!(ρs, r)
+                push!(θs, θ)
+                push!(ϕs, ϕ)
+            end
+        end
+    end
+    qs = zeros(length(ρs))
+    qs .= A / length(qs)
+    return ρs, θs, ϕs, qs
+end
+
+#=
 fig = figure("choose p")
 fig.clear()
 fig.add_subplot(211, xlabel="relative error tolerance", ylabel="resulting error")
@@ -211,3 +290,38 @@ tight_layout()
 
 BSON.@save "choose_p2_d$distance.bson" m l old_ps m_new l_new new_ps
 savefig("choose_p2_d$distance.png")
+=#
+
+fig = figure("choose p multipole")
+fig.clear()
+fig.add_subplot(121, xlabel="relative error tolerance", ylabel="multipole error")
+fig.add_subplot(122, xlabel="relative error tolerance", ylabel="expansion order")
+ax1 = fig.get_axes()[0]
+ax0 = fig.get_axes()[1]
+
+distance = -0.25
+eps = [10.0^x for x in -12:0.005:0]
+
+nx = 5
+A = 1.0
+ρs, θs, ϕs, qs = uniform_sphere(nx, A)
+#ρs, θs, ϕs, qs = uniform_box(nx, A)
+m, l, old_ps, m_new, l_new, new_ps = choose_p(eps, distance; multipole_only=true, ρs, θs, ϕs, qs)
+
+ax0.plot(eps, new_ps)
+ax0.plot(eps, old_ps)
+ax1.plot(eps, m_new)
+ax1.plot(eps, m)
+ax1.plot(eps, eps, "red")
+
+ax0.set_xscale("log")
+ax0.legend(["uniform source", "worst case"])
+
+ax1.set_xscale("log")
+ax1.set_yscale("log")
+ax1.legend(["uniform source", "worst case", "tolerance"])
+
+tight_layout()
+
+#BSON.@save "choose_p2_d$distance.bson" m l old_ps m_new l_new new_ps
+#savefig("choose_p2_d$distance.png")