Mercurial > forge
changeset 11572:cdd6e0bbf196 octave-forge
newton method
| author | cdf |
|---|---|
| date | Mon, 25 Mar 2013 11:04:56 +0000 |
| parents | 10f755d08a29 |
| children | 383535d6d361 |
| files | extra/secs1d/inst/demos/diode_reverse_noscale.m extra/secs1d/inst/demos/diode_tran_noscale.m extra/secs1d/inst/demos/thyristor_tran_noscale.m extra/secs1d/inst/secs1d_dd_newton_noscale.m extra/secs1d/inst/secs1d_tran_dd_newton_noscale.m |
| diffstat | 5 files changed, 601 insertions(+), 89 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/secs1d/inst/demos/diode_reverse_noscale.m Sun Mar 24 15:35:39 2013 +0000 +++ b/extra/secs1d/inst/demos/diode_reverse_noscale.m Mon Mar 25 11:04:56 2013 +0000 @@ -25,8 +25,8 @@ device.tn = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'n'); % initial guess for n, p, V, phin, phip -V_p = 1; -V_n = 0; +V_p = 0; +V_n = -100; Fp = V_p * (device.x <= xm); Fn = Fp; @@ -45,34 +45,50 @@ V = Fn + constants.Vth * log (n ./ device.ni); % tolerances for convergence checks -algorithm.toll = 1e-10; -algorithm.maxit = 100; -algorithm.ptoll = 1e-12; -algorithm.pmaxit = 1000; +algorithm.toll = 1e-10; +algorithm.maxit = 1000; +algorithm.ptoll = 1e-12; +algorithm.pmaxit = 1000; algorithm.colscaling = [10 1e23 1e23]; algorithm.rowscaling = [1e6 1e23 1e23]; -algorithm.maxnpincr = constants.Vth / 10; +algorithm.maxnpincr = constants.Vth; + % solve the problem using the full DD model -[n, p, V, Fn, Fp, Jn, Jp, it, res] = ... - secs1d_dd_gummel_map_noscale (device, material, - constants, algorithm, - V, n, p, Fn, Fp); +%algorithm.maxit = 1; +%[n, p, V, Fn, Fp, Jn, Jp, it, res] = ... +# secs1d_dd_gummel_map_noscale (device, material, +# constants, algorithm, +# V, n, p, Fn, Fp); +# pause -algorithm.maxit = 1000; -[n_newt, p_newt, V_newt, Fn_newt, Fp_newt, Jn, Jp, it, res] = secs1d_dd_newton_noscale ... +[V, n, p, res, niter] = ... + secs1d_nlpoisson_newton_noscale (device, material, + constants, algorithm, + V, n, p, Fn, Fp); + +[ng, pg, Vg, Fng, Fpg, Jn, Jp, it, res] = secs1d_dd_gummel_map_noscale ... + (device, material, constants, algorithm, V, n, p, Fn, Fp); + +[n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_newton_noscale ... (device, material, constants, algorithm, V, n, p, Fn, Fp); + +semilogy (device.x, n, device.x, p, device.x, ng, device.x, pg) +legend ('n', 'p', 'ng', 'pg') +keyboard + % band structure Efn = -Fn; Efp = -Fp; Ec = constants.Vth * log (material.Nc./n) + Efn; Ev = -constants.Vth * log (material.Nv./p) + Efp; -close all + plot (device.x, Efn, device.x, Efp, device.x, Ec, device.x, Ev) + legend ('Efn', 'Efp', 'Ec', 'Ev') axis tight
--- a/extra/secs1d/inst/demos/diode_tran_noscale.m Sun Mar 24 15:35:39 2013 +0000 +++ b/extra/secs1d/inst/demos/diode_tran_noscale.m Mon Mar 25 11:04:56 2013 +0000 @@ -62,33 +62,37 @@ % tolerances for convergence checks algorithm.toll = 1e-6; -algorithm.maxit = 300; +algorithm.maxit = 100; algorithm.ptoll = 1e-12; algorithm.pmaxit = 1000; -algorithm.maxnpincr = constants.Vth / 10; +algorithm.colscaling = [10 1e23 1e23]; +algorithm.rowscaling = [1e6 1e23 1e23]; +algorithm.maxnpincr = constants.Vth; %% 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); +[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, 'x-', device.x, pin, 'x-'); pause %% (pseudo)transient simulation -[n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... - secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm, - Vin, nin, pin, Fnin, Fpin, tspan, vbcs); +% [n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... +% secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm, +% Vin, nin, pin, Fnin, Fpin, tspan, vbcs); + + [n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... + secs1d_tran_dd_newton_noscale (device, material, constants, algorithm, + Vin, nin, pin, Fnin, Fpin, tspan, vbcs); dV = diff (V, [], 1); dx = diff (device.x); E = -dV ./ dx; %% band structure -Efn = -Fn; -Efp = -Fp; -Ec = constants.Vth * log (material.Nc ./ n) + Efn; -Ev = -constants.Vth * log (material.Nv ./ p) + Efp; +%% Efn = -Fn; +%% Efp = -Fp; +%% Ec = constants.Vth * log (material.Nc ./ n) + Efn; +%% Ev = -constants.Vth * log (material.Nv ./ p) + Efp; %## figure (1) %## plot (x, Efn, x, Efp, x, Ec, x, Ev) @@ -102,7 +106,7 @@ ivectora = mean (Jn + Jp, 1); area = 150e-6 * 50e-6; figure (1) -plot (vvector, area*ivector, vvector, area*ivectorn, vvector, area*ivectora) +plotyy (t, vvector, t, area*ivector) legend('J_L','J_0') drawnow
--- a/extra/secs1d/inst/demos/thyristor_tran_noscale.m Sun Mar 24 15:35:39 2013 +0000 +++ b/extra/secs1d/inst/demos/thyristor_tran_noscale.m Mon Mar 25 11:04:56 2013 +0000 @@ -74,11 +74,14 @@ vbcs = {@vbcs_1, @vbcs_2}; % tolerances for convergence checks -algorithm.toll = 1e-6; +algorithm.toll = 1e-4; algorithm.maxit = 100; 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; logplot = @(x) asinh (x/2) / log(10); close all; plot (device.x, logplot (device.D)); @@ -93,12 +96,46 @@ constants, algorithm, V, n, p, Fn, Fp); -close all; semilogy (device.x, nin, 'x-', device.x, pin, 'x-'); pause +close all; semilogy (device.x, nin, 'x-', device.x, pin, 'x-'); +hold on +semilogy (device.x, device.Na, device.x, device.Nd); +hold off +pause + +%% 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 = V(:, end); +nin = n(:, end); +pin = p(:, end); +Fnin = Fn(:, end); +Fpin = Fp(:, end); -%% (pseudo)transient simulation -[n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... - secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm, - Vin, nin, pin, Fn, Fp, tspan, vbcs); +function fn = vbcs_1 (t); + fn = [0; 0]; + fn(2) = -10*t; +endfunction + +function fp = vbcs_2 (t); + fp = [0; 0]; + fp(2) = -10*t; +endfunction + +vbcs = {@vbcs_1, @vbcs_2}; + +%% (pseudo)transient simulation with negative applied volatage +# [n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... +# secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm, +# V(:, end), n(:, end), p(:, end), +# Fn(:, end), Fp(:, end), tspan, vbcs); + +[n, p, V, Fn, Fp, Jnneg, Jpneg, t, it, res] = ... + secs1d_tran_dd_newton_noscale (device, material, constants, algorithm, + Vin, nin, pin, + Fnin, Fpin, tspan, vbcs); +vvectorneg = Fn(end, :); function fn = vbcs_1 (t); fn = [0; 0]; @@ -112,33 +149,20 @@ vbcs = {@vbcs_1, @vbcs_2}; -%% (pseudo)transient simulation -[n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... - secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm, - V(:, end), n(:, end), p(:, end), - Fn(:, end), Fp(:, end), tspan, vbcs); +%% (pseudo)transient simulation with positive applied volatage +# [n, p, V, Fn, Fp, Jn, Jp, t, it, res] = ... +# secs1d_tran_dd_gummel_map_noscale (device, material, constants, algorithm, +# V(:, end), n(:, end), p(:, end), +# Fn(:, end), Fp(:, end), tspan, vbcs); -dV = diff (V, [], 1); -dx = diff (device.x); -E = -dV ./ dx; - -%% band structure -Efn = -Fn(:, end); -Efp = -Fp(:, end); -Ec = constants.Vth * log (material.Nc ./ n(:, end)) + Efn; -Ev = -constants.Vth * log (material.Nv ./ p(:, end)) + Efp; - -figure (1) -plot (device.x, Efn, device.x, Efp, device.x, Ec, device.x, Ev) -legend ('Efn', 'Efp', 'Ec', 'Ev') -axis tight -drawnow +[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); -vvector = Fn(end, :); -ivector = (Jn(end, :) + Jp(end, :)); -ivectorn = (Jn(1, :) + Jp(1, :)); +vvectorpos = Fn(end, :); +ivectorpos = (Jnpos(end, :) + Jppos(end, :)); +ivectorneg = (Jnneg(end, :) + Jpneg(end, :)); -figure (2) -plot (vvector, ivector, vvector, ivectorn) -legend('J_L','J_0') -drawnow +plot (vvectorpos, ivectorpos, vvectorneg, ivectorneg) +drawnow \ No newline at end of file
--- a/extra/secs1d/inst/secs1d_dd_newton_noscale.m Sun Mar 24 15:35:39 2013 +0000 +++ b/extra/secs1d/inst/secs1d_dd_newton_noscale.m Mon Mar 25 11:04:56 2013 +0000 @@ -53,52 +53,64 @@ function [n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_dd_newton_noscale ... (device, material, constants, algorithm, Vin, nin, pin, Fnin, Fpin) - dampcoeff = 2; Nnodes = numel (device.x); - V = Vin; n = nin; p = pin; [res, jac] = residual_jacobian (device, material, constants, algorithm, V, n, p); - normr(1) = norm (res, inf); - normrnew = normr(1); - for it = 1 : algorithm.maxit delta = - jac \ res; - - tk = 1; - dv = norm (delta(1:Nnodes-2), inf); - if (dv > algorithm.maxnpincr) - tk = algorithm.maxnpincr / dv; + + dv = delta(1:Nnodes-2) * algorithm.colscaling(1); + dn = delta((Nnodes-2)+(1:Nnodes-2)) * algorithm.colscaling(2); + dp = delta(2*(Nnodes-2)+(1:Nnodes-2)) * algorithm.colscaling(3); + + ndv = norm (dv, inf); + ndn = norm (dn, inf); + ndp = norm (dp, inf); + + tkv = 1; + if (ndv > algorithm.maxnpincr) + tkv = algorithm.maxnpincr / ndv; endif - Vnew = Vin; - Vnew(2:end-1) = V(2:end-1) + tk * delta(1:Nnodes-2) * algorithm.colscaling(1); - - nnew = nin; - nnew(2:end-1) = n(2:end-1) + tk * delta((Nnodes-2)+(1:Nnodes-2)) * algorithm.colscaling(2); + tkn = tkv; + [howmuch, where] = min (n(2:end-1) + tkn * dn); + if (howmuch <= 0) + tkn = - .9 * n(2:end-1)(where) / dn(where) + endif - pnew = pin; - pnew(2:end-1) = p(2:end-1) + tk * delta(2*(Nnodes-2)+(1:Nnodes-2)) * algorithm.colscaling(3); + tkp = tkn; + [howmuch, where] = min (p(2:end-1) + tkp * dp); + if (howmuch <= 0) + tkp = - .9 * p(2:end-1)(where) / dp(where) + endif - - V = Vnew; n = nnew; p = pnew; + tk = min ([tkv, tkn, tkp]) + + V(2:end-1) += tk * dv; + n(2:end-1) += tk * dn; + p(2:end-1) += tk * dp; + [res, jac] = residual_jacobian (device, material, constants, algorithm, V, n, p); - resvec(it) = reldnorm = + ... - norm (delta(1:Nnodes-2), inf) / norm (V, inf) + ... - norm (delta((Nnodes-2)+(1:Nnodes-2)), inf) / norm (n, inf) + ... - norm (delta(2*(Nnodes-2)+(1:Nnodes-2)), inf) / norm (p, inf); + resvec(it) = reldnorm = ... + ndv / norm (V, inf) + ... + ndn / norm (n, inf) + ... + ndp / norm (p, inf); - plotyy ( - 1:it, log10 (resvec), - 1:it+1, log10 (normr) - ) + figure (2) + semilogy (1:it, resvec) drawnow - + + figure (3) + semilogy (device.x, n, device.x, p) + drawnow + pause + if (reldnorm <= algorithm.toll) break endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/secs1d/inst/secs1d_tran_dd_newton_noscale.m Mon Mar 25 11:04:56 2013 +0000 @@ -0,0 +1,456 @@ +%% Copyright (C) 2004-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/>. + +%% Solve the scaled stationary bipolar DD equation system using Newton's method. +%% +%% [n, p, V, Fn, Fp, Jn, Jp, it, res] = secs1d_tran_dd_newton_noscale ... +%% (device, material, constants, algorithm, Vin, nin, pin, Fnin, Fpin, tspan, vbcs) +%% +%% input: +%% device.x spatial grid +%% tspan = [tmin, tmax] time integration interval +%% vbcs = {fhnbc, fpbc} cell aray of function handles to compute applied potentials as a function of time +%% device.{D,Na,Nd} doping profile +%% Vin initial guess for electrostatic potential +%% nin initial guess for electron concentration +%% pin initial guess for hole concentration +%% Fnin initial guess for electron Fermi potential +%% Fpin initial guess for hole Fermi potential +%% material.{u0n, uminn, vsatn, Nrefn} +%% electron mobility model coefficients +%% material.{u0p, uminp, vsatp, Nrefp} +%% hole mobility model coefficients +%% device.ni intrinsic carrier density +%% material.{tn,tp,Cn,Cp,an,ap,Ecritnin,Ecritpin} +%% generation recombination model parameters +%% algorithm.toll tolerance for Gummel iterarion convergence test +%% algorithm.maxit maximum number of Gummel iterarions +%% +%% output: +%% n electron concentration +%% p hole concentration +%% V electrostatic potential +%% Fn electron Fermi potential +%% Fp hole Fermi potential +%% Jn electron current density +%% Jp hole current density +%% it number of Newton iterations performed +%% res total potential increment at each step + +function [n, p, V, Fn, Fp, Jn, Jp, tout, numit, res] = ... + secs1d_tran_dd_newton_noscale (device, material, constants, algorithm, + Vin, nin, pin, Fnin, Fpin, tspan, vbcs) + + rejected = 0; + + V = Vin; n = nin; p = pin; Fn = Fnin; Fp = Fpin; + + Nnodes = numel (device.x); + Nelements = Nnodes -1; + + [mobilityn, mobilityp] = compute_mobilities ... + (device, material, constants, algorithm, V, n, p); + + [Jn, Jp] = compute_currents ... + (device, material, constants, algorithm, mobilityn, + mobilityp, V, n, p); + + tmax = tspan(2); + tmin = tspan(1); + dt = (tmax - tmin) / 2000; %% initial time step + + tstep = 1; + t = tout(tstep) = tmin; + + resist = 1e7 * 1e-3 ^ 2; + + while (t < tmax) + + try + + t = tout(++tstep) = min (t + dt, tmax); + + n0 = n(:, tstep - 1); + p0 = p(:, tstep - 1); + + if (tstep <= 2) + Fn0 = Fn(:,tstep-1); + Fp0 = Fp(:,tstep-1); + else + Fn0 = Fn(:, tstep-2) .* (tout(tstep-1) - t)/(tout(tstep-1) - tout(tstep-2)) + ... + Fn(:, tstep-1) .* (t - tout(:, tstep-2))/(tout(:, tstep-1) - tout(:, tstep-2)); + Fp0 = Fp(:, tstep-2) .* (tout(tstep-1) - t)/(tout(tstep-1) - tout(tstep-2)) + ... + Fp(:, tstep-1) .* (t - tout(:, tstep-2))/(tout(:, tstep-1) - tout(:, tstep-2)); + endif + + Va = vbcs{1} (t); + Jn(:,tstep) = Jn(:,tstep-1); + Jp(:,tstep) = Jp(:,tstep-1); + + Fn0([1 end]) = Va; + Fp0([1 end]) = Va; + + 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, nrfnp] =... + __one_step__ (device, material, constants, algorithm, + V0, n0, p0, Fn0, Fp0, Jn(:,tstep), Jp(:,tstep), + n(:, tstep - 1), p(:, tstep - 1), dt); + + # [V_gummel, n_gummel, p_gummel, Fn_gummel, ... + # Fp_gummel, Jn_gummel, Jp_gummel, ~, ~, nrfnp] =... + # __one_gummel_step__ (device, material, constants, algorithm, + # V0, n0, p0, Fn0, Fp0, Jn(:,tstep), Jp(:,tstep), + # n(:, tstep - 1), p(:, tstep - 1), dt); + + # figure(1001) + # plot (device.x, n_gummel ./ n(:,end)) + + # figure(1002) + # plot (device.x, p_gummel ./ p(:,end)) + + # figure(1003) + # plot (device.x, Fp_gummel - Fp(:,end)) + + # figure(1004) + # plot (device.x, Fn_gummel - Fn(:,end)) + + %keyboard + + %assert (nrfnp <= algorithm.maxnpincr) + + figure (1) + plot (tout, (mean (Jn(:, 1:tstep) + Jp(:, 1:tstep), 1))) + drawnow + + res(tstep,1:it) = res_; + printf ("dt incremented by %g\n", .95 * sqrt (algorithm.maxnpincr / nrfnp)) + dt *= .95 * sqrt (algorithm.maxnpincr / nrfnp); + numit(tstep) = it; + + catch + + warning (lasterr) + warning ("algorithm.maxnpincr = %17.17g, nrfnp = %17.17g, dt %17.17g reduced to %17.17g", + algorithm.maxnpincr, nrfnp, dt, dt * .95 * sqrt (algorithm.maxnpincr / nrfnp)) + + dt *= .95 * sqrt (algorithm.maxnpincr / nrfnp); + + t = tout(--tstep); + if (dt < 100*eps) + warning ("secs1d_tran_dd_newton_noscale: minimum time step reached.") + return + endif + rejected ++ ; + + end_try_catch + + endwhile + printf ("number of rejected time steps: %d\n", rejected) + +endfunction + + +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); + + Fnref = Fn = Fn0; + Fpref = Fp = Fp0; + n = n0; + p = p0; + V = V0; + +% [V, n, p] = secs1d_nlpoisson_newton_noscale ... +% (device, material, constants, algorithm, V0, n0, p0, Fn0, Fp0); + + + Jn = Jn0; + Jp = Jp0; + + Nnodes = numel (device.x); + + [res, jac] = residual_jacobian (device, material, constants, + algorithm, V, n, p, n_old, p_old, dt); + + for it = 1 : algorithm.maxit + + delta = - jac \ res; + + dv = delta(1:Nnodes-2) * algorithm.colscaling(1); + dn = delta((Nnodes-2)+(1:Nnodes-2)) * algorithm.colscaling(2); + dp = delta(2*(Nnodes-2)+(1:Nnodes-2)) * algorithm.colscaling(3); + + ndv = norm (dv, inf); + ndn = norm (dn, inf); + ndp = norm (dp, inf); + + tkv = 1; + if (ndv > constants.Vth) + tkv = constants.Vth / ndv; + endif + + tkn = tkv; + [howmuch, where] = min (n(2:end-1) + tkn * dn); + if (howmuch <= 0) + tkn = - .9 * n(2:end-1)(where) / dn(where); + endif + + tkp = tkn; + [howmuch, where] = min (p(2:end-1) + tkp * dp); + if (howmuch <= 0) + tkp = - .9 * p(2:end-1)(where) / dp(where); + endif + + tk = min ([tkv, tkn, tkp]) + + V(2:end-1) += tk * dv; + n(2:end-1) += tk * dn; + p(2:end-1) += tk * dp; + + [res, jac] = residual_jacobian (device, material, constants, + algorithm, V, n, p, n_old, + p_old, dt); + + resvec(it) = reldnorm = ... + ndv / norm (V, inf) + ... + ndn / norm (n, inf) + ... + ndp / norm (p, inf); + + figure (3) + semilogy (1:it, resvec) + drawnow + + Fp = V + constants.Vth * log (p ./ device.ni); + Fn = V - constants.Vth * log (n ./ device.ni); + nrfnp = max (norm (Fn - Fnref, inf), norm (Fp - Fpref, inf)); + + if (reldnorm <= algorithm.toll) + %|| (it > 5 && nrfnp > algorithm.maxnpincr)) + break + endif + + endfor + + [mobilityn, mobilityp] = compute_mobilities ... + (device, material, constants, algorithm, V, n, p); + + [Jn, Jp] = compute_currents ... + (device, material, constants, algorithm, mobilityn, + mobilityp, V, n, p); + + res = resvec; + +endfunction + +function [Rn, Rp, Gn, Gp, II] = generation_recombination_model ... + (device, material, constants, algorithm, mobilityn, + mobilityp, V, n, p) + + [Rn_srh, Rp_srh, Gn_srh, Gp_srh] = secs1d_srh_recombination_noscale ... + (device, material, constants, algorithm, n, p); + + [Rn_aug, Rp_aug, G_aug] = secs1d_auger_recombination_noscale ... + (device, material, constants, algorithm, n, p); + + Rn = Rn_srh + Rn_aug; + Rp = Rp_srh + Rp_aug; + + Gp = Gn = Gn_srh + G_aug; + + Fp = V + constants.Vth * log (p ./ device.ni); + Fn = V - constants.Vth * log (n ./ device.ni); + + E = -diff (V) ./ diff (device.x); + + [Jn, Jp] = compute_currents ... + (device, material, constants, algorithm, mobilityn, + mobilityp, V, n, p); + + II = secs1d_impact_ionization_noscale ... + (device, material, constants, algorithm, + E, Jn, Jp, V, n, p, Fn, Fp); + +endfunction + + +function [res, jac] = residual_jacobian (device, material, constants, + algorithm, V, n, p, n0, p0, dt) + + %persistent mobilityn mobilityp Rn Rp Gn Gp II + [mobilityn, mobilityp] = compute_mobilities ... + (device, material, constants, algorithm, V, n, p); + + [Rn, Rp, Gn, Gp, II] = generation_recombination_model ... + (device, material, constants, algorithm, mobilityn, + mobilityp, V, n, p); + + nm = (n(2:end) + n(1:end-1)) / 2; + pm = (p(2:end) + p(1:end-1)) / 2; + + Nnodes = numel (device.x); + Nelements = Nnodes - 1; + + epsilon = material.esi * ones (Nelements, 1) / constants.q; + + A11 = bim1a_laplacian (device.x, epsilon, 1); + A12 = bim1a_reaction (device.x, 1, 1); + A13 = - A12; + R1 = A11 * V + bim1a_rhs (device.x, 1, n - p - device.D); + + A21 = - bim1a_laplacian (device.x, mobilityn .* nm, 1); + A22 = bim1a_advection_diffusion ... + (device.x, mobilityn * constants.Vth, 1, 1, V / constants.Vth) + ... + bim1a_reaction (device.x, 1, Rn + 0/dt); + A23 = bim1a_reaction (device.x, 1, Rp); + R2 = A22 * n + bim1a_rhs (device.x, 1, (Rn + 0/dt) .* n - (Gn + n0 * 0/ dt)); + + A31 = bim1a_laplacian (device.x, mobilityp .* pm, 1); + A32 = bim1a_reaction (device.x, 1, Rn); + A33 = bim1a_advection_diffusion ... + (device.x, mobilityp * constants.Vth , 1, 1, - V / constants.Vth) + ... + bim1a_reaction (device.x, 1, Rp + 0/dt); + R3 = A33 * p + bim1a_rhs (device.x, 1, (Rp + 0/dt) .* p - (Gp + p0 * 0/ dt)); + + N1 = algorithm.rowscaling(1); + N2 = algorithm.rowscaling(2); + N3 = algorithm.rowscaling(3); + + M1 = algorithm.colscaling(1); + M2 = algorithm.colscaling(2); + M3 = algorithm.colscaling(3); + + res = [R1(2:end-1)/N1; R2(2:end-1)/N2; R3(2:end-1)/N2]; + + jac = [[A11(2:end-1, 2:end-1)*M1, A12(2:end-1, 2:end-1)*M2, A13(2:end-1, 2:end-1)*M3]/N1; + [A21(2:end-1, 2:end-1)*M1, A22(2:end-1, 2:end-1)*M2, A23(2:end-1, 2:end-1)*M3]/N2; + [A31(2:end-1, 2:end-1)*M1, A32(2:end-1, 2:end-1)*M2, A33(2:end-1, 2:end-1)*M3]/N3]; + +endfunction + +function [Jn, Jp] = compute_currents (device, material, constants, algorithm, + mobilityn, mobilityp, V, n, p, Fn, Fp) + dV = diff (V); + dx = diff (device.x); + + [Bp, Bm] = bimu_bernoulli (dV ./ constants.Vth); + + Jn = constants.q * constants.Vth * mobilityn .* ... + (n(2:end) .* Bp - n(1:end-1) .* Bm) ./ dx; + + Jp = -constants.q * constants.Vth * mobilityp .* ... + (p(2:end) .* Bm - p(1:end-1) .* Bp) ./ dx; + +endfunction + +function [mobilityn, mobilityp] = compute_mobilities ... + (device, material, constants, algorithm, V, n, p) + + Fp = V + constants.Vth * log (p ./ device.ni); + Fn = V - constants.Vth * log (n ./ device.ni); + + E = -diff (V) ./ diff (device.x); + + mobilityn = secs1d_mobility_model_noscale ... + (device, material, constants, algorithm, E, V, n, p, Fn, Fp, 'n'); + + mobilityp = secs1d_mobility_model_noscale ... + (device, material, constants, algorithm, E, V, n, p, Fn, Fp, 'p'); + +endfunction + + +# function [V, n, p, Fn, Fp, Jn, Jp, res, it, nrfnp] = __one_gummel_step__ ... +# (device, material, constants, algorithm, V0, n0, p0, Fn0, +# Fp0, Jn0, Jp0, n_old, p_old, dt); + +# Fnref = Fn = Fn0; +# Fpref = Fp = Fp0; +# n = n0; +# p = p0; +# V = V0; +# Jn = Jn0; +# Jp = Jp0; + +# dx = diff (device.x); +# Nnodes = numel (device.x); +# Nelements = Nnodes -1; + +# for it = 1 : algorithm.maxit + +# [V, n, p] = secs1d_nlpoisson_newton_noscale ... +# (device, material, constants, algorithm, V, n, p, Fn, Fp); + +# [mobilityn, mobilityp] = compute_mobilities ... +# (device, material, constants, algorithm, V, n, p); + +# [Rn, Rp, Gn, Gp, II] = generation_recombination_model ... +# (device, material, constants, algorithm, mobilityn, +# mobilityp, V, n, p); + + +# An = bim1a_advection_diffusion (device.x, mobilityn, constants.Vth, +# 1, V / constants.Vth); +# An += bim1a_reaction (device.x, 1, Rn + 1/dt); +# R = bim1a_rhs (device.x, 1, Gn + n_old / dt) + ... +# bim1a_rhs (device.x, II, 1); + +# n(2:end-1) = An(2:end-1, 2:end-1) \ (R(2:end-1) - +# An(2:end-1, [1 end]) * n([1 end])); +# Fn = V - constants.Vth * log (n ./ device.ni); + +# Ap = bim1a_advection_diffusion (device.x, mobilityp, constants.Vth, +# 1, -V / constants.Vth); +# Ap += bim1a_reaction (device.x, 1, Rp + 1/dt); +# R = bim1a_rhs (device.x, 1, Gp + p_old / dt) + ... +# bim1a_rhs (device.x, II, 1); + + +# p(2:end-1) = Ap(2:end-1, 2:end-1) \ (R(2:end-1) - +# Ap(2:end-1, [1 end]) * p([1 end])); +# Fp = V + constants.Vth * log (p ./ device.ni); + +# [Jn, Jp] = compute_currents ... +# (device, material, constants, algorithm, mobilityn, mobilityp, +# V, n, p, Fn, Fp); + +# nrfn = norm (Fn - Fn0, inf); +# nrfp = norm (Fp - Fp0, inf); + +# 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 +# || (nrfnp > algorithm.maxnpincr)) +# break; +# endif + +# V0 = V; +# p0 = p ; +# n0 = n; +# Fn0 = Fn; +# Fp0 = Fp; + +# endfor + +# endfunction