octave_packages/secs2d-0.0.8/QDDGOX/QDDGOXnlpoisson.m

   1 function [V,n,p,res,niter] = QDDGOXnlpoisson (mesh,Dsides,Sinodes,SiDnodes,...
   2                                               Sielements,Vin,nin,pin,...
   3                                               Fnin,Fpin,Gin,Gpin,D,l2,l2ox,...
   4                                               toll,maxit,verbose)
   5
   6 %
   7 %   [V,n,p,res,niter] = QDDGOXnlpoisson (mesh,Dsides,Sinodes,SiDnodes,...
   8 %     Sielements,Vin,nin,pin,...
   9 %     Fnin,Fpin,Gin,Gpin,D,l2,l2ox,...
  10 %     toll,maxit,verbose)
  11 %
  12 %  solves $$ -\lambda^2 V'' + (n(V,Fn) - p(V,Fp) -D)$$
  13 %
  14
  15 global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS
  16
  17 %% Set some useful constants
  18 dampit          = 3;
  19 dampcoeff       = 2;
  20
  21 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  22 %% convert input vectors to columns
  23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  24
  25 if Ucolumns(D)>Urows(D)
  26         D=D';
  27 end
  28
  29
  30 if Ucolumns(nin)>Urows(nin)
  31         nin=nin';
  32 end
  33
  34 if Ucolumns(pin)>Urows(pin)
  35         pin=pin';
  36 end
  37
  38 if Ucolumns(Vin)>Urows(Vin)
  39         Vin=Vin';
  40 end
  41
  42 if Ucolumns(Fnin)>Urows(Fnin)
  43         Fnin=Fnin';
  44 end
  45
  46 if Ucolumns(Fpin)>Urows(Fpin)
  47         Fpin=Fpin';
  48 end
  49 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  50 %% setup FEM data structures
  51 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  52
  53 nodes=mesh.p;
  54 elements=mesh.t;
  55 Nnodes = length(nodes);
  56 Nelements = length(elements);
  57
  58 Dedges    =[];
  59
  60 for ii = 1:length(Dsides)
  61         Dedges=[Dedges,find(mesh.e(5,:)==Dsides(ii))];
  62 end
  63
  64 % Set list of nodes with Dirichelet BCs
  65 Dnodes = mesh.e(1:2,Dedges);
  66 Dnodes = [Dnodes(1,:) Dnodes(2,:)];
  67 Dnodes = unique(Dnodes);
  68
  69 % Set values of Dirichelet BCs
  70 Bc     = zeros(length(Dnodes),1);
  71 % Set list of nodes without Dirichelet BCs
  72 Varnodes = setdiff([1:Nnodes],Dnodes);
  73
  74
  75 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  76 %%
  77 %%              initialization:
  78 %%              we're going to solve
  79 %%              $$ - \lambda^2 (\delta V)'' +  (\frac{\partial n}{\partial V} - \frac{\partial p}{\partial V})= -R $$
  80 %%
  81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  82
  83 %% set $$ n_1 = nin $$ and $$ V = Vin $$
  84 V = Vin;
  85 Fn = Fnin;
  86 Fp = Fpin;
  87 G  = Gin;
  88 Gp = Gpin;
  89 n = exp(V(Sinodes)+G-Fn);
  90 p = exp(-V(Sinodes)-Gp+Fp);
  91 n(SiDnodes) = nin(SiDnodes);
  92 p(SiDnodes) = pin(SiDnodes);
  93
  94
  95 %%%
  96 %%% Compute LHS matrices
  97 %%%
  98
  99 %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 100 if (isempty(DDGOXNLPOISSON_LAP))
 101         coeff = l2ox * ones(Nelements,1);
 102         coeff(Sielements)=l2;
 103         DDGOXNLPOISSON_LAP = Ucomplap (mesh,coeff);
 104 end
 105
 106 %% compute $$ Mv = ( n + p)  $$
 107 %% and the (lumped) mass matrix M
 108 if (isempty(DDGOXNLPOISSON_MASS))
 109     Cvect = zeros(Nelements,1);
 110     Cvect(Sielements)=1;
 111         DDGOXNLPOISSON_MASS = Ucompmass2 (mesh,ones(Nnodes,1),Cvect);
 112 end
 113 freecarr=zeros(Nnodes,1);
 114 freecarr(Sinodes)=(n + p);
 115 Mv      =  freecarr;
 116 MV(SiDnodes) = 0;
 117 M       =  DDGOXNLPOISSON_MASS*sparse(diag(Mv));
 118
 119 %%%
 120 %%% Compute RHS vector (-residual)
 121 %%%
 122
 123 %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
 124 if (isempty(DDGOXNLPOISSON_RHS))
 125         DDGOXNLPOISSON_RHS = Ucompconst (mesh,ones(Nnodes,1),ones(Nelements,1));
 126 end
 127 totcharge = zeros(Nnodes,1);
 128 totcharge(Sinodes)=(n - p - D);
 129 Tv0   = totcharge;
 130 T0    = Tv0 .* DDGOXNLPOISSON_RHS;
 131
 132 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 133 A               = DDGOXNLPOISSON_LAP + M;
 134 R               = DDGOXNLPOISSON_LAP * V  + T0;
 135
 136 %% Apply boundary conditions
 137 A (Dnodes,:) = [];
 138 A (:,Dnodes) = [];
 139 R(Dnodes)  = [];
 140
 141 %% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test
 142 normr(1)        =  norm(R,inf);
 143 relresnorm      = 1;
 144 reldVnorm       = 1;
 145 normrnew        = normr(1);
 146 dV              = V*0;
 147
 148 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 149 %%
 150 %% START OF THE NEWTON CYCLE
 151 %%
 152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 153 for newtit=1:maxit
 154         if (verbose>0)
 155                 fprintf(1,'\n***\nNewton iteration: %d, reldVnorm = %e\n***\n',newtit,reldVnorm);
 156
 157         end
 158
 159         dV(Varnodes) =(A)\(-R);
 160         dV(Dnodes)=0;
 161
 162
 163         %%%%%%%%%%%%%%%%%%
 164         %% Start of th damping procedure
 165         %%%%%%%%%%%%%%%%%%
 166         tk = 1;
 167         for dit = 1:dampit
 168                 if (verbose>0)
 169                         fprintf(1,'\ndamping iteration: %d, residual norm = %e\n',dit,normrnew);
 170                 end
 171                 Vnew   = V + tk * dV;
 172
 173                 n = exp(Vnew(Sinodes)+G-Fn);
 174                 p = exp(-Vnew(Sinodes)-Gp+Fp);
 175                 n(SiDnodes) = nin(SiDnodes);
 176                 p(SiDnodes) = pin(SiDnodes);
 177
 178
 179                 %%%
 180                 %%% Compute LHS matrices
 181                 %%%
 182
 183                 %% let's compute  FEM approximation of $$ L = -  \frac{d^2}{x^2} $$
 184                 %L      = Ucomplap (mesh,ones(Nelements,1));
 185
 186                 %% compute $$ Mv =  ( n + p)  $$
 187                 %% and the (lumped) mass matrix M
 188                 freecarr=zeros(Nnodes,1);
 189                 freecarr(Sinodes)=(n + p);
 190                 Mv   =  freecarr;
 191                 M       =  DDGOXNLPOISSON_MASS*sparse(diag(Mv));%M     = Ucompmass (mesh,Mv);
 192
 193                 %%%
 194                 %%% Compute RHS vector (-residual)
 195                 %%%
 196
 197                 %% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$
 198                 totcharge( Sinodes)=(n - p - D);
 199                 Tv0   = totcharge;
 200                 T0    = Tv0 .* DDGOXNLPOISSON_RHS;%T0    = Ucompconst (mesh,Tv0,ones(Nelements,1));
 201
 202                 %% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step
 203                 A               = DDGOXNLPOISSON_LAP + M;
 204                 R               = DDGOXNLPOISSON_LAP * Vnew  + T0;
 205
 206                 %% Apply boundary conditions
 207                 A (Dnodes,:) = [];
 208                 A (:,Dnodes) = [];
 209                 R(Dnodes)  = [];
 210
 211                 %% compute $$ | R_{k+1} | $$ for the convergence test
 212                 normrnew= norm(R,inf);
 213
 214                 % check if more damping is needed
 215                 if (normrnew > normr(newtit))
 216                         tk = tk/dampcoeff;
 217                 else
 218                         if (verbose>0)
 219                           fprintf(1,'\nexiting damping cycle because residual norm = %e \n-----------\n',normrnew);
 220                         end
 221                         break
 222                 end
 223         end
 224
 225         V                           = Vnew;
 226         normr(newtit+1)         = normrnew;
 227         dVnorm              = norm(tk*dV,inf);
 228         pause(.1);
 229         % check if convergence has been reached
 230         reldVnorm           = dVnorm / norm(V,inf);
 231         if (reldVnorm <= toll)
 232           if(verbose>0)
 233                         fprintf(1,'\nexiting newton cycle because reldVnorm= %e \n',reldVnorm);
 234           end
 235           break
 236         end
 237
 238 end
 239
 240 res = normr;
 241 niter = newtit;
 242
 243
 244