# HG changeset patch # User cdf # Date 1364136802 0 # Node ID 35c57f32488aeaac5ea6f51b1efd2fc5d35758f9 # Parent 332c5ee57ec0987f407355c3d0066dd2ef626e9a fix generation recombination and scaling in stationary Newton diff -r 332c5ee57ec0 -r 35c57f32488a extra/secs1d/inst/demos/diode_reverse_noscale.m --- a/extra/secs1d/inst/demos/diode_reverse_noscale.m Fri Mar 22 09:12:39 2013 +0000 +++ b/extra/secs1d/inst/demos/diode_reverse_noscale.m Sun Mar 24 14:53:22 2013 +0000 @@ -25,7 +25,7 @@ device.tn = secs1d_carrier_lifetime_noscale (device.Na, device.Nd, 'n'); % initial guess for n, p, V, phin, phip -V_p = -1; +V_p = 1; V_n = 0; Fp = V_p * (device.x <= xm); @@ -45,8 +45,8 @@ V = Fn + constants.Vth * log (n ./ device.ni); % tolerances for convergence checks -algorithm.toll = 1e-3; -algorithm.maxit = 1000; +algorithm.toll = 1e-10; +algorithm.maxit = 100; algorithm.ptoll = 1e-12; algorithm.pmaxit = 1000; @@ -56,6 +56,7 @@ constants, algorithm, V, n, p, Fn, Fp); +algorithm.maxit = 1000; [n_newt, p_newt, V_newt, Fn_newt, Fp_newt, Jn, Jp, it, res] = secs1d_dd_newton_noscale ... (device, material, constants, algorithm, V, n, p, Fn, Fp); @@ -65,6 +66,9 @@ Ec = constants.Vth * log (material.Nc./n) + Efn; Ev = -constants.Vth * log (material.Nv./p) + Efp; -plot (device.x, Efn, device.x, Efp, device.x, Ec, device.x, Ev) +plot (device.x, Efn, + device.x, Efp, + device.x, Ec, + device.x, Ev) legend ('Efn', 'Efp', 'Ec', 'Ev') axis tight diff -r 332c5ee57ec0 -r 35c57f32488a extra/secs1d/inst/secs1d_dd_newton_noscale.m --- a/extra/secs1d/inst/secs1d_dd_newton_noscale.m Fri Mar 22 09:12:39 2013 +0000 +++ b/extra/secs1d/inst/secs1d_dd_newton_noscale.m Sun Mar 24 14:53:22 2013 +0000 @@ -56,13 +56,11 @@ dampcoeff = 2; Nnodes = numel (device.x); - V = Vin; - n = nin; - p = pin; + V = Vin; n = nin; p = pin; [res, jac] = residual_jacobian (device, material, constants, algorithm, V, n, p); - keyboard + normr(1) = norm (res, inf); normrnew = normr(1); @@ -86,6 +84,7 @@ algorithm, Vnew, nnew, pnew); normrnew = norm (resnew, inf); + if (normrnew > normr(it)) tk = tk / dampcoeff; else @@ -99,7 +98,18 @@ V = Vnew; n = nnew; p = pnew; normr(it+1) = normrnew; - if (normr(it+1) <= algorithm.toll) + 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); + + plotyy ( + 1:it, log10 (resvec), + 1:it+1, log10 (normr) + ) + drawnow + + if (reldnorm <= algorithm.toll) break endif @@ -115,6 +125,8 @@ Fp = V + constants.Vth * log (p ./ device.ni); Fn = V - constants.Vth * log (n ./ device.ni); + res = resvec; + endfunction @@ -152,46 +164,50 @@ function [res, jac] = residual_jacobian (device, material, constants, algorithm, V, n, p) - [mobilityn, mobilityp] = compute_mobilities ... + [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; + 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); - + epsilon = material.esi * ones (Nelements, 1) / constants.q; + A11 = bim1a_laplacian (device.x, epsilon, 1); - A12 = bim1a_reaction (device.x, constants.q, 1); + A12 = bim1a_reaction (device.x, 1, 1); A13 = - A12; - R1 = A11 * V + bim1a_rhs (device.x, constants.q, n - p - device.D); + 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); A23 = bim1a_reaction (device.x, 1, Rp); - R2 = A22 * n + bim1a_rhs (device.x, 1, Gn - Rn .* n); + R2 = A22 * n + bim1a_rhs (device.x, 1, Rn .* n - Gn); 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); - R3 = A33 * p + bim1a_rhs (device.x, 1, Gp - Rp .* p); + R3 = A33 * p + bim1a_rhs (device.x, 1, Rp .* p - Gp); - res = [R1(2:end-1); R2(2:end-1); R3(2:end-1)]; + N1 = norm(R1(2:end-1), inf); + N2 = norm(R2(2:end-1), inf); + N3 = norm(R3(2:end-1), inf); - jac = [A11(2:end-1, 2:end-1), A12(2:end-1, 2:end-1), A13(2:end-1, 2:end-1); - A21(2:end-1, 2:end-1), A22(2:end-1, 2:end-1), A23(2:end-1, 2:end-1); - A31(2:end-1, 2:end-1), A32(2:end-1, 2:end-1), A33(2:end-1, 2:end-1)]; + res = [R1(2:end-1)/N1; R2(2:end-1)/N2; R3(2:end-1)/N2]; + jac = [[A11(2:end-1, 2:end-1), A12(2:end-1, 2:end-1), A13(2:end-1, 2:end-1)]/N1; + [A21(2:end-1, 2:end-1), A22(2:end-1, 2:end-1), A23(2:end-1, 2:end-1)]/N2; + [A31(2:end-1, 2:end-1), A32(2:end-1, 2:end-1), A33(2:end-1, 2:end-1)]/N3]; + endfunction function [Jn, Jp] = compute_currents (device, material, constants, algorithm, diff -r 332c5ee57ec0 -r 35c57f32488a extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m --- a/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m Fri Mar 22 09:12:39 2013 +0000 +++ b/extra/secs1d/inst/secs1d_tran_dd_gummel_map_noscale.m Sun Mar 24 14:53:22 2013 +0000 @@ -157,7 +157,7 @@ end_try_catch endwhile - rejected + printf ("number of rejected time steps: %d\n", rejected) endfunction function [Jn, Jp] = compute_currents (device, material, constants, algorithm,