1function [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 15global DDGOXNLPOISSON_LAP DDGOXNLPOISSON_MASS DDGOXNLPOISSON_RHS 16 17%% Set some useful constants 18dampit = 3; 19dampcoeff = 2; 20 21%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 22%% convert input vectors to columns 23%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 24 25if Ucolumns(D)>Urows(D) 26 D=D'; 27end 28 29 30if Ucolumns(nin)>Urows(nin) 31 nin=nin'; 32end 33 34if Ucolumns(pin)>Urows(pin) 35 pin=pin'; 36end 37 38if Ucolumns(Vin)>Urows(Vin) 39 Vin=Vin'; 40end 41 42if Ucolumns(Fnin)>Urows(Fnin) 43 Fnin=Fnin'; 44end 45 46if Ucolumns(Fpin)>Urows(Fpin) 47 Fpin=Fpin'; 48end 49%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 50%% setup FEM data structures 51%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 52 53nodes=mesh.p; 54elements=mesh.t; 55Nnodes = length(nodes); 56Nelements = length(elements); 57 58Dedges =[]; 59 60for ii = 1:length(Dsides) 61 Dedges=[Dedges,find(mesh.e(5,:)==Dsides(ii))]; 62end 63 64% Set list of nodes with Dirichelet BCs 65Dnodes = mesh.e(1:2,Dedges); 66Dnodes = [Dnodes(1,:) Dnodes(2,:)]; 67Dnodes = unique(Dnodes); 68 69% Set values of Dirichelet BCs 70Bc = zeros(length(Dnodes),1); 71% Set list of nodes without Dirichelet BCs 72Varnodes = 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 $$ 84V = Vin; 85Fn = Fnin; 86Fp = Fpin; 87G = Gin; 88Gp = Gpin; 89n = exp(V(Sinodes)+G-Fn); 90p = exp(-V(Sinodes)-Gp+Fp); 91n(SiDnodes) = nin(SiDnodes); 92p(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} $$ 100if (isempty(DDGOXNLPOISSON_LAP)) 101 coeff = l2ox * ones(Nelements,1); 102 coeff(Sielements)=l2; 103 DDGOXNLPOISSON_LAP = Ucomplap (mesh,coeff); 104end 105 106%% compute $$ Mv = ( n + p) $$ 107%% and the (lumped) mass matrix M 108if (isempty(DDGOXNLPOISSON_MASS)) 109 Cvect = zeros(Nelements,1); 110 Cvect(Sielements)=1; 111 DDGOXNLPOISSON_MASS = Ucompmass2 (mesh,ones(Nnodes,1),Cvect); 112end 113freecarr=zeros(Nnodes,1); 114freecarr(Sinodes)=(n + p); 115Mv = freecarr; 116MV(SiDnodes) = 0; 117M = DDGOXNLPOISSON_MASS*sparse(diag(Mv)); 118 119%%% 120%%% Compute RHS vector (-residual) 121%%% 122 123%% now compute $$ T0 = \frac{q}{\epsilon} (n - p - D) $$ 124if (isempty(DDGOXNLPOISSON_RHS)) 125 DDGOXNLPOISSON_RHS = Ucompconst (mesh,ones(Nnodes,1),ones(Nelements,1)); 126end 127totcharge = zeros(Nnodes,1); 128totcharge(Sinodes)=(n - p - D); 129Tv0 = totcharge; 130T0 = Tv0 .* DDGOXNLPOISSON_RHS; 131 132%% now we're ready to build LHS matrix and RHS of the linear system for 1st Newton step 133A = DDGOXNLPOISSON_LAP + M; 134R = DDGOXNLPOISSON_LAP * V + T0; 135 136%% Apply boundary conditions 137A (Dnodes,:) = []; 138A (:,Dnodes) = []; 139R(Dnodes) = []; 140 141%% we need $$ \norm{R_1} $$ and $$ \norm{R_k} $$ for the convergence test 142normr(1) = norm(R,inf); 143relresnorm = 1; 144reldVnorm = 1; 145normrnew = normr(1); 146dV = V*0; 147 148%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 149%% 150%% START OF THE NEWTON CYCLE 151%% 152%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 153for 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 238end 239 240res = normr; 241niter = newtit; 242 243 244 245