Comparação da Numérica com a Analíticca

clc
clear
%===============Spacial and Temporal mesh creation=========================
nSpacialPoints=500;
nElements=nSpacialPoints-1;
L=1;
T=1;
nTemporalPoints=100;
[x_nodes,x_centers,temporalMesh,dx,dt]=malha(L,nSpacialPoints,T,nTemporalPoints);

%=================Boundary and initialConditions===========================
pressureIn=1;
pressureOut=0;
%
concentrationIn=1;
concentrationOut=0;
%=====================Physical Properties ================================
permeability=ones([1, nElements]);
compressibility=0;
porosity=1;
waterViscosity=1e0;
difusionCoeff = 1*ones([1, nElements]);
%====================Numerical Solution definiton==========================
numericalConcentrationPredictor = zeros(nElements,1);
numericalConcentrationPredictor(1:end) = concentrationOut;
numericalConcentrationPredictorMatrix = zeros(nElements,nTemporalPoints);
numericalConcentrationPredictorMatrix(:,1) = numericalConcentrationPredictor;
%
numericalConcentrationCorrector = zeros(nElements,1);
numericalConcentrationCorrector(1:end) = concentrationOut;
numericalConcentrationCorrectorMatrix = zeros(nElements,nTemporalPoints);
%
numericalConcentrationFiniteVolumeArray = zeros(nSpacialPoints,1);
numericalConcentrationFiniteVolumeMatrix = zeros(nSpacialPoints,nTemporalPoints);
numericalConcentrationFiniteVolumeMatrix(:,1) = zeros(nSpacialPoints,1);
%
numericalPressurePrevious = zeros(nElements,1);
numericalPressurePrevious(:) = pressureOut;
%=======================Permeability array definition======================
condutivity=condutividade(permeability,waterViscosity,x_centers);
%===========================Hydrodynamics solver===========================
[A_hidro,b_hidro]=A_global_hidro(x_nodes,nSpacialPoints,nElements,condutivity,pressureIn,pressureOut,...
    compressibility,dx,dt,numericalPressurePrevious);
solucao_hidro=A_hidro\b_hidro;
A_hidro(:,:)=0.0;
b_hidro(:)=0.0;
pressure=solucao_hidro(nSpacialPoints+1:nSpacialPoints+nElements);
velocity=solucao_hidro(1:nSpacialPoints);
pressureprevious=pressure;
%============================Temporal Loop=================================
dtMicro=0;
for j=2:1:nTemporalPoints
    %========================CFL CHECK + MICRO MESH========================
    dt=temporalMesh(j)-temporalMesh(j-1);
    micro=false;
    if((max(velocity)*max(dt))/(min(porosity)*min(dx))>=1)
        dtMicro = ((min(porosity)*min(dx))/velocity(1))-1e-40;
        micro=true;
    end
    numericalConcentrationFiniteVolumeArray = volumes_finitos(porosity,velocity(1),difusionCoeff(1),...
    nSpacialPoints,dx,dt,numericalConcentrationFiniteVolumeArray,concentrationIn,concentrationOut);
    %==========================PREDICTOR STEP==============================
    [numericalConcentrationPredictor]=transporte_concentracao_explicito_vetorial(numericalConcentrationCorrector,...
    dx,dt,velocity,j,porosity,x_centers,temporalMesh,micro,dtMicro,concentrationIn);
    numericalConcentrationPredictorMatrix(:,j) = numericalConcentrationPredictor;
    numericalConcentrationFiniteVolumeMatrix(:,j) = numericalConcentrationFiniteVolumeArray;
    %==========================CORRECTOR STEP==============================
    [A_hidro,b_hidro]=A_global_hidro(x_nodes,nSpacialPoints,nElements,difusionCoeff,concentrationIn,concentrationOut,...
    porosity,dx,dt,numericalConcentrationPredictor);
    solucao_hidro=A_hidro\b_hidro;
    A_hidro(:,:)=0.0;
    b_hidro(:)=0.0;
    numericalConcentrationCorrector=solucao_hidro(nSpacialPoints+1:nSpacialPoints+nElements);
    flux=solucao_hidro(1:nSpacialPoints);
    numericalConcentrationCorrectorMatrix(:,j) = numericalConcentrationCorrector;

end

% Parâmetros do problema
U0 = concentrationIn; % Condição de contorno em x = 0
UL = concentrationOut; % Condição inicial
phi = porosity; % Porosidade
vel = velocity(1); % Velocidade
D = difusionCoeff(1); % Coeficiente de difusão
L = 1; % Comprimento característico
x = x_centers;
t = temporalMesh;

% Calcula a solução analítica
solucao_analitica = analitica_conveccao_difusao_cartesiana(U0,UL,phi,vel,D,L,x,t);

% Calcula a solução numérica
solucao_numerica1 = zeros(nElements, nTemporalPoints);
solucao_numerica2 = zeros(nElements, nTemporalPoints);

for i = 1:nTemporalPoints
    solucao_numerica1(:, i) = numericalConcentrationCorrectorMatrix(:, i);
    solucao_numerica2(:, i) = numericalConcentrationPredictorMatrix(:, i);
end

% Plotagem das soluções
for i = 1:nTemporalPoints
    clf;
    plot(x, solucao_numerica1(:, i), 'b', 'LineWidth', 2);
    hold on;
    plot(x, solucao_numerica2(:, i), 'm', 'LineWidth', 2);
    hold on;
    plot(x, solucao_analitica(:, i), 'r--', 'LineWidth', 2);
    hold on;
    xlabel('Posição (m)');
    ylabel('Concentração');
    title(['Comparação entre solução numérica e analítica (t = ', num2str(t(i)), ')']);
    legend('Corretor', 'Upwind implícito','Analítica');
    grid on;

    pause(0.1);
end