%-----------------------------------------------------------------------------------
% anomalous_diff3D_mancvg_test.m
% 04/2018, David Shirokoff
%
% Description: Performs a manufactured convergence test in 3D 
% for ImEx coefficients and ImEx splitting parameters (delta, sigma).
% Test is for the a nonlinear coefficient diffusion equation 
% u_t = Au + Bu + f(x,t), 
% Au ~ sigma d^2u/dx^2,  Bu ~ d/dx(  ( u(x)^gamma - 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 anomalous_diff3D_mancvg_test
    format long;    
                                        %----------------------------------------------
	option = 2;                         % Specify different manufactured solution test
                                        % options (default = 2) 
	gam = 5/3;
                            
                                        %----------------------------------------------
                                        % Variable Initializations
    n = 64;                             % number grid points.
    h = 1/n;                            % grid size, domain = [0, 1]
    x = (0:n-1)*h;                      % set grid
    y = (0:n-1)*h;    
    z = (0:n-1)*h;
    [X, Y, Z] = meshgrid(x, y, z);    
    
    tf = .3;                            % final integration time.
    del_t_vect = (2.^(-5:-1:-7)');      % set different time steps
                                    
    num_points = size(del_t_vect, 1);   % save number of integration runs.          
    
    
    %-------------------------------------------------------------------------------------
    r_max = 5;                          % Change this to change the max order being tested
                                        % NOTE: changing r_max may change
                                        % the values of (delta, sigma) --
                                        % which will change the error
                                        % constant observed.
    order_vect = 1:r_max;               % Array of the orders to loop over in cvg test. 
                                        % (should all be less than r_max).
    num_order = size(order_vect, 2);
        
    % Get ImEx parameter and splitting parameter.
    [delta sigma] = GetImExParameters(option, r_max, gam);
    
    % Error array -- saves errors for table output at end of code.
    err_vect = zeros(num_points, num_order);    
    
    
    %---------------------------------------------------------------------
    %Loop over different orders
    for ord_counter = 1:num_order
        
        % get order and ImEx coefficients, for the given order.
        order = order_vect(ord_counter);
        [a b c] = GetImExCoefficients(order, delta);
        
        %Loop over different dt values
        for m = 1:num_points                
            
            % Get time step del_t
            del_t = del_t_vect(m);
            num_steps = round(tf/del_t);
            true_tf = num_steps*del_t;
            
            % Output loop information to screen.
            fprintf('\nOrder = %f, \t dt = %f \n', ord_counter, del_t)
            
            % Call ImEx time integrator. Time step to tf.             
            [Uf] = ImEx_TimeIntegrate(option, del_t, num_steps, X, Y, Z, a, b, c, n, sigma, gam); 
            
            % Get exact solution and compute the error.
            Uexact = ExactSolution(option, X, Y, Z, true_tf);                
            err_vect(m, ord_counter) = max(max( max(abs(Uf - Uexact))));            
        end
    end

    % Compute instantaneous convergence rates
    cvgrate = zeros(num_points, num_order);
    cvgrate(2:end, :) = log( err_vect(1:end-1, :)./err_vect(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 counter = 1:num_points
        fprintf('%f \t & ', del_t_vect(counter));
        for jj = 1:num_order
            fprintf('\t  %1.1e \t & \t %1.1f &', err_vect(counter, jj), cvgrate(counter, jj));
        end
        fprintf('\n');        
    end
    

end

%------------------------------------------------------------------
% Function: GetImExParameters
%
% Description: Computes the ImEx parameter delta, and splitting
% parameter sigma for a given manufactured solution option number.
%
% Inputs: option -- specify which test is being run (default 2)
% gam -- adiabatic constant (default 5/3).
% r_max -- output coefficients (delta,sigma) stabilize schemes with 
% orders 1 <= r <= r_max. 
%
% Outputs: (delta, sigma) parameters. 
%------------------------------------------------------------------
function [delta sigma] = GetImExParameters(option, r_max, gam)        
        switch option 
            case {1, 2}          
                % `Gap parameter' -- closer to 0 gives a better
                % error constant, but is less robust.
                eta = 0.1;          
                
                % Estimate umax and umin for the solution over all x, t.
                umax = 2*exp(1) + exp(1);
                umin = 2*exp(1) - exp(1);
                
                % Compute max and min diffusion coefficients.
                dmin = umin^gam;
                dmax = umax^gam;
    
                % Temporary variables
                kap = dmin/dmax*(1-eta);
                cr = cos(pi/r_max)^(-r_max);
                
                % Use formulas to compute (delta, sigma).
                sigma = dmin*(1-1/2*eta)*( 1 + cr )/(1 + kap*cr);
                delta = min(2 - 2*( (1-kap)/(1+kap*cr ) )^(1/r_max), 1);
        end
end

%------------------------------------------------------------------
% Function: ImEx_TimeIntegrate
%
% Description: Uses ImEx multistep methods to integrate to the final 
% time.
%
% Outputs: Uf -- solution u at grid points x and final time 
%------------------------------------------------------------------
function [Uf] = ImEx_TimeIntegrate(option, del_t, num_steps, X, Y, Z, a, b, c, n, sigma, gam)
                
        % Build Fourier wavenumbers. 
        L = 1;      %domain length 1
        ikx = 1i*(mod((1:n) - ceil(n/2+1),n) - floor(n/2))*(2*pi)/L;
        ikx(n/2 + 1) = 1i*(n/2)*(2*pi)/L;   ikx(1) = 0;
        iky = 1i*(mod((1:n) - ceil(n/2+1),n) - floor(n/2))*(2*pi)/L;
        iky(n/2 + 1) = 1i*(n/2)*(2*pi)/L;    iky(1) = 0;
        ikz = 1i*(mod((1:n) - ceil(n/2+1),n) - floor(n/2))*(2*pi)/L;
        ikz(n/2 + 1) = 1i*(n/2)*(2*pi)/L;    ikz(1) = 0;
    
        kxyz = zeros(n, n, n);      %Build Laplacian multipliers
        ikx_mat = zeros(n, n, n);
        iky_mat = zeros(n, n, n);
        ikz_mat = zeros(n, n, n);
        for jj = 1:n
            for kk = 1:n
                for ll = 1:n
                    kxyz(jj, kk, ll) = iky(jj)^2 + ikx(kk)^2 + ikz(ll)^2;
                    ikx_mat(jj, kk, ll) = ikx(kk);
                    iky_mat(jj, kk, ll) = iky(jj);
                    ikz_mat(jj, kk, ll) = ikz(ll);
                end
            end
        end  
        order = size(a, 2)-1;
                
        % Compute initial data: u_n, Au_n and Bu_n, for n = 1, 2.., r-1. 
        % We use the exact initial data for the manufactured solution.
        u_old = zeros(n, n, n, order);
        bu_old = zeros(n, n, n, order);
        au_old = zeros(n, n, n, order);
        f_old = zeros(n, n, n, order);
        for jj = 1:order
            t = -(jj-1)*del_t;
            u_old(:,:, :, jj) = ExactSolution(option, X, Y, Z, t);
            [diff_coeffb] = GetDiffusionCoeff(option, X, Y, Z, u_old(:, :, :, jj), sigma, gam);
            
            % Compute the Bu terms
            utemp_hat = fftn(u_old(:, :, :, jj) );
            ux_temp = ifftn(ikx_mat.*utemp_hat);
            uy_temp = ifftn(iky_mat.*utemp_hat);
            uz_temp = ifftn(ikz_mat.*utemp_hat);
            vx = ux_temp.*diff_coeffb; 
            vy = uy_temp.*diff_coeffb;
            vz = uz_temp.*diff_coeffb;
            divhat = fftn(vx).*ikx_mat + fftn(vy).*iky_mat + fftn(vz).*ikz_mat;
            bu_old(:, :, :, jj) = real(ifftn(divhat)); 
            
            % Compute the Au terms
            au_old(:, :, :, jj) = sigma*real(ifftn(kxyz.*utemp_hat));   
            f_old(:, :, :, jj) = Forcing(option, X, Y, Z, -(jj-1)*del_t, gam);
        end           
        imp_part = a(1)*ones(n, n, n) - sigma*del_t*c(1)*kxyz;
        
        % Main time loop -- integrate num_steps.
        for counter = 1:num_steps            

            % compute new right hand side.
            rhs = zeros(n,n,n);
            for ind = 1:order
                rhs = rhs - u_old(:, :, :, ind)*a(1+ind) + del_t*( au_old(:,:, :, ind)*c(1+ind) + ...
                    bu_old(:, :, :, ind)*b(1+ind) ) + del_t*f_old(:, :, :, ind)*b(1+ind);            
            end

            % solve for new u_{n+1} via Fourier transform.
            u = real(ifftn(fftn(rhs)./imp_part));

            % update stored vectors, one time step.
            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
            
            % For efficiency, compute and store A u_{n+1}, B u_{n+1}, and
            % f(t_{n+1})
            [diff_coeffb] = GetDiffusionCoeff(option, X, Y, Z, u, sigma, gam);
            utemp_hat = fftn(u);
            ux_temp = ifftn(ikx_mat.*utemp_hat);
            uy_temp = ifftn(iky_mat.*utemp_hat);
            uz_temp = ifftn(ikz_mat.*utemp_hat);
            vx = ux_temp.*diff_coeffb; 
            vy = uy_temp.*diff_coeffb;
            vz = uz_temp.*diff_coeffb;
            divhat = fftn(vx).*ikx_mat + fftn(vy).*iky_mat + fftn(vz).*ikz_mat;
            bu_1 = real(ifftn(divhat)); 
            au_1 = sigma*real(ifftn(kxyz.*utemp_hat));           
            u_old(:, :, :, 1) = u;
            bu_old(:, :, :, 1) = bu_1;
            au_old(:, :, :, 1) = au_1;
            f_old(:, :, :, 1) =  Forcing(option,  X, Y, Z, counter*del_t, gam);
        end
        % return solution u at final time.
        Uf =  u;        
end

%--------------------------------------------------------------------------
% Function: GetDiffusionCoeff
%
% Description: Computes the nonlinear diffusion coefficient.
%
% Inputs: option -- manufactured solution option (default 2)
% (X,Y,Z) -- allows for variable spatial dependence to be incorporated 
% u_in -- function value
% gam -- adiabatic constant (default 5/3)
% sigma -- splitting parameter
%
% Outputs: diff_coeffb -- diffusion coefficient used in the B (explicit)
% term of the ImEx scheme.
%--------------------------------------------------------------------------
function [diff_coeffb] = GetDiffusionCoeff(option, X, Y, Z, u_in, sigma, gam)
    switch option
        case {1, 2}            
        	diff_coeffb = u_in.^gam - sigma;   
    end            
end

%--------------------------------------------------------------------------
% Function: ExactSolution
%
% Description: Returns the exact solution for the manufacture solution test
% at time t. 
%
% Inputs: option -- manufactured solution option (default 2)
% (X,Y,Z) -- allows for variable spatial dependence to be incorporated 
% t -- time to be evaluated at.
%
% Outputs: Uexact -- exact solution (as a 3d array at locations (x,y,z))
%--------------------------------------------------------------------------
function [Uexact] = ExactSolution(option, X, Y, Z, t)
    switch option
        case 1
            C = 2*exp(1);
            w = 1;
            Uexact = C + exp(sin(4*pi*X))*cos(w*t);
        case 2
            C = 2*exp(1);
            w = 1;
            Uexact = C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(w*t);            
    end
end

%--------------------------------------------------------------------------
% Function: Forcing
%
% Description: Computes the forcing required for the manufactured solution
% option at time t. 
%
% Inputs: option -- manufactured solution option (default 2)
% (X,Y,Z) -- allows for variable spatial dependence to be incorporated 
% gam -- adiabatic constant (default 5/3)
%
% Outputs: forcing -- forcing (in form of 3d array at (x,y,z)) at time t. 
%--------------------------------------------------------------------------
function [forcing] = Forcing(option, X, Y, Z, t, gam)
    switch option
        case 1
            C = 2*exp(1);
            w = 1;
            ut = -w*exp(sin(4*pi*X))*sin(t*w);
            diffterm = 16*pi^2*exp(sin(4*pi*X)).*cos(4*pi*X).^2*cos(t*w).*(C + exp(sin(4*pi*X)).*cos(t*w)).^gam -...
                16*pi^2*exp(sin(4*pi*X)).*sin(4*pi*X).*cos(t*w).*(C + exp(sin(4*pi*X)).*cos(t*w)).^gam + ...
                16*gam*pi^2*exp(2*sin(4*pi*X)).*cos(4*pi*X).^2.*cos(t*w).^2.*(C + exp(sin(4*pi*X)).*cos(t*w)).^(gam - 1);
            
            forcing = ut - diffterm;            
        case 2
            C = 2*exp(1);
            w = 1;
            ut = -w*exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*sin(t*w);
            u1 = 16*pi^2*exp(sin(4*pi*X)).*cos(4*pi*X).^2.*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w).*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^gam -...
                16*pi^2*exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*sin(4*pi*X).*cos(t*w).*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^gam +...
                16*gam*pi^2*exp(2*sin(4*pi*X)).*cos(4*pi*X).^2.*cos(2*pi*Y).^2.*cos(2*pi*Z).^2.*cos(t*w).^2.*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^(gam - 1);
            
            u2 = 4*gam*pi^2*exp(2*sin(4*pi*X)).*cos(2*pi*Z).^2.*sin(2*pi*Y).^2.*cos(t*w).^2.*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^(gam - 1) - ...
                4*pi^2*exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w).*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^gam;
            
            u3 = 4*gam*pi^2*exp(2*sin(4*pi*X)).*cos(2*pi*Y).^2.*sin(2*pi*Z).^2.*cos(t*w).^2.*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^(gam - 1) - ...
            4*pi^2*exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w).*(C + exp(sin(4*pi*X)).*cos(2*pi*Y).*cos(2*pi*Z).*cos(t*w)).^gam;
            
            forcing = ut - (u1 + u2 + u3);    
    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