1 function [V,n,p,res,niter] = ThDDGOXnlpoisson (mesh,Dsides,Sinodes,SiDnodes,...
2 Sielements,Vin,Vthn,Vthp,...
4 Fnin,Fpin,D,l2,l2ox,...
8 %% [V,n,p,res,niter] = DDGOXnlpoisson (mesh,Dsides,Sinodes,Vin,nin,pin,...
9 %% Fnin,Fpin,D,l2,l2ox,toll,maxit,verbose)
11 %% solves $$ -\lambda^2 V'' + (n(V,Fn,Tn) - p(V,Fp,Tp) -D) = 0$$
15 %% This file is part of
17 %% SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
18 %% -------------------------------------------------------------------
19 %% Copyright (C) 2004-2006 Carlo de Falco
23 %% SECS2D is free software; you can redistribute it and/or modify
24 %% it under the terms of the GNU General Public License as published by
25 %% the Free Software Foundation; either version 2 of the License, or
26 %% (at your option) any later version.
28 %% SECS2D is distributed in the hope that it will be useful,
29 %% but WITHOUT ANY WARRANTY; without even the implied warranty of
30 %% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 %% GNU General Public License for more details.
33 %% You should have received a copy of the GNU General Public License
34 %% along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
36 global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS
38 %% Set some useful constants
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44 %% setup FEM data structures
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 Nnodes = length(nodes);
50 Nelements = length(elements);
52 % Set list of nodes with Dirichelet BCs
53 Dnodes = Unodesonside(mesh,Dsides);
55 % Set values of Dirichelet BCs
56 Bc = zeros(length(Dnodes),1);
57 % Set list of nodes without Dirichelet BCs
58 Varnodes = setdiff([1:Nnodes],Dnodes);
61 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 %% we're going to solve
65 %% $$ - \lambda^2 (\delta V)''
66 %% + (\frac{\partial n}{\partial V}
67 %% - \frac{\partial p}{\partial V})= -R $$
69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71 %% set $$ n_1 = nin $$ and $$ V = Vin $$
75 n = exp((V(Sinodes)-Fn)./Vthn);
76 p = exp((-V(Sinodes)+Fp)./Vthp);
77 n(SiDnodes) = nin(SiDnodes);
78 p(SiDnodes) = pin(SiDnodes);
82 %%% Compute LHS matrices
85 %% let's compute FEM approximation of $$ L = - \frac{d^2}{x^2} $$
86 if (isempty(DDGOXNLPOISSON_LAP))
87 coeff = l2ox * ones(Nelements,1);
89 DDGOXNLPOISSON_LAP = Ucomplap (mesh,coeff);
92 %% compute $$ Mv = ( n + p) $$
93 %% and the (lumped) mass matrix M
94 if (isempty(DDGOXNLPOISSON_MASS))
95 coeffe = zeros(Nelements,1);
97 DDGOXNLPOISSON_MASS = Ucompmass2(mesh,ones(Nnodes,1),coeffe);
99 freecarr=zeros(Nnodes,1);
100 freecarr(Sinodes)=(n./Vthn + p./Vthp);
102 M = DDGOXNLPOISSON_MASS*spdiag(Mv);
105 %%% Compute RHS vector (-residual)
108 %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
109 if (isempty(DDGOXNLPOISSON_RHS))
110 coeffe = zeros(Nelements,1);
111 coeffe(Sielements)=1;
112 DDGOXNLPOISSON_RHS = Ucompconst (mesh,ones(Nnodes,1),coeffe);
114 totcharge = zeros(Nnodes,1);
115 totcharge(Sinodes)=(n - p - D);
117 T0 = Tv0 .* DDGOXNLPOISSON_RHS;
119 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
120 A = DDGOXNLPOISSON_LAP + M;
121 R = DDGOXNLPOISSON_LAP * V + T0;
123 %% Apply boundary conditions
128 %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
129 normr(1) = norm(R,inf);
135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137 %% START OF THE NEWTON CYCLE
139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 fprintf(1,'\n***\nNewton iteration: %d, reldVnorm = %e\n***\n',newtit,reldVnorm);
146 dV(Varnodes) =(A)\(-R);
151 %% Start of th damping procedure
156 fprintf(1,'\ndamping iteration: %d, residual norm = %e\n',dit,normrnew);
160 n = exp((Vnew(Sinodes)-Fn)./Vthn);
161 p = exp((-Vnew(Sinodes)+Fp)./Vthp);
162 n(SiDnodes) = nin(SiDnodes);
163 p(SiDnodes) = pin(SiDnodes);
167 %%% Compute LHS matrices
170 %% let's compute FEM approximation of $$ L = - \frac{d^2}{x^2} $$
171 %L = Ucomplap (mesh,ones(Nelements,1));
173 %% compute $$ Mv = ( n + p) $$
174 %% and the (lumped) mass matrix M
175 freecarr=zeros(Nnodes,1);
176 freecarr(Sinodes)=(n./Vthn + p./Vthp);
178 M = DDGOXNLPOISSON_MASS*spdiag(Mv);
181 %%% Compute RHS vector (-residual)
184 %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
185 totcharge( Sinodes)=(n - p - D);
187 T0 = Tv0 .* DDGOXNLPOISSON_RHS;%T0 = Ucompconst (mesh,Tv0,ones(Nelements,1));
189 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
190 A = DDGOXNLPOISSON_LAP + M;
191 R = DDGOXNLPOISSON_LAP * Vnew + T0;
193 %% Apply boundary conditions
198 %% compute $$ | R_{k+1} | $$ for the convergence test
199 normrnew= norm(R,inf);
201 % check if more damping is needed
202 if (normrnew > normr(newtit))
206 fprintf(1,'\nexiting damping cycle because residual norm = %e \n-----------\n',normrnew);
213 normr(newtit+1) = normrnew;
214 dVnorm = norm(tk*dV,inf);
216 % check if convergence has been reached
217 reldVnorm = dVnorm / norm(V,inf);
218 if (reldVnorm <= toll)
220 fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);