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

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["NDSolveFEM"]
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]

• What are the boundary condition w.r.t x? – zhk May 3 '17 at 14:45
• @zhk you're quite right I hadn't defined boundary conditions wrt x, but the numerical solution still works without them just fine. I got the answer already, thanks anyway! – Dan May 4 '17 at 12:01

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]


• Awesome! Thanks @user21 ! Saved me days of work. I had not fully looked into the NeumannValue function as at first glance I was confused by the documentation. Now that I have read up a bit on different types of boundary conditions I can see mine looks like a Robin boundary which can also be defined through the NeumannValue as you have done. – Dan May 4 '17 at 11:57
• I am glad I could help. Would mind explaining to me what it was that confused you in the documentation - perhaps it could be improved but for that I'd need to know how. What would have helped you? – user21 May 4 '17 at 13:10
• Maybe confused is the wrong word, I guess overwhelmed. I was reading link and it was unfamiliar + thrown off by the symbols and length of the article. It all clicked when I saw a table on wikipedia of the different possible types of boundary conditions with notation more familiar to me link. That allowed me to "label" my pde and see it was actually covered quite well in the documentation once i knew where to look. – Dan May 4 '17 at 22:29
• @Dan, thanks. I'll keep the idea for a table of boundary conditions in mind. – user21 May 4 '17 at 23:18