--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/secs1d/inst/demos/diode_refine_mesh.m	Thu Mar 14 16:21:08 2013 +0000
@@ -0,0 +1,85 @@
+% physical constants and parameters
+constants = secs1d_physical_constants_fun ();
+material  = secs1d_silicon_material_properties_fun (constants);
+
+% geometry
+Nelements = 100;
+L  = 50e-6;          % [m] 
+xm = L/2;
+device.W = 150e-6 * 50e-6;
+
+device.x  = linspace (0, L, Nelements+1)';
+device.sinodes = [1:length(device.x)];
+
+converged = false;
+iters = 1;
+maxiters = 10;
+while (! converged)
+
+  %% doping profile [m^{-3}]
+  device.Na = 1e23 * exp (- (device.x / 2e-6) .^ 2);
+  device.Nd = 1e25 * exp (- ((device.x-L) / 2.4e-6) .^ 2) + 1e19;
+
+  %% avoid zero doping
+  device.D  = device.Nd - device.Na;  
+
+  Fn = Fp = zeros (size (device.x));
+
+  %% carrier lifetime
+  device.tp = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'p');
+  device.tn = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'n');
+
+  %% 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 = 100;
+  algorithm.maxnpincr = constants.Vth;
+
+  secs1d_logplot (device.x, device.D);
+  drawnow
+
+  %% initial guess via stationary simulation
+  [Vin, nin, pin, res, niter] = ...
+      secs1d_nlpoisson_newton_noscale (device, material, 
+                                       constants, algorithm, 
+                                       V, n, p, Fn, Fp);
+  
+  els_to_refine = find ((abs (diff (log10 (nin))) > .02) 
+                        | (abs (diff (log10 (pin))) > .02) 
+                        | (abs (diff (log10 (abs (device.D)))) > .05));
+
+  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)];
+
+endwhile
+
+save diode_mesh x L xm
\ No newline at end of file

--- a/extra/secs1d/inst/demos/diode_tran_noscale.m	Thu Mar 14 13:04:01 2013 +0000
+++ b/extra/secs1d/inst/demos/diode_tran_noscale.m	Thu Mar 14 16:21:08 2013 +0000
@@ -1,14 +1,16 @@
+
 % physical constants and parameters
 constants = secs1d_physical_constants_fun ();
 material  = secs1d_silicon_material_properties_fun (constants);
 
 % geometry
-L  = 50e-6;          % [m] 
-xm = L/2;
+diode_refine_mesh
+load diode_mesh L xm x
+# Nelements = 1000;
+# L  = 50e-6;          % [m] 
+# xm = L/2;
 device.W = 150e-6 * 50e-6;
-
-Nelements = 1000;
-device.x  = linspace (0, L, Nelements+1)';
+# device.x  = linspace (0, L, Nelements+1)';
 device.sinodes = [1:length(device.x)];
 
 % doping profile [m^{-3}]
@@ -20,7 +22,7 @@
 
 % time span for simulation
 tmin = 0;
-tmax = 769;
+tmax = 3/500;
 tspan = [tmin, tmax];
 
 Fn = Fp = zeros (size (device.x));
@@ -48,22 +50,22 @@
 
 function fn = vbcs_1 (t);
   fn = [0; 0];
-  fn(2) = t;
+  fn(2) = 1 * sin (500 * pi * t);
 endfunction
 
 function fp = vbcs_2 (t);
   fp = [0; 0];
-  fp(2) = t;
+  fp(2) = 1 * sin (500 * pi * t);
 endfunction
 
 vbcs = {@vbcs_1, @vbcs_2};
 
 % tolerances for convergence checks
-algorithm.toll  = 1e-10;
-algorithm.maxit = 10;
+algorithm.toll  = 1e-6;
+algorithm.maxit = 300;
 algorithm.ptoll = 1e-12;
 algorithm.pmaxit = 1000;
-algorithm.maxnpincr = constants.Vth;
+algorithm.maxnpincr = constants.Vth / 10;
 
 %% initial guess via stationary simulation
 [nin, pin, Vin, Fnin, Fpin, Jn, Jp, it, res] = ...

--- a/extra/secs1d/inst/demos/thyristor_refine_mesh.m	Thu Mar 14 13:04:01 2013 +0000
+++ b/extra/secs1d/inst/demos/thyristor_refine_mesh.m	Thu Mar 14 16:21:08 2013 +0000
@@ -7,10 +7,10 @@
 material  = secs1d_silicon_material_properties_fun (constants);
 
 % geometry
-L  = 400e-6;          % [m] 
+L  = 80e-6;          % [m] 
 xm = L/2;
 
-Nelements      = 1000;
+Nelements      = 10;
 device.x       = linspace (0, L, Nelements+1)';
 device.sinodes = [1:length(device.x)];
 
@@ -54,11 +54,7 @@
   algorithm.pmaxit = 100;
   algorithm.maxnpincr = constants.Vth;
 
-  logplot = @(x) asinh (x/2) / log(10);
-  close all; plot (device.x, logplot (device.D)); 
-  set (gca, "yticklabel",
-       arrayfun (@(x) sprintf ("%s10^{%g}", ifelse (x >=0, "+", "-"), abs(x)),
-                 get (gca, "ytick"), "uniformoutput", false));
+  secs1d_logplot (device.x, device.D);
   drawnow
 
   %% initial guess via stationary simulation
@@ -67,7 +63,10 @@
                                        constants, algorithm, 
                                        V, n, p, Fn, Fp);
   
-  els_to_refine = find ((abs (diff (log10 (nin))) > .1) | (abs (diff (log10 (pin))) > .1) | (abs (diff (log10 (abs (device.D)))) > .1));
+  els_to_refine = find ((abs (diff (log10 (nin))) > .05) 
+                        | (abs (diff (log10 (pin))) > .05) 
+                        | (abs (diff (log10 (abs (device.D)))) > .1));
+
   if isempty (els_to_refine)
     converged = true;
   endif
@@ -80,4 +79,6 @@
   newx = (x(els_to_refine) + x(els_to_refine+1)) / 2;
   device.x = sort (vertcat (x, newx));
   device.sinodes = [1:length(device.x)];
-endwhile
\ No newline at end of file
+endwhile
+
+save thyristor_mesh x L
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/secs1d/inst/secs1d_logplot.m	Thu Mar 14 16:21:08 2013 +0000
@@ -0,0 +1,29 @@
+%% Copyright (C) 2013  Carlo de Falco
+%%
+%% This file is part of 
+%% SECS1D - A 1-D Drift--Diffusion Semiconductor Device Simulator
+%%
+%% SECS1D 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 3 of the License, or
+%% (at your option) any later version.
+%%
+%% SECS1D 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 SECS1D; If not, see <http://www.gnu.org/licenses/>.
+
+function secs1d_logplot (x, y)
+
+  plot (x, __logplot__ (y)); 
+  set (gca, "yticklabel",
+       arrayfun (@(x) sprintf ("%s10^{%g}", ifelse (x >=0, "+", "-"), abs(x)),
+                 get (gca, "ytick"), "uniformoutput", false));
+endfunction
+
+function logplot = __logplot__ (x)
+  logplot = asinh (x/2) / log(10);
+endfunction
\ No newline at end of file

--- a/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m	Thu Mar 14 13:04:01 2013 +0000
+++ b/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m	Thu Mar 14 16:21:08 2013 +0000
@@ -61,6 +61,8 @@
       secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm,
                                          Vin, nin, pin, Fnin, Fpin, tspan, vbcs)
 
+  rejected = 0;
+
   V  = Vin;
   n  = nin;
   p  = pin;
@@ -120,12 +122,13 @@
 
       V0 = Fn0 + constants.Vth * log (n0 ./ device.ni);
       
+      nrfnp = 4 * algorithm.maxnpincr; %% initialize so it is not undefined in case of exception
       [V(:,tstep), n(:,tstep), p(:,tstep), Fn(:,tstep), ...
-       Fp(:,tstep), Jn(:,tstep), Jp(:,tstep), res_, it] =...
+       Fp(:,tstep), Jn(:,tstep), Jp(:,tstep), res_, it, nrfnp] =...
           __one_step__ (device, material, constants, algorithm, 
                         V0, n0, p0, Fn0, Fp0, Jn(:,tstep), Jp(:,tstep),
-                        n(:, tstep - 1), p(:, tstep - 1), dt);      
-        
+                        n(:, tstep - 1), p(:, tstep - 1), dt);  
+      assert (nrfnp <= algorithm.maxnpincr)
       figure (1)
       plot (tout, (mean (Jn(:, 1:tstep) + Jp(:, 1:tstep), 1)))
       drawnow
@@ -134,29 +137,34 @@
       %if (tstep > 2)
       %  dt *= (algorithm.maxnpincr - res(tstep, 1)) / (res(tstep, 1) - res(tstep-1, 1));
       %endif
-      dt *= 2;
+      printf ("dt incremented by %g\n", .5 * sqrt (algorithm.maxnpincr / nrfnp))
+      dt *= .5 * sqrt (algorithm.maxnpincr / nrfnp);
       numit(tstep) = it;      
       
     catch
       
       warning (lasterr)
-      dt /= 2;
+      warning ("algorithm.maxnpincr = %17.17g, nrfnp = %17.17g, dt reduced by %17.17g", 
+               algorithm.maxnpincr, nrfnp, .5 * sqrt (algorithm.maxnpincr / nrfnp))
+      dt *=  .5 * sqrt (algorithm.maxnpincr / nrfnp);
       t = tout(--tstep);
       if (dt < 100*eps)
         warning ("secs1d_tran_dd_gummel_map_noscale: minimum time step reached.")
         return
       endif
+      rejected ++ ;
 
     end_try_catch
     
   endwhile  
-  
+  rejected
 endfunction
 
 function [Jn, Jp] = compute_currents (device, material, constants, algorithm,
                                       mobilityn, mobilityp, V, n, p, Fn, Fp)
-  dV = diff (V);
+  dV  = diff (V);
   dx = diff (device.x);
+
   [Bp, Bm] = bimu_bernoulli (dV ./ constants.Vth);
 
   Jn =  constants.q * constants.Vth * mobilityn .* ...
@@ -212,7 +220,7 @@
 endfunction
 
 
-function [V, n, p, Fn, Fp, Jn, Jp, res, it] = __one_step__ ...
+function [V, n, p, Fn, Fp, Jn, Jp, res, it, nrfnp] = __one_step__ ...
       (device, material, constants, algorithm, V0, n0, p0, Fn0,
        Fp0, Jn0, Jp0, n_old, p_old, dt); 
 
@@ -229,14 +237,14 @@
   Nelements = Nnodes -1;
   
   for it = 1 : algorithm.maxit
-    
+
     [V, n, p] = secs1d_nlpoisson_newton_noscale ...
         (device, material, constants, algorithm, V, n, p, Fn, Fp); 
       
     dV = diff (V);
     E  = - dV ./ dx;
 
-    [Bp, Bm] = bimu_bernoulli (dV);
+    %% [Bp, Bm] = bimu_bernoulli (dV);
     
     [Rn, Rp, Gn, Gp, II] = generation_recombination_model ...
         (device, material, constants, algorithm, E, Jn, Jp, 
@@ -271,16 +279,16 @@
          V, n, p, Fn, Fp);
     
     nrfn   = norm (Fn - Fn0, inf);
-    assert (norm (Fn - Fnref, inf) <= algorithm.maxnpincr);
-    
     nrfp   = norm (Fp - Fp0, inf);
-    assert (norm (Fp - Fpref, inf) <= algorithm.maxnpincr);
+
+    nrfnp = max (norm (Fn - Fnref, inf), norm (Fp - Fpref, inf));
     
     nrv    = norm (V  - V0,  inf);
     res(it) = max  ([nrfn; nrfp; nrv]);
     
-    if (res(it) < algorithm.toll)
-      break
+    if (res(it) < algorithm.toll
+        || (nrfnp > algorithm.maxnpincr))
+      break;
     endif
     
     V0  =  V;