octave_packages/secs1d-0.0.8/DDN/DDNnewtonmap.m

   1 function [odata,it,res] = DDNnewtonmap (x,idata,toll,maxit,verbose)
   2
   3 %
   4 % [odata,it,res] = DDNnewtonmap (x,idata,toll,maxit,verbose)
   5 %
   6 %     Solves the scaled stationary bipolar DD equation system
   7 %     using a coupled Newton algorithm
   8 %
   9 %     input: x              spatial grid
  10 %            idata.D        doping profile
  11 %            idata.p        initial guess for hole concentration
  12 %            idata.n        initial guess for electron concentration
  13 %            idata.V        initial guess for electrostatic potential
  14 %            idata.Fn       initial guess for electron Fermi potential
  15 %            idata.Fp       initial guess for hole Fermi potential
  16 %            idata.l2       scaled electric permittivity (diffusion coefficient in Poisson equation)
  17 %            idata.un       scaled electron mobility
  18 %            idata.up       scaled electron mobility
  19 %            idata.nis      scaled intrinsic carrier density
  20 %            idata.tn       scaled electron lifetime
  21 %            idata.tp       scaled hole lifetime
  22 %            toll           tolerance for Newton iterarion convergence test
  23 %            maxit          maximum number of Newton iterarions
  24 %            verbose        verbosity level: 0,1,2
  25 %
  26 %     output: odata.n     electron concentration
  27 %             odata.p     hole concentration
  28 %             odata.V     electrostatic potential
  29 %             odata.Fn    electron Fermi potential
  30 %             odata.Fp    hole Fermi potential
  31 %             it          number of Newton iterations performed
  32 %             res         residual at each step
  33
  34 ## This file is part of
  35 ##
  36 ## SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
  37 ## -------------------------------------------------------------------
  38 ## Copyright (C) 2004-2007  Carlo de Falco
  39 ##
  40 ##
  41 ##
  42 ##  SECS1D is free software; you can redistribute it and/or modify
  43 ##  it under the terms of the GNU General Public License as published by
  44 ##  the Free Software Foundation; either version 2 of the License, or
  45 ##  (at your option) any later version.
  46 ##
  47 ##  SECS1D is distributed in the hope that it will be useful,
  48 ##  but WITHOUT ANY WARRANTY; without even the implied warranty of
  49 ##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  50 ##  GNU General Public License for more details.
  51 ##
  52 ##  You should have received a copy of the GNU General Public License
  53 ##  along with SECS1D; If not, see <http://www.gnu.org/licenses/>.
  54
  55 odata = idata;
  56
  57 Nnodes=rows(x);
  58
  59 Nelements=Nnodes-1;
  60 elements=[1:Nnodes-1;2:Nnodes]';
  61 BCnodesp = [1,Nnodes];
  62 BCnodes = [BCnodesp,BCnodesp+Nnodes,BCnodesp+2*Nnodes];
  63 totaldofs= Nnodes-2;
  64 dampcoef = 2;
  65 maxdamp  = 2;
  66
  67 V = idata.V;
  68 n = idata.n;
  69 p = idata.p;
  70 D = idata.D;
  71
  72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  73 %% create the complete unknown vector
  74 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  75 u = [V; n; p];
  76
  77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  78 %% build fem matrices
  79 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  80
  81 L = Ucomplap (x,Nnodes,elements,Nelements,idata.l2*ones(Nelements,1));
  82 M = Ucompmass (x,Nnodes,elements,Nelements,ones(Nnodes,1),ones(Nelements,1));
  83 DDn = Uscharfettergummel(x,Nnodes,elements,Nelements,idata.un,1,V);
  84 DDp = Uscharfettergummel(x,Nnodes,elements,Nelements,idata.up,1,-V);
  85
  86 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  87 %% Initialise RHS
  88 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  89 r1  = L * V + M * (n - p - D);
  90 r2  = DDn * n;
  91 r3  = DDp * p;
  92
  93 RHS =- [...
  94 r1;
  95 r2;
  96 r3
  97 ];
  98
  99 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 100 %  Apply BCs
 101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 102 RHS(BCnodes,:)= [];
 103
 104 nrm = norm(RHS,inf);
 105 res(1) = nrm;
 106
 107 %%%%%%%%%%%%%%%%%%%%%%%
 108 %% Begin Newton Cycle
 109 %%%%%%%%%%%%%%%%%%%%%%%
 110 for count = 1: maxit
 111   if verbose
 112     fprintf (1,'\n\n\nNewton Iteration Number:%d\n',count);
 113   end
 114     Ln = Ucomplap (x,Nnodes,elements,Nelements,Umediaarmonica(idata.un*n));
 115     Lp = Ucomplap (x,Nnodes,elements,Nelements,Umediaarmonica(idata.up*p));
 116     Z  = sparse(zeros(Nnodes));
 117     DDn = Uscharfettergummel(x,Nnodes,elements,Nelements,idata.un,1,V);
 118     DDp = Uscharfettergummel(x,Nnodes,elements,Nelements,idata.up,1,-V);
 119
 120     A   = L;
 121     B   = M;
 122     C   = -M;
 123     DDD = -Ln;
 124     E   = DDn;
 125     F   = Z;
 126     G   = Lp;
 127     H   = Z;
 128     I   = DDp;
 129
 130     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 131     %% Build LHS
 132     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 133     LHS =sparse([
 134         [A      B C];
 135         [DDD    E F];
 136         [G      H I];
 137     ]);
 138
 139     %Apply BCs
 140     LHS(BCnodes,:)=[];
 141     LHS(:,BCnodes)=[];
 142
 143     %%%%%%%%%%%%%%%%%%%%%%%
 144     % Solve the linearised system
 145     %%%%%%%%%%%%%%%%%%%%%%%
 146     dutmp = (LHS) \ (RHS);
 147     dv    = dutmp(1:totaldofs);
 148     dn    = dutmp(totaldofs+1:2*totaldofs);
 149     dp    = dutmp(2*totaldofs+1:3*totaldofs);
 150     du    = [0;dv;0;0;dn;0;0;dp;0];
 151
 152     %%%%%%%%%%%%%%%%%%%%%%%
 153     %% Check Convergence
 154     %%%%%%%%%%%%%%%%%%%%%%%
 155
 156     nrm_u = norm(u,inf);
 157
 158     nrm_du = norm(du,inf);
 159
 160     ratio = nrm_du/nrm_u;
 161     if verbose
 162       fprintf (1,'ratio = %e\n', ratio);
 163     end
 164
 165     if (ratio <= toll)
 166         V           = u(1:Nnodes);
 167         n           = u(Nnodes+1:2*Nnodes);
 168         p           = u(2*Nnodes+1:end);
 169         res(count)  = nrm;
 170         break;
 171     end
 172
 173
 174     %%%%%%%%%%%%%%%%%%%%%%
 175     %% begin damping cycle
 176     %%%%%%%%%%%%%%%%%%%%%%
 177     tj = 1;
 178
 179
 180     for cc = 1:maxdamp
 181       if verbose
 182         fprintf (1,'\ndamping iteration number:%d\n',cc);
 183         fprintf (1,'reference residual norm:%e\n',nrm);
 184       end
 185         %%%%%%%%%%%%%%%%%%%%%%%%%
 186         %       Update the unknown vector
 187         %%%%%%%%%%%%%%%%%%%%%%%%%
 188         utmp    = u + tj*du;
 189         Vnew        = utmp(1:Nnodes);
 190         nnew        = utmp(Nnodes+1:2*Nnodes);
 191         pnew        = utmp(2*Nnodes+1:end);
 192
 193         %%%%%%%%%%%%%%%%%
 194         %% try a new RHS
 195         %%%%%%%%%%%%%%%%
 196         DDn = Uscharfettergummel(x,Nnodes,elements,Nelements,idata.un,1,Vnew);
 197         DDp = Uscharfettergummel(x,Nnodes,elements,Nelements,idata.up,1,-Vnew);
 198
 199         r1  = L * V + M * (nnew - pnew - D);
 200         r2  = DDn * nnew;
 201         r3  = DDp * pnew;
 202
 203         RHS =- [...
 204         r1;
 205         r2;
 206         r3
 207         ];
 208
 209         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 210         %  Apply BCs
 211         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 212         RHS(BCnodes,:)= [];
 213
 214         nrmtmp=norm(RHS,inf);
 215
 216         %%%%%%%%%%%%%%%%%%%%%%%%%
 217         %% Update the damping coefficient
 218         %%%%%%%%%%%%%%%%%%%%%%%%
 219         if verbose
 220           fprintf(1,'residual norm:%e\n\n', nrmtmp);
 221         end
 222
 223                 if (nrmtmp>nrm)
 224                         tj = tj/(dampcoef*cc);
 225                         if verbose
 226                           fprintf (1,'\ndamping coefficients = %e',tj);
 227                         end
 228         else
 229                         break;
 230         end
 231     end
 232
 233     nrm_du = norm(tj*du,inf);
 234
 235     u   = utmp;
 236
 237     if (count>1)
 238         ratio = nrm_du/nrm_du_old;
 239         if (ratio<.005)
 240             V       = u(1:Nnodes);
 241             n       = u(Nnodes+1:2*Nnodes);
 242             p       = u(2*Nnodes+1:end);
 243             res(count)  = nrm;
 244             break;
 245         end
 246     end
 247
 248     nrm = nrmtmp;
 249
 250     res(count)  = nrm;
 251     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 252     %% convert result vector into distinct output vectors
 253     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 254
 255     V       = u(1:Nnodes);
 256     n       = u(Nnodes+1:2*Nnodes);
 257     p       = u(2*Nnodes+1:end);
 258
 259     nrm_du_old = nrm_du;
 260 end
 261
 262 odata.V = V;
 263 odata.n = n;
 264 odata.p = p;
 265 Fn   = V - log(n);
 266 Fp   = V + log(p);
 267
 268 it   = count;
 269
 270
 271
 272