Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / secs2d-0.0.8 / ThDDGOX / ThDDGOXthermaliteration.m

diff --git a/octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXthermaliteration.m b/octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXthermaliteration.m
new file mode 100644 (file)
index 0000000..abb8ba8
--- /dev/null
+++ b/octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXthermaliteration.m
@@ -0,0 +1,156 @@
+function [thermdata,nrm] = ThDDGOXthermaliteration(imesh,Dsides,...
+                                                  Simesh,Sinodes,Sielements,...
+                                                  SiDsides,thermdata,toll,...
+                                                  maxit,verbose)
+  
+
+  %%  [thermdata,innrm] = ThDDGOXthermaliteration(imesh,Dsides,...
+  %%                                             Simesh,Sinodes,Sielements,...
+  %%                                              SiDsides,thermdata,toll,...
+  %%                                              maxit,verbose)
+
+  %%%%%%%%%%%%%%%
+  %% RRE param %%
+  RREnnit      = [10,2];
+  RRErank      = maxit;  
+  RREpattern   = URREcyclingpattern(RREnnit,RRErank,maxit);
+  %%%%%%%%%%%%%%%
+
+  %% Set list of nodes with Dirichlet BCs
+  Dnodes = Unodesonside(imesh,Dsides);
+  SiDnodes = Unodesonside(Simesh,SiDsides);
+
+  Tl  = thermdata.Tl;
+  Tn  = thermdata.Tn;
+  Tp  = thermdata.Tp;
+  
+  tldampcoef = 1;
+  tndampcoef = 10;
+  tpdampcoef = 10;
+  
+  mobn0 = thermdata.mobn0(imesh,Simesh,Sinodes,Sielements,thermdata);
+  mobp0 = thermdata.mobp0(imesh,Simesh,Sinodes,Sielements,thermdata);
+  mobn1 = thermdata.mobn1(imesh,Simesh,Sinodes,Sielements,thermdata);
+  mobp1 = thermdata.mobp1(imesh,Simesh,Sinodes,Sielements,thermdata);
+  twn0  = thermdata.twn0 (imesh,Simesh,Sinodes,Sielements,thermdata);
+  twp0  = thermdata.twp0 (imesh,Simesh,Sinodes,Sielements,thermdata);
+  twn1  = thermdata.twn1 (imesh,Simesh,Sinodes,Sielements,thermdata);
+  twp1  = thermdata.twp1 (imesh,Simesh,Sinodes,Sielements,thermdata);
+  [Ex,Ey] = Updegrad(Simesh,-thermdata.V(Sinodes));
+  E =  [Ex;Ey];  
+
+  [jnx,jny] = Ufvsgcurrent3(Simesh,thermdata.n,...
+                           mobn0,mobn1,Tn,thermdata.V(Sinodes)-Tn);
+  [jpx,jpy] = Ufvsgcurrent3(Simesh,thermdata.p,...
+                           -mobp0,mobp1,Tp,-thermdata.V(Sinodes)-Tp);
+
+  Jn = [jnx;jny];
+  Jp = [jpx;jpy];
+  
+  for ith=1:maxit
+    
+    if (verbose>=1)
+      fprintf(1,"*** start of inner iteration number: %d\n",ith);
+    end
+
+    if (verbose>=1)
+      fprintf(1,'\t***updating electron temperature\n');
+    end
+    
+    Tn =  ThDDGOXupdateelectron_temp(Simesh,SiDnodes,thermdata.Tn,...
+                                    thermdata.n,thermdata.p,...
+                                    thermdata.Tl,Jn,E,mobn0,...
+                                     twn0,twn1,thermdata.tn,thermdata.tp,...
+                                    thermdata.ni,thermdata.ni);
+
+    ##Tn(Tn<thermdata.Tl) = thermdata.Tl(Tn<thermdata.Tl);
+    
+    dtn     = norm(Tn-thermdata.Tn,inf);
+    if (dtn>0) 
+      tndampfact_n = log(1+tndampcoef*dtn)/(tndampcoef*dtn);
+      Tn  = tndampfact_n * Tn + (1-tndampfact_n) * thermdata.Tn;
+    end
+    
+    if (verbose>=1)
+      fprintf(1,'\t***updating hole temperature\n');
+    end
+    
+
+    Tp = ThDDGOXupdatehole_temp(Simesh,SiDnodes,thermdata.Tp,...
+                               thermdata.n,thermdata.p,...
+                               thermdata.Tl,Jp,E,mobp0,...
+                                twp0,twp1,thermdata.tn,thermdata.tp,...
+                               thermdata.ni,thermdata.ni);
+
+    ##Tp(Tp<thermdata.Tl) = thermdata.Tl(Tp<thermdata.Tl);
+
+    dtp     = norm(Tp-thermdata.Tp,inf);
+    if (dtp>0) 
+      tpdampfact_p = log(1+tpdampcoef*dtp)/(tpdampcoef*dtp);
+      Tp  = tpdampfact_p * Tp + (1-tpdampfact_p) * thermdata.Tp;
+    end
+    
+    if (verbose>=1)
+      fprintf(1,'\t***updating lattice temperature\n');
+    end
+    
+    ## store result for RRE
+    if RREpattern(ith)>0
+      TEMPstore(:,RREpattern(ith)) = [Tn;Tp;Tl];
+      if RREpattern(ith+1)==0 % Apply RRE extrapolation
+       if (verbose>=1)         
+         fprintf(1,"\n\t**********\n\tRRE EXTRAPOLATION STEP\n\t**********\n\n");
+       end
+       TEMP = Urrextrapolation(TEMPstore);
+       Tn = TEMP(1:rows(Tn));
+       Tp = TEMP(rows(Tn)+1:rows(Tn)+rows(Tp));
+        Tl = TEMP(rows(Tn)+rows(Tp)+1:end);
+    end
+  end
+
+  Tl = ThDDGOXupdatelattice_temp(Simesh,SiDnodes,thermdata.Tl,...
+                                Tn,Tp,thermdata.n,...
+                                thermdata.p,thermdata.kappa,thermdata.Egap,...
+                                thermdata.tn,thermdata.tp,twn0,...
+                                twp0,twn1,twp1,...
+                                thermdata.ni,thermdata.ni);
+  
+  ##Tl(Tl<thermdata.Tl) = thermdata.Tl(Tl<thermdata.Tl);
+
+  dtl = norm(Tl-thermdata.Tl,inf);
+  if (dtl > 0)
+    tldampfact = log(1+tldampcoef*dtl)/(tldampcoef*dtl);
+    Tl    = tldampfact * Tl + (1-tldampfact) * thermdata.Tl; 
+  end
+  
+  
+  if (verbose>=1)
+    fprintf(1,"\t*** checking for convergence:\n ");
+  end
+  
+  nrm(ith) = max([dtl,dtn,dtp]);       
+  
+  if (verbose>=1)
+    fprintf (1,"\t\t|dTL|= %g\n",dtl);
+    fprintf (1,"\t\t|dTn|= %g\n",dtn);
+    fprintf (1,"\t\t|dTp|= %g\n",dtp);
+  end
+  
+  thermdata.Tl = Tl;
+  thermdata.Tn = Tn;
+  thermdata.Tp = Tp;
+  if (verbose>1)
+    subplot(1,3,2);
+    title("max(|dTl|,|dTn|,|dTn|)")
+    semilogy(nrm)
+    pause(.1)
+  end
+  if nrm(ith)< toll
+    if (verbose>=1)
+      fprintf(1,"\n***\n***\texit from thermal iteration \n");
+    end
+    
+    break
+  end
+  
+end
\ No newline at end of file