X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fsecs2d-0.0.8%2FThDDGOX%2FThDDGOXnlpoisson.m;fp=octave_packages%2Fsecs2d-0.0.8%2FThDDGOX%2FThDDGOXnlpoisson.m;h=a6d2bafd534d84d93033129aca10b3700a6cd335;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXnlpoisson.m b/octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXnlpoisson.m new file mode 100644 index 0000000..a6d2baf --- /dev/null +++ b/octave_packages/secs2d-0.0.8/ThDDGOX/ThDDGOXnlpoisson.m @@ -0,0 +1,229 @@ +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 . + + 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; +