%%
%% Listing 7.15 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

% Należy odkomentować linię z oznaczeniem wybranej metody 

SOLVER = 'CVODES';	% MATLAB/OCTAVE, o ile zainstalowano sundialsTB
% SOLVER = 'ODE15S';	% tylko MATLAB
% SOLVER = 'ODE5R';	% tylko Octave

if ~exist('OCTAVE_VERSION')
% MATLAB
	stderr = 2;
	stb = '/home/przykry/Matlab/sundialsTB';
else
% Octave
	stb = '/home/przykry/Octave/sundialsTB';
end

addpath(stb);
startup_STB(stb); % inicjalizacja sundialsTB

global Lap2D;
global beta;

%------------------------
% PARAMETRY i PODSTAWOWE OBIEKTY
%------------------------

% parametry zadania

a = 1; b = 1; % rozmiary obszaru $\Omega$
Tmin = 0.0; Tmax = 15.0; % długość przedziału czasowego
beta = 5e-4; % współczynnik dyfuzji

% parametry dyskretyzacji przestrzennej

Nx = 266; Ny = 267; % wewnętrzne węzły siatki
NEQ = Nx*Ny; % rozmiar zadania
hx = a/(Nx+1); hy = b/(Ny+1); % rozmiar oczka siatki
xx = hx*[1:Nx]'; yy = hy*[1:Ny]'; % współrzędne węzłów dyskretyzacji 

% dyskretny operator dyfuzji

Lap1Dx = lap1ddset(Nx, hx); Ix = speye(Nx, Nx);
Lap1Dy = lap1ddset(Ny, hy); Iy = speye(Ny, Ny);
Lap2D = kron(Iy, Lap1Dx) + kron(Lap1Dy, Ix);

clear Lap1Dx;
clear Lap1Dy;

%------------------------
% DANE POCZĄTKOWE
%------------------------

u0 = zeros(Nx,Ny);

% dane początkowe wczytujemy z pliku w formacie row-major-order (jak w OpenDX)
u0datNx = 256; u0datNy = u0datNx;
u0dat = load(['cvrd-N', num2str(u0datNx), '-step0.txt']);
u0dat = permute(u0dat,[2 1]); % indeks wiersza ma odpowiadać współrzędnej OX
mNx = min(Nx,u0datNx); mNy = min(Ny,u0datNy);
u0(1:mNx,1:mNy) = u0dat(1:mNx,1:mNy);

% musimy "rozprostować" wektor tak, by mogły na nim działać solvery ODE
u0 = u0(:);

%------------------------
% PARAMETRY PRACY SOLVERA
%------------------------

% toleracja błędu
rtol = 1.0e-5; atol = 1.0e-6;

switch SOLVER

case 'CVODES'
	data.P = []; % przestrzeń robocza
	options = CVodeSetOptions('UserData', data, 'RelTol', rtol, 'AbsTol', atol, 'LinearSolver', 'GMRES');
	CVodeInit(@cvrdrhs, 'BDF', 'Newton', Tmin, u0, options);
case {'ODE15S', 'ODE5R'}
	options = odeset('AbsTol', atol, 'RelTol', rtol, 'JPattern', (Lap2D ~= 0), 'Stats', 'on'); 
otherwise
	error('Wybrano nieznana metode');
end

%------------------------
% ROZWIĄZANIE RÓWNANIA RÓŻNICZKOWEGO
%------------------------

nout = 100;
T = linspace(Tmin, Tmax, nout+1); % chwile, w których chcemy poznać rozwiązanie

U = zeros(nout+1, NEQ); % tablica na zapisanie wartości rozwiązania w chwilach T
U(1,:) = u0'; % profil początkowy

fprintf(stderr,'Start symulacji\n'); tic;
switch SOLVER
	case 'CVODES'
		for i = 2:nout+1
			tout = T(i);
			[status,t,u] = CVode(tout,'Normal');
			U(i,:) = u'; T(i) = t;
			si = CVodeGetStats;
			fprintf(stderr, 'status = %d   t = %.2e   nst = %d   q = %d   h = %.2e\n', ...
			 status, t, si.nst, si.qlast, si.hlast);
		end
		CVodeFree;
	case 'ODE15S'
		[T, U] = ode15s(@cvrdrhs, T, u0, options); 
	case 'ODE5R'
		[T, U] = ode5r(@cvrdrhs, [Tmin, Tmax], u0, options);
		nout = length(T)-1;
end
wct = toc; % pomiar czasu działania solvera

%------------------------
% ZAPIS lub WIZUALIZACJA WYNIKÓW
%------------------------

fprintf(stderr,'Czas obliczen: %g (s)\n', wct);
filename = [lower(SOLVER),'-',num2str(Nx),'x',num2str(Ny),'-step'];
for i = 1:nout+1
% animacja rozwiązania w czasie
	t = T(i); u = reshape(U(i,:),Nx,Ny);
	imagesc(xx,yy,u);
	axis('xy', 'square');	
	title(['t = ',num2str(t)]);
	if (i == 1) || (i == nout+1) 			
		colorbar('EastOutside');
		print('-depsc',[filename,num2str(i-1),'.eps']);
	end
	pause(0.1);
end
