# NDSolveValue for Laplace equation not converging to analytic solution

I'm solving Laplace equation $$\nabla^2 \phi = 0$$ with BC's $$\phi_x(x=\pm 1) = 0,\, \phi_y(y=-h) = 0$$ with a specified BC along the circular arc $$x^2 + (-1 + y)^2 = 4$$, which I call $$\Gamma$$ (so the domain $$\Omega$$ is like a rectangular box with the circular arc $$\Gamma$$ as the lid).

What I've done is taken the analytic solution $$\phi_A = \cos(\pi m (x + 1)/2) \cosh(\pi m (y + h)/2) : m \in \mathbb{N}$$ and evaluated its normal derivative on $$\Gamma$$, and then set this as a Neumann boundary condition for the top lid. Since this does not give a unique solution (no Dirichlet BC) I also give the PDE the BC $$\phi = \phi_A(1,\Gamma(x=1))$$. However, the numerical solution does not converge to the analytic solution. Code is supplied below, which shows the domain $$\Omega$$ for clarity.

h = 2;
r = 2;
yp = 1;
m = 1;
Ω = ImplicitRegion[-1 <= x <= 1 && -h <= y <= 2 && ! (x^2 + (y - yp)^2 <= r^2), {x, y}];
Region[Ω]
ϕA = Cos[(π m)/2 (x + 1)] Cosh[(π m)/2 (y + h)];
g = ϕA /. y -> -Sqrt[r^2 - x^2] + yp /. x -> 1;
dnϕA = Grad[ϕA, {x, y}]. ({x, y}/Sqrt[x^2 + y^2]);
dnϕΓA = dnϕA /. y -> -Sqrt[r^2 - x^2] + yp;
op = Laplacian[ϕD[x, y], {x, y}];
ΓNV = NeumannValue[dnϕA, x^2 + (y - yp)^2 == r^2];
Γ = {DirichletCondition[ϕD[x, y] == g, x == 1 && y == -Sqrt[r^2 - 1^2] + yp]};
Needs["NDSolveFEM"]
mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.01];
mesh["Wireframe"]
ϕ = NDSolveValue[{op == ΓNV, Γ},
ϕD, {x, y} ∈ mesh, InterpolationOrder -> 3];
Plot3D[{ϕ[x, y], ϕA}, {x, y} ∈ Ω,
AxesLabel -> {"x", "y", "ϕ"}]

• I think your Neumann value is wrong: since the circle lid is not centered around (0,0), the unit normal in dn\[Phi]A is not (x,y)/norm(...) but (x,y-1)/norm(...), so dn\[Phi]A = Grad[\[Phi]A, {x, y}].(Normalize[{x, y - 1}]); – egwene sedai Jan 14 '19 at 22:22
• @egwenesedai More precisely dn\[Phi]A = Grad[\[Phi]A, {x, y}].({x, y - yp}/Sqrt[x^2 + (y - yp)^2]); – Alex Trounev Jan 15 '19 at 13:50
• Thank you egwene sedai and Alex Trounev. I'm new here and unsure how to upvote your comments, but I agree with you both. (Also, a little unclear if I should be thanking each of you...) – Josh McCraney Jan 15 '19 at 21:13
• @JoshMcCraney: something that AFAIK is very welcome on this site in cases where noone answers but people provide useful information in comments is to answer your own question with the information given by others in comments (which also is a good place to mention their help/input)... – Albert Retey Jan 15 '19 at 22:28

Edit: my normal derivative to the circular arc $$\Gamma$$ assumed the center was the origin. After recognizing the center was intentionally above the origin, the normal vector should be changed to
({x, y - yp}/Sqrt[x^2 + (y - yp)^2])