Mercurial > forge
changeset 11745:a6f6a1deefbb octave-forge
cleanup mods
| author | cdf |
|---|---|
| date | Wed, 05 Jun 2013 05:28:04 +0000 |
| parents | d995950329f1 |
| children | 94f687098b2e |
| files | extra/secs1d/inst/demos/simple_diode_tran_noscale.m extra/secs1d/inst/demos/simple_diode_tran_noscale_res.m extra/secs1d/inst/demos/thyristor_tran_noscale.m extra/secs1d/inst/secs1d_impact_ionization_noscale.m extra/secs1d/inst/secs1d_newton_res.m |
| diffstat | 5 files changed, 149 insertions(+), 50 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/secs1d/inst/demos/simple_diode_tran_noscale.m Wed Jun 05 05:28:04 2013 +0000 @@ -0,0 +1,91 @@ +% physical constants and parameters +constants = secs1d_physical_constants_fun (); +material = secs1d_silicon_material_properties_fun (constants); + +% geometry +Nelements = 1000; +L = 1e-6; % [m] +xm = L/2; +device.W = 1e-6 * 1e-6; +device.x = linspace (0, L, Nelements+1)'; +device.sinodes = [1:length(device.x)]; + +% doping profile [m^{-3}] +device.Na = 1e23 * (device.x <= xm); +device.Nd = 1e23 * (device.x > xm); +device.R = [0; 0]; + +% avoid zero doping +device.D = device.Nd - device.Na; + +% time span for simulation +tmin = 0; +tmax = 100; +tspan = [tmin, tmax]; + +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); + +function fn = vbcs_1 (t); + fn = [0; 0]; + fn(2) = -t; +endfunction + +function fp = vbcs_2 (t); + fp = [0; 0]; + fp(2) = -t; +endfunction + +vbcs = {@vbcs_1, @vbcs_2}; + +% tolerances for convergence checks +algorithm.toll = 1e-5; +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 1e20 1e20 1]; +algorithm.maxnpincr = 1e-2; + +%% 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); + +close all; semilogy (device.x, nin .* pin, 'x-'); pause + +%% (pseudo)transient simulation +[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs1d_newton (device, material, constants, algorithm, + Vin, nin, pin, tspan, @vbcs_1); + +dV = diff (V, [], 1); +dx = diff (device.x); +E = -dV ./ dx; + +vvector = Fn(end, :); +ivector = Itot (2, :); + +plotyy (tout, vvector, tout, ivector) +drawnow +
--- a/extra/secs1d/inst/demos/simple_diode_tran_noscale_res.m Wed Jun 05 05:23:23 2013 +0000 +++ b/extra/secs1d/inst/demos/simple_diode_tran_noscale_res.m Wed Jun 05 05:28:04 2013 +0000 @@ -3,7 +3,7 @@ material = secs1d_silicon_material_properties_fun (constants); % geometry -Nelements = 19; +Nelements = 3000; L = 1e-6; % [m] xm = L/2; device.W = 1e-6 * 1e-6; @@ -11,15 +11,15 @@ device.sinodes = [1:length(device.x)]; % doping profile [m^{-3}] -device.Na = 1e21 * (device.x <= xm); -device.Nd = 1e21 * (device.x > xm); +device.Na = 1e23 * exp ( - .5 * device.x.^2 / xm^2); +device.Nd = 1e23 * exp ( - .5 * (device.x - L) .^2 / xm^2); % avoid zero doping device.D = device.Nd - device.Na; % time span for simulation tmin = 0; -tmax = 100; +tmax = 10; tspan = [tmin, tmax]; Fn = Fp = zeros (size (device.x)); @@ -46,30 +46,30 @@ V = Fn + constants.Vth * log (n ./ device.ni); function [g, j] = vbcs (t, dt); - g = -[1; 1]; - j = [t 0]; + g = -[10; 10]; + j = [t*10 0]; endfunction % tolerances for convergence checks -algorithm.toll = 1e-5; +algorithm.toll = 1e-7; algorithm.ltol = 1e-10; -algorithm.maxit = 100; +algorithm.maxit = 5; algorithm.lmaxit = 100; algorithm.ptoll = 1e-12; algorithm.pmaxit = 1000; -algorithm.colscaling = [10 1e21 1e21 1]; -algorithm.rowscaling = [1 1e20 1e20 1]; -algorithm.maxnpincr = 1e-2; +algorithm.colscaling = [10 1e21 1e21 .1]; +algorithm.rowscaling = [1 1e-7 1e-7 .1]; +algorithm.maxnpincr = 1e-3; %% 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); -close all; semilogy (device.x, nin .* pin, 'x-'); pause +close all; secs1d_logplot (device.x, device.D, 'x-'); pause %% (pseudo)transient simulation -[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs1d_newton (device, material, constants, algorithm, - Vin, nin, pin, tspan, @vbcs); +[V, n, p, Fn, Fp, Jn, Jp, Itot, tout] = secs1d_newton_res (device, material, constants, algorithm, + Vin, nin, pin, tspan, @vbcs); dV = diff (V, [], 1); dx = diff (device.x);
--- a/extra/secs1d/inst/demos/thyristor_tran_noscale.m Wed Jun 05 05:23:23 2013 +0000 +++ b/extra/secs1d/inst/demos/thyristor_tran_noscale.m Wed Jun 05 05:28:04 2013 +0000 @@ -79,9 +79,9 @@ algorithm.ptoll = 1e-12; algorithm.pmaxit = 100; algorithm.maxnpincr = constants.Vth; -algorithm.colscaling = [10 1e23 1e23]; -algorithm.rowscaling = [1e6 1e23 1e23]; -algorithm.maxnpincr = constants.Vth; +algorithm.colscaling = [1 1e23 1e23]; +algorithm.rowscaling = [1 1e7 1e7]; +algorithm.maxnpincr = constants.Vth / 10; logplot = @(x) asinh (x/2) / log(10); close all; plot (device.x, logplot (device.D)); @@ -105,7 +105,7 @@ %% initial guess via (pseudo)transient simulation [n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... secs1d_tran_dd_newton_noscale (device, material, constants, algorithm, - Vin, nin, pin, Fn, Fp, tspan, vbcs); + Vin, nin, pin, Fn, Fp, tspan, vbcs); Vin = V(:, end); nin = n(:, end); @@ -115,12 +115,12 @@ function fn = vbcs_1 (t); fn = [0; 0]; - fn(2) = -10*t; + fn(2) = -10*2*t/7; endfunction function fp = vbcs_2 (t); fp = [0; 0]; - fp(2) = -10*t; + fp(2) = -10*2*t/7; endfunction vbcs = {@vbcs_1, @vbcs_2}; @@ -157,8 +157,8 @@ [n, p, V, Fn, Fp, Jnpos, Jppos, t, it, res] = ... secs1d_tran_dd_newton_noscale (device, material, constants, algorithm, - Vin, nin, pin, - Fnin, Fpin, tspan, vbcs); + Vin, nin, pin, + Fnin, Fpin, tspan, vbcs); vvectorpos = Fn(end, :); ivectorpos = (Jnpos(end, :) + Jppos(end, :));
--- a/extra/secs1d/inst/secs1d_impact_ionization_noscale.m Wed Jun 05 05:23:23 2013 +0000 +++ b/extra/secs1d/inst/secs1d_impact_ionization_noscale.m Wed Jun 05 05:28:04 2013 +0000 @@ -16,7 +16,8 @@ %% You should have received a copy of the GNU General Public License %% along with SECS1D; If not, see <http://www.gnu.org/licenses/>. -%% II = secs1d_impact_ionization_noscale (E, Jn, Jp) +%% II = secs1d_impact_ionization_noscale ... +%% (device, material, constants, algorithm, E, Jn, Jp, V, n, p, Fn, Fp) function II = secs1d_impact_ionization_noscale ... (device, material, constants, algorithm, E, Jn, Jp, V, n, p, Fn, Fp)
--- a/extra/secs1d/inst/secs1d_newton_res.m Wed Jun 05 05:23:23 2013 +0000 +++ b/extra/secs1d/inst/secs1d_newton_res.m Wed Jun 05 05:28:04 2013 +0000 @@ -4,9 +4,11 @@ Nnodes = numel (device.x); dt = (tspan(2) - tspan(1)) / 20; t(tstep = 1) = tspan (1); - [V, n, p] = deal (Vin, nin, pin); + [V, n, p] = deal (Vin, nin, pin); + F = V([1 end], 1) - constants.Vth * log (n([1 end], 1) ... + ./ device.ni([1 end], :)); M = bim1a_reaction (device.x, 1, 1); - + while (t < tspan(2)) reject = false; @@ -15,8 +17,8 @@ [gi, ji] = va (t, dt); [V0, n0, p0, F0] = predict (device, material, constants, ... - algorithm, V, n, p, tstep, tout, gi, ji); - + algorithm, V, n, p, F, tstep, tout, gi, ji); + [V2, n2, p2, F2] = deal (V0, n0, p0, F0); for in = 1 : algorithm.maxit @@ -58,8 +60,8 @@ ## else J = [([algorithm.colscaling(1)*jac{1,1}, ... - algorithm.colscaling(2)*jac{1,2}(2:end-1, 2:end-1), ... - algorithm.colscaling(3)*jac{1,3}(2:end-1, 2:end-1), ... + algorithm.colscaling(2)*jac{1,2}(:, 2:end-1), ... + algorithm.colscaling(3)*jac{1,3}(:, 2:end-1), ... algorithm.colscaling(4)*jac{1,4}]/algorithm.rowscaling(1)); ... ## ([algorithm.colscaling(1)*jac{2,1}(2:end-1, :), ... @@ -85,10 +87,10 @@ delta = - J \ f; - dV = delta * algorithm.colscaling(1); + dV = delta (1:Nnodes) * algorithm.colscaling(1); dn = delta (Nnodes+1:2*(Nnodes)-2) * algorithm.colscaling(2); dp = delta (2*(Nnodes)-1:3*(Nnodes)-4) * algorithm.colscaling(3); - dF = delta (3*(Nnodes)-4:end) * algorithm.colscaling(4); + dF = delta (3*(Nnodes)-3:end) * algorithm.colscaling(4); ## endif @@ -110,7 +112,7 @@ tkp = .9 * min (p1(2:end-1) ./ abs (dp)) endif - tk = min ([tkv, tkn, tkp]) + tk = min ([tkv, tkn, tkp]); if (tk <= 0) error ("relaxation parameter too small") endif @@ -173,8 +175,7 @@ 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, F); + [V(:, tstep), n(:, tstep), p(:, tstep), F(:, tstep)] = deal (V2, n2, p2, F2); [mobilityn, mobilityp] = compute_mobilities ... (device, material, constants, algorithm, V2, n2, p2); @@ -227,8 +228,8 @@ A11 = bim1a_laplacian (device.x, epsilon, 1); res{1} = A11 * V + bim1a_rhs (device.x, 1, constants.q * (n - p - device.D)); - res{1}(1) = A11(1,1) (V(1) + constants.Vth*log(p(1)./device.ni(1))-F(1)); - res{1}(end) = A11(end,end) (V(end)+constants.Vth*log(p(end)./device.ni(end))-F(end)); + res{1}(1) = A11(1,1) * (V(1) + constants.Vth*log(p(1)./device.ni(1))-F(1)); + res{1}(end) = A11(end,end) * (V(end)+constants.Vth*log(p(end)./device.ni(end))-F(end)); A22 = bim1a_advection_diffusion ... (device.x, mobilityn * constants.Vth, 1, 1, V / constants.Vth); @@ -240,17 +241,20 @@ I1 = - constants.q * A22(1,:) * n + constants.q * A33(1,:) * p; - I1 += constants.e0 * material.esir * ... - ((V(2, tstep) - V(1, tstep)) - - (V(2, tstep-1) - V(1, tstep-1))) / ... - (device.x(2) - device.x(1)); - I1 *= device.W; + I2 = - constants.q * A22(end,:) * n + constants.q * A33(end,:) * p; + + if (columns (V) >= 2) + I1 += 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 * ... + ((V(end, end) - V(end-1, end)) - + (V(end, end-1) - V(end-1, end-1))) / ... + (device.x(end) - device.x(end-1)); + endif - I2 = - constants.q * A22(end,:) * n + constants.q * A33(end,:) * p; - I2 += 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)); + I1 *= device.W; I2 *= device.W; res{4} = [gi(1)*F(1)+ji(1)+I1; @@ -377,9 +381,9 @@ endfunction -function [V0, n0, p0] = predict (device, material, constants, ... - algorithm, V, n, p, F, tstep, ... - tout, gi, ji) +function [V0, n0, p0, F0] = predict (device, material, constants, ... + algorithm, V, n, p, F, tstep, ... + tout, gi, ji) if (tstep > 2) @@ -391,12 +395,15 @@ Fn(:, 2) = V(:, it(2)) - constants.Vth * log (n(:, it(2)) ./ device.ni); dFndt = diff (Fn, 1, 2) ./ diff (tout(it)); + dFdt = diff (F(:, it(1:2)), 1, 2) ./ diff (tout(it)); Fn0 = Fn(:, 2) + dFndt * dt; + F0 = F(:, it(2)) + dFdt * dt; else Fn0 = V(:, tstep-1) - constants.Vth * log (n(:, tstep-1) ./ device.ni); + F0 = F(:, tstep-1); endif @@ -406,5 +413,5 @@ %Fn0([1 end]) = va (tout (tstep)); V0 = Fn0 + constants.Vth * log (n0 ./ device.ni); - + endfunction