Mercurial > forge

changeset 11595:344a62cf0446 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix transient
author cdf
date Wed, 03 Apr 2013 16:37:41 +0000
parents 7495ebb646aa
children 493ef9d1c25d
files extra/secs1d/inst/demos/diode_tran_noscale.m extra/secs1d/inst/secs1d_newton.m
diffstat 2 files changed, 13 insertions(+), 7 deletions(-) [+]
line wrap: on
line diff
--- a/extra/secs1d/inst/demos/diode_tran_noscale.m	Wed Apr 03 10:31:34 2013 +0000
+++ b/extra/secs1d/inst/demos/diode_tran_noscale.m	Wed Apr 03 16:37:41 2013 +0000
@@ -22,7 +22,7 @@
 
 % time span for simulation
 tmin = 0;
-tmax = 3/5e3;
+tmax = 2/5e3;
 tspan = [tmin, tmax];
 
 Fn = Fp = zeros (size (device.x));
@@ -50,12 +50,12 @@
 
 function fn = vbcs_1 (t);
   fn = [0; 0];
-  fn(2) = - 30 * sin (5e3 * pi * t);
+  fn(2) = -1 * (t > 1e-4);
 endfunction
 
 function fp = vbcs_2 (t);
   fp = [0; 0];
-  fp(2) = - 30 * sin (5e3 * pi * t);
+  fp(2) = -1 * (t > 1e-4);
 endfunction
 
 vbcs = {@vbcs_1, @vbcs_2};
--- a/extra/secs1d/inst/secs1d_newton.m	Wed Apr 03 10:31:34 2013 +0000
+++ b/extra/secs1d/inst/secs1d_newton.m	Wed Apr 03 16:37:41 2013 +0000
@@ -20,8 +20,8 @@
 
       [V1, n1, p1] = deal (V2, n2, p2);
 
-      res = compute_residual (device, material, constants, algorithm, V2, n2, p2, n2, p2, dt);
-      jac = compute_jacobian (device, material, constants, algorithm, V2, n2, p2, n2, p2, dt);
+      res = compute_residual (device, material, constants, algorithm, V2, n2, p2, n(:, tstep-1), p(:, tstep-1), dt);
+      jac = compute_jacobian (device, material, constants, algorithm, V2, n2, p2, n(:, tstep-1), p(:, tstep-1), dt);
       dn = dp = dV = zeros(rows (n) - 2, 1);
 
       direct = true;
@@ -96,6 +96,12 @@
       incr1p = norm (log (p2./p1), inf) / (norm (log (p0), inf) + log (algorithm.colscaling(3)));
 
       resnlin(in) = incr1 = max ([incr1v, incr1n, incr1p]);
+      if (in > 3 && resnlin(in) > resnlin(in-3))
+        printf ("newton step is diverging\n")
+        tstep
+        reject = true;
+        break;
+      endif
 
       incr0v = norm (V2 - V0, inf) / (norm (V0, inf) + algorithm.colscaling(1));
       incr0n = norm (log (n2./n0), inf) / (norm (log (n0), inf) + log (algorithm.colscaling(2)));
@@ -115,7 +121,7 @@
       endif
 
       if (incr1 < algorithm.toll)
-        printf ("iteration %d time step %d  convergence reached incr = %g\n", in, tstep, incr1)
+        printf ("iteration %d, time step %d, model time %g: convergence reached incr = %g\n", in, tstep, t, incr1)
         break;
       endif
       
@@ -160,7 +166,7 @@
       plot (tout, Itot);
       drawnow
 
-      dt *= .5 * sqrt (algorithm.maxnpincr / incr0)
+      dt *= .75 * sqrt (algorithm.maxnpincr / incr0)
 
     endif