Mercurial > forge

changeset 11751:5bba72d50632 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
more derivatives
author cdf
date Thu, 06 Jun 2013 06:58:13 +0000
parents 956d6a9e060e
children 5fa40e1d044c
files extra/secs1d/inst/demos/nplus_n_nnplus.m extra/secs1d/inst/demos/nplus_n_nnplus_mesh.m extra/secs1d/inst/demos/simple_diode_tran_noscale_res.m extra/secs1d/inst/secs1d_newton_res.m
diffstat 4 files changed, 216 insertions(+), 14 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/secs1d/inst/demos/nplus_n_nnplus.m	Thu Jun 06 06:58:13 2013 +0000
@@ -0,0 +1,100 @@
+pkg load secs1d
+clear all
+close all
+
+                                % physical constants and parameters
+constants = secs1d_physical_constants_fun ();
+material  = secs1d_silicon_material_properties_fun (constants);
+
+                                % geometry
+
+loaded_mesh    = load ("nplus_n_nnplus_mesh");
+device.x       = loaded_mesh.x;
+L              = loaded_mesh.L;          % [m] 
+device.D       = loaded_mesh.D;          % [m^-3] 
+device.Na      = loaded_mesh.Na;         % [m^-3] 
+device.Nd      = loaded_mesh.Nd;         % [m^-3] 
+device.W       = 1e-12;                  % [m^2]
+tspan          = [0 60];
+
+device.sinodes = [1:length(device.x)];
+
+Fn = Fp = zeros (size (device.x));
+
+%% bandgap narrowing correction
+device.ni = (material.ni) * exp (secs1d_bandgap_narrowing_model
+                                 (device.Na, device.Nd) / constants.Vth); 
+
+%% carrier lifetime
+device.tp = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'p');
+device.tn = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'n');
+
+%% initial guess for n, p, V, phin, phip
+p = ((abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+     (device.D <= 0)) / 2 + 2 * device.ni.^2 ./ ...
+                            (abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+                            (device.D > 0);
+
+n = ((abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+     (device.D > 0)) / 2 + 2 * device.ni.^2 ./ ...
+                           (abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+                           (device.D <= 0);
+
+V = Fn + constants.Vth * log (n ./ device.ni);
+
+
+% tolerances for convergence checks
+algorithm.toll       = 1e-6;
+algorithm.ltol       = 1e-10;
+algorithm.maxit      = 100;
+algorithm.lmaxit     = 100;
+algorithm.ptoll      = 1e-12;
+algorithm.pmaxit     = 1000;
+algorithm.colscaling = [10 1e21 1e21 1];
+algorithm.rowscaling = [1  1e-7 1e-7 1];
+algorithm.maxnpincr  = 1e-3;
+
+secs1d_logplot (device.x, device.D);
+drawnow
+
+%% initial guess via stationary simulation
+[nin, pin, Vin, Fnin, Fpin, Jn, Jp, it, res] = secs1d_dd_gummel_map_noscale ...
+    (device, material, constants, algorithm, V, n, p, Fn, Fp);  
+
+
+function [g, j, r] = vbcs (t, dt);
+  g =  [1; 1];
+  j = -[t  0];
+  r =  [0  0];
+endfunction
+
+%% (pseudo)transient simulation
+[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = ...
+secs1d_newton_res (device, material, constants, algorithm,
+                   Vin(:,end), nin(:,end), pin(:,end), 
+                   tspan, @vbcs);
+
+dV   = diff (V, [], 1);
+dx   = diff (device.x);
+E    = -dV ./ dx;
+
+%% band structure
+Efn  = -Fn(:, end);
+Efp  = -Fp(:, end);
+Ec   =  constants.Vth * log (material.Nc ./ n(:, end)) + Efn;
+Ev   = -constants.Vth * log (material.Nv ./ p(:, end)) + Efp;
+
+figure (1)
+plot (device.x, Efn, device.x, Efp, device.x, Ec, device.x, Ev)
+legend ('Efn', 'Efp', 'Ec', 'Ev')
+axis tight
+drawnow
+
+vvector  = Fn(end, :);
+ivector  = (Jn(end, :) + Jp(end, :));
+ivectorn = (Jn(1, :)   + Jp(1, :));
+
+figure (2) 
+plot (vvector, ivector, vvector, ivectorn)
+legend('J_L','J_0')
+drawnow
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/secs1d/inst/demos/nplus_n_nnplus_mesh.m	Thu Jun 06 06:58:13 2013 +0000
@@ -0,0 +1,101 @@
+pkg load secs1d
+clear all
+close all
+
+% physical constants and parameters
+constants = secs1d_physical_constants_fun ();
+material  = secs1d_silicon_material_properties_fun (constants);
+
+% geometry
+L  = 600e-9;          % [m] 
+xm = L/2;
+
+Nelements      = 10;
+device.x       = linspace (0, L, Nelements+1)';
+device.sinodes = [1:length(device.x)];
+
+converged = false;
+iters = 1;
+maxiters = 100;
+while (! converged)
+
+  %% doping profile [m^{-3}]
+  device.Nd = 1e22 + 4.9e23 * (device.x < 100e-9) + 4.9e23 * (device.x > 500e-9);
+  device.Na = 2;
+
+  %% avoid zero doping
+  device.D  = device.Nd - device.Na;  
+
+
+  Fn = Fp = zeros (size (device.x));
+
+  %% bandgap narrowing correction
+  device.ni = (material.ni) * exp (secs1d_bandgap_narrowing_model
+                                   (device.Na, device.Nd) / constants.Vth); 
+
+  %% initial guess for n, p, V, phin, phip
+  p = ((abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+       (device.D <= 0)) / 2 + 2 * device.ni.^2 ./ ...
+      (abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+      (device.D > 0);
+
+  n = ((abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+       (device.D > 0)) / 2 + 2 * device.ni.^2 ./ ...
+      (abs(device.D) + sqrt (abs(device.D) .^ 2 + 4 * device.ni .^2)) .* ...
+      (device.D <= 0);
+
+  V = Fn + constants.Vth * log (n ./ device.ni);
+
+
+  %% tolerances for convergence checks
+  algorithm.toll  = 1e-6;
+  algorithm.maxit = 100;
+  algorithm.ptoll = 1e-12;
+  algorithm.pmaxit = 1000;
+  algorithm.maxnpincr = constants.Vth;
+
+  figure(1)
+  secs1d_logplot (device.x, device.D);
+  drawnow
+
+  assert (~ any (diff (device.x < 0)))
+
+  %% initial guess via stationary simulation
+  [Vin, nin, pin, res, niter] = ...
+      secs1d_nlpoisson_newton_noscale (device, material, 
+                                       constants, algorithm, 
+                                       V, n, p, Fn, Fp);
+  figure(2)
+  secs1d_logplot (device.x, nin);
+  drawnow
+
+  figure(3)
+  secs1d_logplot (device.x, pin);
+  drawnow
+
+
+  assert (~ any (isnan (Vin)))
+  assert (~ any (isnan (nin)))
+  assert (~ any (isnan (pin)))
+
+  els_to_refine = find ((abs (diff (log10 (nin))) > .4) 
+                        | (abs (diff (log10 (pin))) > .4));
+
+  if isempty (els_to_refine)
+    converged = true;
+  endif
+
+  if (converged || (++iters > maxiters))
+    break
+  endif
+
+  x = device.x;
+  newx = (x(els_to_refine) + x(els_to_refine+1)) / 2;
+  device.x = sort (vertcat (x, newx));
+  device.sinodes = [1:length(device.x)];
+  D = device.D;
+  Na = device.Na;
+  Nd = device.Nd;
+endwhile
+
+save nplus_n_nnplus_mesh x L D Na Nd
--- a/extra/secs1d/inst/demos/simple_diode_tran_noscale_res.m	Wed Jun 05 16:46:45 2013 +0000
+++ b/extra/secs1d/inst/demos/simple_diode_tran_noscale_res.m	Thu Jun 06 06:58:13 2013 +0000
@@ -19,7 +19,7 @@
 
 % time span for simulation
 tmin  = 0;
-tmax  = 1;
+tmax  = 10;
 tspan = [tmin, tmax];
 
 Fn = Fp = zeros (size (device.x));
@@ -48,11 +48,11 @@
 function [g, j, r] = vbcs (t, dt);
   g = [1; 1];
   j = -[0 -t];
-  r = [1e3 0];
+  r = [1e4 0e4];
 endfunction
 
 % tolerances for convergence checks
-algorithm.toll       = 1e-6;
+algorithm.toll       = 1e-3;
 algorithm.ltol       = 1e-10;
 algorithm.maxit      = 100;
 algorithm.lmaxit     = 100;
@@ -60,7 +60,7 @@
 algorithm.pmaxit     = 1000;
 algorithm.colscaling = [10 1e21 1e21 .1];
 algorithm.rowscaling = [1  1e-7 1e-7 .1];
-algorithm.maxnpincr  = 1e-1;
+algorithm.maxnpincr  = 1e-2;
 
 %% initial guess via stationary simulation
 [nin, pin, Vin, Fnin, Fpin, Jn, Jp, it, res] = secs1d_dd_gummel_map_noscale ...
--- a/extra/secs1d/inst/secs1d_newton_res.m	Wed Jun 05 16:46:45 2013 +0000
+++ b/extra/secs1d/inst/secs1d_newton_res.m	Thu Jun 06 06:58:13 2013 +0000
@@ -189,12 +189,12 @@
 
       Itot(:, tstep)  = Jn([1 end], tstep) + Jp([1 end], tstep);
 
-      Itot(1, tstep) += constants.e0 * material.esir * ...
+      Itot(1, tstep) += (1/dt) * constants.e0 * material.esir * ...
           ((V(2, tstep) - V(1, tstep)) -
            (V(2, tstep-1) - V(1, tstep-1))) / ...
           (device.x(2) - device.x(1));
 
-      Itot(2, tstep) += constants.e0 * material.esir * ...
+      Itot(2, tstep) += (1/dt) * constants.e0 * material.esir * ...
           ((V(end, tstep) - V(end-1, tstep)) -
            (V(end, tstep-1) - V(end-1, tstep-1))) / ...
           (device.x(end) - device.x(end-1));
@@ -248,11 +248,11 @@
   I2  = - constants.q * (A22(end,:) * n -  A33(end,:) * p);
   
   if (columns (V) >= 2)
-    I1 += constants.e0 * material.esir * ...
+    I1 += (1/dt) * constants.e0 * material.esir * ...
           ((V(2, end) - V(1, end)) -
            (V(2, end-1) - V(1, end-1))) / ...
           (device.x(2) - device.x(1));
-    I2 += constants.e0 * material.esir * ...
+    I2 += (1/dt) * constants.e0 * material.esir * ...
           ((V(end, end) - V(end-1, end)) -
            (V(end, end-1) - V(end-1, end-1))) / ...
           (device.x(end) - device.x(end-1));
@@ -316,13 +316,14 @@
       bim1a_reaction (device.x, 1, Rp + 1/dt);
   jac{3,4} = sparse (Nnodes, 2);
 
+
+  A41 = bim1a_laplacian (device.x, 
+                         (-mobilityn .* nm + mobilityp .* pm), 1);
   
-  jac{4,1} = sparse(2, Nnodes);
-  jac{4,2}(1,:) = - ri(1) * device.W * constants.q * A22(1, 2:end-1);
-  jac{4,2}(2,:) = - ri(2) * device.W * constants.q * A22(end, 2:end-1);
-  jac{4,3}(1,:) =   ri(1) * device.W * constants.q * A33(1, 2:end-1);
-  jac{4,3}(2,:) =   ri(2) * device.W * constants.q * A33(end, 2:end-1);
-  jac{4,4} = spdiags (gi(:), 0, 2, 2);
+  jac{4,1} =   diag (ri(:)) * device.W * constants.q * A41 ([1, end], :) ;
+  jac{4,2} = - diag (ri(:)) * device.W * constants.q * A22([1 end], 2:end-1);
+  jac{4,3} =   diag (ri(:)) * device.W * constants.q * A33([1 end], 2:end-1);
+  jac{4,4} =   spdiags (gi(:), 0, 2, 2);
 
 endfunction