octave_packages/secs2d-0.0.8/DDGOX/DDGOXnlpoisson.m

   1 function [V,n,p,res,niter] = DDGOXnlpoisson (mesh,Dsides,Sinodes,SiDnodes,...
   2                                              Sielements,Vin,nin,pin,...
   3                                              Fnin,Fpin,D,l2,l2ox,...
   4                                              toll,maxit,verbose)
   5
   6 %
   7 %   [V,n,p,res,niter] = DDGOXnlpoisson (mesh,Dsides,Sinodes,Vin,nin,pin,...
   8 %                                             Fnin,Fpin,D,l2,l2ox,toll,maxit,verbose)
   9 %
  10 %  solves $$ -\lambda^2 V'' + (n(V,Fn) - p(V,Fp) -D)$$
  11 %
  12
  13
  14 % This file is part of
  15 %
  16 %            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  17 %         -------------------------------------------------------------------
  18 %            Copyright (C) 2004-2006  Carlo de Falco
  19 %
  20 %
  21 %
  22 %  SECS2D is free software; you can redistribute it and/or modify
  23 %  it under the terms of the GNU General Public License as published by
  24 %  the Free Software Foundation; either version 2 of the License, or
  25 %  (at your option) any later version.
  26 %
  27 %  SECS2D is distributed in the hope that it will be useful,
  28 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
  29 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  30 %  GNU General Public License for more details.
  31 %
  32 %  You should have received a copy of the GNU General Public License
  33 %  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
  34
  35 global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS
  36
  37 %% Set some useful constants
  38 dampit          = 10;
  39 dampcoeff       = 2;
  40
  41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  42 %% convert input vectors to columns
  43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  44
  45 if Ucolumns(D)>Urows(D)
  46         D=D';
  47 end
  48
  49
  50 if Ucolumns(nin)>Urows(nin)
  51         nin=nin';
  52 end
  53
  54 if Ucolumns(pin)>Urows(pin)
  55         pin=pin';
  56 end
  57
  58 if Ucolumns(Vin)>Urows(Vin)
  59         Vin=Vin';
  60 end
  61
  62 if Ucolumns(Fnin)>Urows(Fnin)
  63         Fnin=Fnin';
  64 end
  65
  66 if Ucolumns(Fpin)>Urows(Fpin)
  67         Fpin=Fpin';
  68 end
  69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  70 %% setup FEM data structures
  71 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  72
  73 nodes=mesh.p;
  74 elements=mesh.t;
  75 Nnodes = length(nodes);
  76 Nelements = length(elements);
  77
  78 % Set list of nodes with Dirichelet BCs
  79 Dnodes = Unodesonside(mesh,Dsides);
  80
  81 % Set values of Dirichelet BCs
  82 Bc     = zeros(length(Dnodes),1);
  83 % Set list of nodes without Dirichelet BCs
  84 Varnodes = setdiff([1:Nnodes],Dnodes);
  85
  86
  87 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  88 %%
  89 %%              initialization:
  90 %%              we're going to solve
  91 %%              $$ - \lambda^2 (\delta V)'' +  (\frac{\partial n}{\partial V} - \frac{\partial p}{\partial V})= -R $$
  92 %%
  93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  94
  95 %% set $$ n_1 = nin $$ and $$ V = Vin $$
  96 V = Vin;
  97 Fn = Fnin;
  98 Fp = Fpin;
  99 n = exp(V(Sinodes)-Fn);
 100 p = exp(-V(Sinodes)+Fp);
 101 n(SiDnodes) = nin(SiDnodes);
 102 p(SiDnodes) = pin(SiDnodes);
 103
 104
 105 %%%
 106 %%% Compute LHS matrices
 107 %%%
 108
 109 %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 110 if (isempty(DDGOXNLPOISSON_LAP))
 111         coeff = l2ox * ones(Nelements,1);
 112         coeff(Sielements)=l2;
 113         DDGOXNLPOISSON_LAP = Ucomplap (mesh,coeff);
 114 end
 115
 116 %% compute $$ Mv = ( n + p)  $$
 117 %% and the (lumped) mass matrix M
 118 if (isempty(DDGOXNLPOISSON_MASS))
 119     coeffe = zeros(Nelements,1);
 120     coeffe(Sielements)=1;
 121         DDGOXNLPOISSON_MASS = Ucompmass2(mesh,ones(Nnodes,1),coeffe);
 122 end
 123 freecarr=zeros(Nnodes,1);
 124 freecarr(Sinodes)=(n + p);
 125 Mv      =  freecarr;
 126 M       =  DDGOXNLPOISSON_MASS*spdiag(Mv);
 127
 128 %%%
 129 %%% Compute RHS vector (-residual)
 130 %%%
 131
 132 %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
 133 if (isempty(DDGOXNLPOISSON_RHS))
 134     coeffe = zeros(Nelements,1);
 135     coeffe(Sielements)=1;
 136         DDGOXNLPOISSON_RHS = Ucompconst (mesh,ones(Nnodes,1),coeffe);
 137 end
 138 totcharge = zeros(Nnodes,1);
 139 totcharge(Sinodes)=(n - p - D);
 140 Tv0   = totcharge;
 141 T0    = Tv0 .* DDGOXNLPOISSON_RHS;
 142
 143 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 144 A               = DDGOXNLPOISSON_LAP + M;
 145 R               = DDGOXNLPOISSON_LAP * V  + T0;
 146
 147 %% Apply boundary conditions
 148 A (Dnodes,:) = [];
 149 A (:,Dnodes) = [];
 150 R(Dnodes)  = [];
 151
 152 %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
 153 normr(1)        =  norm(R,inf);
 154 relresnorm      = 1;
 155 reldVnorm   = 1;
 156 normrnew        = normr(1);
 157 dV          = V*0;
 158
 159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 160 %%
 161 %% START OF THE NEWTON CYCLE
 162 %%
 163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 164 for newtit=1:maxit
 165   if (verbose>0)
 166     fprintf(1,'\n***\nNewton iteration: %d, reldVnorm = %e\n***\n',newtit,reldVnorm);
 167   end
 168
 169 %  A(1,end)=realmin;
 170   dV(Varnodes) =(A)\(-R);
 171   dV(Dnodes)=0;
 172
 173
 174   %%%%%%%%%%%%%%%%%%
 175         %% Start of th damping procedure
 176         %%%%%%%%%%%%%%%%%%
 177         tk = 1;
 178         for dit = 1:dampit
 179                 if (verbose>0)
 180                         fprintf(1,'\ndamping iteration: %d, residual norm = %e\n',dit,normrnew);
 181                 end
 182                 Vnew   = V + tk * dV;
 183
 184                 n = exp(Vnew(Sinodes)-Fn);
 185                 p = exp(-Vnew(Sinodes)+Fp);
 186                 n(SiDnodes) = nin(SiDnodes);
 187                 p(SiDnodes) = pin(SiDnodes);
 188
 189
 190                 %%%
 191                 %%% Compute LHS matrices
 192                 %%%
 193
 194                 %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 195                 %L      = Ucomplap (mesh,ones(Nelements,1));
 196
 197                 %% compute $$ Mv =  ( n + p)  $$
 198                 %% and the (lumped) mass matrix M
 199                 freecarr=zeros(Nnodes,1);
 200                 freecarr(Sinodes)=(n + p);
 201                 Mv   =  freecarr;
 202                 M       =  DDGOXNLPOISSON_MASS*spdiag(Mv);
 203
 204                 %%%
 205                 %%% Compute RHS vector (-residual)
 206                 %%%
 207
 208                 %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
 209                 totcharge( Sinodes)=(n - p - D);
 210                 Tv0   = totcharge;
 211                 T0    = Tv0 .* DDGOXNLPOISSON_RHS;%T0    = Ucompconst (mesh,Tv0,ones(Nelements,1));
 212
 213                 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 214                 A               = DDGOXNLPOISSON_LAP + M;
 215                 R               = DDGOXNLPOISSON_LAP * Vnew  + T0;
 216
 217                 %% Apply boundary conditions
 218                 A (Dnodes,:) = [];
 219                 A (:,Dnodes) = [];
 220                 R(Dnodes)  = [];
 221
 222                 %% compute $$ | R_{k+1} | $$ for the convergence test
 223                 normrnew= norm(R,inf);
 224
 225                 % check if more damping is needed
 226                 if (normrnew > normr(newtit))
 227                         tk = tk/dampcoeff;
 228                 else
 229                         if (verbose>0)
 230                                 fprintf(1,'\nexiting damping cycle because residual norm = %e \n-----------\n',normrnew);
 231                         end
 232                         break
 233                 end
 234         end
 235
 236         V                           = Vnew;
 237         normr(newtit+1)         = normrnew;
 238         dVnorm              = norm(tk*dV,inf);
 239         pause(.1);
 240     % check if convergence has been reached
 241         reldVnorm           = dVnorm / norm(V,inf);
 242         if (reldVnorm <= toll)
 243                 if(verbose>0)
 244                         fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
 245                 end
 246                 break
 247         end
 248
 249 end
 250
 251 res = normr;
 252 niter = newtit;
 253