%Scalar advection demo
clf;close all;
for jj=1:3
    %set number of unknowns, domain size
    Nx   = 5*2^jj;      Ny   = 5*2^jj;
    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);

    %Decide to plot solution or not
    ifplot=1;

    %% 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);
    % P(Nbcs+1:end) = exp(-10*((XP-a/2+a/8*0).^2+(YP-b/2+b/4).^2).^.5);
    u(Nbcs+1:end) = a;%bcsDu(1);%-YU/a+1/2;
    v(Nbcs+1:end) = b;%bcsDv(2);% XV/b-1/2;
    % u(Nbcs+1:end) = -2*pi*(YU-b/2);
    % v(Nbcs+1:end) =  2*pi*(XV-a/2);
    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
    %% 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 50 times
        if ~mod(i,Nt/50)
            surf(xp,yp,P(NodeP(2:Ny+1,2:Nx+1)),'LineStyle','none') ,colorbar
            xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
            TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
            TT(4)=title(sprintf('Numerical Solution',Nx,Ny));
            formatplot([],gca,TT);caxis([-1,1])%view(2)
            drawnow
        end
    end
    %Calculate error
    ERR(jj) = max(abs((Pi-P)));
    %Plot errors
    if ifplot
        figure,subplot(131)
        surf(xp,yp,P(NodeP(2:Ny+1,2:Nx+1)),'LineStyle','none') ,colorbar
        xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
        TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
        TT(4)=title(sprintf('Numerical solution',Nx,Ny));
        formatplot([],gca,TT);view(2)
        subplot(132)
        surf(xp,yp,Pi(NodeP(2:Ny+1,2:Nx+1)),'LineStyle','none') ,colorbar
        xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
        TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
        TT(4)=title(sprintf('Analytical Solution',Nx,Ny));
        formatplot([],gca,TT);view(2)
        subplot(133)
        surf(xp,yp,Pi(NodeP(2:Ny+1,2:Nx+1))-P(NodeP(2:Ny+1,2:Nx+1)),'LineStyle','none') ,colorbar
        xlim([min(xp) max(xp)]);ylim([min(yp) max(yp)]);
        TT(1)=xlabel('x(m)');TT(2)=ylabel('y(m)');TT(3)=zlabel('T(^oC)');
        TT(4)=title(sprintf('Analytical - Numerical',Nx,Ny));
        formatplot([],gca,TT);view(2)
        figure(1)
    end
end
fprintf('Convergence rate is in infinity norm is %g\n',log(ERR(1:end-1)./ERR(2:end))/log(2))
