2D PDE solving continuity eq: NDsolve stops without output or error

Question

My aim: Model distribution of electrons in semiconductor in 2D, using steady-state continuity equation. i.e. model $p(x,y)$ based on a partial differential equation and some boundary conditions with derivatives. Mathematica v10.0

The issue: NDsolve stops by itself, no output or error is given and then memory is cleared (or at least all variables cleared and font goes blue etc). Any advice would be greatly appreciated:) this is my first attempt doing anything like this in mathematica. From my tests at the end, I'm thinking it's something to do with my boundaries, or maybe this problem is too large for my laptop?

Code that's not working

Needs["NDSolve`FEM`"]
Dp = 50;    (* diffusivity *)
G = 10^10;            (* generation rate  *)
S = 100000; (* Surface Recombination Velocity *)
tau = 0.01; 
w = 0.5; (*  WIDTH  *)
wd = w; (*  depth  *)

bmesh = ToBoundaryMesh[Rectangle[{0, 0}, {w, wd}]];
mesh = ToElementMesh[bmesh,MeshRefinementFunction ->Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices]; If[x < 0.05*w || x > 0.95*w || y < 0.05*wd || y > 0.95*wd, area > 0.00001, area > 0.0001]]] ];
mesh["Wireframe"]

equation =Dp*(D[p[x, y], {x, 2}] + D[p[x, y], {y, 2}]) - p[x, y]/tau+G == 0;
cond3 = Dp*(D[p[x, y], y] /. y -> 0) == S p[x, 0];
cond4 = Dp*(D[p[x, y], y] /. y -> wd) == -S p[x, wd];
sol2D = NDSolve[{equation, cond3, cond4}, p, {x, 0, w}, {y, 0, wd}];
p2d[x_, y_] = p[x, y] /. sol2D[[1]];
Plot3D[Evaluate[p2d[x, y]], {x, 0, w}, {y, 0, wd}, PlotRange -> Full]

In plain english:

(1) I setup a 2D region which is a square x and y between (0,w). The square represents a semiconductor of uniform quality ("$\tau$"), but electrons can recombine (die) at the 4 edges of the square. I found specifying a mesh was necessary to avoid spikes in the results.

(2) I setup the continuity equation in steady state. For each area element, the electrons are conserved. This is a balance of generation of electrons, recombination of electrons, and diffusion of electrons in/out of the area. i.e. $D( \frac{\delta p^2(x,y)}{\delta x^2} + \frac{\delta p^2(x,y)}{\delta y^2})+G-\frac{p(x,y)}{\tau}=0$

(3) I setup the boundary conditions related to the edges. This appears to be the hurdle in the NDsolve (since it works without this boundary condition or for simpler ones). In the e.g. here, I apply the boundary condition to the top and bottom edge only (i.e. for y=0 and y=w boundary of square) as $D \frac{\delta p(x,0)}{\delta y}=S*p(x,0)$ and $D \frac{\delta p(x,w)}{\delta y}=-S*p(x,w)$. D and S are constants.

Code that works:

The code works if I set a simpler boundary condition such as $p(x,0)=0$ and $p(x,w)=0$. Replace above code 3rd paragraph with:

equation =Dp*(D[p[x, y], {x, 2}] + D[p[x, y], {y, 2}]) - p[x, y]/tau+G == 0;
cond1 = p[x, 0] == 0;
cond2 = p[x, wd] == 0;
sol2D = NDSolve[{equation, cond1, cond2}, p, {x, 0, w}, {y, 0, wd}];
p2d[x_, y_] = p[x, y] /. sol2D[[1]];
Plot3D[Evaluate[p2d[x, y]], {x, 0, w}, {y, 0, wd}, PlotRange -> Full]

And it also works for a similar case in 1D, even with the boundary condition with derivatives:

Dp = 50;    (* diffusivity *)
G = 10^10;            (* generation rate  *)
S = 100000; (* Surface Recombination Velocity *)
tau = 0.01; 
w = 0.5; (*  WIDTH  *)
wd = w; (*  depth  *)

sol1D = NDSolve[{Dp*p''[z] - p[z]/tau   + G == 0, 
Dp*p'[0] == S p[0], Dp*p'[w] == -S p[w]}, p[z], z];
p1d[z_] = p[z] /. sol1D[[1]];
LogPlot[{p1d[x]}, {x, 0, w}, PlotRange -> All]

Answer 1

Something like this get's you started (tried in version 11.1):

Dp = 50;(*diffusivity*)
G = 
 10^10;(*generation rate*)
S = 100000;(*Surface Recombination \
Velocity*)
tau = 0.01;
w = 0.5;(*WIDTH*)
wd = w;(*depth*)
equation = 
 Dp*(D[p[x, y], {x, 2}] + D[p[x, y], {y, 2}]) - p[x, y]/tau + G == 
  NeumannValue[S p[x, y], y == 0] + 
   NeumannValue[-S p[x, y], y == wd];
sol2D = NDSolveValue[{equation}, p, {x, 0, w}, {y, 0, wd}];
Plot3D[sol2D[x, y], {x, 0, w}, {y, 0, wd}, PlotRange -> Full]