%%  Finite differences for the advection-diffusion-reaction equation
%
%     u_t   =   v u_x    + D u_xx    + A (1-u) exp(-E/(u-p))
%
%    on the domain x=[0,1],
%    with boundary conditions    u(0)=1,   u(1)=0.
% 
%    This equation models combustion (reactions) unless A=0.
%
%    The equation is discretized spatially by finite differences.
%    
%    When A=0 (no reaction terms), the linear system at steady state 
%    is solved directly, via "Backslash". It is also sovled via an ODE
%    solver (Method 1). 
%
%    When A>0 the nonlinear system at steady state is solved via
%    an ODE solver (Method 1) and via two iterative methods:
%
%    Method 1: The resulting system of ODEs is solved by MATLAB ode15s to reach steady state.
%
%    Method 2: Newton iteration 
%    
%    Method 3: Fixed point iteration
%              
%


clear
clc
close all

%% parameters
t=2;
D=0.1; v=10; % diffusion and advection parameters
% reaction parameters: A, E, p
% experiment with parameter A=0 or A>0 
A=2e6  
%A=0 % A=0 so no reactions
E=8; p=-0.1;
params.A=A; params.E=E; params.p=p;

% Number of points in discretization
xL=0; xR=1;
N=20+1; % N=20+1;
h=(xR-xL)/(N-1); % spatial discretization size
x=xL:h:xR; x=x';

u0=1-erf(8*(x-0.7)); u0=0.5*u0; % initial values in time-dependent case
% u0=1-x; % try experimenting with different initial values
u0(1)=1; u0(end)=0; % boundary conditions


% diffusion matrix
e = ones(N,1);
K = spdiags([e -2*e e], -1:1, N, N);
%K=full(K)  % full matrix instead of sparse matrix

% first order upwind difference matrix
DelPlus = spdiags([-e e], 0:1, N, N);
%DelPlus=full(DelPlus) % full matrix instead of sparse matrix

vconst=true;
% we may change the advection term so that it depends on space, x
% only do this if you are interested; it is NOT part of the homework.
% the following line replaces constant parameter v, with v(x)=v*x
% xd=spdiags(x, 0, N, N); DelPlus=xd*DelPlus; vconst=false; 

DelPlus=v*DelPlus/h;
K=D*K/h^2;
AD=DelPlus+K;            % advection-diffusion matrix
AD(1,:)=0; AD(end,:)=0;  % boundary conditions

%% ODE solver
tol=1e-8;
initial_state=u0;
options = odeset('RelTol',tol,'AbsTol',tol*ones(size(initial_state)));
disp('Solving with the MATLAB code ode15s...')
tic;
t1=min(0.1,t/10); t2=min(0.15,t/5); tspan=[0 t1 t2 t];
[T,U] = ode15s(@(t, u) adr_dt(t,u,params,AD), tspan, initial_state, options);
ODEsolverTime=toc
ut=U(end,:);

if (A==0)
    %% Linear solve
    % if A==0, we can solve for the steady state, u_ss, with backslash:
    b=zeros(size(x));
    AD(1,1)=1; AD(end,end)=1; b(1)=1; % boundary conditions
    tic;
    u_ss=AD\b;
    LinearSolverTime=toc
else
    %% Nonlinear solve 
    % if A is nonzero we need a nonlinear solver
    % try experimenting with different initial guesses, u_guess
    u_guess=1-x;  
    % u_guess=u0;
    sstol=1e-10;  % tolerance for iterative solver for steady state
    tic
    [u_ssN errorsN iterationsN]= NewtonIteration(u_guess, AD, params, sstol);
    NewtonIterationTime=toc
    
    % fixed point iteration
    tic
    [u_ssFP errorsFP iterationsFP]=FPIteration(u_guess, AD, params, sstol);
    FixedPointIterationTime=toc
    
    % Newton iteration, with initial guess from fixed point iteration
    tic
    [u_ssFPcrude errorsFPcrude iterationsFPcrude]=FPIteration(u_guess, AD, params, 20);
    [u_ssNFP errorsNFP iterationsNFP]= NewtonIteration(u_ssFPcrude, AD, params, sstol);   
    NewtonIterationPreCondTime=toc
    
    % choose the most accurate of the three approximations to the steady
    % state from the three nonlinear solvers
    es=[errorsN(end) errorsFP(end) errorsNFP(end)]; [emin  imin]=min(es);
    if imin==1
        u_ss=u_ssN;
    else if imin==2
            u_ss=u_ssFP;
        else
            u_ss=u_ssNFP;
        end
    end
end


%% plot results
disp('Plotting solution...')
clf('reset')
FS='FontSize';
figure(1)
plot(x,U, x, u_ss, '--', 'LineWidth', 3)
if (A==0)
    if (vconst)
        title({['u_t  =  D u_{xx}   + v u_x  '],  ['   D=', num2str(D), '   v=', num2str(v), '         h=', num2str(h)]}, FS, 16)
    else
        title({['u_t  =  D u_{xx}   + v x u_x  '],  ['   D=', num2str(D), '   v=', num2str(v), '         h=', num2str(h)]}, FS, 16)
    end
    legend('initial condition  u(0)', ['u(t=', num2str(t1), ')'], ['u(t=', num2str(t2), ')'], ['u(t=', num2str(t), ')'], 'Steady state, u_{ss}, via linear solver ')
    xlabel('x', FS, 16)
    ylabel('u', FS, 16)
else
    if (vconst)
        title({['u_t  =  D u_{xx}   + v u_x  +   A (1-u) exp(-E/(u-p))'],  ['   D=', num2str(D), '   v=', num2str(v), '         h=', num2str(h)]}, FS, 16)
    else
        title({['u_t  =  D u_{xx}   + v x u_x  +   A (1-u) exp(-E/(u-p)) '],  ['   D=', num2str(D), '   v=', num2str(v), '         h=', num2str(h)]}, FS, 16)
    end
    legend('initial condition  u(0)', ['u(t=', num2str(t1), ')'], ['u(t=', num2str(t2), ')'], ['u(t=', num2str(t), ')'], 'Steady state, u_{ss}, via nonlinear solver ')
    xlabel('x', FS, 16)
    ylabel('u', FS, 16)
    
    
    figure(2)
    semilogy(1:iterationsFP,errorsFP, 'x-', 1:iterationsN,errorsN,'*--', 1:iterationsNFP,errorsNFP,'o--', 'LineWidth', 2)
    xlabel('Iteration', FS, 16)
    ylabel('Error', FS, 16)
    legend('Fixed point iteration', 'Newton iteration', 'Newton iteration, with initial guess from fixed point iteration')
    if (vconst)
        title({['0  =  D u_{xx}   + v u_x  +   A (1-u) exp(-E/(u-p))'],  ['   D=', num2str(D), '   v=', num2str(v), '         h=', num2str(h)]}, FS, 16)
    else
        title({['0  =  D u_{xx}   + v x u_x  +   A (1-u) exp(-E/(u-p)) '],  ['   D=', num2str(D), '   v=', num2str(v), '         h=', num2str(h)]}, FS, 16)
    end    
end


disp('Finished')

