octave-secs2d 0.0.8-5 / usr / share / octave / packages / secs2d-0.0.8 / ThDDGOX / 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;