octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXnlpoisson.m

   1 function [V,n,p,res,niter] = ThDDGOXnlpoisson (mesh,Dsides,Sinodes,SiDnodes,...
   2                                                Sielements,Vin,Vthn,Vthp,...
   3                                                nin,pin,...
   4                                                Fnin,Fpin,D,l2,l2ox,...
   5                                                toll,maxit,verbose)
   6
   7   %%
   8   %%   [V,n,p,res,niter] = DDGOXnlpoisson (mesh,Dsides,Sinodes,Vin,nin,pin,...
   9   %%                                       Fnin,Fpin,D,l2,l2ox,toll,maxit,verbose)
  10   %%
  11   %%  solves $$ -\lambda^2 V'' + (n(V,Fn,Tn) - p(V,Fp,Tp) -D) = 0$$
  12   %%
  13
  14
  15   %% This file is part of
  16   %%
  17   %%            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  18   %%         -------------------------------------------------------------------
  19   %%            Copyright (C) 2004-2006  Carlo de Falco
  20   %%
  21   %%
  22   %%
  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.
  27   %%
  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.
  32   %%
  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/>.
  35
  36         global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS
  37
  38         %% Set some useful constants
  39         dampit          = 3;
  40         dampcoeff       = 2;
  41
  42
  43         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  44         %% setup FEM data structures
  45         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  46
  47         nodes=mesh.p;
  48         elements=mesh.t;
  49         Nnodes = length(nodes);
  50         Nelements = length(elements);
  51
  52         % Set list of nodes with Dirichelet BCs
  53         Dnodes = Unodesonside(mesh,Dsides);
  54
  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);
  59
  60
  61         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  62         %%
  63         %%              initialization:
  64         %%              we're going to solve
  65         %%              $$ - \lambda^2 (\delta V)''
  66         %%                              + (\frac{\partial n}{\partial V}
  67         %%                              - \frac{\partial p}{\partial V})= -R $$
  68         %%
  69         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  70
  71         %% set $$ n_1 = nin $$ and $$ V = Vin $$
  72         V = Vin;
  73         Fn = Fnin;
  74         Fp = Fpin;
  75         n = exp((V(Sinodes)-Fn)./Vthn);
  76         p = exp((-V(Sinodes)+Fp)./Vthp);
  77         n(SiDnodes) = nin(SiDnodes);
  78         p(SiDnodes) = pin(SiDnodes);
  79
  80
  81         %%%
  82         %%% Compute LHS matrices
  83         %%%
  84
  85         %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
  86         if (isempty(DDGOXNLPOISSON_LAP))
  87           coeff = l2ox * ones(Nelements,1);
  88           coeff(Sielements)=l2;
  89           DDGOXNLPOISSON_LAP = Ucomplap (mesh,coeff);
  90         end
  91
  92         %% compute $$ Mv = ( n + p)  $$
  93         %% and the (lumped) mass matrix M
  94         if (isempty(DDGOXNLPOISSON_MASS))
  95           coeffe = zeros(Nelements,1);
  96           coeffe(Sielements)=1;
  97           DDGOXNLPOISSON_MASS = Ucompmass2(mesh,ones(Nnodes,1),coeffe);
  98         end
  99         freecarr=zeros(Nnodes,1);
 100         freecarr(Sinodes)=(n./Vthn + p./Vthp);
 101         Mv      =  freecarr;
 102         M       =  DDGOXNLPOISSON_MASS*spdiag(Mv);
 103
 104         %%%
 105         %%% Compute RHS vector (-residual)
 106         %%%
 107
 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);
 113         end
 114         totcharge = zeros(Nnodes,1);
 115         totcharge(Sinodes)=(n - p - D);
 116         Tv0   = totcharge;
 117         T0    = Tv0 .* DDGOXNLPOISSON_RHS;
 118
 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;
 122
 123           %% Apply boundary conditions
 124           A (Dnodes,:) = [];
 125           A (:,Dnodes) = [];
 126           R(Dnodes)  = [];
 127
 128           %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
 129             normr(1)    =  norm(R,inf);
 130             relresnorm  = 1;
 131             reldVnorm   = 1;
 132             normrnew    = normr(1);
 133             dV          = V*0;
 134
 135             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 136             %%
 137             %% START OF THE NEWTON CYCLE
 138             %%
 139             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 140             for newtit=1:maxit
 141               if (verbose>2)
 142                 fprintf(1,'\n***\nNewton iteration: %d, reldVnorm = %e\n***\n',newtit,reldVnorm);
 143               end
 144
 145               %  A(1,end)=realmin;
 146             dV(Varnodes) =(A)\(-R);
 147             dV(Dnodes)=0;
 148
 149
 150             %%%%%%%%%%%%%%%%%%
 151             %% Start of th damping procedure
 152             %%%%%%%%%%%%%%%%%%
 153             tk = 1;
 154             for dit = 1:dampit
 155               if (verbose>2)
 156                 fprintf(1,'\ndamping iteration: %d, residual norm = %e\n',dit,normrnew);
 157               end
 158               Vnew   = V + tk * dV;
 159
 160               n = exp((Vnew(Sinodes)-Fn)./Vthn);
 161               p = exp((-Vnew(Sinodes)+Fp)./Vthp);
 162               n(SiDnodes) = nin(SiDnodes);
 163               p(SiDnodes) = pin(SiDnodes);
 164
 165
 166               %%%
 167               %%% Compute LHS matrices
 168               %%%
 169
 170               %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 171               %L      = Ucomplap (mesh,ones(Nelements,1));
 172
 173               %% compute $$ Mv =  ( n + p)  $$
 174               %% and the (lumped) mass matrix M
 175               freecarr=zeros(Nnodes,1);
 176               freecarr(Sinodes)=(n./Vthn + p./Vthp);
 177               Mv   =  freecarr;
 178               M    =  DDGOXNLPOISSON_MASS*spdiag(Mv);
 179
 180               %%%
 181               %%% Compute RHS vector (-residual)
 182               %%%
 183
 184               %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
 185               totcharge( Sinodes)=(n - p - D);
 186               Tv0   = totcharge;
 187               T0    = Tv0 .* DDGOXNLPOISSON_RHS;%T0    = Ucompconst (mesh,Tv0,ones(Nelements,1));
 188
 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;
 192
 193                 %% Apply boundary conditions
 194                 A (Dnodes,:) = [];
 195                 A (:,Dnodes) = [];
 196                 R(Dnodes)  = [];
 197
 198                 %% compute $$ | R_{k+1} | $$ for the convergence test
 199                   normrnew= norm(R,inf);
 200
 201                   % check if more damping is needed
 202                     if (normrnew > normr(newtit))
 203                       tk = tk/dampcoeff;
 204                     else
 205                       if (verbose>2)
 206                         fprintf(1,'\nexiting damping cycle because residual norm = %e \n-----------\n',normrnew);
 207                       end
 208                       break
 209                     end
 210                   end
 211
 212                   V                         = Vnew;
 213                   normr(newtit+1)       = normrnew;
 214                   dVnorm              = norm(tk*dV,inf);
 215                   pause(.1);
 216                   % check if convergence has been reached
 217                     reldVnorm           = dVnorm / norm(V,inf);
 218                     if (reldVnorm <= toll)
 219                       if(verbose>2)
 220                         fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
 221                       end
 222                       break
 223                     end
 224
 225                   end
 226
 227                   res = normr;
 228                   niter = newtit;
 229