%Script prova metodo FDM
%Conduzione calore - non stazionario - 2D
%
%Scritto da Pavan Andrea - 19/11/2021
clear;
clc;

%% dati input
dx = 0.005;      %discretizzazione lungo x
dy = 0.005;      %discretizzazione lungo y
x = 0:dx:0.1;       %griglia
x = x(:);
y = 0:dy:0.1;       %griglia
k = 0.1;        %conducibilità termica [W/(m*K)]
rho = 8.0;      %densità [kg/m3]
c = 3.8;        %calore specifico [J/(kg*K)]
Q = @(x,y,t) -exp(2*x+y).*(rho*c*sin(t)+5*k*cos(t));       %sorgenti calore [W/m3]
u0 = @(x,y) exp(2*x+y);     %condizioni iniziali [K]
uB = @(x,y,t) exp(2*x+y)*cos(t);        %condizioni contorno [K]
t0 = 0;     %istante iniziale [s]
tf = 1;     %istante finale [s]
ue = @(x,y,t) exp(2*x+y).*cos(t);       %soluzione esatta [K]


%% evoluzione temporale
dt = 0.001;      %passo temporale iniziale
u = zeros(length(x),length(y),2);       %soluzione approssimata
err = 0;      %errore massimo
t(1) = t0;
u(:,:,1) = u0(x,y);
l = 2;

while t(l-1)<tf
    t(l) = t(l-1)+dt;
    u(:,:,l) = uB(x,y,t(l));        %condizioni al contorno

    %punti interni
    for i=2:length(x)-1
        for j=2:length(y)-1
            %du_dt = F(t,u)
            F = (k/(rho*c))*((u(i+1,j,l-1)-2*u(i,j,l-1)+u(i-1,j,l-1))/(dx^2) + ((u(i,j+1,l-1)-2*u(i,j,l-1)+u(i,j-1,l-1))/(dy^2))) + Q(x(i),y(j),t(l-1))/(rho*c);
            u(i,j,l) = u(i,j,l-1) + dt*F;       %Eulero Esplicito
        end
    end

    %grafico evoluzione
    [X1,Y1] = meshgrid(x,y);
    Z1 = u(:,:,l);
    figure(1);
    [~,h] = contourf(X1,Y1,Z1,10*length(x));
    colormap(hot);
    set(h,'LineColor','none')
    colorbar
    title(['Soluzione numerica (t=' num2str(t(l)) 's)']);
    xlabel('x');
    ylabel('y');
    axis square;
    axis equal;

    err(l) = max(max(abs(u(:,:,l)-ue(x,y,t(l)))));
    l = l+1;
end

%grafico comparazione a tf
[X,Y] = meshgrid(x,y);
figure(1);
tiledlayout(1,2);
nexttile;
Z1 = u(:,:,end);
[~,h] = contourf(X,Y,Z1,10*length(x));
colormap(hot);
set(h,'LineColor','none')
colorbar
title(['Soluzione numerica (t=' num2str(t(end)) 's)']);
xlabel('x');
ylabel('y');
axis square;
axis equal;
nexttile;
Z2 = ue(x,y,t(end));
[~,h] = contourf(X,Y,Z2,10*length(x));
colormap(hot);
set(h,'LineColor','none')
colorbar
title(['Soluzione esatta (t=' num2str(t(end)) 's)']);
xlabel('x');
ylabel('y');
axis square;
axis equal;

%grafico errore
figure(3);
semilogy(t,err)
title('Andamento errore nel tempo');
xlabel('t');
ylabel('err');

