Mercurial > forge

changeset 11784:b9165047d2b4 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge diffs
author cdf
date Wed, 12 Jun 2013 14:32:26 +0000
parents cb452966068b
children 6da7dad71918
files extra/secs1d/inst/coupled_circuit_coeff.m extra/secs1d/inst/demos/coupled_circuit_coeff.m extra/secs1d/inst/demos/decomposematrixes.m extra/secs1d/inst/demos/simple_resistor_tran_noscale_res.m extra/secs1d/inst/secs1d_newton_res.m
diffstat 5 files changed, 85 insertions(+), 150 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/secs1d/inst/coupled_circuit_coeff.m	Wed Jun 12 14:32:26 2013 +0000
@@ -0,0 +1,60 @@
+## Copyright (C) 2013 davide
+## 
+## This program 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.
+## 
+## This program 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 Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## [ g, j, r ] = coupled_circuit_coeff (A, B, C, dt, x)
+## 
+## Compute coefficients of circuit-device coupling of the form:
+##
+## g * F + j(x) + r * I = 0;
+## 
+## 
+## INPUT
+## A B C : descriptor system 
+## dt : the time step in the backward Euler discretization;
+## x : state variables.
+## 
+## OUTPUT
+## g : coefficient of the patch voltage
+## j : constant bias term
+## r : coefficient of the current in the node
+## 
+## Author: davide <davide@davide-K53SV>
+## Created: 2013-06-10
+
+function [g, j, r] = coupled_circuit_coeff (A, B, C, dt, x)
+
+  freq = 1/dt;
+  
+  a12 = A(1,2:end);
+  a22 = A(2:end,2:end);
+
+  b11 = B(1,1);
+  b12 = B(1,2:end);
+  b21 = B(2:end,1);
+  b22 = B(2:end,2:end);
+
+  e11 = b11;
+  e12 = a12 / dt + b12;
+  e21 = b21;
+  e22 = a22 / dt + b22;
+  f1 = C(1);
+  f2 = C(2:end);
+
+  g = e11 - e12 * (e22 \ e21);
+  j = f1 - e12 * (e22 \ f2) + (-a12 + e12 * (e22 \ a22)) * x / dt;
+  r = 1;
+
+endfunction
--- a/extra/secs1d/inst/demos/coupled_circuit_coeff.m	Wed Jun 12 14:28:22 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,63 +0,0 @@
-## Copyright (C) 2013 davide
-## 
-## This program 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.
-## 
-## This program 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 Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## [ g, j, r ] = coupled_circuit_coeff (A, B, C, dt, x)
-## 
-## Compute coefficients of circuit-device coupling of the form:
-##
-## g * F + j(x) + r * I = 0;
-## 
-## 
-## INPUT
-## A B C : descriptor system 
-## dt : the time step in the backward Euler discretization;
-## x : state variables.
-## 
-## OUTPUT
-## g : coefficient of the patch voltage
-## j : constant bias term
-## r : coefficient of the current in the node
-## 
-## Author: davide <davide@davide-K53SV>
-## Created: 2013-06-10
-
-function [ g, j, r ] = coupled_circuit_coeff (A, B, C, dt, x)
-
-freq = 1/dt;
-
-%a{1,1} = A(1,1); % 0 by design
-a{1,2} = A(1,2:end);
-%a{2,1} = A(2:end,1); % 0 by design
-a{2,2} = A(2:end,2:end);
-
-b{1,1} = B(1,1);
-b{1,2} = B(1,2:end);
-b{2,1} = B(2:end,1);
-b{2,2} = B(2:end,2:end);
-
-e{1,1} = b{1,1};
-e{1,2} = a{1,2}*freq+b{1,2};
-e{2,1} = b{2,1};
-e{2,2} = a{2,2}*freq+b{2,2};
-f{1} = C(1);
-f{2} = C(2:end);
-
-g = e{1,1} - e{1,2} * (e{2,2} \ e{2,1});
-j = f{1} - e{1,2} * (e{2,2} \ f{2}) + freq*( -a{1,2} + e{1,2} * (e{2,2} \ a{2,2})) * x;
-r = 1;
-
-return
-endfunction
--- a/extra/secs1d/inst/demos/decomposematrixes.m	Wed Jun 12 14:28:22 2013 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,63 +0,0 @@
-## Copyright (C) 2013 davide
-## 
-## This program 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.
-## 
-## This program 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 Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## [ g, j, r ] = decomposematrixes (A, B, C, dt, x)
-## 
-## Compute coefficients of circuit-device coupling of the form:
-##
-## g * F + j(x) + r * I = 0;
-## 
-## 
-## INPUT
-## A B C : descriptor system 
-## dt : the time step in the backward Euler discretization;
-## x : state variables.
-## 
-## OUTPUT
-## g : coefficient of the patch voltage
-## j : constant bias term
-## r : coefficient of the current in the node
-## 
-## Author: davide <davide@davide-K53SV>
-## Created: 2013-06-10
-
-function [ g, j, r ] = decomposematrixes (A, B, C, dt, x)
-
-freq = 1/dt;
-
-%a{1,1} = A(1,1); % 0 by design
-a{1,2} = A(1,2:end);
-%a{2,1} = A(2:end,1); % 0 by design
-a{2,2} = A(2:end,2:end);
-
-b{1,1} = B(1,1);
-b{1,2} = B(1,2:end);
-b{2,1} = B(2:end,1);
-b{2,2} = B(2:end,2:end);
-
-e{1,1} = b{1,1};
-e{1,2} = a{1,2}*freq+b{1,2};
-e{2,1} = b{2,1};
-e{2,2} = a{2,2}*freq+b{2,2};
-f{1} = C(1);
-f{2} = C(2:end);
-
-g = e{1,1} - e{1,2} * (e{2,2} \ e{2,1});
-j = f{1} - e{1,2} * (e{2,2} \ f{2}) + freq*( -a{1,2} + e{1,2} * (e{2,2} \ a{2,2})) * x;
-r = 1;
-
-return
-endfunction
--- a/extra/secs1d/inst/demos/simple_resistor_tran_noscale_res.m	Wed Jun 12 14:28:22 2013 +0000
+++ b/extra/secs1d/inst/demos/simple_resistor_tran_noscale_res.m	Wed Jun 12 14:32:26 2013 +0000
@@ -3,7 +3,7 @@
 material  = secs1d_silicon_material_properties_fun (constants);
 
 % geometry
-Nelements = 100;
+Nelements = 10;
 L  = 100e-6;          % [m] 
 xm = L/2;
 device.W = 1e-6 * 1e-6;
@@ -19,7 +19,7 @@
 
 % time span for simulation
 tmin  = 0;
-tmax  = 20;
+tmax  = 16;
 tspan = [tmin, tmax];
 
 Fn = Fp = zeros (size (device.x));
@@ -51,33 +51,30 @@
   if (isempty (A1))
     load ("resistor_circuit_matrices")
   endif
-  C1(5) = -min(t,1);
+  C1(4) = -min(t,1);
   [g(1) j(1) r(1)] = coupled_circuit_coeff (A1, B1, C1, dt, x1);
   g(2) = 1;   j(2) = 0;   r(2) = 0;
 endfunction
 
+
 function [x1, x2] = update_states (t, dt, F1)
-   persistent A1 B1 C1 x1 
-   %in this case  it's not necessary for A2 B2 C2 x2
+   persistent A1 B1 C1 x1 %A2 B2 C2 x2
    if (isempty (A1))
      load ("resistor_circuit_matrices")
-   endif
+   else
+     C1(4) = -min(t,1);
+     a22 = A1(2:end,2:end);
+     b21 = B1(2:end,1);
+     b22 = B1(2:end,2:end);
+     e22 = a22/dt+b22;
+     f2 = C1(2:end);
 
-   C1(5) = -min(t,1);
-   freq = 1 / dt;
-   
-   a{2,2} = A1(2:end,2:end);
-   b{2,1} = B1(2:end,1);
-   b{2,2} = B1(2:end,2:end);
-   e{2,2} = a{2,2} * freq + b{2,2};
-   f{2} = C1(2:end);
-
-   x1 = e{2,2} \ (freq*a{2,2}*x1 - b{2,1}*F1 - f{2});
-   x2=[];
+     x1 = e22 \ (a22 * x1 / dt - b21 * F1 - f2);
+     x2 = [];
 endfunction
 
 % tolerances for convergence checks
-algorithm.toll       = 1e-8;
+algorithm.toll       = 1e-6;
 algorithm.ltol       = 1e-10;
 algorithm.maxit      = 100;
 algorithm.lmaxit     = 100;
@@ -85,7 +82,7 @@
 algorithm.pmaxit     = 1000;
 algorithm.colscaling = [10 1e21 1e21 1];
 algorithm.rowscaling = [1e7 1e-7 1e-7 1];
-algorithm.maxnpincr  = 1e-5;
+algorithm.maxnpincr  = 5e-2;
 
 %% compute resistance
 u = secs1d_mobility_model_noscale ...
@@ -102,8 +99,9 @@
 close all; secs1d_logplot (device.x, device.D, 'x-'); pause
 
 %% (pseudo)transient simulation
-[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs1d_newton_res (device, material, constants, algorithm,
-                                                           Vin, nin, pin, tspan, @vbcs, @update_states);
+[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs1d_newton_res ...
+                                          (device, material, constants, algorithm,
+                                           Vin, nin, pin, tspan, @vbcs, @update_states);
 
 dV   = diff (V, [], 1);
 dx   = diff (device.x);
--- a/extra/secs1d/inst/secs1d_newton_res.m	Wed Jun 12 14:28:22 2013 +0000
+++ b/extra/secs1d/inst/secs1d_newton_res.m	Wed Jun 12 14:32:26 2013 +0000
@@ -7,8 +7,10 @@
   dt = (tspan(2) - tspan(1)) / 20;  
   t(tstep = 1) = tspan (1);
   [V, n, p] = deal (Vin, nin, pin);  
-  F = V([1 end], 1) - constants.Vth * log (n([1 end], 1) ...
-                                           ./ device.ni([1 end], :));
+  F = V([1 end], 1) - constants.Vth * ...
+                      log (n([1 end], 1) ./ device.ni([1 end], :));
+  [x1, x2] = upd (t, dt, F(1), F(2));
+
   M = bim1a_reaction (device.x, 1, 1);
 
   [x1, x2] = upd(t, dt, F(1), F(2));
@@ -18,7 +20,7 @@
     reject = false;
     t = tout(++tstep) = min (t + dt, tspan(2)); 
     incr0 = 4 * algorithm.maxnpincr;
-
+    
     [gi, ji, ri] = va (t, dt, x1, x2);
     [V0, n0, p0, F0] = predict (device, material, constants, ...
                                 algorithm, V, n, p, F, tstep, ...
@@ -182,6 +184,7 @@
 
     else
 
+      [x1, x2] = upd (t, dt, F(1), F(2));
       Fp(:, tstep) = V2 + constants.Vth * log (p2 ./ device.ni);
       Fn(:, tstep) = V2 - constants.Vth * log (n2 ./ device.ni);        
       [V(:, tstep), n(:, tstep), p(:, tstep), F(:, tstep)] = deal (V2, n2, p2, F2);