octave_packages/secs2d-0.0.8/QDDGOX/QDDGOXgummelmap.m

   1 function [odata,it,res] = QDDGOXgummelmap (imesh,Dsides,...
   2     Simesh,Sinodes,Sielements,SiDsides,Intsides,...
   3     idata,toll,maxit,ptoll,pmaxit,stoll,smaxit,verbose,options)
   4
   5 % [odata,it,res] = QDDGOXgummelmap (imesh,Dsides,...
   6 %     Simesh,Sinodes,Sielements,SiDsides,Intsides,...
   7 %     idata,toll,maxit,ptoll,pmaxit,stoll,smaxit,verbose,options);
   8
   9
  10 clear global
  11 global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS DDG_RHS DDG_MASS
  12
  13 %%%%%%%%%%%%%%%
  14 %% RRE param %%
  15 RREnnit      = [1,0];
  16 RRErank      = 5;
  17 RREpattern   = URREcyclingpattern(RREnnit,RRErank,maxit);
  18
  19 RREnnit2      = [1,0];
  20 RRErank2      = 4;
  21 RREpattern2    = URREcyclingpattern(RREnnit2,RRErank2,smaxit);
  22 %%%%%%%%%%%%%%%
  23
  24 Nnodes      = max(size(imesh.p));
  25 Nelements   = max(size(imesh.t));
  26 SiNnodes    = max(size(Simesh.p));
  27 SiNelements = max(size(Simesh.t));
  28
  29 if (options.SRH==1)
  30     tn= idata.tn;tp=idata.tp;
  31 else
  32     tn=Inf;tp=Inf;
  33 end
  34
  35 V (:,1) = idata.V;
  36 p (:,1) = idata.p;
  37 n (:,1) = idata.n;
  38 Fn(:,1) = idata.Fn;
  39 Fp(:,1) = idata.Fp;
  40
  41 D       = idata.D;
  42
  43 Dedges    =[];
  44
  45 for ii = 1:length(Dsides)
  46     Dedges=[Dedges,find(imesh.e(5,:)==Dsides(ii))];
  47 end
  48
  49 % Set list of nodes with Dirichelet BCs
  50 Dnodes = imesh.e(1:2,Dedges);
  51 Dnodes = [Dnodes(1,:) Dnodes(2,:)];
  52 Dnodes = unique(Dnodes);
  53
  54 Intedges    =[];
  55
  56 for ii = 1:length(Intsides)
  57     Intedges=[Intedges,find(Simesh.e(5,:)==Intsides(ii))];
  58 end
  59
  60 % Set list of Interface nodes
  61 Intnodes = Simesh.e(1:2,Intedges);
  62 Intnodes = [Intnodes(1,:) Intnodes(2,:)];
  63 Intnodes = unique(Intnodes);
  64
  65 SiDedges    =[];
  66
  67 for ii = 1:length(SiDsides)
  68     SiDedges=[SiDedges,find(Simesh.e(5,:)==SiDsides(ii))];
  69 end
  70
  71 % Set list of nodes with Dirichelet BCs
  72 SiDnodes = Simesh.e(1:2,SiDedges);
  73 SiDnodes = [SiDnodes(1,:) SiDnodes(2,:)];
  74 SiDnodes = unique(SiDnodes);
  75
  76 if (options.FD==1)
  77     FDn  = idata.FDn;
  78     FDp  = idata.FDp;
  79 else
  80     FDn  = zeros(SiNnodes,1);
  81     FDp  = zeros(SiNnodes,1);
  82 end
  83
  84 G (:,1) =  Fn(:,1) - V(Sinodes,1) - FDn + log(n(:,1));
  85 if (options.holes==1)
  86     Gp (:,1) = Fp(:,1) - V(Sinodes,1) - FDp - log(p(:,1));
  87 else
  88     Gp (:,1) = G(:,1)*0;
  89 end
  90
  91
  92
  93 nrm=1;
  94 for i=1:1:maxit
  95     if (verbose>=1)
  96         fprintf(1,'*****************************************************************\n');
  97         fprintf(1,'****    start of gummel iteration number: %d\n',i);
  98         fprintf(1,'*****************************************************************\n');
  99     end
 100
 101     V(:,2)= V(:,1);
 102     G(:,2)= G(:,1);
 103     Gp(:,2)= Gp(:,1);
 104     n(:,2)= n(:,1);
 105     bohmdeltav=inf;
 106     for j=1:smaxit
 107         if (verbose>=1)
 108             fprintf(1,'*---------------------------------------------------------------*\n');
 109             fprintf(1,'****    start of Poisson-Bohm iteration number: %d (bohmdeltav=%g)\n',j,bohmdeltav);
 110             fprintf(1,'*---------------------------------------------------------------*\n');
 111         end
 112
 113
 114         if (verbose>1)
 115             fprintf(1,'solving non linear poisson equation\n\n');
 116         end
 117
 118         [V(:,3),n(:,2),p(:,2)] =...
 119             QDDGOXnlpoisson (imesh,Dsides,Sinodes,[SiDnodes,Intnodes] ,Sielements,...
 120             V(:,2),n(:,1),p(:,1),Fn(:,1),Fp(:,1),G(:,2)+FDn,Gp(:,2)+FDp,D,...
 121             idata.l2,idata.l2ox,ptoll,pmaxit,verbose-1);
 122
 123         n([SiDnodes,Intnodes],2) = idata.n([SiDnodes,Intnodes]);
 124         p([SiDnodes,Intnodes],2) = idata.p([SiDnodes,Intnodes]);
 125         V(Dnodes,3) = idata.V(Dnodes);
 126
 127         if (verbose>1)
 128             fprintf(1,'solving non linear Bohm equation for electrons\n\n');
 129         end
 130         n(Intnodes,2) = idata.n(Intnodes);
 131
 132
 133         w       = QDDGOXcompdens(Simesh,[SiDsides,Intsides],sqrt(n(:,2)),V(Sinodes,3) + FDn,Fn(:,1),idata.dn2,stoll,smaxit,verbose-1);
 134         n(:,2)  = w.^2;
 135         n([SiDnodes,Intnodes],2) = idata.n([SiDnodes,Intnodes]);
 136         G(:,3)  = Fn(:,1) - V(Sinodes,3) - FDn + log(n(:,2));
 137         if (verbose>1)
 138             fprintf(1,'solving non linear Bohm equation for holes\n\n');
 139         end
 140
 141         if (options.holes==1)
 142
 143             p(Intnodes,2) = idata.p(Intnodes);
 144             wp      = QDDGOXcompdens(Simesh,[SiDsides,Intsides],sqrt(p(:,2)),-V(Sinodes,3) - FDp,...
 145                 -Fp(:,1),idata.dp2,ptoll,pmaxit,verbose-1);
 146             p(:,2)  = wp.^2;
 147             p([SiDnodes,Intnodes],2) = idata.p([SiDnodes,Intnodes]);
 148             Gp(:,3) = Fp(:,1) - V(Sinodes,3) - FDp - log(p(:,2));
 149
 150         else
 151             Gp(:,3)=G(:,3)*0;
 152         end
 153
 154
 155         if (options.FD==1)
 156             fprintf(1,'\n*** APPLYING FD STATISTICS ***\n')
 157             n(:,2) = idata.Nc*Ufermidirac(V(Sinodes,3)+G(:,3)-Fn(:,1)-log(idata.Nc),1/2);
 158             n(SiDnodes,2) = idata.n(SiDnodes);
 159             nMBtmp = exp(V(Sinodes,3)+G(:,3)-Fn(:,1));
 160             FDn    = log(n(:,2)./ nMBtmp);
 161             FDn(SiDnodes) = idata.FDn(SiDnodes);
 162
 163             p(:,2) = idata.Nv*Ufermidirac(-V(Sinodes,3)-Gp(:,3)+Fp(:,1)-log(idata.Nv),1/2);
 164             p([SiDnodes,Intnodes],2) = idata.p([SiDnodes,Intnodes]);
 165             pMBtmp = exp(-V(Sinodes,3)-Gp(:,3)+Fp(:,1));
 166             FDp    = -log(p(:,2)./ pMBtmp);
 167             FDp(SiDnodes) = idata.FDp(SiDnodes);
 168         end
 169
 170         bohmdeltav = norm(G(:,3)-G(:,2),inf) +...
 171             norm(Gp(:,3)-Gp(:,2),inf) +...
 172             norm(V(:,3)-V(:,2),inf);
 173
 174
 175         %%%% store result for RRE
 176         if RREpattern2(j)>0
 177             Gstore(:,RREpattern2(j)) = G(:,3);
 178             if RREpattern2(j+1)==0 % Apply RRE extrapolation
 179                 G(:,3) = Urrextrapolation(Gstore);
 180             end
 181         end
 182
 183         G(:,2)=G(:,3);
 184         Gp(:,2)=Gp(:,3);
 185
 186
 187         V(:,2)=V(:,3);
 188
 189
 190         if (bohmdeltav<=stoll)
 191             if (verbose>1)
 192                 fprintf(1,'Exiting poisson-bohm iteration because bohmdeltav=%g\n\n',bohmdeltav);
 193             end
 194             break;
 195         end
 196     end
 197
 198     if (verbose>1)
 199         fprintf (1,'\n\nupdating electron qfl\n\n');
 200     end
 201
 202     mob = Ufielddepmob(Simesh,idata.un,Fn(:,1), ...
 203         idata.vsatn,idata.mubn);
 204
 205
 206     n(:,3) = DDGOXelectron_driftdiffusion(Simesh,SiDsides,n(:,2),p(:,2),...
 207                                           V(Sinodes,3)+G(:,3)+FDn,mob,...
 208                                           tn,tp,idata.n0,idata.p0);
 209
 210
 211     Fn(:,2)                  = V(Sinodes,3) + G(:,3) + FDn - log(n(:,3));
 212     Fn(SiDnodes,2)           = idata.Fn(SiDnodes);
 213     n([SiDnodes,Intnodes],3) = idata.n([SiDnodes,Intnodes]);
 214
 215     %%%% store result for RRE
 216     if RREpattern(i)>0
 217         Fnstore(:,RREpattern(i)) = Fn(:,2);
 218         if RREpattern(i+1)==0 % Apply RRE extrapolation
 219             Fn(:,2) = Urrextrapolation(Fnstore);
 220         end
 221     end
 222
 223     if (verbose>1)
 224         fprintf(1,'updating hole qfl\n\n');
 225     end
 226
 227     mob = Ufielddepmob(Simesh,idata.up,Fp(:,1),idata.vsatp,idata.mubp);
 228     p(:,3) =DDGOXhole_driftdiffusion(Simesh,SiDsides,n(:,3),p(:,2),...
 229         V(Sinodes,3)+Gp(:,3)+FDp,mob,...
 230         tn,tp,idata.n0,idata.p0);
 231
 232
 233     if (options.holes==1)
 234         Fp(:,2)=V(Sinodes,3) + Gp(:,3) + FDp + log(p(:,3));
 235         p([SiDnodes,Intnodes],3) = idata.p([SiDnodes,Intnodes]);
 236     else
 237         Fp(:,2)=Fn(:,2) + 2 * log(idata.ni);
 238         p(:,3) = exp(Fp(:,2)-V(Sinodes,3)-FDp);
 239         p([SiDnodes],3) = idata.p([SiDnodes]);
 240     end
 241     Fp(SiDnodes,2) = idata.Fp(SiDnodes);
 242
 243     if (verbose>1)
 244         fprintf(1,'checking for convergence\n\n');
 245     end
 246
 247     nrfn= norm(Fn(:,2)-Fn(:,1),inf);
 248     nrfp= norm (Fp(:,2)-Fp(:,1),inf);
 249     nrv = norm (V(:,3)-V(:,1),inf);
 250     nrg = norm (G(:,3)-G(:,1),inf);
 251     nrgp = norm (Gp(:,3)-Gp(:,1),inf);
 252     nrm(i) = max([nrfn;nrfp;nrv;nrg;nrgp]);
 253
 254     figure(2)
 255     semilogy(nrm)
 256     pause(.1)
 257
 258
 259     if (verbose>1)
 260         fprintf (1,' max(|phin_(k+1)-phinn_(k)| , |phip_(k+1)-phip_(k)| , |v_(k+1)-v_(k)|  |g_(k+1)-g_(k)|)= %d\n',nrm(i));
 261     end
 262     if (nrm(i)<toll)
 263         break
 264     end
 265
 266     V(:,1) = V(:,end);
 267     G(:,1) = G(:,end);
 268     Gp(:,1) = Gp(:,end);
 269     n(:,1) = n(:,end);
 270     p(:,1) = p(:,end);
 271     Fn(:,1)= Fn(:,end);
 272     Fp(:,1)= Fp(:,end);
 273
 274
 275 end
 276
 277 it = i;
 278 res = nrm
 279
 280 if (verbose>0)
 281     fprintf(1,'\n\nDD: # of Gummel iterations = %d\n\n',it);
 282 end
 283
 284 odata = idata;
 285
 286 odata.n  = n(:,end);
 287 odata.p  = p(:,end);
 288 odata.V  = V(:,end);
 289 odata.Fn = Fn(:,end);
 290 odata.Fp = Fp(:,end);
 291 odata.G  = G(:,end);
 292 odata.Gp  = Gp(:,end);
 293
 294