[HELP] Mur's ABC for 2D FDTD

Dear all,

I was working with the Mur's first order and second order ABCs for the 2D FDTD programming.
There is a very strange problem happening right now.
For BOTH the first order AND the second order Mur's ABC in my code, they only work for the upper and right boundaries, not for the top and bottom boundaries.
The top and bottom boundaries cannot absorb the wave somehow!

Could anyone help me with this?
I've been working on this problem for days and cannot figure out what is happening.

Thank you in advance!

---------- Post added at 23:06 ---------- Previous post was at 23:04 ----------

well, cannot upload file, says file invalid somehow. I'll just paste my file here.

clear all;

%Initialize matrices
Nx=201;
Ny=201;
N=500;
c=3e8;
dx=0.0005;
dy=0.0005;
dt=1/(c*sqrt((1/dx)^2+(1/dy)^2));
f=30e9;
epo=8.85419*10^-12;
epr=1;
mu=4*pi*10^-7;

Ez=zeros(Nx,Ny,N);
Hx=zeros(Nx,Ny,N);
Hy=zeros(Nx,Ny,N);

A=(c*dt-dx)/(c*dt+dx);
B=(2*dx)/(c*dt+dx);
C=((c*dt)^2)*dx/(2*(c*dt+dx)*dy^2);
R=(c*dt - dx)/(c*dt + dx);

for n=2:N
for i=2:Nx-1
for j=2:Ny-1
if ((i==30)&&(j==30))
Ez(i,j,n+1)=sin(2*pi*f*n*dt); %source
else
Hx(i,j,n+1)=Hx(i,j,n)-dt/(mu*dy)*(Ez(i,j+1,n)-Ez(i,j,n));
Hy(i,j,n+1)=Hy(i,j,n)+dt/(mu*dx)*(Ez(i+1,j,n)-Ez(i,j,n));
Ez(i,j,n+1)=Ez(i,j,n)+dt/(epo*epr*dx)*(Hy(i,j,n+1)-Hy(i-1,j,n+1))-dt/(epo*epr*dy)*(Hx(i,j,n+1)-Hx(i,j-1,n+1));
end
end
end

clear i; clear j;

%Boundary conditions
for j=2:Ny-1
Ez(Nx,j,n+1) = Ez(Nx-1,j,n) + R*(Ez(Nx-1,j,n+1) - Ez(Nx,j,n)); %right boundary
%Ez(Nx,j,n+1)=-Ez(Nx-1,j,n-1)+A*(Ez(Nx,j,n-1)+Ez(Nx-1,j,n+1))+B*(Ez(Nx,j,n)+Ez(Nx-1,j,n))+C*(Ez(Nx,j+1,n)-2*Ez(Nx,j,n)+Ez(Nx,j-1,n)+Ez(Nx-1,j+1,n)-2*Ez(Nx-1,j,n)+Ez(Nx-1,j-1,n));
Ez(1,j,n+1) = Ez(2,j,n) + R*(Ez(2,j,n+1)-Ez(1,j,n)); %left boundary
%Ez(1,j,n+1)=-Ez(2,j,n-1)+A*(Ez(1,j,n-1)+Ez(2,j,n+1))+B*(Ez(1,j,n)+Ez(2,j,n))+C*(Ez(1,j+1,n)-2*Ez(1,j,n)+Ez(1,j-1,n)+Ez(2,j+1,n)-2*Ez(2,j,n)+Ez(2,j-1,n));
end
for i=2:Nx-1
Ez(i,Ny,n+1) = Ez(i,Ny-1,n) + R*(Ez(i,Ny-1,n+1) - Ez(i,Ny,n)); %top boundary
%Ez(i,Ny,n+1)=-Ez(i,Ny-1,n-1)+A*(Ez(i,Ny,n-1)+Ez(i,Ny-1,n+1))+B*(Ez(i,Ny,n)+Ez(i,Ny-1,n))+C*(Ez(i+1,Ny,n)-2*Ez(i,Ny,n)+Ez(i-1,Ny,n)+Ez(i+1,Ny-1,n)-2*Ez(i,Ny-1,n)+Ez(i-1,Ny-1,n));
Ez(i,1,n+1) = Ez(i,2,n) + R*(Ez(i,2,n+1) - Ez(i,1,n)); %bottom boundary
%Ez(i,1,n+1)=-Ez(i,2,n-1)+A*(Ez(i,1,n-1)+Ez(i,2,n+1))+B*(Ez(i,1,n)+Ez(i,2,n))+C*(Ez(i+1,1,n)-2*Ez(i,1,n)+Ez(i-1,1,n)+Ez(i+1,2,n)-2*Ez(i,2,n)+Ez(i-1,2,n));
end
% Boundary conditions for four corner nodes
Ez(1,1,n+1)=.5*(Ez(1,2,n+1)+Ez(2,1,n+1));
Ez(Nx,1,n+1)=.5*(Ez(Nx-1,1,n+1)+Ez(Nx,2,n+1));
Ez(1,Ny,n+1)=.5*(Ez(1,Ny-1,n+1)+Ez(2,Ny,n+1));
Ez(Nx,Ny,n+1)=.5*(Ez(Nx-1,Ny,n+1)+Ez(Nx,Ny-1,n+1));

x=(0:0.05:10);
y=(0:0.05:10);
mesh(x,y,Ez,:,n+1));
axis([0 10 0 10 -1 1]);
pause(0.001);
end

thanks again

Hi,

I am able to run your code in matlab.....I am able to see the mesh plot with wave propagation......I will like to know what you are expecting as output......I am uploading the last Mesh Plot file for your reference...... It is possible to debug your code .....But need more info.....I think the boundary condition that you are applying for top to bottom with index of interation i and from left to right by j.....Please conform the request
with regards,

milind

Here is image of the last mesh plot......I could not able to upload it in last communication

Hi Milind,

Thank you very much for your reply.
Actually I was thinking to email you directly but afraid to bother you so I just post this.
I know the code can be run, but the only problem is that on the left and bottom boundaries it is not absorbing.
My goal is to have the all four boundaries absorbing the wave, not reflecting.
If you plot the contour use contour(Ez,:,n+1)), you'll see there are ripples around, not smoothy curves.

and I think the index is correct, but don't know where is wrong.

---------- Post added at 01:36 ---------- Previous post was at 01:17 ----------


here is attached file for clearance. If the left and bottom boundaries are absorbing, it should be all circles around the source, not some small circles (which is due to the reflection) near the source.

Hi zmcddn,

The Problem that you could not able to feel that right and top boundaries are not absorbing because....I think your source is shifted to 30 30 position in mesh........I had attempted to center source condtion by making line no 26 i.e. if ((i==101)&&(j==101)) then I can see effect of right and top bounday condition too.....



Please check that.....I think the boundary condition that you are applying are right.......only it is a problem of the source position

Hope your problem will get solved after this fix......

regards,

Milind

Dear Milind,

Thank you for your reply.
My goal is to have all four boundaries absoring the wave, regardless of the source position.
If the boundary conditions are correct, then why it is a problem of the source position?
If I want to achieve my goal, should I do any modifications?

Just as a reference, the attached polot below is a correct output (i.e. all boundaries absorbing regardless of the source position)


Thank you again for your help

zmcddn said:

Dear Milind,

Thank you for your reply.
My goal is to have all four boundaries absoring the wave, regardless of the source position.
If the boundary conditions are correct, then why it is a problem of the source position?
If I want to achieve my goal, should I do any modifications?

Just as a reference, the attached polot below is a correct output (i.e. all boundaries absorbing regardless of the source position)


Thank you again for your help

Click to expand...

How do you got the 'correct' plot (analytically, different program perhaps different type of ABC?)

In any case ABCs like 2nd order Mur make a number of simplifying assumptions. The closer a source is to the boundary the greater the error these assumptions produce. PML is much better in this regard but still has problems in some situations.

Essentially, Mur's problems get bigger the smaller the incidence angle. If your beams come in almost parallel to the boundary there is little absorption.

I used a different scheme, but the same ABC. So in the above program, I calculated H and then E. To plot the 'correct' plot, I manually eleminated the H part and use only E field, then discretize it and program it. However, I use exactly the same ABC and no matter where I put the source, the boundaries can absorb the waves very well. So I was so confused about where I was wrong in the above program.

Thank you iyami.

Help!~help!!

zmcddn,

I am still not sure where exactly your problem is (I do not have time to look at your code). Are the 'good' and 'bad' snapshots taken at identical moments?

Also try all four corners (same distance) and compare the output. If one corner has 'better' behaviour I would use a debugger or dumps of all field values near the corners for some steps to check where the difference in behaviour is coming from.

Sorry, i am quite busy now but will try to check this thread over the next few days.

Hi iyami

Thank you for your enthusiasm.

The problem is that the left and bottom boundaries are not absorbing but reflecting. What I want is all four boundaries absorbing.
The first plot shows that the left and bottom boundaries are reflecting (can be seen from the irregular circules). The second plot (the 'good' one) shows that all boundaries are absorbing.
I have tried moving the point source to a mid point that is very close to the left boundary (i.e. i=100, j=30), it clearly shows that the left boundary is reflecting somehow.
I think there is only one small little thing that is wrong but I cannot find it. I'll keep debuging but still looking for help at the same time.

Thank you very much

zmcddn said:

Hi iyami

Thank you for your enthusiasm.

The problem is that the left and bottom boundaries are not absorbing but reflecting. What I want is all four boundaries absorbing.

Click to expand...

First some general hints:
perfect reflection =>
usually some nonzero constant may be mistakenly set to 0
reflection and increase in amplitude (nearly doubling) =>
usually a sign error

But let me repeat the earlier suggestions.
Record the field values when the excitation arrives for the first time. Do it once at the top and once at the bottom. Choose identical distances from the boundary. Compare the field values (sign changes are possible) and look at the first value which is different. The update equations for that field are wrong. Usually that lets you find the error.

Repeat with right and left boundaries.