# Laplace's equation in spherical coordinates with Neumann b.c

I am trying to find the temperature field in a semi-infinite solid on whose surface there is an isotherm spherical cap sunken by a length p. For example:

R = 10;
p = 3;
SphericalPlot3D[
1/2 (Sqrt[2] Sqrt[Cos[2 θ] (R - p)^2 + 2 p R + R^2 - p^2] +
2 Cos[θ] (R - p)), {θ, Pi/2, 3/2 Pi}, {ϕ, 0,
2 Pi}]


The rest of the surface is adiabatic.

In a spherical coordinate system (physical convention) centered at the "center" of the cap, the PDE is: $$\frac{\partial }{\partial r}\left(r^2\frac{\partial T}{\partial r}\right)+\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial T}{\partial \theta }\right)=0$$ With boundary conditions in dimensionless form given by:

$$T=0$$, $$r\to \infty$$, this sets $$T$$ equal to the initial value far from the cap,

$$T=1$$, $$r=\frac{1}{2} \left(\sqrt{2} \sqrt{-p^2+(1-p)^2 \cos (2 \theta )+2 p+1}+2 (1-p) \cos (\theta )\right)$$, this imposes $$T$$ on the cap

$$\frac{\partial T}{\partial \theta}\bigg| _{\theta=\pi/2}=0$$, adiabatic condition

$$\frac{\partial T}{\partial \theta}\bigg| _{\theta=\pi}=0$$, symmetry condition.

I tried this:

p = 0.2;

boundaries = {-r + 1/2 (Sqrt[2] Sqrt[Cos[2 θ] (1 - p)^2 + 2 p + 1 - p^2] +
2 Cos[θ] (1 - p)), r - 100, -θ + Pi/2, θ - Pi}
Ω = ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {r, θ}];
RegionPlot[Ω, PlotRange -> {{0, 3}, {1, 5}}]
NDSolveValue[{r^2 D[T[r, θ], {r, 2}] +
2 r D[T[r, θ], r] + D[T[r, θ], {θ, 2}] +
Cot[θ] D[T[r, θ], θ] == {NeumannValue[0.,
boundaries[[3]] == 0],
NeumannValue[0., boundaries[[4]] == 0]}, {DirichletCondition[
T[r, θ] == 1., boundaries[[1]] == 0.],
DirichletCondition[T[r, θ] == 0.,
boundaries[[2]] ==
0.]}}, T, {r, θ} ∈ Ω]


But it does not seem to work. I have two Dirichlet conditions and two Neumann conditions, but I don't know if I inserted them in NDSolve in the right way.

• 1. You need to express the b.c. involving derivative with NeumannValue, related: mathematica.stackexchange.com/q/224812/1871 2. How can $\theta=3\pi/2$ in spherical coordinates? What convention do you follow? – xzczd Apr 27 at 11:01
• Using physical convention $\theta \leq \pi$, you're right. The question was corrected according to your comment. Thank you. – umby Apr 28 at 11:17

Two issues here.

First of all, you've chosen 100 to approximate Infinity, which is way too large in this case. Something like 5 is OK:

p = 0.2; inf= 5;
boundaries = {-r + 1/2 (Sqrt[2] Sqrt[Cos[2 θ] (1 - p)^2 + 2 p + 1 - p^2] +
2 Cos[θ] (1 - p)), r - inf, -θ + Pi/2, θ - Pi};
<< NDSolveFEM
Ω = ToElementMesh@ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {r, θ}];


ToElementMesh is added to help NDSolve analyzing the domain properly, otherwise the femcbtd warning will pop up, at least in v12.2. (Alternatively, DiscretizeRegion can be used in place of ToElementMesh, which is a bit slower. )

The next issue is, you haven't set NeumannValue correctly. If you read the Details section of NeumannValue and the FEM document carefully, you might notice NeumannValue is actually defined based on the formal form of a PDE. How can we check the formal form of a PDE? The new-in-12.2 NDSolveFEMGetInactivePDE does the work. (If you're not yet in v12.2, try the function in this post. )

seq = Sequence[{r^2 D[T[r, θ], {r, 2}] + 2 r D[T[r, θ], r] +
D[T[r, θ], {θ, 2}] + Cot[θ] D[T[r, θ], θ] ==
0, {DirichletCondition[T[r, θ] == 1., boundaries[[1]] == 0.],
DirichletCondition[T[r, θ] == 0., boundaries[[2]] == 0.]}},
T, {r, θ} ∈ Ω]

{state} = NDSolveProcessEquations@seq

GetInactivePDE@state


exprInBlueBox = -{{-r^2, 0}, {0, -1}} . Inactive[Grad][T[r, θ], {r, θ}];


The normal vector $$\overset{\rightharpoonup }{n}=(0,1)$$ at $$\theta=\pi$$, so the left hand side (LHS) of Neumann b.c. at $$\theta=\pi$$ is:

normalVector = {0, 1};
- exprInBlueBox . normalVector // Activate
(* - Derivative[0, 1][T][r, θ] *)


The normal vector $$\overset{\rightharpoonup }{n}=(0,-1)$$ at $$\theta=\pi/2$$, so the LHS of Neumann b.c. at $$\theta=\pi/2$$ is:

normalVector = {0, -1};
- exprInBlueBox . normalVector // Activate
(* Derivative[0, 1][T][r, θ] *)


Thus $$\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi/2}=0$$ and $$\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi}=0$$ are equivalent to zero NeumannValue in this case. Once again, according to the Details section of NeumannValue:

…not specifying a boundary condition at all is equivalent to specifying a Neumann 0 condition.

In other words, we don't need to explicitly set NeumannValue for your problem. So the problem can be solved with:

sol = NDSolveValue@seq

DensityPlot[sol[Sqrt[x^2 + z^2], ArcTan[z, Abs@x]], {x, -5, 5}, {z, -5, 0},
AspectRatio -> Automatic, PlotPoints -> 100, PlotRange -> All]


You may further adjust value of inf and MaxCellMeasure option of ToElementMesh to see how the solution varies.

• Great. Small suggestion, you could use NDSolveFEMToElementMesh in place of DiscretizeRegion as that would return a second order mesh with curved elements and should approximate the region and solution better. But this is probably a minor thing. – user21 Apr 29 at 11:46
• @user21 The mesh quality seems to be similar in this case, but ToElementMesh turns out to be faster. Edited. Thx for the suggestion. – xzczd Apr 29 at 12:02
• @xzczd Many thx. Knowing the formal form of a PDE is fundamental to set Neumann BC. If I understand, according to the documentation (can u check pls?), in my problem $a=f=0$, $\alpha=\gamma=0$, $\beta=\left\{c,cot\theta\right\}$ and $c=\left\{-r^{2},-1\right\}$. Then to set $\frac{\partial T}{\partial \theta}=0$ on the boundary3, I have to write: == NeumannValue[-1, boundaries[[3]] == 0.]. While, to impose the same b.c. on both boundary 3 and 4 implies: == NeumannValue[-1, boundaries[[3]] == 0.]+NeumannValue[-1, boundaries[[4]] == 0.]. Unfortunately, I use v11 so no NDSolveFEMGetInactivePDE. – umby Apr 30 at 10:44
• @umby You can create a .m of course, but it's not necessary, just copy and execute the code in a notebook is OK. In principle, $\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi/2}=0$ should be interpreted to NeumannValue[1, boundaries[[3]] == 0.], because the normal vector is {0,-1} here. – xzczd Apr 30 at 13:13
• @umby Yes, but as mentioned in the comment above, NeumannValue[-1, boundaries[[4]] == 0.] is ill-posed. – xzczd Apr 30 at 14:30

In a previous answer 240190, I showed how one could use anisotropic meshing to add a DirichletCondition at "infinity" for a 1D problem. In this answer, I shall extend the technique to a 2D problem.

# Geometry description

In many FEM software packages, problems with spherical symmetry can be posed as an axisymmetric problem. Since it is easier for me to think in these terms, I will recast the problem.

As I understood the system, a spherical cap is embedded in a semi-infinite domain, as I have sketched below. The y-axis is the symmetry axis.

# Helper functions

## Mesh helper functions

I use some of the following helper functions to construct an anisotropic mesh based on connecting and extending edge segments. A structured Quad mesh can then be easily constructed using RegionProduct

(*Import required FEM package*)
Needs["NDSolveFEM"];
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] :=
Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] :=
Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] :=
Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=(*Quiet@*)
Abs@FindRoot[
firstElmHeight[x0, xf, n, r] - fElm, {r, 0.00000001,
100000000/fElm}, Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] :=
N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] :=
Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]]
reflectRight[pts_] :=
With[{rt = ReflectionTransform[{1}, {Last@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] :=
With[{rt = ReflectionTransform[{-1}, {First@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]


## Model specific helper functions

Typically, the structured Quad mesh is used on rectangular domains. I use the following helper functions to map a square UV space mesh onto the curved domain.

Clear[β, γ, rcap, rinf, rl, capMesh]
β[R_, h_] := ArcCos[(R - h)/R]
γ[R_, h_, ρ_] := ArcCos[(R - h)/ρ]
rcap[R_, h_][u_] := Module[{angle = β[R, h], r = R},
R {Sin[angle u], -Cos[angle u]}]
rinf[R_, h_, ρ_][u_] :=
Module[{angle = γ[R, h, ρ], r = ρ},
r {Sin[angle u], -Cos[angle u]}]
rl[R_, h_, ρ_][u_, v_] :=
Module[{rc = rcap[R, h][u], ri = rinf[R, h, ρ][u]},
(ri - rc) v + rc]
capMesh[R_, h_, ρ_][rh_, rv_] :=
Module[{sqr, crd, inc, msh, mean, mrkrs, bmrkrs, pEle, pe, pm, pcrd,
sdf, n, leIds, bcEle, z = {0, 0}, ex = {1, 0}, ey = {0, 1}, f, g},
sqr = RegionProduct[rh, rv];
crd = MeshCoordinates[sqr];
inc = ( Delete[0] /@ MeshCells[sqr, 2]);
mean = Mean /@ GetElementCoordinates[crd, #] & /@ {inc} // First;
mrkrs = If[#2 > rf/ρ, 2, 1] & @@@ mean;
msh = ToElementMesh["Coordinates" -> crd,
pm = epm[msh] /. {0 -> 4};
pe = epi[msh];
pcrd = crd[[Flatten@pe]];
sdf = Flatten@
Position[SignedRegionDistance[#, pcrd], _?(Abs[#] < 0.0000001 &),
1] &;
g = (pm[[#1]] = First@#2) &;
MapIndexed[g,
sdf /@ Table[
TransformedRegion[Line[{z, ex}],
RotationTransform[i 90 °, 1/2 (ex + ey)]], {i, 0, 3}]];
pEle = {PointElement[pe, pm]};
bmrkrs = ebm[msh];
n = ebn[msh];
leIds = Range@Length@n;
f = Function[{d}, Flatten@Position[n, _?(0.9999 < d . # &), 1]];
g = (bmrkrs[[#1]] = First@#2) &;
MapIndexed[g,
f /@ Table[RotationTransform[i 90 °][-ey], {i, 0, 3}]];
bcEle = {LineElement[ebi[msh], bmrkrs]};
crd = rl[R, h, ρ][#1, #2] & @@@ crd;
inc = inc /. {{i_, j_, k_, l_} :> {l, k, j, i}};
ToElementMesh["Coordinates" -> crd,
"BoundaryElements" -> bcEle, "PointElements" -> pEle]
]


# Mesh construction

The following workflow constructs a mesh with an angular resolution of 1°. Radially, there are two segments. There is a fine mesh in the region of interest (defined as 5X the radius) and an infinite segment that extends 10,000X the region of interest.

(*Define geometric and meshing parameters*)
R = 1; h = 1/5; rf = 5 R; Rinf =
10000 rf; nelmr = 80; nelminf = 40; nelmang = 90;
Print["Angular discretization segment"]
segu = Subdivide[0, 1, nelmang];
ru = pointsToMesh@segu
Print["Mesh segment in the radial region of interest"]
segr = leftSegmentGrowth[rf, nelmr, rf/100];
pointsToMesh@segr
seginf = meshGrowthByElm0[Rinf - rf, nelminf, Last@segr - segr[[-2]]];
reginf = pointsToMesh@seginf
rv = pointsToMesh@(#/Last[#] &@extendMesh[segr, seginf])
mesh = capMesh[R, h, Rinf][ru, rv];
Print["Full domain"]
Show[mesh["Wireframe"], Axes -> True]
Print["Zoomed region"]
Show[mesh["Wireframe"], PlotRange -> {{0, 2}, {-R + h, -2}},
Axes -> True]


# PDE set up and solution

In the Heat Transfer Verification Manual there are some helper functions to create a well-formed operator for axisymmetric heat transfer problems. The code is reproduced here:

Clear[HeatTransferModelAxisymmetric, TimeHeatTransferModelAxisymmetric]
HeatTransferModelAxisymmetric[T_, {r_, z_}, k_, ρ_, Cp_,
Velocity_, Source_] :=
Module[{V, Q},
V = If[Velocity === "NoFlow",
0, ρ*Cp*Velocity . Inactive[Grad][T, {r, z}]];
Q = If[Source === "NoSource", 0, Source];
1/r*D[-k*r*D[T, r], r] + D[-k*D[T, z], z] + V - Q]
TimeHeatTransferModelAxisymmetric[T_, TimeVar_, {r_, z_}, k_, ρ_,
Cp_, Velocity_, Source_] := ρ*Cp*D[T, {TimeVar, 1}] +
HeatTransferModelAxisymmetric[T, {r, z}, k, ρ, Cp, Velocity,
Source]


After all the heavy lifting has been done to create the mesh, the construction and solution of the PDE are straightforward.

parms = {k -> 1, ρ -> 1, Cp -> 1, hc -> 10, Ta -> 0};
Γhot =
DirichletCondition[θ[r, z] == 1, ElementMarker == 1];
Γcold =
DirichletCondition[θ[r, z] == 0, ElementMarker == 3];
Γconv = 0;
parmop = HeatTransferModelAxisymmetric[θ[r, z], {r, z},
k, ρ, Cp, "NoFlow", "NoSource"];
op = parmop /. parms;
pde = {op == Γconv, Γhot, \
Γcold};
Tfun = NDSolveValue[pde, θ, {r, z} ∈ mesh];


Now, we can construct some plots:

Plot[{Tfun[0, -z], Tfun[z, -R + h]}, {z, 0, 5}, PlotPoints -> 100,
PlotRange -> {0, 1.0}, PlotLegends -> "Expressions",
PlotLabel -> "Temperature along symmetry edges"]
uRange = MinMax[Tfun["ValuesOnGrid"]];
legendBar =
BarLegend[{"TemperatureMap", uRange}, 50,
LegendLabel -> Style["[°C]", Opacity[0.6]]];
options = {PlotRange -> {{-2, 2}, {-2.8, -R + h}, uRange},
ColorFunction -> ColorData[{"TemperatureMap", uRange}],
ContourStyle -> Opacity[0.5], ColorFunctionScaling -> False,
Contours -> 10, AspectRatio -> Automatic,
PlotPoints -> 100, FrameLabel -> {"r", "z"},
PlotLabel -> Style["Temperature Field: θ(r,z)", 18],
ImageSize -> 650};
Legended[ContourPlot[Tfun[Abs[r], z], {r, -2, 2}, {z, -2.8, 0},
Evaluate[options]], legendBar]


As you can see in the first plot, at a distance of 5, the temperature has decayed about 90%.

A benefit of anisotropic meshing is that one can pose some stringent questions to the model with minimal computational cost. For example, suppose you had a requirement that you needed to know the distance where temperature decayed 99.99% of the spherical cap value. One can easily find that this occurs at a distance of about 4000 as shown below:

FindRoot[Tfun[0, -z] - 0.0001, {z, 100}]
(* {z -> 4024.02} *)


# Convectively cooled top surface

It is straightforward to create a convectively cooled top surface (ElementMarker==2) using a Robin-type condition. From the previously defined parms, I defined a convective heat transfer coefficient of 10 and an ambient fluid temperature of 0°. To set up, we simply need to modify the NeumannValue.

Γconv =
NeumannValue[hc (Ta - θ[r, z]), ElementMarker == 2] /. parms;
pde = {op == Γconv, Γhot, \
Γcold};
Tfun = NDSolveValue[pde, θ, {r, z} ∈ mesh];


We can plot the solution as before:

Plot[{Tfun[0, -z], Tfun[z, -R + h]}, {z, 0, 5}, PlotPoints -> 100,
PlotRange -> {0, 1.0}, PlotLegends -> "Expressions",
PlotLabel -> "Temperature along symmetry edges"]
uRange = MinMax[Tfun["ValuesOnGrid"]];
legendBar =
BarLegend[{"TemperatureMap", uRange}, 50,
LegendLabel -> Style["[°C]", Opacity[0.6]]];
options = {PlotRange -> {{-2, 2}, {-2.8, -R + h}, uRange},
ColorFunction -> ColorData[{"TemperatureMap", uRange}],
ContourStyle -> Opacity[0.5], ColorFunctionScaling -> False,
Contours -> 10, AspectRatio -> Automatic,
PlotPoints -> 100, FrameLabel -> {"r", "z"},
PlotLabel -> Style["Temperature Field: θ(r,z)", 18],
ImageSize -> 650};
Legended[ContourPlot[Tfun[Abs[r], z], {r, -2, 2}, {z, -2.8, 0},
Evaluate[options]], legendBar]


# Comparison with another code

When possible, it is often conducive to compare the Mathematica results with another simulation code. To simulate boundary conditions at infinity, the FEM software COMSOL introduces an Infinite Element Domain (IED) concept below.

A large scaling factor (e.g., 1000) is applied to the equations in the IED.

As shown below, there is an excellent agreement between the Mathematica and COMSOL simulations. That should give one more confidence in the validity of this approach to solve the infinite domain problem.

• Thanks for the comparative analysis @Tim Laska! Does the IED concept imply using a special elements or it is just increasing of element length far from heat source? – Oleksii Semenov May 5 at 13:06
• @OleksiiSemenov You are welcome! In COMSOL, the IED is either mapped in 2D (quads) or swept in 3D (wedges or hexas) in the infinite direction. I do not think you can apply an IED to an unstructured mesh. I am not quite sure what COMSOL is doing under the hood. It is possible that they are scaling the elements invisible to the user or possibly they are scaling the thermal conductivity in a piecewise fashion to achieve the same effect. – Tim Laska May 5 at 13:51
• @OleksiiSemenov From the [COMSOL documentation] (doc.comsol.com/5.5/doc/com.comsol.help.comsol/…) applies coordinates stretching of the IED not dissimilar to what I have done. In my approach, the first element of the IED matches the width of the outer layer of the finite domain and stretches the elements to provide a gradual transition out to "infinity". – Tim Laska May 5 at 19:45
• Thanks for the link @TimLaska. I also found some papers devoted to this aspect: Zienkiewicz, O. C., C. Emson, and P. Bettess. "A novel boundary infinite element." International Journal for Numerical Methods in Engineering 19.3 (1983): 393-404. and Marques, J. M. M. C., and D. R. J. Owen. "Infinite elements in quasi-static materially nonlinear problems." Computers & structures 18.4 (1984): 739-751. – Oleksii Semenov May 6 at 11:42

One can also consider the 3D statement of the problem. Solution of a such linear problem is not so time consuming nowadays. For mesh generation let's take advantage of OpenCascadeLink procedures which are very useful for tessellation of domains with complex geometry. Let's $$r$$ is a radius of a cap and $$R$$ is "infinity" radius. The computational region is defined by a difference between spherical wedge of radius $$R$$ and a ball of radius $$r$$. Whereas the spherical wedge can be defined as intersection of rectangular cuboid and a ball of radius $$R$$.

Definition of computation domain and FE mesh generation

Needs["OpenCascadeLink"]
Needs["NDSolveFEM"]

r = 10; (*radius of a cap*)
p = 0.2;

shape1 = OpenCascadeShape[Ball[{0, 0, 0}, R]];
Hexahedron[{{0, -R, -R}, {R, -R, -R}, {Sqrt[2] R,
Sqrt[2] R, -R}, {0, R, -R}, {0, -R, R}, {R, -R, R}, {Sqrt[2] R,
Sqrt[2] R, R}, {0, R, R}}]];
shape3 = OpenCascadeShape[Ball[{p - r, 0, 0}, r]];

difference];(*boundary mesh geteration*)

mesh = ToElementMesh[bmesh];(*FE mesh generation*)

groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;

Show[mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]],
Axes -> True, AxesLabel -> {"x", "y", "z"},
AxesStyle -> RGBColor[0, 0, 0], BaseStyle -> 14]


Numerical solution

In the code below we solve the problem by means of FEM low level routines. Boundary elements with ElementMarkers=4;5belong to the surface of cap whereas on "infinity" surface ElementMarkers=1;3. This is taken into account when implementing Dirichlet BC. The rest surface of computational domain is adiabatic.

nr = ToNumericalRegion[mesh];
vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u}, {x, y, z}}];
sd = NDSolveSolutionData["Space" -> nr];

pded = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {{-IdentityMatrix[3]}}];
bcd = InitializeBoundaryConditions[vd,
sd, {{DirichletCondition[u[x, y, z] == 1,
ElementMarker == 4 || ElementMarker == 5],
DirichletCondition[u[x, y, z] == 0,
ElementMarker == 1 || ElementMarker == 3]}}]
md = InitializePDEMethodData[vd, sd];

dpde = DiscretizePDE[pded, md, sd];
dbc = DiscretizeBoundaryConditions[bcd, md, sd];

ufun = ElementMeshInterpolation[{mesh}, res];


Postprocessing

Temperature distributions along axes $$y$$ and $$z$$ should be the same

Show[

Plot[ufun[x, 0, 0], {x, p, R}, PlotRange -> All,
PlotLegends -> {"along x axis"},
PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.005]}]
,

Plot[{ufun[0, x, 0], ufun[0, 0, x]}, {x, Sqrt[r^2 - (p - r)^2], R},
PlotRange -> All, PlotLegends -> {"along y axis", "along z axis"},
PlotStyle -> {{RGBColor[0, 1, 0],
Thickness[0.005]}, {RGBColor[0, 0, 1], Thickness[0.005]}}],
Frame -> True , FrameLabel -> {"Distance", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 18, ImageSize -> 600,
GridLines -> {{p, Sqrt[r^2 - (p - r)^2]}, None},
GridLinesStyle -> {Dashed, RGBColor[0, 0, 0]}
]