% This script produces the TVD report
%
%

clear all; close all;

N_points = 100;

% Linear Advection with N=100, dx=1/100, 1st Order Finite Difference in Space
% Stepped to final time t=1/8
dt_fe = 1/N_points;



%  dudx_template = @SSP_Tools.Discretizers.FiniteDifference;
dudx_template = @(varargin) SSP_Tools.Discretizers.FiniteDifference();
problem_template = @(varargin) SSP_Tools.TestProblems.Advection('a', 1, ...
                                                                'initial_condition', 'stepfunction', ...
                                                                varargin{:} );
                                                                               

% First we're going to catalog all of the methods we're testing by their
% number of stages (s), steps (k), and order (p). By hasing the numeric
% array containing the values [s,k,p] we can easily index methods in
% a matlab containers.Map object.

test_methods = containers.Map;

% For instance, RK3 s=3, k=0, p=3, so we'll
% first generate the hash...

rk3_hash = SSP_Tools.utils.DataHash([3,1,3]);

% Then use it as an intest for the RK3 object

test_methods(rk3_hash) = SSP_Tools.Integrators.RK3();

% Now we'll do the same for RK4 and SSP(10,4)

rk4_hash = SSP_Tools.utils.DataHash([4,1,4]);
test_methods(rk4_hash) = SSP_Tools.Integrators.RK4();

ssp104_hash = SSP_Tools.utils.DataHash([10, 1, 4]);
test_methods(ssp104_hash) = SSP_Tools.Integrators.SSP104();

% Now we're going to do the same for the MSRK methods which
% all use the same MATLAB class, but are differentiated by
% different coefficient files that are passed to that as
% at initialization

msrk_dir = [pwd, '/+SSP_Tools/Method Coefficients/Multistep Multistage (MSRK)'];
file_list = dir([msrk_dir, '/*.mat']);
file_list = { file_list.name };
% Now that we have a list of all the coefficient files, we're going to
% go through them one at a time, extract the [s,k,p] values from the
% filename, hash them, and initialize the appropriate MSRK object.

for i=1:numel(file_list)
	[~, file_name, file_ext] = fileparts(file_list{i});
	match_str = '(.*)s(.*)k(.*)p(?:T2|T3|T|$)';
	skp = regexp(file_name, match_str, 'tokens');
	skp = cellfun(@str2num, skp{1});
	hash = SSP_Tools.utils.DataHash(skp);
	test_methods(hash) = SSP_Tools.Integrators.MSRK('coefficients', [msrk_dir, '/', file_list{i}], ...
	                                                'initial_integrator', 'use-exact' );
end


% Now for the test.

test_results_tvd = containers.Map;
test_results_pos = containers.Map;

methods_to_test_skp = { [2,3,3], [6,3,3], [3,1,3], [7,3,3], [10,1,4], [4,1,4], [2,3,4], ...
                    [3,4,4], [7,3,4], [3,3,5], [6,3,5], [3,4,5], [3,5,5], ...
                    [5,3,6], [9,3,6], [4,4,6], [3,5,6], [6,5,6], [7,3,7], ...
                    [7,4,7], [4,5,7], [8,3,8], [9,5,8], [9,4,9],  ...
                    [20, 3, 10] };
                    
methods_to_test_hash = cellfun(@(id) SSP_Tools.utils.DataHash(id), methods_to_test_skp, 'UniformOutput', false);

diary('tvdpos_advection_log.txt')


for i=1:numel(methods_to_test_skp)

		methods_to_test_skp{i}
		dudt = test_methods(methods_to_test_hash{i});

		dudt_tvd = dudt.copy();
		dudx_tvd = dudx_template();
		problem_tvd = problem_template('integrator', dudt_tvd, 'discretizer', dudx_tvd);
		test_tvd = SSP_Tools.Tests.SSP('problem', problem_tvd, 'N', N_points);
		results_tvd = test_tvd.run();
		
		dudt_pos = dudt.copy();
		dudx_pos = dudx_template();
		problem_pos = problem_template('integrator', dudt_pos, 'discretizer', dudx_pos);
		test_pos = SSP_Tools.Tests.Positivity('problem', problem_pos, 'N', N_points);	
		results_pos = test_pos.run()
		
		test_results_tvd(methods_to_test_hash{i}) = results_tvd;
		test_results_pos(methods_to_test_hash{i}) = results_pos;
		
end
	

% Print the results
header_format = '%11s & %11s & %11s & %11s & %11s & %11s & %11s \n\n';
fprintf(header_format, 'Method', 'dt_tvd/dx', 'dt_tvd/dx/s', 'r', 'r/s', 'dt_pos/dx', 'dt_pos/dx/s');

for p=1:10
	for s=1:20
		for k=0:10
			hash = SSP_Tools.utils.DataHash([s,k,p]);
			
			if ~test_results_tvd.isKey(hash) | ~test_results_pos.isKey(hash)
				continue;
			end			
			
			if test_results_tvd(hash).r/N_points - test_results_tvd(hash).max_dt > 1e-6
				flagged = '*';
			else
				flagged = ' ';
			end
			
			line_format = '%8s & %6.3f & %6.3f & %6.3f & %6.3f & %6.3f & %6.3f  \\\\ \n';
			fprintf(line_format, sprintf('(%d,%d,%d)', s,k,p), ...
			                     test_results_tvd(hash).max_dt*N_points, ...
			                     test_results_tvd(hash).max_dt*N_points / s, ...
			                     test_results_tvd(hash).r, ...
			                     test_results_tvd(hash).r /s, ...
										test_results_pos(hash).max_dt*N_points, ...
			                     test_results_pos(hash).max_dt*N_points / s);
		end
	end
end

diary('off')