%TEMPLE8200_AP_METHODS
%   Test of two numerical methods for the system
%   (P_1 model of radiation transport)
%      u_t + v_x = 0
%      v_t + u_x = -sigma*v
%   The hyperbolic part is approximated (A) via Godunov's
%   method, and (B) via a staggered grid method. In both
%   cases, the relaxation term if approximated via a backward
%   Euler step. Both methods are tested in the underresolved
%   regime when h and dt are both larger than 1/sigma. The
%   result are that method (B) is asymptotic preserving,
%   while method (A) is not.
%
%   (C) 2024/11/17 by Benjamin Seibold
%   http://www.math.temple.edu/~seibold
%   
%   This code may be used and modified for teaching,
%   learning, and research, given its orgin is credited.

% Parameters
ax = [0,1]; % spatial domain
tf = 50; % final time
h = 1e-2; % spatial grid resolution
dt = .9e-2; % time step size
sigma = 1e4; % decay parameter of v component
u0 = @(x) exp(-10*cos(pi*x).^2); v0 = @(x) x*0; % initial conditions

% Initialization
nt = ceil(tf/dt); dt = tf/nt; % actual time step
L = diff(ax); % domain size
n = ceil(L/h); h = L/n; % actual mesh size
x1 = linspace(ax(1)+h/2,ax(2)-h/2,n); % cell center grid
x2 = x1+h/2; % right-stagged grid
U = cell(1,2); V = cell(1,2);
name{1} = 'Godunov method';
U{1} = u0(x1); V{1} = v0(x1); % evaluate i.c. for Godunov method
name{2} = 'Staggered grid method';
U{2} = u0(x1); V{2} = v0(x2); % evaluate i.c. for staggered grid method
U_d = u0(x1); % evaluate i.c. for diffusion approximation

% Run simulation
for n = 0:nt % time loop
    if ceil((n+1)/100)>ceil(n/100)||n==nt % plot only every few steps
        % Plotting
        clf
        for s = 1:length(U) % Loop over methods
            subplot(length(U),1,s)
            plot(x1,U_d,'LineWidth',3,'Color',[.7 .7 1]), hold on
            plot(x1,U{s},'LineWidth',1,'Color',[0 0 1]), hold off
            axis([ax,-.1,1])
            legend('reference (diffusion approx.)','direct method')
            xlabel('x'), ylabel('u')
            title(sprintf('%s at t = %0.2f',name{s},n*dt))
        end
        drawnow
    end

    % Update step Godunov method
    F = (U{1}-U{1}([2:end,1]))/2 + (V{1}+V{1}([2:end,1]))/2;
    G = (V{1}-V{1}([2:end,1]))/2 + (U{1}+U{1}([2:end,1]))/2;
    U{1} = U{1}-dt*(F-F([end,1:end-1]))/h;
    V{1} = V{1}-dt*(G-G([end,1:end-1]))/h-dt*sigma*V{1};
    V{1} = V{1}/(1+dt*sigma); % implicit relaxation

    % Update step staggered grid method
    U{2} = U{2}-dt*(V{2}-V{2}([end,1:end-1]))/h;
    V{2} = V{2}-dt*(U{2}([2:end,1])-U{2})/h;
    V{2} = V{2}/(1+dt*sigma); % implicit relaxation

    % Diffusion approximation
    mu = dt/h.^2/sigma;
    U_d = mu*U_d([end,1:end-1])+(1-2*mu)*U_d+mu*U_d([2:end,1]);
end
