changeset 12117:d830ec7b1ca6 octave-forge

added circuit variables output to secs3d_coupled_circuit_newton.m
author dvd7587
date Fri, 25 Oct 2013 10:33:29 +0000
parents 88c10b8baba2
children ecb65e433faf
files extra/secs1d/inst/demos/coupled_3d_diode.m extra/secs1d/inst/demos/coupled_3d_resistor_and_circuit.m extra/secs1d/inst/secs3d_coupled_circuit_newton.m
diffstat 3 files changed, 50 insertions(+), 34 deletions(-) [+]
line wrap: on
line diff
--- a/extra/secs1d/inst/demos/coupled_3d_diode.m	Wed Oct 23 17:35:22 2013 +0000
+++ b/extra/secs1d/inst/demos/coupled_3d_diode.m	Fri Oct 25 10:33:29 2013 +0000
@@ -122,7 +122,7 @@
 
 save -binary -z datafile_rlc_circuit.octbin.gz  device material constants algorithm tspan Vin nin pin  % A B C x r 
 
-[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs3d_coupled_circuit_newton ...
+[V, n, p, F, Fn, Fp, Jn, Jp, Itot, tout] = secs3d_coupled_circuit_newton ...
                                            (device, material, constants, algorithm,
                                             Vin, nin, pin, tspan, @vbcs);
 
--- a/extra/secs1d/inst/demos/coupled_3d_resistor_and_circuit.m	Wed Oct 23 17:35:22 2013 +0000
+++ b/extra/secs1d/inst/demos/coupled_3d_resistor_and_circuit.m	Fri Oct 25 10:33:29 2013 +0000
@@ -1,3 +1,4 @@
+#!/usr/bin/octave
 %%
 clear all;
 close all;
@@ -119,10 +120,11 @@
 pause
 
 %% (pseudo)transient simulation
-
-save -binary -z datafile_rlc_circuit.octbin.gz  device material constants algorithm tspan Vin nin pin  % A B C x r 
-
-[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs3d_coupled_circuit_newton ...
+[~, ~, ~, ~, Fin, ~] = vbcs (0);
+device.msh = bim3c_mesh_properties (device.msh);
+save -binary -z datafile_rlc_circuit.octbin.gz  device material constants algorithm tspan Vin nin pin Fin  % A B C x r 
+return
+[V, n, p, F, Fn, Fp, Jn, Jp, Itot, tout] = secs3d_coupled_circuit_newton ...
                                            (device, material, constants, algorithm,
                                             Vin, nin, pin, tspan, @vbcs);
 
--- a/extra/secs1d/inst/secs3d_coupled_circuit_newton.m	Wed Oct 23 17:35:22 2013 +0000
+++ b/extra/secs1d/inst/secs3d_coupled_circuit_newton.m	Fri Oct 25 10:33:29 2013 +0000
@@ -1,7 +1,13 @@
-function [V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = ...
+function [V, n, p, F, Fn, Fp, Jn, Jp, Itot, tout] = ...
          secs3d_coupled_circuit_newton ...
            (device, material, constants, algorithm, 
             Vin, nin, pin, tspan, va)    
+
+  if (isfield (algorithm, 'plotsOff'));
+    plotsOff = algorithm.plotsOff;
+  else
+    plotsOff = 0;
+  endif
   
   tags = {'V', 'n', 'p', 'F'};
   rejected = 0;
@@ -21,9 +27,9 @@
     this_dnodes = bim3c_unknowns_on_faces (
                    device.msh, device.contacts(iii));
     if (! sum(this_dnodes));
-      error((["circuit pin #%d is not connected",
-              " to the device. Check mesh bound"
-              "ary labeling."])'(:)', iii)
+      error((['circuit pin #%d is not connected',
+              ' to the device. Check mesh bound'
+              'ary labeling.'])'(:)', iii)
     endif
     inodes = setdiff (inodes, this_dnodes);
     dnodes(iii, this_dnodes) = 1;
@@ -64,19 +70,23 @@
                                         dt, A, B, C, r, pins,
                                         indexing, dnodes, inodes); 
       
-      jac =  __secs3d_newton_jacobian__ (device, material, constants, 
-                                         algorithm, V2, n2, p2, F2, 
-                                         dt, A, B, C, r, pins,
-                                         indexing, dnodes, inodes);
+      jac = __secs3d_newton_jacobian__ (device, material, constants, 
+                                        algorithm, V2, n2, p2, F2, 
+                                        dt, A, B, C, r, pins,
+                                        indexing, dnodes, inodes);
 
       J = sparse(rows(jac{1, 1}), columns(jac{1, 1}));
 
+
       for iii = 1:4
         for jjj = 1:4
           J += jac{iii, jjj};
         end
       end
 
+      save -binary -z octave_jacobian_and_residual.octbin.gz J res
+      return
+      
       delta = - J \ res;
 
       dn = dp = dV = zeros(rows (n), 1);
@@ -102,7 +112,7 @@
 
       tk = min ([tkv, tkn, tkp]);
       if (tk <= 0)
-        error ("relaxation parameter too small, die!")
+        error ('relaxation parameter too small, die!')
       endif
       V2 += tk * dV;
       n2 += tk * dn;
@@ -127,8 +137,8 @@
 
       [incr0, whichone] = max ([incr0v, incr0n, incr0p, incr0F]);
       if (incr0 > algorithm.maxnpincr)
-        printf ("at time step %d, fixed point iteration %d, the ", tstep, in);
-        printf ("increment in %s has grown too large\n", tags{whichone});
+        printf ('at time step %d, fixed point iteration %d, the ', tstep, in);
+        printf ('increment in %s has grown too large\n', tags{whichone});
         reject = true;
         break;
       endif
@@ -148,22 +158,24 @@
       [incr1, whichone] = max ([incr1v, incr1n, incr1p, incr1F]);
       resnlin(in) = incr1;
       if (in > 3 && resnlin(in) > resnlin(in-3))
-        printf ("at time step %d, fixed point iteration %d, ", tstep, in);
-        printf ("the Newton algorithm is diverging: ");
-        printf ("the increment in %s is not decreasing\n", tags{whichone});
+        printf ('at time step %d, fixed point iteration %d, ', tstep, in);
+        printf ('the Newton algorithm is diverging: ');
+        printf ('the increment in %s is not decreasing\n', tags{whichone});
         reject = true;
         break;
       endif
 
-      figure (1)
-      semilogy (1:in, resnlin(1:in),'bo-');
-      xlim([1,15]);
-      ylim([5e-9,5e-2]);
-      drawnow        
+      if (! plotsOff) 
+        figure (1)
+        semilogy (1:in, resnlin(1:in),'bo-');
+        xlim([1,15]);
+        ylim([5e-9,5e-2]);
+        drawnow
+      endif         
 
       if (incr1 < algorithm.toll)
-        printf ("fixed point iteration %d, time step %d, ", in, tstep);
-        printf ("model time %g: convergence reached incr = %g ", t, incr1);
+        printf ('fixed point iteration %d, time step %d, ', in, tstep);
+        printf ('model time %g: convergence reached incr = %g ', t, incr1);
         break;
       endif
 
@@ -175,10 +187,10 @@
       t = tout (--tstep);
       dt /= 2;
 
-      printf ("reducing time step\n");
-      printf ("\ttime step #%d, ", tstep);
-      printf ("model time %g s\n", t);
-      printf ("\tnew dt %g s\n", dt);
+      printf ('reducing time step\n');
+      printf ('\ttime step #%d, ', tstep);
+      printf ('model time %g s\n', t);
+      printf ('\tnew dt %g s\n', dt);
 
     else
 
@@ -212,12 +224,14 @@
       endfor
       Itot(:, tstep) = (Jn_pins + Jp_pins + Jd_pins)(:);
       
-      figure (2)
-      plotyy (tout, Itot(2, :), tout, F(pins(2), :) - F(pins(1), :));
-      drawnow
+      if (! plotsOff) 
+        figure (2)
+        plotyy (tout, Itot(2, :), tout, F(pins(2), :) - F(pins(1), :));
+        drawnow
+      endif
     
       dt *= .8 * sqrt (algorithm.maxnpincr / incr0);
-      printf ("\nestimate for next time step size: dt = %g \n", dt);
+      printf ('\nestimate for next time step size: dt = %g \n', dt);
     endif
 
   endwhile %% time step