2D advection boundary conditions

I'm trying to produce a simple simulation of a two-dimensional advection equation, but am having trouble with applying periodic boundary conditions. I have the following code:

clear
clc
close all

Xmin = -1; % starting x value
Xmax = 1;% ending x value
xgridpoints = 31;
dx = (Xmax-Xmin)/xgridpoints;
x = Xmin:dx:Xmax;

Ymin = -1; % starting y value
Ymax = 1;% ending y value
ygridpoints = 31;
dy = (Ymax-Ymin)/ygridpoints;
y = Ymin:dy:Ymax;

lambda = 0.4; % can change to desired lambda value
Tmin = 0;
Tmax = 2;
dt = dx*lambda;
t = Tmin:dt:Tmax;

a = 1; %x-wave speed
b = 1; %y-wave speed

%Creating 2D grid for u=u(x,y)
un = zeros(length(x), length(y));

%Initial conditions for quantity u within the xy-grid
for i = 1:length(x)
    for j = 1:length(y)
       if abs(x(i)) <= (1/3) && abs(y(j)) <= (1/3)
           un(i,j) = cos(x(i))*sin(y(j));
       end
    end
end

i = 2:length(x) %Creating i index (starting from 2 due to i-1 term in FTBS)
j = 2:length(y) %Creating j index (starting from 2 due to j-1 term in FTBS)
for n = 1:length(t)
    h=surf(x,y,un);
    axis ([-1 1 -1 1 -0.2 0.2]);
    xlabel('x')
    ylabel('y')
    drawnow; 
    %pause(0.1)
    un(i,j) = un(i,j) - a*lambda*(un(i,j) - un(i-1,j)) - b*lambda*(un(i,j) - un(i,j-1));
    %FIX BOUNDARY CONDITIONS!
    un(1,:) = un(length(x),:); %periodic BC applied explicitly
    un(:,1) = un(:,length(y)); %periodic BC applied explicitly       
end

This produces a moving wave that initially looks as follows:

where the wave speeds of a=1 and b=1 force the wave to move towards the (x,y)=(1,1) corner. All is well up to this point, but I'd like the periodic boundary conditions to make the wave "reappear" at the (x,y)=(-1,-1) corner. Currently, this is unfortunately not happening, and I'm hoping someone could scrutinize my boundary condition implementation to see where I'm going wrong. Thank you!

EDIT: Essentially, the wave moves along the black line seen in the image below:

[Torsten]

I think the main problem with your code is that by setting

un(i,j) = un(i,j) - a*lambda*(un(i,j) - un(i-1,j)) - b*lambda*(un(i,j) - un(i,j-1));

you overwrite the old value for "un(i,j)", but this value is still needed later on in your loop (e.g. in the computation of un(i+1,j)).

Use two arrays ("un_old" and "un", e.g.) to avoid this mixture of old and updated values.

Best wishes

Torsten.

Answer 1

[JeffR1992]

Unfortunately this does not seem to fix the boundary conditions, as they still act in the same strange way.

[Torsten]

Your loop over n must be the outer loop, not the inner loop.

for n = 1:length(t)
  un(1,:) = ...;
  un(:,1) = ...;
  for i = 2:length(x)
    for j = 2:length(y)
      ...
    end
  end 
end

Best wishes

Torsten.

Answer 2

[JeffR1992]

Hi Torsten, thanks for your help, if you run the code and see the surface animation you'll see that those methods unfortunately don't work, as the solution goes "haywire" when the wave reaches the top right corner of the domain. There is something fundamentally wrong with how I coded the two lines of code that dictate the boundary conditions. It's difficult for me to explain it in words, but if you run the attached code and see the animation you'll understand.