1 function [V,n,p,res,niter] = DDGnlpoisson (x,sinodes,Vin,nin,...
2 pin,Fnin,Fpin,D,l2,toll,maxit,verbose)
4 % [V,n,p,res,niter] = DDGnlpoisson (x,sinodes,Vin,nin,...
5 % pin,Fnin,Fpin,D,l2,toll,maxit,verbose)
6 % Solves the non linear Poisson equation
7 % $$ - lamda^2 *V'' + (n(V,Fn) - p(V,Fp) -D)=0 $$
8 % input: x spatial grid
9 % sinodes index of the nodes of the grid which are in the
10 % semiconductor subdomain
11 % (remaining nodes are assumed to be in the oxide subdomain)
12 % Vin initial guess for the electrostatic potential
13 % nin initial guess for electron concentration
14 % pin initial guess for hole concentration
15 % Fnin initial guess for electron Fermi potential
16 % Fpin initial guess for hole Fermi potential
18 % l2 scaled electric permittivity (diffusion coefficient)
19 % toll tolerance for convergence test
20 % maxit maximum number of Newton iterations
21 % verbose verbosity level: 0,1,2
22 % output: V electrostatic potential
23 % n electron concentration
24 % p hole concentration
25 % res residual norm at each step
26 % niter number of Newton iterations
28 ## This file is part of
30 ## SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
31 ## -------------------------------------------------------------------
32 ## Copyright (C) 2004-2007 Carlo de Falco
36 ## SECS1D is free software; you can redistribute it and/or modify
37 ## it under the terms of the GNU General Public License as published by
38 ## the Free Software Foundation; either version 2 of the License, or
39 ## (at your option) any later version.
41 ## SECS1D is distributed in the hope that it will be useful,
42 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
43 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
44 ## GNU General Public License for more details.
46 ## You should have received a copy of the GNU General Public License
47 ## along with SECS1D; If not, see <http://www.gnu.org/licenses/>.
50 %% Set some useful constants
55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 %% convert grid info to FEM form
57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 sielements = sinodes(1:end-1);
63 Nnodes = length(nodes);
66 totdofs = Nnodes - Ndiricheletnodes ;
68 elements = [[1:Nnodes-1]' , [2:Nnodes]'];
69 Nelements = size(elements,1);
73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76 %% we're going to solve
77 %% $$ - lamda^2*(\delta V)'' + (\frac{\partial n}{\partial V} -
78 %% \frac{\partial p}{\partial V})= -R $$
80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 %% set $$ n_1 = nin $$ and $$ V = Vin $$
86 n = DDGphin2n(V(sinodes),Fn);
87 p = DDGphip2p(V(sinodes),Fp);
92 if (sinodes(end)==Nnodes)
98 %%% Compute LHS matrices
101 %% let's compute FEM approximation of $$ L = - \frac{d^2}{x^2} $$
102 L = Ucomplap (nodes,Nnodes,elements,Nelements,l2.*ones(Nelements,1));
104 %% compute $$ Mv = ( n + p) $$
105 %% and the (lumped) mass matrix M
106 Mv = zeros(Nnodes,1);
107 Mv(sinodes) = (n + p);
108 Cv = zeros(Nelements,1);
110 M = Ucompmass (nodes,Nnodes,elements,Nelements,Mv,Cv);
113 %%% Compute RHS vector (-residual)
116 %% now compute $$ T0 = (n - p - D) $$
117 Tv0 = zeros(Nnodes,1);
118 Tv0(sinodes) = (n - p - D);
119 Cv = zeros(Nelements,1);
121 T0 = Ucompconst (nodes,Nnodes,elements,Nelements,Tv0,Cv);
123 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
127 %% Apply boundary conditions
133 %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
134 normr(1) = norm(R,inf);
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 %% START OF THE NEWTON CYCLE
144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147 fprintf(1,'\n newton iteration: %d, reldVnorm = %e',newtit,reldVnorm);
156 %% Start of the damping procedure
161 fprintf(1,'\n damping iteration: %d, residual norm = %e',dit,normrnew);
165 n = DDGphin2n(Vnew(sinodes),Fn);
166 p = DDGphip2p(Vnew(sinodes),Fp);
171 if (sinodes(end)==Nnodes)
176 %%% Compute LHS matrices
179 %% let's compute FEM approximation of $$ L = - \frac{d^2}{x^2} $$
180 %L = Ucomplap (nodes,Nnodes,elements,Nelements,ones(Nelements,1));
182 %% compute $$ Mv = ( n + p) $$
183 %% and the (lumped) mass matrix M
184 Mv = zeros(Nnodes,1);
185 Mv(sinodes) = (n + p);
186 Cv = zeros(Nelements,1);
188 M = Ucompmass (nodes,Nnodes,elements,Nelements,Mv,Cv);
191 %%% Compute RHS vector (-residual)
194 %% now compute $$ T0 = (n - p - D) $$
195 Tv0 = zeros(Nnodes,1);
196 Tv0(sinodes) = (n - p - D);
197 Cv = zeros(Nelements,1);
199 T0 = Ucompconst (nodes,Nnodes,elements,Nelements,Tv0,Cv);
201 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
203 Rnew = L * Vnew + T0;
205 %% Apply boundary conditions
206 Anew(BCnodes,:) = [];
207 Anew(:,BCnodes) = [];
211 if ((dit>1)&(norm(Rnew,inf)>=norm(R,inf)))
213 fprintf(1,'\nexiting damping cycle \n');
221 %% compute $$ | R_{k+1} | $$ for the convergence test
222 normrnew= norm(R,inf);
224 % check if more damping is needed
225 if (normrnew > normr(newtit))
229 fprintf(1,'\nexiting damping cycle because residual norm = %e \n',normrnew);
236 normr(newtit+1) = normrnew;
237 dVnorm = norm(tk*dV,inf);
238 % check if convergence has been reached
239 reldVnorm = dVnorm / norm(V,inf);
240 if (reldVnorm <= toll)
242 fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
253 % $Date: 2008-02-04 16:26:27 +0100 (man, 04 feb 2008) $