octave_packages/secs1d-0.0.8/DDG/DDGnlpoisson.m

   1 function [V,n,p,res,niter] = DDGnlpoisson (x,sinodes,Vin,nin,...
   2                                            pin,Fnin,Fpin,D,l2,toll,maxit,verbose)
   3
   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
  17 %             D       doping profile
  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
  27
  28 ## This file is part of
  29 ##
  30 ## SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
  31 ## -------------------------------------------------------------------
  32 ## Copyright (C) 2004-2007  Carlo de Falco
  33 ##
  34 ##
  35 ##
  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.
  40 ##
  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.
  45 ##
  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/>.
  48
  49
  50 %% Set some useful constants
  51 dampit          = 10;
  52 dampcoeff       = 2;
  53
  54
  55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  56 %% convert grid info to FEM form
  57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  58
  59 Ndiricheletnodes        = 2;
  60 nodes                   = x;
  61 sielements          = sinodes(1:end-1);
  62
  63 Nnodes                  = length(nodes);
  64
  65
  66 totdofs             = Nnodes - Ndiricheletnodes ;
  67
  68 elements            = [[1:Nnodes-1]' , [2:Nnodes]'];
  69 Nelements           = size(elements,1);
  70
  71 BCnodes             = [1;Nnodes];
  72
  73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  74 %%
  75 %%              initialization:
  76 %%              we're going to solve
  77 %%              $$ -  lamda^2*(\delta V)'' +  (\frac{\partial n}{\partial V} -
  78 %%            \frac{\partial p}{\partial V})= -R $$
  79 %%
  80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  81
  82 %% set $$ n_1 = nin $$ and $$ V = Vin $$
  83 V = Vin;
  84 Fn = Fnin;
  85 Fp = Fpin;
  86 n = DDGphin2n(V(sinodes),Fn);
  87 p = DDGphip2p(V(sinodes),Fp);
  88 if (sinodes(1)==1)
  89     n(1)=nin(1);
  90     p(1)=pin(1);
  91 end
  92 if (sinodes(end)==Nnodes)
  93     n(end)=nin(end);
  94     p(end)=pin(end);
  95 end
  96
  97 %%%
  98 %%% Compute LHS matrices
  99 %%%
 100
 101 %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 102 L      = Ucomplap (nodes,Nnodes,elements,Nelements,l2.*ones(Nelements,1));
 103
 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);
 109 Cv(sielements)=  1;
 110 M             =  Ucompmass (nodes,Nnodes,elements,Nelements,Mv,Cv);
 111
 112 %%%
 113 %%% Compute RHS vector (-residual)
 114 %%%
 115
 116 %% now compute $$ T0 =  (n - p - D) $$
 117 Tv0            =  zeros(Nnodes,1);
 118 Tv0(sinodes)   = (n - p - D);
 119 Cv            =  zeros(Nelements,1);
 120 Cv(sielements)=  1;
 121 T0             =  Ucompconst (nodes,Nnodes,elements,Nelements,Tv0,Cv);
 122
 123 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 124 A               = L + M;
 125 R               = L * V  + T0;
 126
 127 %% Apply boundary conditions
 128 A(BCnodes,:)    = [];
 129 A(:,BCnodes)    = [];
 130
 131 R(BCnodes)      = [];
 132
 133 %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
 134 normr(1)                =  norm(R,inf);
 135 relresnorm      = 1;
 136 reldVnorm   = 1;
 137 normrnew        = normr(1);
 138
 139
 140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 141 %%
 142 %% START OF THE NEWTON CYCLE
 143 %%
 144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 145 for newtit=1:maxit
 146     if verbose
 147         fprintf(1,'\n newton iteration: %d, reldVnorm = %e',newtit,reldVnorm);
 148
 149     end
 150
 151     dV =[0;(A)\(-R);0];
 152
 153
 154
 155     %%%%%%%%%%%%%%%%%%
 156     %% Start of the damping procedure
 157     %%%%%%%%%%%%%%%%%%
 158     tk = 1;
 159     for dit = 1:dampit
 160         if verbose
 161             fprintf(1,'\n damping iteration: %d, residual norm = %e',dit,normrnew);
 162         end
 163         Vnew   = V + tk * dV;
 164
 165         n = DDGphin2n(Vnew(sinodes),Fn);
 166         p = DDGphip2p(Vnew(sinodes),Fp);
 167         if (sinodes(1)==1)
 168             n(1)=nin(1);
 169             p(1)=pin(1);
 170         end
 171         if (sinodes(end)==Nnodes)
 172             n(end)=nin(end);
 173             p(end)=pin(end);
 174         end
 175         %%%
 176         %%% Compute LHS matrices
 177         %%%
 178
 179         %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 180         %L      = Ucomplap (nodes,Nnodes,elements,Nelements,ones(Nelements,1));
 181
 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);
 187         Cv(sielements)=  1;
 188         M    = Ucompmass (nodes,Nnodes,elements,Nelements,Mv,Cv);
 189
 190         %%%
 191         %%% Compute RHS vector (-residual)
 192         %%%
 193
 194         %% now compute $$ T0 =  (n - p - D) $$
 195         Tv0            =  zeros(Nnodes,1);
 196         Tv0(sinodes)   =  (n - p - D);
 197         Cv            =  zeros(Nelements,1);
 198         Cv(sielements)=  1;
 199         T0     = Ucompconst (nodes,Nnodes,elements,Nelements,Tv0,Cv);
 200
 201         %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 202         Anew            = L + M;
 203         Rnew            = L * Vnew  + T0;
 204
 205         %% Apply boundary conditions
 206         Anew(BCnodes,:) = [];
 207         Anew(:,BCnodes) = [];
 208
 209         Rnew(BCnodes)   = [];
 210
 211         if ((dit>1)&(norm(Rnew,inf)>=norm(R,inf)))
 212             if verbose
 213                 fprintf(1,'\nexiting damping cycle \n');
 214             end
 215             break
 216         else
 217             A = Anew;
 218             R = Rnew;
 219         end
 220
 221         %% compute $$ | R_{k+1} | $$ for the convergence test
 222         normrnew= norm(R,inf);
 223
 224         % check if more damping is needed
 225         if (normrnew > normr(newtit))
 226             tk = tk/dampcoeff;
 227         else
 228             if verbose
 229                 fprintf(1,'\nexiting damping cycle because residual norm = %e \n',normrnew);
 230             end
 231             break
 232         end
 233     end
 234
 235     V                       = Vnew;
 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)
 241         if(verbose)
 242             fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
 243         end
 244         break
 245     end
 246 end
 247
 248 res = normr;
 249 niter = newtit;
 250
 251 % Last Revision:
 252 % $Author: adb014 $
 253 % $Date: 2008-02-04 16:26:27 +0100 (man, 04 feb 2008) $
 254
 255