clear all;close all;clc;

%Grid
Nx = 1024; Nz = 256;
Lx = 4.0; Lz = 1.0;
hx = Lx/Nx; hz = Lz/Nz;
x = -Lx/2:hx:(Lx/2)-hx;   z = -Lz/2:hz:(Lz/2)-hz;
[xx,zz]= meshgrid(x,z);

% Model parameters
Ra = 1000;
var_lnk=1; %variance of log(natural) permeability
corr_lenx=10*hx; %correlation length in x
corr_lenz=10*hz; %correlation length in z

% Initial condition
c = zeros(Nz,Nx);
c0 = 1; % Prescribed concentration at the top
c(1,:) = c0; % Prescribed concentration at the top
% c(1,:) = 1e-03*rand(1,Nx); % Small perturbation if needed
% Visualize the initial concentration field
% surf(xx,zz,c,'facecolor','interp','edgecolor','none','facelighting','phong'); view([0,0,1]); 
% axis equal tight; colorbar

% Generate multigaussian random permeability field
[perm,var_lnk_actual]=random_perm(var_lnk,corr_lenx,corr_lenz,Nx,Nz,Lx,Lz);

% Visualize log permeability field
surf(xx,zz,log(perm),'facecolor','interp','edgecolor','none','facelighting','phong');view([0,0,1]); 
axis equal tight off; colorbar; %caxis([2 7]);

% Save the permeability field to use in future runs. 
% Since random number generator is used in the function below, every call generates a new field with presricbed variance and correlation lengths.
% dlmwrite(['perm_'  num2str(var_lnk) '_' num2str(corr_lenx/hx) '_' num2str(corr_lenz/hz) '.txt'],perm);
save perm.mat perm

t = 0; istep= 0; 
dt = 0.0001;
% dtadv = min(hz,hx);
% dtdiff = Ra/2*dtadv.^2;
% dt = 1/2*min(dtdiff,dtadv);

while (t<1)
    istep= istep+1;
	
    % Solve for pressure and velocities implicitly using TPFA/MPFA 
    [ux,uz]= solve_vel(c,perm,hx,hz);
    
	% Solve for concentration using RK explicit time integration 
    for istage= 1:3
        Res= residual_conc(c,ux,uz,Ra,hx,hz,c0); % Residual of the concentration eqn.
        if istage==1
            c= c + dt*( (8/15)*Res );
            Res1= Res;
        elseif istage==2
            c= c + dt*( (5/12)*Res + (-17/60)*Res1 );
            Res1= Res;
        else
            c= c + dt*( (3/4)*Res + (-5/12)*Res1 );
        end
    end
    t = t+dt;
    dt = min(hx,hz)/(4*max(max(ux)));
	
    %%% save output %%%%%%%%%%%%%%%%%%        
	surf(xx,zz,c,'facecolor','interp','edgecolor','none','facelighting','phong'); view([0,0,1]); 
	axis equal tight; colorbar;
	drawnow    
end

