%--------------------------------------------------------------------------
% variable_diffusion1D_mancvg_test.m
% 04/2018, David Shirokoff
%
% Description: Performs a manufactured convergence test in 1D 
% for ImEx coefficients and ImEx splitting parameters (delta, sigma).
% Test is for the variable coefficient diffusion equation 
% u_t = Au + Bu + f(x,t), 
% Au ~ sigma d^2u/dx^2,  Bu ~ d/dx(  (d(x) - sigma )du/dx)
% Au - implicit, Bu - explicit.
%
% Discretizations: Fourier spectral method.
%
% The code cycles through schemes with orders defined by the variable
% order_vect (default 1:5), and runs convergence tests varying the
% time-step del_t through the values defined by del_t_vect.
%
% No inputs. Output: convergence table to screen.
%--------------------------------------------------------------------------
function variable_diffusion1D_mancvg_test
    format long;    
                            %----------------------------------------------
    option = 1;             % choice of test solution (default option = 1)
                            % modify the following routines below to add
                            % other manufactured solution tests 
                            % (i) GetDiffusionCoeff, (ii) ExactSolution
                            % (iii) Forcing
                            %----------------------------------------------
    
                            %----------------------------------------------
                            % Variable Initializations
    n = 64;                 % number grid points.
    h = 1/n;                % grid size, domain = [0, 1]
    x = (0:n-1)*h;          % set grid
    tf = 5;                 % set final integration time
    order_vect = 1:5;       % convergence orders (default 1-5), change for different orders.
    num_order = size(order_vect, 2);
    
    delta = 0.1732;         % set ImEx deltaeter delta
       
    del_t_vect = 2.^(-9:-1:-11)';   % set different time steps
    num_dt_points = size(del_t_vect, 1);
                                
    err_array = zeros(num_dt_points, num_order);     % store errors to output at end.
    
    %Loop over schemes with different orders 
    for ord_counter = 1:num_order
        
        % get ImEx coefficients for given order
        order = order_vect(ord_counter);        
        [a_coeff b_coeff c_coeff] = GetImExCoefficients(order, delta);  
        
        % Loop over different del_t values
        for m = 1:num_dt_points
            del_t = del_t_vect(m);          
            num_steps = round(tf/del_t);
            true_tf = num_steps*del_t;
            fprintf('\nOrder = %f, \t dt = %f \n', ord_counter, del_t)
            
            % Call the ImEx integration: Integrate to time true_tf
            Uf = ImEx_TimeIntegrate(option, del_t, num_steps, x, a_coeff, b_coeff, c_coeff, n, h);        
            
            % Compute error with respect to exact solution
            Uexact = ExactSolution(option, x, true_tf);    
            err_array(m, ord_counter) = max(abs(Uf' - Uexact));            
        end
    end 
    
    % Compute instantaneous convergence rate
    cvgrate = zeros(num_dt_points, num_order);
    cvgrate(2:end, :) = log( err_array(1:end-1, :)./err_array(2:end, :) )/log(2);
    
    % Output errors and convergence to screen
    fprintf('\ndt\t\t\t\t')
    for jj = 1:num_order
        fprintf('err ord %d \t cvg rate \t', order_vect(jj) );
    end
    fprintf('\n')
    for ll = 1:num_dt_points
        fprintf('%2f \t & ', del_t_vect(ll));
        for jj = 1:num_order
            fprintf('\t  %1.1e \t & \t %1.1f &', err_array(ll, jj), cvgrate(ll, jj));
        end
        fprintf('\n');        
    end
  
end


%------------------------------------------------------------------
% Function: ImEx_TimeIntegrate
%
% Description: Uses ImEx multistep methods to integrate to the final 
% time.
%------------------------------------------------------------------
function Uf = ImEx_TimeIntegrate(option, del_t, num_steps, x, a, b, c, n, h)
        % Build Fourier wavenumbers. 
        ikx = zeros(n, 1);
        for ll = 1:n/2-1
            ikx(ll) = (ll-1)*i;   
            ikx(n - ll + 1) = -ll*i;
        end
        ikx(n/2) = (n/2 - 1)*i;
        ikx(n/2 + 1) = n/2*i;
        dw = 2*pi/1;            %Multiply by 2pi/domain length
        ikx = ikx*dw;         
        lap = ikx.^2;           %Laplacian -- Fourier space multipliers
        
        % Get the implicit and explicit diffusion coefficients. 
        [diff_coeffa, diff_coeffb] = GetDiffusionCoeff(option, x);
      
        % Initialize the ImEx scheme with the exact initial conditions.
        order = size(a, 2)-1;
        u_old = zeros(n, order);
        bu_old = zeros(n, order);
        au_old = zeros(n, order);
        f_old = zeros(n, order);
        for jj = 1:order
            u_old(:, jj) =  ExactSolution(option, x, -(jj-1)*del_t);
            bu_old(:,jj) = real(ifft(ikx.*fft(diff_coeffb.*ifft(ikx.*fft(u_old(:,jj)))) ));
            au_old(:,jj) = diff_coeffa*real(ifft(ikx.*fft(ifft(ikx.*fft(u_old(:,jj)))) ));
            f_old(:, jj) = Forcing(option, x, -(jj-1)*del_t);
        end    
        
        %Build the implicit vector. The implicit term is diagonal in
        %Fourier space with this vector being the diagonal entries.
        imp_part = a(1)*ones(n, 1) - diff_coeffa*del_t*c(1)*lap;
          
        %Main time-stepping loop. Loop over the total number of steps.
        for counter = 1:num_steps            
            
            % Build multistep right hand size and compute u_(n+1) from
            % previous u_n's.
            rhs = -u_old*a(2:end)' + del_t*( au_old*c(2:end)' + bu_old*b(2:end)' ) + del_t*f_old*b(2:end)';
            u = real(ifft(fft(rhs)./imp_part));
      
            % Update all temporary variables.
            for jj = 1:order-1
                u_old(:, order - jj + 1) = u_old(:, order - jj);
                bu_old(:, order - jj+1) = bu_old(:,order - jj);
                au_old(:, order - jj+1) = au_old(:,order - jj);
                f_old(:, order - jj+1) = f_old(:, order - jj);
            end
            bu_1 = real(ifft(ikx.*fft(diff_coeffb.*ifft(ikx.*fft(u))) ));
            au_1 = diff_coeffa*real(ifft(ikx.*fft(ifft(ikx.*fft(u))) ));
            
            u_old(:, 1) = u;
            bu_old(:, 1) = bu_1;
            au_old(:, 1) = au_1;
            f_old(:, 1) =  Forcing(option, x, counter*del_t);
        end        
        %return u at final time.
        Uf =  u;        
end


%------------------------------------------------------------------
% Function: GetDiffusionCoeff
%
% Description: Inputs the option number (option), and the spatial
% grid point vector x, and outputs both the implicit and explicit
% diffusion coefficients. 
%
% Inputs: option (default 1), x 
% Outputs: diff_coeffa -- Implicit diffusion coefficient (constant), 
% diff_coeffb -- Explicit diffusion coefficient (variable)
%------------------------------------------------------------------
function [diff_coeffa, diff_coeffb] = GetDiffusionCoeff(option, x)
    if (size(x, 2) > size(x, 1))
        x = x';
    end
    switch option
        case 1
            sigma = 2.69;           %Sigma is the ImEx splitting parameter
            diff_coeffa = sigma;
            diff_coeffb = 4 + 3*cos(2*pi*x) - sigma;
    end            
end

%------------------------------------------------------------------
% Function: ExactSolution
%
% Description: Inputs the option number (option), the spatial
% grid point vector x, and time t.  Outputs the exact solution vector 
% at grid points x and time t. 
%
% Inputs: option (default 1), x, t 
% Outputs: Uexact -- exact solution at grid points x and time t.
%------------------------------------------------------------------
function [Uexact] = ExactSolution(option, x, t)
    switch option
        case 1
            w = 20;
            Uexact = sin(w*t)*exp(sin(2*pi*x));
    end
end

%------------------------------------------------------------------
% Function: Forcing
%
% Description: Inputs the option number (option), the spatial
% grid point vector x, and time t.  Outputs the forcing required to 
% obtain the manufactured solution. 
%
% Inputs: option (default 1), x, t 
% Outputs: forcing -- forcing at grid points x and time t, required
% for the exact solution with the same option number.
%------------------------------------------------------------------
function [forcing] = Forcing(option, x, t)
    switch option
        case 1
            w = 20;
            forcing = w*cos(w*t)*exp(sin(2*pi*x)) - ...
                sin(w*t)*( 4*pi^2*exp(sin(2*pi*x)).*cos(2*pi*x).^2.*(3*cos(2*pi*x) + 4) - ...
                12*pi^2*exp(sin(2*pi*x)).*cos(2*pi*x).*sin(2*pi*x) - ...
                4*pi^2*exp(sin(2*pi*x)).*sin(2*pi*x).*(3*cos(2*pi*x) + 4) );       
    end    
end


%------------------------------------------------------------------
% Function: GetImExCoefficients
%
% Description: Inputs the scheme order and ImEx parameter delta, and 
% returns the ImEx time-stepping coefficients. For order r, there are
% r + 1 coefficients (for each, a_coeff, b_coeff, and c_coeff). 
%
% Inputs: order (integer 1-5), delta (float, 0 < delta <= 1)
% Outputs: ImEx coefficients, a_coeff, b_coeff, c_coeff
%
%------------------------------------------------------------------
function [a_coeff b_coeff c_coeff] = GetImExCoefficients(order, delta)
    switch order
        
        case 1
            a_coeff = [delta -delta];
            b_coeff = [0 delta];
            c_coeff = [1 delta-1];

        case 2
            a_coeff = [(2*delta - delta^2/2) (2*delta^2 - 4*delta) (2*delta-(3*delta^2)/2)];
            b_coeff = [0 2*delta (delta - 1)^2 - 1];
            c_coeff = [1 (2*delta - 2) (delta - 1)^2];
            
        case 3
            a_coeff =[(delta^3/3 - (3*delta^2)/2 + 3*delta), -((3*delta^3)/2 - (15*delta^2)/2 + 9*delta), ...
                (3*delta^3 - (21*delta^2)/2 + 9*delta), (-3*delta + (9*delta^2)/2 - (11*delta^3)/6) ];
            b_coeff =[ 0, 3*delta,  ((delta - 1)^2 + (2*delta - 2)*(delta - 1) - 3), (delta - 1)^3 + 1];
            c_coeff = [1, (3*delta - 3), ((delta - 1)^2 + (2*delta - 2)*(delta - 1)), (delta - 1)^3];
            
        case 4 
            a_coeff = [(- delta^4/4 + (4*delta^3)/3 - 3*delta^2 + 4*delta), -(-(4*delta^4)/3 + (22*delta^3)/3 - ...
                18*delta^2 + 16*delta), (- 3*delta^4 + 18*delta^3 - 36*delta^2 + 24*delta), ...
                -(- 4*delta^4 + (58*delta^3)/3 - 30*delta^2 + 16*delta), 4*delta - 9*delta^2 + (22*delta^3)/3 - (25*delta^4)/12 ];
            b_coeff = [ 0, 4*delta, ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) - 6),  ...
                ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + 4), (delta - 1)^4  - 1];
            c_coeff = [1, (4*delta - 4), ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1)), ...
                ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3), (delta - 1)^4]; 
            
        case 5
           a_coeff =[(delta^5/5 - (5*delta^4)/4 + (10*delta^3)/3 - 5*delta^2 + 5*delta), -((5*delta^5)/4 - ...
               (95*delta^4)/12 + (65*delta^3)/3 - 35*delta^2 + 25*delta), ((10*delta^5)/3 - (65*delta^4)/3 + (190*delta^3)/3 ...
               - 90*delta^2 + 50*delta), -(5*delta^5 - 35*delta^4 + (280*delta^3)/3 - 110*delta^2 + 50*delta), ...
               (5*delta^5 - (365*delta^4)/12 + (200*delta^3)/3 - 65*delta^2 + 25*delta), 15*delta^2 - (55*delta^3)/3 ...
               + (125*delta^4)/12 - (137*delta^5)/60 - 5*delta];
           b_coeff = [0, 5*delta, ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) + ...
               (4*delta - 4)*(delta - 1) - 10),  ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + ...
               (delta - 1)^3 + (delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1)) + 10), ...
               ((delta - 1)^4 + ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3)*(delta - 1) - 5), 1 + (delta - 1)^5];
           c_coeff = [1, (5*delta - 5), ((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1) + ...
               (4*delta - 4)*(delta - 1)), ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3 + ...
               (delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1) + (3*delta - 3)*(delta - 1))), ((delta - 1)^4 + ...
               ((delta - 1)*((delta - 1)^2 + (2*delta - 2)*(delta - 1)) + (delta - 1)^3)*(delta - 1)), (delta - 1)^5];           
    end
    
end