Mercurial > forge

changeset 11570:35c57f32488a octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix generation recombination and scaling in stationary Newton
author cdf
date Sun, 24 Mar 2013 14:53:22 +0000
parents 332c5ee57ec0
children 10f755d08a29
files extra/secs1d/inst/demos/diode_reverse_noscale.m extra/secs1d/inst/secs1d_dd_newton_noscale.m extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m
diffstat 3 files changed, 44 insertions(+), 24 deletions(-) [+]
line wrap: on
line diff
--- a/extra/secs1d/inst/demos/diode_reverse_noscale.m	Fri Mar 22 09:12:39 2013 +0000
+++ b/extra/secs1d/inst/demos/diode_reverse_noscale.m	Sun Mar 24 14:53:22 2013 +0000
@@ -25,7 +25,7 @@
 device.tn = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'n');
 
 % initial guess for n, p, V, phin, phip
-V_p = -1;
+V_p =  1;
 V_n =  0;
 
 Fp = V_p * (device.x <= xm);
@@ -45,8 +45,8 @@
 V = Fn + constants.Vth * log (n ./ device.ni);
 
 % tolerances for convergence checks
-algorithm.toll  = 1e-3;
-algorithm.maxit = 1000;
+algorithm.toll  = 1e-10;
+algorithm.maxit = 100;
 algorithm.ptoll = 1e-12;
 algorithm.pmaxit = 1000;
 
@@ -56,6 +56,7 @@
                                     constants, algorithm,
                                     V, n, p, Fn, Fp);  
 
+algorithm.maxit = 1000;
 [n_newt, p_newt, V_newt, Fn_newt, Fp_newt, Jn, Jp, it, res] = secs1d_dd_newton_noscale ...  
     (device, material, constants, algorithm, V, n, p, Fn, Fp);
 
@@ -65,6 +66,9 @@
 Ec   = constants.Vth * log (material.Nc./n) + Efn;
 Ev   = -constants.Vth * log (material.Nv./p) + Efp;
 
-plot (device.x, Efn, device.x, Efp, device.x, Ec, device.x, Ev)
+plot (device.x, Efn, 
+      device.x, Efp, 
+      device.x, Ec, 
+      device.x, Ev)
 legend ('Efn', 'Efp', 'Ec', 'Ev')
 axis tight
--- a/extra/secs1d/inst/secs1d_dd_newton_noscale.m	Fri Mar 22 09:12:39 2013 +0000
+++ b/extra/secs1d/inst/secs1d_dd_newton_noscale.m	Sun Mar 24 14:53:22 2013 +0000
@@ -56,13 +56,11 @@
   dampcoeff  = 2;
   Nnodes     = numel (device.x);
 
-  V  = Vin;
-  n  = nin;
-  p  = pin;
+  V  = Vin; n = nin; p = pin;
 
   [res, jac] = residual_jacobian (device, material, constants, 
                                   algorithm, V, n, p);
-  keyboard
+
   normr(1)   = norm (res, inf);
   normrnew   = normr(1);
 
@@ -86,6 +84,7 @@
                                             algorithm, Vnew, nnew, pnew);
 
       normrnew = norm (resnew, inf);
+
       if (normrnew > normr(it))
         tk = tk / dampcoeff;
       else
@@ -99,7 +98,18 @@
     V = Vnew; n = nnew; p = pnew;
     normr(it+1) = normrnew;
 
-    if (normr(it+1) <= algorithm.toll)
+    resvec(it) = reldnorm    = + ...
+        norm (delta(1:Nnodes-2), inf) / norm (V, inf) + ...
+        norm (delta((Nnodes-2)+(1:Nnodes-2)), inf) / norm (n, inf) + ...
+        norm (delta(2*(Nnodes-2)+(1:Nnodes-2)), inf) / norm (p, inf);
+      
+    plotyy (
+            1:it, log10 (resvec), 
+            1:it+1, log10 (normr)
+            )
+    drawnow
+      
+    if (reldnorm <= algorithm.toll)
       break
     endif
 
@@ -115,6 +125,8 @@
   Fp = V + constants.Vth * log (p ./ device.ni);
   Fn = V - constants.Vth * log (n ./ device.ni);
 
+  res = resvec;
+
 endfunction
 
 
@@ -152,46 +164,50 @@
 function [res, jac] = residual_jacobian (device, material, constants, 
                                          algorithm, V, n, p)
 
- [mobilityn, mobilityp] = compute_mobilities ...
+  [mobilityn, mobilityp] = compute_mobilities ...
       (device, material, constants, algorithm, V, n, p);
-
+  
   [Rn, Rp, Gn, Gp, II] = generation_recombination_model ...
       (device, material, constants, algorithm, mobilityn, 
        mobilityp, V, n, p);
 
-  nm         = (n(2:end) + n(1:end-1))/2;
-  pm         = (p(2:end) + p(1:end-1))/2;
+  nm         = (n(2:end) + n(1:end-1)) / 2;
+  pm         = (p(2:end) + p(1:end-1)) / 2;
 
   Nnodes     = numel (device.x);
   Nelements  = Nnodes - 1;
 
-  epsilon    = material.esi * ones (Nelements, 1);
-
+  epsilon    = material.esi * ones (Nelements, 1) / constants.q;
+  
   A11 = bim1a_laplacian (device.x, epsilon, 1);
-  A12 = bim1a_reaction (device.x, constants.q, 1);
+  A12 = bim1a_reaction (device.x, 1, 1);
   A13 = - A12;
-  R1  = A11 * V + bim1a_rhs (device.x, constants.q, n - p - device.D);
+  R1  = A11 * V + bim1a_rhs (device.x, 1, n - p - device.D);
 
   A21 = - bim1a_laplacian (device.x, mobilityn .* nm, 1);
   A22 = bim1a_advection_diffusion ...
       (device.x, mobilityn * constants.Vth, 1, 1, V / constants.Vth) + ...
       bim1a_reaction (device.x, 1, Rn);
   A23 = bim1a_reaction (device.x, 1, Rp);
-  R2  = A22 * n + bim1a_rhs (device.x, 1, Gn - Rn .* n);
+  R2  = A22 * n + bim1a_rhs (device.x, 1, Rn .* n - Gn);
 
   A31 = bim1a_laplacian (device.x, mobilityp .* pm, 1);
   A32 = bim1a_reaction (device.x, 1, Rn);
   A33 = bim1a_advection_diffusion ...
       (device.x, mobilityp * constants.Vth , 1, 1, - V / constants.Vth) + ...
       bim1a_reaction (device.x, 1, Rp);
-  R3  = A33 * p + bim1a_rhs (device.x, 1, Gp - Rp .* p);
+  R3  = A33 * p + bim1a_rhs (device.x, 1, Rp .* p - Gp);
 
-  res = [R1(2:end-1); R2(2:end-1); R3(2:end-1)];
+  N1 = norm(R1(2:end-1), inf);
+  N2 = norm(R2(2:end-1), inf);
+  N3 = norm(R3(2:end-1), inf);
 
-  jac = [A11(2:end-1, 2:end-1), A12(2:end-1, 2:end-1), A13(2:end-1, 2:end-1);
-         A21(2:end-1, 2:end-1), A22(2:end-1, 2:end-1), A23(2:end-1, 2:end-1);
-         A31(2:end-1, 2:end-1), A32(2:end-1, 2:end-1), A33(2:end-1, 2:end-1)];
+  res = [R1(2:end-1)/N1; R2(2:end-1)/N2; R3(2:end-1)/N2];
 
+  jac = [[A11(2:end-1, 2:end-1), A12(2:end-1, 2:end-1), A13(2:end-1, 2:end-1)]/N1;
+         [A21(2:end-1, 2:end-1), A22(2:end-1, 2:end-1), A23(2:end-1, 2:end-1)]/N2;
+         [A31(2:end-1, 2:end-1), A32(2:end-1, 2:end-1), A33(2:end-1, 2:end-1)]/N3];
+  
 endfunction
 
 function [Jn, Jp] = compute_currents (device, material, constants, algorithm,
--- a/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m	Fri Mar 22 09:12:39 2013 +0000
+++ b/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m	Sun Mar 24 14:53:22 2013 +0000
@@ -157,7 +157,7 @@
     end_try_catch
     
   endwhile  
-  rejected
+  printf ("number of rejected time steps: %d\n", rejected)
 endfunction
 
 function [Jn, Jp] = compute_currents (device, material, constants, algorithm,