--- /dev/null
+function [V,n,p,res,niter] = DDGnlpoisson (x,sinodes,Vin,nin,...
+ pin,Fnin,Fpin,D,l2,toll,maxit,verbose)
+
+% [V,n,p,res,niter] = DDGnlpoisson (x,sinodes,Vin,nin,...
+% pin,Fnin,Fpin,D,l2,toll,maxit,verbose)
+% Solves the non linear Poisson equation
+% $$ - lamda^2 *V'' + (n(V,Fn) - p(V,Fp) -D)=0 $$
+% input: x spatial grid
+% sinodes index of the nodes of the grid which are in the
+% semiconductor subdomain
+% (remaining nodes are assumed to be in the oxide subdomain)
+% Vin initial guess for the electrostatic potential
+% nin initial guess for electron concentration
+% pin initial guess for hole concentration
+% Fnin initial guess for electron Fermi potential
+% Fpin initial guess for hole Fermi potential
+% D doping profile
+% l2 scaled electric permittivity (diffusion coefficient)
+% toll tolerance for convergence test
+% maxit maximum number of Newton iterations
+% verbose verbosity level: 0,1,2
+% output: V electrostatic potential
+% n electron concentration
+% p hole concentration
+% res residual norm at each step
+% niter number of Newton iterations
+
+## This file is part of
+##
+## SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
+## -------------------------------------------------------------------
+## Copyright (C) 2004-2007 Carlo de Falco
+##
+##
+##
+## SECS1D is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## SECS1D is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with SECS1D; If not, see <http://www.gnu.org/licenses/>.
+
+
+%% Set some useful constants
+dampit = 10;
+dampcoeff = 2;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% convert grid info to FEM form
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Ndiricheletnodes = 2;
+nodes = x;
+sielements = sinodes(1:end-1);
+
+Nnodes = length(nodes);
+
+
+totdofs = Nnodes - Ndiricheletnodes ;
+
+elements = [[1:Nnodes-1]' , [2:Nnodes]'];
+Nelements = size(elements,1);
+
+BCnodes = [1;Nnodes];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%
+%% initialization:
+%% we're going to solve
+%% $$ - lamda^2*(\delta V)'' + (\frac{\partial n}{\partial V} -
+%% \frac{\partial p}{\partial V})= -R $$
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% set $$ n_1 = nin $$ and $$ V = Vin $$
+V = Vin;
+Fn = Fnin;
+Fp = Fpin;
+n = DDGphin2n(V(sinodes),Fn);
+p = DDGphip2p(V(sinodes),Fp);
+if (sinodes(1)==1)
+ n(1)=nin(1);
+ p(1)=pin(1);
+end
+if (sinodes(end)==Nnodes)
+ n(end)=nin(end);
+ p(end)=pin(end);
+end
+
+%%%
+%%% Compute LHS matrices
+%%%
+
+%% let's compute FEM approximation of $$ L = - \frac{d^2}{x^2} $$
+L = Ucomplap (nodes,Nnodes,elements,Nelements,l2.*ones(Nelements,1));
+
+%% compute $$ Mv = ( n + p) $$
+%% and the (lumped) mass matrix M
+Mv = zeros(Nnodes,1);
+Mv(sinodes) = (n + p);
+Cv = zeros(Nelements,1);
+Cv(sielements)= 1;
+M = Ucompmass (nodes,Nnodes,elements,Nelements,Mv,Cv);
+
+%%%
+%%% Compute RHS vector (-residual)
+%%%
+
+%% now compute $$ T0 = (n - p - D) $$
+Tv0 = zeros(Nnodes,1);
+Tv0(sinodes) = (n - p - D);
+Cv = zeros(Nelements,1);
+Cv(sielements)= 1;
+T0 = Ucompconst (nodes,Nnodes,elements,Nelements,Tv0,Cv);
+
+%% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
+A = L + M;
+R = L * V + T0;
+
+%% Apply boundary conditions
+A(BCnodes,:) = [];
+A(:,BCnodes) = [];
+
+R(BCnodes) = [];
+
+%% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
+normr(1) = norm(R,inf);
+relresnorm = 1;
+reldVnorm = 1;
+normrnew = normr(1);
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%
+%% START OF THE NEWTON CYCLE
+%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+for newtit=1:maxit
+ if verbose
+ fprintf(1,'\n newton iteration: %d, reldVnorm = %e',newtit,reldVnorm);
+
+ end
+
+ dV =[0;(A)\(-R);0];
+
+
+
+ %%%%%%%%%%%%%%%%%%
+ %% Start of the damping procedure
+ %%%%%%%%%%%%%%%%%%
+ tk = 1;
+ for dit = 1:dampit
+ if verbose
+ fprintf(1,'\n damping iteration: %d, residual norm = %e',dit,normrnew);
+ end
+ Vnew = V + tk * dV;
+
+ n = DDGphin2n(Vnew(sinodes),Fn);
+ p = DDGphip2p(Vnew(sinodes),Fp);
+ if (sinodes(1)==1)
+ n(1)=nin(1);
+ p(1)=pin(1);
+ end
+ if (sinodes(end)==Nnodes)
+ n(end)=nin(end);
+ p(end)=pin(end);
+ end
+ %%%
+ %%% Compute LHS matrices
+ %%%
+
+ %% let's compute FEM approximation of $$ L = - \frac{d^2}{x^2} $$
+ %L = Ucomplap (nodes,Nnodes,elements,Nelements,ones(Nelements,1));
+
+ %% compute $$ Mv = ( n + p) $$
+ %% and the (lumped) mass matrix M
+ Mv = zeros(Nnodes,1);
+ Mv(sinodes) = (n + p);
+ Cv = zeros(Nelements,1);
+ Cv(sielements)= 1;
+ M = Ucompmass (nodes,Nnodes,elements,Nelements,Mv,Cv);
+
+ %%%
+ %%% Compute RHS vector (-residual)
+ %%%
+
+ %% now compute $$ T0 = (n - p - D) $$
+ Tv0 = zeros(Nnodes,1);
+ Tv0(sinodes) = (n - p - D);
+ Cv = zeros(Nelements,1);
+ Cv(sielements)= 1;
+ T0 = Ucompconst (nodes,Nnodes,elements,Nelements,Tv0,Cv);
+
+ %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
+ Anew = L + M;
+ Rnew = L * Vnew + T0;
+
+ %% Apply boundary conditions
+ Anew(BCnodes,:) = [];
+ Anew(:,BCnodes) = [];
+
+ Rnew(BCnodes) = [];
+
+ if ((dit>1)&(norm(Rnew,inf)>=norm(R,inf)))
+ if verbose
+ fprintf(1,'\nexiting damping cycle \n');
+ end
+ break
+ else
+ A = Anew;
+ R = Rnew;
+ end
+
+ %% compute $$ | R_{k+1} | $$ for the convergence test
+ normrnew= norm(R,inf);
+
+ % check if more damping is needed
+ if (normrnew > normr(newtit))
+ tk = tk/dampcoeff;
+ else
+ if verbose
+ fprintf(1,'\nexiting damping cycle because residual norm = %e \n',normrnew);
+ end
+ break
+ end
+ end
+
+ V = Vnew;
+ normr(newtit+1) = normrnew;
+ dVnorm = norm(tk*dV,inf);
+ % check if convergence has been reached
+ reldVnorm = dVnorm / norm(V,inf);
+ if (reldVnorm <= toll)
+ if(verbose)
+ fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
+ end
+ break
+ end
+end
+
+res = normr;
+niter = newtit;
+
+% Last Revision:
+% $Author: adb014 $
+% $Date: 2008-02-04 16:26:27 +0100 (man, 04 feb 2008) $
+
+
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff