%Scalar advection demo
%set domain size
a    = 1;           b    = 1;

%Calculate dx and dy
dx   = a./(Nx);     dy   = b./(Ny);

%Set periodic boundary conditions
bcsD = [];          bcsN = [];
bcsDu= [];          bcsNu=[];
bcsDv= [];          bcsNv=[];

%Set time integration parameters
T    = 1;          Nt   = 1000;       dt   = 1/(Nt);
if isfield(app,'T')
    T = app.T;
end
if isfield(app,'dt')
    dt = app.dt;
end
Nt = ceil(T/dt);

%% create geometry mask
periodic=1;
bnd=[0 0 0 0];
Mask = zeros(Ny, Nx);
Nbcs =length(bcsD)+length(bcsN);

%create Node matrices
[NodeP, Nodeu, Nodev, idsP, idsu, idsv] = NodePad(Mask,bnd,periodic);

%% Initialize vectors
%Create x vectors for intialization
xp = linspace(dx/2, a-dx/2, Nx  );yp = fliplr(linspace(dy/2, b-dy/2, Ny  ));
xu = linspace( dx   ,a-dx , Nx-1);yu = fliplr(linspace(dy/2, b-dy/2, Ny  ));
xv = linspace(dx/2, a-dx/2, Nx  );yv = fliplr(linspace(dy    , b-dy, Ny-1));
[XP YP] = meshgrid(xp,yp);[XU YU] = meshgrid(xu,yu);[XV YV] = meshgrid(xv,yv);

%initialize vectors
P = [bcsD' ; bcsN'*NaN; zeros( Nx    * Ny    ,1)];
u = [bcsDu'; bcsNu'*NaN; zeros((Nx-1) * Ny    ,1)];
v = [bcsDv'; bcsNv'*NaN; zeros( Nx    * (Ny-1),1)];
P(Nbcs+1:end) = sin(2*pi.*XP/a).*sin(2*pi.*YP/b);

u(Nbcs+1:end) = a;
v(Nbcs+1:end) = b;

Pi=P;


%% Reformat for plotting
if ~periodic
    xp = linspace(-dx/2, a+dx/2, Nx+2);
    yp = linspace(-dy/2, b+dy/2, Ny+2);
    xu = linspace( 0   , a     , Nx+1);
    yu = linspace(-dy/2, b+dy/2, Ny+2);
    xv = linspace(-dx/2, a+dx/2, Nx+2);
    yv = linspace(0    , b     , Ny+1);
    [XP YP] = meshgrid(xp,yp);
    [XU YU] = meshgrid(xu,yu);
    [XV YV] = meshgrid(xv,yv);
end

surf(xp,yp,P(NodeP(2:Ny+1,2:Nx+1)),'LineStyle','none') ,colorbar
xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);zlim([-1 1]);
TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
TT(4)=title(sprintf('Numerical Solution Nx=%g, Ny=%g',Nx,Ny));
formatplot([],gca,TT);caxis([-1,1])%view(2)
display('Paused at initial condition. Hit any key to continue');drawnow;pause

%% Integrate in time
for i=1:Nt
    P(Nbcs+1:end,1) = P(Nbcs+1:end,1)...
        + dt*SCAadvect(Nx, Ny, dx, dy, P, u, v, NodeP, Nodeu, Nodev, idsP);
    %Plot solution 20 times
    if ~mod(i,Nt/20)
        surf(xp,yp,P(NodeP(2:Ny+1,2:Nx+1)),'LineStyle','none') ,colorbar
        xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);zlim([-1 1]);
        TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
        TT(4)=title(sprintf('Numerical Solution Nx=%g, Ny=%g',Nx,Ny));
        formatplot([],gca,TT);caxis([-1,1])%view(2)
        drawnow
    end
end

