Kulikovsky_SG_ISG0

%% Example demonstrating Scharfetter-Gummel (SG) and Improved SG methods 
%% for 1D convection-diffusion equation from paper J. Comp. Phys. 119, pp.
%% 149-155 (1995) by A. A. Kulikovsky + Method of Lines (MOL) on
%% non-uniform grid
function Kulikovsky_SG_ISG0
close all; clc;
global x

% Number of non-uniform grid points
N = 200; 

% Define non-uniform refined at x = 0 smooth grid
xis = 1.0;
x = 0 + (1 - 0)*(exp(linspace(0.0, 1.0, N)*xis)-1)./(exp(xis)-1);

% Initial conditions (see (35) or (36))
% u0 = 1.0e+2 + 0.5*(1 + tanh((x-0.8)/0.02))*1.0e+12;
u0 = 1.0e+2 + 1.0e+12*exp(-(x-0.6).^2/0.004);

% Uniform time interval
t = linspace(0.0, 4.0e-5, 100);

%% Solve MOL ODEs for convection-diffusio equation with SG and ISG-0 fluxes
disp('Computing...');
tic; 
[ZZ, u1] = ode45(@MOL_eqs_SG, t, u0);
[ZZ, u2] = ode45(@MOL_eqs_ISG, t, u0); 
toc
disp('Complete!');

%% Plot solution in dynamics
for I=1:100
    semilogy(x, u0, 'b-', x, u1(I,:), 'r-', x, u2(I,:), 'g-', 'LineWidth', 2);
    legend('Initial', 'MOL + SG', 'MOL + ISG-0', 'Location', 'NorthWest');
    xlim([0.0 1.0]);
    grid on;
    pause(0.01)
    drawnow;
end

%% MOL ODEs system with SG flux for equation 
%% convection-diffusion equation in 1D (32)
function [res] = MOL_eqs_SG(t, u)
global x
N = length(u); 
% Constants from original paper
mue = 1.0; De = 1.0; A = 1.0e+4;
% Linear field
E = A*x;

% Define spatial grid steps for non-uniform grid
h(1:N-1) = x(2:N)-x(1:(N-1)); h(N) = h(N-1); h(N+1) = h(N);

for k=1:(N-1)
    % Restrict averaged E_{k+1/2} to prevent run-time errors
    if 0.5*abs(E(k+1)+E(k))> eps
        alpha = mue*h(k)/De*0.5*(E(k+1)+E(k));
        I0 = ( exp(alpha) - 1.0 )/alpha;
        j(k+1) = De/h(k)*(u(k) - exp(alpha)*u(k+1))/I0;
    else
        j(k+1) = De/h(k)*(u(k) - u(k+1));
    end
end
% Apply simple BCs for fluxes at ghost nodes
j(1) = -A*x(1)*u(1);   j(N+1) = -A*x(N)*u(N);

% Form MOL ODEs
k = 1:N;
res(k) = -(j(k+1)-j(k))./(0.5*h(k)+0.5*h(k+1));
res = res';

%% MOL ODEs system with ISG-0 flux for equation 
%% convection-diffusion equation in 1D (32)
function [res] = MOL_eqs_ISG(t,u)
global x
N = length(u); 
% Constants from original paper
mue = 1.0; De = 1.0; A = 1.0e+4; epsilon = 0.01;
% Linear field
E = A*x;

h(1:N-1) = x(2:N)-x(1:(N-1)); h(N) = h(N-1); h(N+1) = h(N);

for k=1:(N-1)
    % Distance between virtual nodes h_v (see (21))
    h_v = sqrt(2.0*epsilon*De*h(k)/mue/abs(E(k+1)-E(k)));
    % If h_v>=h then use SG scheme
    if (h_v >= h(k))        
        if 0.5*abs(E(k+1)+E(k))> eps
            alpha = mue*h(k)/De*0.5*(E(k+1)+E(k));
            I0 = ( exp(alpha) - 1.0 )/alpha;
            j(k+1) = De/h(k)*(u(k) - exp(alpha)*u(k+1))/I0;
        else
            j(k+1) = De/h(k)*(u(k) - u(k+1));
        end
    % If h_v<h then use ISG-0 (24)-(25)   
    else
        if 0.5*abs(E(k+1)+E(k))> eps
            x_L = (x(k)+x(k+1)-h_v)/2; 
            x_R = (x(k)+x(k+1)+h_v)/2;
            alpha_v = h_v*mue/De*0.5*(E(k+1)+E(k));
            I0 = ( exp(alpha_v) - 1.0 )/alpha_v;
            j(k+1) = De/h_v*(interp1(x, u, x_L, 'pchip') - exp(alpha_v)*interp1(x, u, x_R, 'pchip'))/I0;
        else
            j(k+1) = De/h_v*(interp1(x, u, x_L, 'pchip') - interp1(x, u, x_R, 'pchip'));            
        end
    end
end
% Apply simple BCs for fluxes at ghost nodes
j(1) = -A*x(1)*u(1);   j(N+1) = -A*x(N)*u(N);

% Form MOL ODEs
k = 1:N;
res(k) = -(j(k+1)-j(k))./(0.5*h(k)+0.5*h(k+1));
res = res';