Mercurial > forge
changeset 11590:3a8c8c109b36 octave-forge
fixed bug in newton solver
| author | cdf |
|---|---|
| date | Tue, 02 Apr 2013 12:38:55 +0000 |
| parents | ccd2ae2dc974 |
| children | 4599bf7580aa |
| files | extra/secs1d/inst/demos/diode_tran_noscale.m extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m extra/secs1d/inst/secs1d_tran_dd_newton_noscale.m |
| diffstat | 3 files changed, 117 insertions(+), 132 deletions(-) [+] |
line wrap: on
line diff
--- a/extra/secs1d/inst/demos/diode_tran_noscale.m Mon Apr 01 16:59:36 2013 +0000 +++ b/extra/secs1d/inst/demos/diode_tran_noscale.m Tue Apr 02 12:38:55 2013 +0000 @@ -61,13 +61,15 @@ vbcs = {@vbcs_1, @vbcs_2}; % tolerances for convergence checks -algorithm.toll = 1e-6; -algorithm.maxit = 100; -algorithm.ptoll = 1e-12; +algorithm.toll = 1e-10; +algorithm.ltol = 1e-10; +algorithm.maxit = 100; +algorithm.lmaxit = 100; +algorithm.ptoll = 1e-12; algorithm.pmaxit = 1000; -algorithm.colscaling = [10 1e23 1e23]; -algorithm.rowscaling = [1e6 1e23 1e23]; -algorithm.maxnpincr = constants.Vth; +algorithm.colscaling = [1 1e23 1e25]; +algorithm.rowscaling = [1 1e7 1e7]; +algorithm.maxnpincr = constants.Vth; %% initial guess via stationary simulation [nin, pin, Vin, Fnin, Fpin, Jn, Jp, it, res] = secs1d_dd_gummel_map_noscale ... @@ -76,13 +78,9 @@ %% 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_newton_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); @@ -102,11 +100,7 @@ vvector = Fn(end, :); ivector = (Jn(end, :) + Jp(end, :)); -ivectorn = (Jn(1, :) + Jp(1, :)); -ivectora = mean (Jn + Jp, 1); -area = 150e-6 * 50e-6; -figure (1) -plotyy (t, vvector, t, area*ivector) -legend('J_L','J_0') + +plotyy (t, vvector, t, device.W*ivector) drawnow
--- a/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m Mon Apr 01 16:59:36 2013 +0000 +++ b/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m Tue Apr 02 12:38:55 2013 +0000 @@ -291,6 +291,10 @@ break; endif + figure (3) + semilogy (1:it, res) + drawnow + V0 = V; p0 = p ; n0 = n;
--- a/extra/secs1d/inst/secs1d_tran_dd_newton_noscale.m Mon Apr 01 16:59:36 2013 +0000 +++ b/extra/secs1d/inst/secs1d_tran_dd_newton_noscale.m Tue Apr 02 12:38:55 2013 +0000 @@ -38,8 +38,8 @@ %% 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 +%% algorithm.toll tolerance for Newton iterarion convergence test +%% algorithm.maxit maximum number of Newton iterarions %% %% output: %% n electron concentration @@ -59,7 +59,6 @@ rejected = 0; V = Vin; n = nin; p = pin; Fn = Fnin; Fp = Fpin; - Nnodes = numel (device.x); Nelements = Nnodes -1; @@ -72,7 +71,7 @@ tmax = tspan(2); tmin = tspan(1); - dt = (tmax - tmin) / 2000; %% initial time step + dt = (tmax - tmin) / 200; %% initial time step tstep = 1; t = tout(tstep) = tmin; @@ -109,55 +108,42 @@ nrfnp = 4 * algorithm.maxnpincr; %% initialize so it is not undefined in case of exception + # [V_gummel, n_gummel, p_gummel, Fn_gummel, ... + # Fp_gummel, Jn_gummel, Jp_gummel, Rn, Rp, Gn, Gp, II, mun, mup, ~, ~, nrfnp] =... + # __one_gummel_step__ (device, material, constants, algorithm, + # V0, n0, p0, Fn0, Fp0, Jn(:,tstep), Jp(:,tstep), + # n(:, tstep - 1), p(:, tstep - 1), dt, 1); + + + [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) + + 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); + res(tstep,1:it) = res_; + printf ("dt incremented by %g\n", .5 * sqrt (algorithm.maxnpincr / nrfnp)) + dt *= .5 * 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)) + algorithm.maxnpincr, nrfnp, dt, dt * .5 * sqrt (algorithm.maxnpincr / nrfnp)) - dt *= .95 * sqrt (algorithm.maxnpincr / nrfnp); + dt *= .5 * sqrt (algorithm.maxnpincr / nrfnp); t = tout(--tstep); if (dt < 100*eps) warning ("secs1d_tran_dd_newton_noscale: minimum time step reached.") - return endif rejected ++ ; @@ -179,17 +165,14 @@ 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); + algorithm, V, n, p, n_old, p_old, + dt); for it = 1 : algorithm.maxit @@ -203,32 +186,33 @@ ndn = norm (dn, inf); ndp = norm (dp, inf); - tkv = 1; + tkv = 1; if (ndv > constants.Vth) tkv = constants.Vth / ndv; endif tkn = tkv; - [howmuch, where] = min (n(2:end-1) + tkn * dn); + [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); + [howmuch, where] = min (p(2:end-1) + tkp * dp); dp(where) if (howmuch <= 0) tkp = - .9 * p(2:end-1)(where) / dp(where); endif - tk = min ([tkv, tkn, tkp]) - + 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); + p_old, % Rn, Rp, Gn, Gp, II, mun, mup, + dt); resvec(it) = reldnorm = ... ndv / norm (V, inf) + ... @@ -236,7 +220,7 @@ ndp / norm (p, inf); figure (3) - semilogy (1:it, resvec) + semilogy (1:it, resvec); drawnow Fp = V + constants.Vth * log (p ./ device.ni); @@ -251,7 +235,7 @@ endfor [mobilityn, mobilityp] = compute_mobilities ... - (device, material, constants, algorithm, V, n, p); + (device, material, constants, algorithm, V, n, p); [Jn, Jp] = compute_currents ... (device, material, constants, algorithm, mobilityn, @@ -293,7 +277,8 @@ function [res, jac] = residual_jacobian (device, material, constants, - algorithm, V, n, p, n0, p0, dt) + algorithm, V, n, p, n0, p0, % Rn, Rp, Gn, Gp, II, mobilityn, mobilityp, + dt) %persistent mobilityn mobilityp Rn Rp Gn Gp II [mobilityn, mobilityp] = compute_mobilities ... @@ -303,8 +288,8 @@ (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; + nm = 2./(1./n(2:end) + 1./n(1:end-1)); + pm = 2./(1./p(2:end) + 1./p(1:end-1)); Nnodes = numel (device.x); Nelements = Nnodes - 1; @@ -317,18 +302,20 @@ 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); + A22_stiff = bim1a_advection_diffusion ... + (device.x, mobilityn * constants.Vth, 1, 1, V / constants.Vth); + A22 = A22_stiff + ... + bim1a_reaction (device.x, 1, Rn + 1/dt); A23 = bim1a_reaction (device.x, 1, Rp); - R2 = A22 * n + bim1a_rhs (device.x, 1, (Rn + 0/dt) .* n - (Gn + n0 * 0/ dt)); + R2 = A22_stiff * n + bim1a_rhs (device.x, 1, (Rn + 1/dt) .* n - (Gn + n0 * 1/ 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)); + A33_stiff = bim1a_advection_diffusion ... + (device.x, mobilityp * constants.Vth , 1, 1, - V / constants.Vth); + A33 = A33_stiff + ... + bim1a_reaction (device.x, 1, Rp + 1/dt); + R3 = A33_stiff * p + bim1a_rhs (device.x, 1, (Rp + 1/dt) .* p - (Gp + p0 * 1/ dt)); N1 = algorithm.rowscaling(1); N2 = algorithm.rowscaling(2); @@ -378,79 +365,79 @@ 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); +function [V, n, p, Fn, Fp, Jn, Jp, Rn, Rp, Gn, Gp, II, mobilityn, mobilityp, res, it, nrfnp] = __one_gummel_step__ ... + (device, material, constants, algorithm, V0, n0, p0, Fn0, + Fp0, Jn0, Jp0, n_old, p_old, dt, NIT); -# Fnref = Fn = Fn0; -# Fpref = Fp = Fp0; -# n = n0; -# p = p0; -# V = V0; -# Jn = Jn0; -# Jp = Jp0; + 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; + dx = diff (device.x); + Nnodes = numel (device.x); + Nelements = Nnodes -1; -# for it = 1 : algorithm.maxit + for it = 1 : NIT -# [V, n, p] = secs1d_nlpoisson_newton_noscale ... -# (device, material, constants, algorithm, V, n, p, Fn, Fp); + [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); + [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); + [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); + 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); + 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); + 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); + 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); + [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); + nrfn = norm (Fn - Fn0, inf); + nrfp = norm (Fp - Fp0, inf); -# nrfnp = max (norm (Fn - Fnref, inf), norm (Fp - Fpref, inf)); + nrfnp = max (norm (Fn - Fnref, inf), norm (Fp - Fpref, inf)); -# nrv = norm (V - V0, inf); -# res(it) = max ([nrfn; nrfp; nrv]); + nrv = norm (V - V0, inf); + res(it) = max ([nrfn; nrfp; nrv]); -# if (res(it) < algorithm.toll -# || (nrfnp > algorithm.maxnpincr)) -# break; -# endif + if (res(it) < algorithm.toll + || (nrfnp > algorithm.maxnpincr)) + break; + endif -# V0 = V; -# p0 = p ; -# n0 = n; -# Fn0 = Fn; -# Fp0 = Fp; + V0 = V; + p0 = p ; + n0 = n; + Fn0 = Fn; + Fp0 = Fp; -# endfor + endfor -# endfunction +endfunction