Find the largest subdomain of a disk with negative solution to the Dirichlet Poisson equation

Question

Question

I am looking for a way to find the largest subdomain $U$ of the unit disk $D(0,1)$ for which the solution to $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in U\\ u&=0,\quad (x,y)\in\partial U \end{align} $$ is negative in $U$.

Some additional details

We start by looking at the problem in the disk. Solving the Dirichlet problem $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in D(0,1)\\ u&=0,\quad (x,y)\in\partial D(0,1) \end{align} $$ using

usol1[x_, y_] =
NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0.3 - y, 
 DirichletCondition[u[x, y] == 0, True]}, 
 u[x, y], 
 Element[{x, y},Disk[]]
]

gives a solution that is negative in some part (blue) and positive in some part (pink):

ContourPlot[Evaluate@usol1[x, y], Element[{x, y},Disk[]], 
 Contours -> {0, -0.05, -0.1}, 
 ContourStyle -> {Black, Black, {Thick, Black}}, 
 ContourShading -> {Darker[Blue], Blue, Lighter[Blue], Pink}, 
 PlotPoints -> 100, AspectRatio -> 1]

If we instead try a smaller domain $$ \Omega=D(0,1)\cap \{(x,y)\in D(0,1)~|~y<0.3\} $$ defined in Mathematica as

\[CapitalOmega] := 
  RegionDifference[Disk[], Rectangle[{-2, 0.3}, {2, 2}]]

We solve $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in \Omega\\ u&=0,\quad (x,y)\in\partial \Omega \end{align} $$ with

usol2[x_, y_] = 
NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0.3 - y, 
 DirichletCondition[u[x, y] == 0, True]}, 
 u[x, y], 
 Element[{x, y}, \[CapitalOmega]]
]

and show the corresponding contour plot with

ContourPlot[Evaluate@usol2[x, y], 
 Element[{x, y}, \[CapitalOmega]], 
 Contours -> {0, -0.05, -0.1}, 
 ContourStyle -> {Black, Black, {Thick, Black}}, 
 ContourShading -> {Darker[Blue], Blue, Lighter[Blue], Pink}, 
 PlotPoints -> 100, AspectRatio -> 1.3/2]

We find a strictly negative solution.

Now, clearly the domain $\Omega$ is not the largest possible, so I restate my question,

Is there some way, using Mathematica, probably with some numerical method, to find the largest subdomain $U\subset D(0,1)$ such that the Poisson problem $$ \begin{align} \Delta u&=0.3-y,\quad (x,y)\in U\\ u&=0,\quad (x,y)\in\partial U \end{align} $$ has a strictly negative solution in $U$?

Update

Just to clarify, I do not want to find the area of the part where my usol1 is negative.

Answer 1

Define the unit disc centred on the origin. This is used only to highlight the unit disc in the contour plot below.

disk = Disk[{0, 0}, 1];

Define a grid of points on the boundary of the disk. The parameter n controls the fineness of the grid.

n = 20;
angles = Range[0, 2Pi-Pi/n, Pi/n];
boundarypts = AngleVector /@ angles;

From here on down one iteration of the algorithm is being described. Iterations 2, 3, ... should begin at this point. The solution converges quickly, so I have left the iterations as being manual, rather than wrap the code up in a loop (or whatever).

Define the corresponding boundary mesh region. Updated versions of this mesh region are used to define the region over which the partial differential equation is to be solved.

region = BoundaryMeshRegion[boundarypts, Line[Append[Range[2 n], 1]]];

Solve the PDE using the current boundary mesh region.

soln = NDSolve[{Laplacian[u[x, y], {x, y}] == 0.3 - y, DirichletCondition[u[x, y] == 0, True]}, u, {x, y} \[Element] region];

Substitute this solution into a function for ease of use.

u0[{x_, y_}] = u[x, y] /. soln[[1]];

Generate a contour plot of the solution. The boundary mesh is highlighted in red, and the original unit disc is highlighted in blue.

plot = ContourPlot[u0[{x, y}], {x, y} \[Element] disk, Contours -> {0, -0.05, -0.1}, BoundaryStyle -> Red, Epilog -> {Blue, Circle[]}, PlotRange -> {{-1, 1}, {-1, 1}}, PlotRangePadding -> 0.1]

Display the boundary mesh region - i.e. the solution to the posed problem.

region

The above plot shows the solution after 20 update iterations, using the update procedure described below.

Compute the sign of the solution (the "signature") just inside each boundary point. Where the signature is negative/positive the boundary should be moved outwards/inwards (respectively).

signatures = u0 /@ (0.99 boundarypts);
maxsig = Max@Abs@signatures;

Update the boundary points:

1. If a signature is negative then move the corresponding boundary point outwards a bit, and truncate it back to the unit circle if it jumps outside the circle.
2. If the signature is positive then move the corresponding boundary point inwards a bit.

boundarypts = MapThread[ Switch[Sign[#1], -1, If[Norm[#] > 1, Normalize[#], #] &[(1 - #1) #2], 1, #2/(1 + #1), _, #2] &, {0.3/maxsig signatures, boundarypts}];

This algorithm for moving the boundary points is very crude, and with some experimentation it could no doubt be improved.