Prev: SQL type of data operators needed in Matlab
Next: how to plot line between any two pixels?
From: hasan on 11 Jun 2010 08:22
Dear Members of Mathworks Forum,
I am working on FDTD 3D simulation in Matlab. But the performance is slower than fortran. I have three for loops which are updating electrical and magnetic fields with looking material type array. I have to make them in vector form but unfortunately i don't. Is there anybody to help me about vectorization.
Thank you.


PS:I don't want to MEX.




for K=2:NZ1,
KK=K;
for J=2:NY1,
JJ=J;
for I=1:NX1,

% DETERMINE MATERIAL TYPE

if (IDONE(I,J,K) == 0)

% FREE SPACE

EXS(I,J,K)=EXS(I,J,K)+(HZS(I,J,K)-HZS(I,J-1,K))*DTEDY-(...
HYS(I,J,K)-HYS(I,J,K-1))*DTEDZ;
elseif (IDONE(I,J,K) == 1)

% PERFECT CONDUCTOR

II = I;
EXS(I,J,K)=-EXI(II,JJ,KK);
else

% LOSSY DIELECTRIC

II = I;
EXS(I,J,K)=EXS(I,J,K)*ESCTC(IDONE(I,J,K))-...
EINCC(IDONE(I,J,K))*exi(II,JJ,KK)-...
EDEVCN(IDONE(I,J,K))*dexi(II,JJ,KK)+(HZS(I,J,K)-...
HZS(I,J-1,K))*ECRLY(IDONE(I,J,K))-(HYS(I,J,K)-...
HYS(I,J,K-1))*ECRLZ(IDONE(I,J,K));
end
end
end
end
From: koubaa kamal on 14 Jun 2010 05:52
try the code bellow, it may not work, because i don't have all the variable defined, but it may need only small changes,

NZ1=10;
NY1=10;
NX1=10;

EXS(1:NX1,1:NY1,1:NZ1)=0;

HZS1=HZS(:,1:end-1,2:end);
HYS1=HYS(:,2:end,1:end-1);
HZS2=HZS(:,2:end,2:end);
HYS2=HYS(:,2:end,2:end);
IDONE2=IDONE(:,2:end,2:end);
EXS2=EXS(:,2:end,2:end);
% DETERMINE MATERIAL TYPE
% FREE SPACE
EXS2(IDONE2 == 0)=EXS2(IDONE2==0)+...
(HZS2(IDONE2 == 0) -HZS1(IDONE2 == 0))*DTEDY-...
(HYS2(IDONE2 == 0,J,K)-HYS1(IDONE2 == 0))*DTEDZ;

% PERFECT CONDUCTOR
EXS2(IDONE2 == 1)=-EXS2(IDONE2== 1);

% LOSSY DIELECTRIC
EXS2((IDONE2~= 1)&(IDONE2~= 0))=...
EXS2((IDONE2~= 1)&(IDONE2~= 0)).*ESCTC((IDONE2~= 1)&(IDONE2~= 0))- ...
EINCC(IDONE2((IDONE2~= 1)&(IDONE2~= 0))).*exi((IDONE2~= 1)&(IDONE2~= 0))- ...
EDEVCN(IDONE2((IDONE2~= 1)&(IDONE2~= 0))).*dexi((IDONE2~= 1)&(IDONE2~= 0))+ ...
(HZS2((IDONE2~= 1)&(IDONE2~= 0))-HZS1((IDONE2~= 1)&(IDONE2~= 0))).*ECRLY(IDONE((IDONE2~= 1)&(IDONE2~= 0)))-...
(HYS2((IDONE2~= 1)&(IDONE2~= 0))-HYS1((IDONE2~= 1)&(IDONE2~= 0))).*ECRLZ(IDONE((IDONE2~= 1)&(IDONE2~= 0)));

EXS(:,2:end,2:end)=EXS2;


"hasan " <albada(a)gmail.com> wrote in message <hut9pd$86v$1(a)fred.mathworks.com>...
> Dear Members of Mathworks Forum,
> I am working on FDTD 3D simulation in Matlab. But the performance is slower than fortran. I have three for loops which are updating electrical and magnetic fields with looking material type array. I have to make them in vector form but unfortunately i don't. Is there anybody to help me about vectorization.
> Thank you.
>
>
> PS:I don't want to MEX.
>
>
>
>
> for K=2:NZ1,
> KK=K;
> for J=2:NY1,
> JJ=J;
> for I=1:NX1,
>
> % DETERMINE MATERIAL TYPE
>
> if (IDONE(I,J,K) == 0)
>
> % FREE SPACE
>
> EXS(I,J,K)=EXS(I,J,K)+(HZS(I,J,K)-HZS(I,J-1,K))*DTEDY-(...
> HYS(I,J,K)-HYS(I,J,K-1))*DTEDZ;
> elseif (IDONE(I,J,K) == 1)
>
> % PERFECT CONDUCTOR
>
> II = I;
> EXS(I,J,K)=-EXI(II,JJ,KK);
> else
>
> % LOSSY DIELECTRIC
>
> II = I;
> EXS(I,J,K)=EXS(I,J,K)*ESCTC(IDONE(I,J,K))-...
> EINCC(IDONE(I,J,K))*exi(II,JJ,KK)-...
> EDEVCN(IDONE(I,J,K))*dexi(II,JJ,KK)+(HZS(I,J,K)-...
> HZS(I,J-1,K))*ECRLY(IDONE(I,J,K))-(HYS(I,J,K)-...
> HYS(I,J,K-1))*ECRLZ(IDONE(I,J,K));
> end
> end
> end
> end
From: hasan on 14 Jun 2010 10:49
Dear Kamal,
Thank you so much. Your response is very useful for me. I need your help again.
In Lossy Dielectric part, the exi is a sub function. When the interpreter calls the exi((IDONE2==2)) it gives an error like

??? Input argument "J" is undefined.

Error in ==> fdtda>exi at 978
DIST=((I-1)*DELX+0.5*DELX*OFF)*XDISP+((J-1)*DELY)*YDISP+((K-1)*DELZ)...

Error in ==> fdtda at 313
EXS2(IDONE2==2)=EXS2(IDONE2==2)*ESCTC(2)-EINCC(2).*...

So,
I need a way that finding I,J,K indices when the IDONE2==2 and send them exi and dexi subfunctions.
My Best Regards,
Hasan


Details of the Program;

I reapplied your answer as the program.
HZS1=HZS(:,1:end-1,2:end);
HYS1=HYS(:,2:end,1:end-1);
HZS2=HZS(:,2:end,2:end);
HYS2=HYS(:,2:end,2:end);
IDONE2=IDONE(:,2:end,2:end);
EXS2=EXS(:,2:end,2:end);

EXS2(IDONE2==0)=EXS2(IDONE2==0)+...(HZS2(IDONE2==0)-HZS1(IDONE2==0))*DTEDY-...
(HYS2(IDONE2==0)-HYS1(IDONE2==0))*DTEDZ;

EXS2(IDONE2 == 1)=-EXS2(IDONE2== 1);

EXS2(IDONE2==2)=EXS2(IDONE2==2)*ESCTC(2)-EINCC(2).*...
exi((IDONE2==2))-EDEVCN(2)*dexi(IDONE2==2)+(HZS2(IDONE2==2)-...
HZS1(IDONE2==2))*ECRLY(2)-(HYS2(IDONE2==2)-HYS1(IDONE2==2)).*...
ECRLZ(2);

EXS(:,2:end,2:end)=EXS2;


Initial Constants and Arrays:
NX=50;
NY=50;
NZ=50;
NX1=NX-1;
NY1=NY-1;
NZ1=NZ-1;
NTEST=4;
NSTOP=40;
DELX=1./17;
DELY=1./17;
DELZ=1./17;
THINC=180;
PHINC=180;
ETHINC=1;
EPHINC=0;
AMP=1000;
BETA=64;
EPS0=8.854E-12;
XMU0=1.2566306E-6;
IDONE=zeros(NX,NY,NZ);
IDTWO=zeros(NX,NY,NZ);
IDTHRE=zeros(NX,NY,NZ);
EXS=zeros(NX,NY,NZ);
EYS=zeros(NX,NY,NZ);
EZS=zeros(NX,NY,NZ);
HXS=zeros(NX,NY,NZ);
HYS=zeros(NX,NY,NZ);
HZS=zeros(NX,NY,NZ);
EYSX1(1:4,1:NY1,1:NZ1)=zeros(4,33,33);
EYSX2(1:4,1:NY1,1:NZ1)=zeros(4,33,33);
EZSX1(1:4,1:NY1,1:NZ1)=zeros(4,33,33);
EZSX2(1:4,1:NY1,1:NZ1)=zeros(4,33,33);
EXSY1(1:NX1,1:4,1:NZ1)=zeros(33,4,33);
EXSY2(1:NX1,1:4,1:NZ1)=zeros(33,4,33);
EZSY1(1:NX1,1:4,1:NZ1)=zeros(33,4,33);
EZSY2(1:NX1,1:4,1:NZ1)=zeros(33,4,33);
EXSZ1(1:NX1,1:NY1,1:4)=zeros(33,33,4);
EXSZ2(1:NX1,1:NY1,1:4)=zeros(33,33,4);
EYSZ1(1:NX1,1:NY1,1:4)=zeros(33,33,4);
EYSZ2(1:NX1,1:NY1,1:4)=zeros(33,33,4);

Subfunctions:
function exi = exi(I,J,K)
DIST=((I-1)*DELX+0.5*DELX*OFF)*XDISP+((J-1)*DELY)*YDISP+((K-1)*DELZ)...
*ZDISP + DELAY;
exi=AMPX*source(DIST);
end
function source = source(DIST)
source=0;
TAU=T-DIST/C;
if (TAU < 0)
elseif(TAU > PERIOD)
else
source=exp(-ALPHA*((TAU-BETADT)^2));
end
end

function dexi = dexi(I,J,K)
DIST=((I-1)*DELX+0.5*DELX*OFF)*XDISP+((J-1)*DELY)*YDISP+((K-1)*DELZ)...
*ZDISP + DELAY;
dexi=AMPX*dsrce(DIST);
end
function source = source(DIST)
source=0;
TAU=T-DIST/C;
if (TAU < 0)
elseif(TAU > PERIOD)
else
source=exp(-ALPHA*((TAU-BETADT)^2));
end
end
function dsrce = dsrce(DIST)
dsrce=0;
TAU=T-DIST/C;
if (TAU < 0)
elseif(TAU > PERIOD)
else
dsrce=exp(-ALPHA*((TAU-BETADT)^2))*(-2*ALPHA*(TAU-BETADT));
end
end
From: koubaa kamal on 14 Jun 2010 12:52
Hello,
Can you write down the code you run with my modifications ?



"hasan " <albada(a)gmail.com> wrote in message <hut9pd$86v$1(a)fred.mathworks.com>...
> Dear Members of Mathworks Forum,
> I am working on FDTD 3D simulation in Matlab. But the performance is slower than fortran. I have three for loops which are updating electrical and magnetic fields with looking material type array. I have to make them in vector form but unfortunately i don't. Is there anybody to help me about vectorization.
> Thank you.
>
>
> PS:I don't want to MEX.
>
>
>
>
> for K=2:NZ1,
> KK=K;
> for J=2:NY1,
> JJ=J;
> for I=1:NX1,
>
> % DETERMINE MATERIAL TYPE
>
> if (IDONE(I,J,K) == 0)
>
> % FREE SPACE
>
> EXS(I,J,K)=EXS(I,J,K)+(HZS(I,J,K)-HZS(I,J-1,K))*DTEDY-(...
> HYS(I,J,K)-HYS(I,J,K-1))*DTEDZ;
> elseif (IDONE(I,J,K) == 1)
>
> % PERFECT CONDUCTOR
>
> II = I;
> EXS(I,J,K)=-EXI(II,JJ,KK);
> else
>
> % LOSSY DIELECTRIC
>
> II = I;
> EXS(I,J,K)=EXS(I,J,K)*ESCTC(IDONE(I,J,K))-...
> EINCC(IDONE(I,J,K))*exi(II,JJ,KK)-...
> EDEVCN(IDONE(I,J,K))*dexi(II,JJ,KK)+(HZS(I,J,K)-...
> HZS(I,J-1,K))*ECRLY(IDONE(I,J,K))-(HYS(I,J,K)-...
> HYS(I,J,K-1))*ECRLZ(IDONE(I,J,K));
> end
> end
> end
> end
From: hasan on 14 Jun 2010 13:26
Of course,
Code is too long for this area. I send it the web. You can see it in below url.
http://www.zumodrive.com/share/5TcPNWU5OT
 |  Next  |  Last
Pages: 1 2
Prev: SQL type of data operators needed in Matlab
Next: how to plot line between any two pixels?