% Buckley Leverett with N=100, dx=1/100, Central Scheme with Koren Limiter in Space
% Stepped to final time t=1/8
%

clear all; close all;

data_file = 'tvdpos_report_msrk';

N_points = 100;

dt_fe = 1/N_points/4;



%  dudx_template = @SSP_Tools.Discretizers.FiniteDifference;
dudx_template = @(varargin) SSP_Tools.Discretizers.KorenLimiter();
problem_template = @(varargin) SSP_Tools.TestProblems.BuckleyLeverett( 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', '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], ...
%                      [8,3,7], [7,4,7], [4,5,7], [8,3,8], [6,4,8], [6,5,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);

methods_to_test_hash = test_methods.keys();

if exist([data_file, '.mat'], 'file')
	load([data_file, '.mat']);
else
	test_results_tvd = containers.Map;
	test_results_pos = containers.Map;
end

diary([data_file, '.txt'])

for i=1:numel(methods_to_test_hash)

		% If we've already done this method and saved the results, skip it.
		if test_results_tvd.isKey(methods_to_test_hash{i}) & test_results_pos.isKey(methods_to_test_hash{i})
			continue;
		end
		
		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

save([data_file, '.mat'], 'test_results_tvd', 'test_results_pos');
	


% 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 0.25*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 & %1s  \\\\ \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, ...
			                     0.25*test_results_tvd(hash).r, ...
			                     0.25*test_results_tvd(hash).r /s, ...
										test_results_pos(hash).max_dt*N_points, ...
			                     test_results_pos(hash).max_dt*N_points / s, ...
			                     flagged );
		end
	end
end

diary('off');


s_values = containers.Map;

line_styles = { '-', '--', '-.' };
line_markers = { '*', 'o', '^', 'v' };
line_colors = { 'b', 'r', 'g'};

for p=[3,4]
	
	figure();
	
	k=[3,4,5];
	legend_titles = {};

	for ki=1:numel(k)
	
		obs_dt = [];
		thr_dt = [];
		
		
		sb=2:10;
		s = [];
		
		for si=1:numel(s)
			[sb(si),k(ki),p];
			hash = SSP_Tools.utils.DataHash([sb(si),k(ki),p]);
			if test_results_tvd.isKey(hash)
				tvd_tesults = test_results_tvd(hash);
			else
				continue;
			end
			
			s(end+1) = sb(si);
			obs_dt(end+1) = test_results_tvd(hash).max_dt;
			thr_dt(end+1) = 0.0025*test_results_tvd(hash).r;
		end
		
		plot(s, obs_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '-', 'LineWidth', 2);
		legend_titles{end+1} = sprintf('k=%i Order (Observed)', k(ki));
		
		hold all
		plot(s, thr_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '--', 'LineWidth', 2);
		legend_titles{end+1} = sprintf('k=%i Order (Theoretical)', k(ki));
	end
	
	legend(legend_titles, 'Location', 'SouthEast');
	filename = sprintf('Plot-%dp.eps', p)
	print(filename, '-depsc2');
end

for p=[5,6]
	
	figure();
	
	k=[3,4,5];
	legend_titles = {};

	for ki=1:numel(k)
	
		obs_dt = [];
		thr_dt = [];
		
		
		sb=3:10;
		s = [];
		
		for si=1:numel(s)
			[sb(si),k(ki),p];
			hash = SSP_Tools.utils.DataHash([sb(si),k(ki),p]);
			if test_results_tvd.isKey(hash)
				tvd_tesults = test_results_tvd(hash);
			else
				continue;
			end
			s(end+1) = sb(si);
			obs_dt(end+1) = test_results_tvd(hash).max_dt;
			thr_dt(end+1) = 0.0025*test_results_tvd(hash).r;
		end
		
		plot(s, obs_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '-', 'LineWidth', 2);
		legend_titles{end+1} = sprintf('k=%i Order (Observed)', k(ki));
		
		hold all
		plot(s, thr_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '--', 'LineWidth', 2);
		legend_titles{end+1} = sprintf('k=%i Order (Theoretical)', k(ki));
	end
	
	legend(legend_titles, 'Location', 'SouthEast');
	filename = sprintf('Plot-%dp.eps', p)
	print(filename, '-depsc2');
end

p=7;
figure();

k=[3,4];
legend_titles = {};

for ki=1:numel(k)

	obs_dt = [];
	thr_dt = [];
	
	
	sb=4:10;
	s = [];
	
	for si=1:numel(s)
		[sb(si),k(ki),p];
		hash = SSP_Tools.utils.DataHash([sb(si),k(ki),p]);
		if test_results_tvd.isKey(hash)
			tvd_tesults = test_results_tvd(hash);
		else
			continue;
		end
		s(end+1) = sb(si);
		obs_dt(end+1) = test_results_tvd(hash).max_dt;
		thr_dt(end+1) = 0.0025*test_results_tvd(hash).r;
	end
	
	plot(s, obs_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '-', 'LineWidth', 2);
	legend_titles{end+1} = sprintf('k=%i Order (Observed)', k(ki));
	
	hold all
	plot(s, thr_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '--', 'LineWidth', 2);
	legend_titles{end+1} = sprintf('k=%i Order (Theoretical)', k(ki));
end

k=5;

obs_dt = [];
thr_dt = [];

sb=3:10;
s = [];
for si=1:numel(s)
	[sb(si),k,p];
	hash = SSP_Tools.utils.DataHash([sb(si),k,p]);
	if test_results_tvd.isKey(hash)
		tvd_tesults = test_results_tvd(hash);
	else
		continue;
	end
	s(end+1) = sb(si);
	obs_dt(end+1) = test_results_tvd(hash).max_dt;
	thr_dt(end+1) = 0.0025*test_results_tvd(hash).r;
end

plot(s, obs_dt, 'Marker', line_markers{3}, 'Color', line_colors{3}, 'LineStyle', '-', 'LineWidth', 2);
legend_titles{end+1} = sprintf('k=%i Order (Observed)', k);

hold all
plot(s, thr_dt, 'Marker', line_markers{3}, 'Color', line_colors{3}, 'LineStyle', '--', 'LineWidth', 2);
legend_titles{end+1} = sprintf('k=%i Order (Theoretical)', k);



legend(legend_titles, 'Location', 'SouthEast');
filename = sprintf('Plot-%dp.eps', p)
print(filename, '-depsc2');


for p=[8]
	
	figure();
	
	k=[3,4,5];
	legend_titles = {};

	for ki=1:numel(k)
	
		obs_dt = [];
		thr_dt = [];
		
		sb=5:10;
		s = [];
		for si=1:numel(sb)
			[sb(si),k(ki),p];
			hash = SSP_Tools.utils.DataHash([sb(si),k(ki),p]);
			
			if test_results_tvd.isKey(hash)
				tvd_tesults = test_results_tvd(hash);
			else
				continue;
			end
			
			s(end+1) = sb(si);
			obs_dt(end+1) = test_results_tvd(hash).max_dt;
			thr_dt(end+1) = 0.0025*test_results_tvd(hash).r;
		end
		
		plot(s, obs_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '-', 'LineWidth', 2);
		legend_titles{end+1} = sprintf('k=%i Order (Observed)', k(ki));
		
		hold all
		plot(s, thr_dt, 'Marker', line_markers{ki}, 'Color', line_colors{ki}, 'LineStyle', '--', 'LineWidth', 2);
		legend_titles{end+1} = sprintf('k=%i Order (Theoretical)', k(ki));
	end
	
	legend(legend_titles, 'Location', 'SouthEast');
	filename = sprintf('Plot-%dp.eps', p)
	print(filename, '-depsc2');
end