Logo Search packages:      
Sourcecode: octave-secs2d version File versions  Download package

ThDDGOXnlpoisson.m

function [V,n,p,res,niter] = ThDDGOXnlpoisson (mesh,Dsides,Sinodes,SiDnodes,...
                                               Sielements,Vin,Vthn,Vthp,...
                                     nin,pin,...
                                               Fnin,Fpin,D,l2,l2ox,...
                                               toll,maxit,verbose)

  %%  
  %%   [V,n,p,res,niter] = DDGOXnlpoisson (mesh,Dsides,Sinodes,Vin,nin,pin,...
  %%                                       Fnin,Fpin,D,l2,l2ox,toll,maxit,verbose)
  %%
  %%  solves $$ -\lambda^2 V'' + (n(V,Fn,Tn) - p(V,Fp,Tp) -D) = 0$$
  %%


  %% This file is part of 
  %%
  %%            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  %%         -------------------------------------------------------------------
  %%            Copyright (C) 2004-2006  Carlo de Falco
  %%
  %%
  %%
  %%  SECS2D is free software; you can redistribute it and/or modify
  %%  it under the terms of the GNU General Public License as published by
  %%  the Free Software Foundation; either version 2 of the License, or
  %%  (at your option) any later version.
  %%
  %%  SECS2D is distributed in the hope that it will be useful,
  %%  but WITHOUT ANY WARRANTY; without even the implied warranty of
  %%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  %%  GNU General Public License for more details.
  %%
  %%  You should have received a copy of the GNU General Public License
  %%  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.

      global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS 

      %% Set some useful constants
      dampit            = 3;
      dampcoeff   = 2;


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %% setup FEM data structures
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      nodes=mesh.p;
      elements=mesh.t;
      Nnodes = length(nodes);
      Nelements = length(elements);

      % Set list of nodes with Dirichelet BCs
      Dnodes = Unodesonside(mesh,Dsides);

      % Set values of Dirichelet BCs
      Bc     = zeros(length(Dnodes),1);
      % Set list of nodes without Dirichelet BCs
      Varnodes = setdiff([1:Nnodes],Dnodes);


      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%
      %%          initialization:
      %%          we're going to solve
      %%          $$ - \lambda^2 (\delta V)'' 
      %%                              + (\frac{\partial n}{\partial V} 
      %%                              - \frac{\partial p}{\partial V})= -R $$
      %%
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      %% set $$ n_1 = nin $$ and $$ V = Vin $$
      V = Vin;
      Fn = Fnin;
      Fp = Fpin;
      n = exp((V(Sinodes)-Fn)./Vthn);
      p = exp((-V(Sinodes)+Fp)./Vthp);
      n(SiDnodes) = nin(SiDnodes);
      p(SiDnodes) = pin(SiDnodes);


      %%%
      %%% Compute LHS matrices
      %%%

      %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
      if (isempty(DDGOXNLPOISSON_LAP))
        coeff = l2ox * ones(Nelements,1);
        coeff(Sielements)=l2;
        DDGOXNLPOISSON_LAP = Ucomplap (mesh,coeff);
      end

      %% compute $$ Mv = ( n + p)  $$
      %% and the (lumped) mass matrix M
      if (isempty(DDGOXNLPOISSON_MASS))
        coeffe = zeros(Nelements,1);
        coeffe(Sielements)=1;
        DDGOXNLPOISSON_MASS = Ucompmass2(mesh,ones(Nnodes,1),coeffe);
      end
      freecarr=zeros(Nnodes,1);
      freecarr(Sinodes)=(n./Vthn + p./Vthp);
      Mv      =  freecarr;
      M       =  DDGOXNLPOISSON_MASS*spdiag(Mv);

      %%%
      %%% Compute RHS vector (-residual)
      %%%

      %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
      if (isempty(DDGOXNLPOISSON_RHS))
        coeffe = zeros(Nelements,1);
        coeffe(Sielements)=1;
        DDGOXNLPOISSON_RHS = Ucompconst (mesh,ones(Nnodes,1),coeffe);
      end
      totcharge = zeros(Nnodes,1);
      totcharge(Sinodes)=(n - p - D);
      Tv0   = totcharge;
      T0    = Tv0 .* DDGOXNLPOISSON_RHS;

      %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
        A         = DDGOXNLPOISSON_LAP + M;
        R         = DDGOXNLPOISSON_LAP * V  + T0; 

        %% Apply boundary conditions
        A (Dnodes,:) = [];
        A (:,Dnodes) = [];
        R(Dnodes)  = [];

        %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test   
          normr(1)      =  norm(R,inf);
          relresnorm    = 1;
          reldVnorm   = 1;
          normrnew      = normr(1);
          dV          = V*0;

          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          %%
          %% START OF THE NEWTON CYCLE
          %%
          %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
          for newtit=1:maxit
            if (verbose>2)
            fprintf(1,'\n***\nNewton iteration: %d, reldVnorm = %e\n***\n',newtit,reldVnorm);
            end

            %  A(1,end)=realmin;
          dV(Varnodes) =(A)\(-R);
          dV(Dnodes)=0;
          
          
          %%%%%%%%%%%%%%%%%%
          %% Start of th damping procedure
          %%%%%%%%%%%%%%%%%%
          tk = 1;
          for dit = 1:dampit
            if (verbose>2)
            fprintf(1,'\ndamping iteration: %d, residual norm = %e\n',dit,normrnew);
            end
            Vnew   = V + tk * dV;
            
            n = exp((Vnew(Sinodes)-Fn)./Vthn);
            p = exp((-Vnew(Sinodes)+Fp)./Vthp);
            n(SiDnodes) = nin(SiDnodes);
            p(SiDnodes) = pin(SiDnodes);
            
            
            %%%
            %%% Compute LHS matrices
            %%%
            
            %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
            %L      = Ucomplap (mesh,ones(Nelements,1));
            
            %% compute $$ Mv =  ( n + p)  $$
            %% and the (lumped) mass matrix M
            freecarr=zeros(Nnodes,1);
            freecarr(Sinodes)=(n./Vthn + p./Vthp);  
            Mv   =  freecarr;
            M    =  DDGOXNLPOISSON_MASS*spdiag(Mv);
            
            %%%
            %%% Compute RHS vector (-residual)
            %%%
            
            %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
            totcharge( Sinodes)=(n - p - D);
            Tv0   = totcharge;
            T0    = Tv0 .* DDGOXNLPOISSON_RHS;%T0    = Ucompconst (mesh,Tv0,ones(Nelements,1));
            
            %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
            A           = DDGOXNLPOISSON_LAP + M;
            R           = DDGOXNLPOISSON_LAP * Vnew  + T0; 
            
            %% Apply boundary conditions
            A (Dnodes,:) = [];
            A (:,Dnodes) = [];
            R(Dnodes)  = [];
            
            %% compute $$ | R_{k+1} | $$ for the convergence test
              normrnew= norm(R,inf);
              
              % check if more damping is needed
                if (normrnew > normr(newtit))
                  tk = tk/dampcoeff;
                else
                  if (verbose>2)
                  fprintf(1,'\nexiting damping cycle because residual norm = %e \n-----------\n',normrnew);
                  end         
                  break
                end     
              end
              
              V                     = Vnew;     
              normr(newtit+1)       = normrnew;
              dVnorm              = norm(tk*dV,inf);
              pause(.1);
              % check if convergence has been reached
                reldVnorm           = dVnorm / norm(V,inf);
                if (reldVnorm <= toll)
                  if(verbose>2)
                  fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
                  end
                  break
                end

              end

              res = normr;
              niter = newtit;


Generated by  Doxygen 1.6.0   Back to index