Error in this code.

Dear all, I have a code as attaches file. This code will be worked OK, but when i modified the coordinate of xc, yc, it mean when coordinate locate on the Domain (boundary) or outside the Domain, the code has error. Someone know that, please help me.

Thanks all.

%======================================================================== %domain size and physical variables Lx=1.0;Ly=1.0;gx=0.0;gy=-100.0; rho1=1.0; rho2=2.0; m0=0.01; rro=rho1; unorth=0;usouth=0;veast=0;vwest=0;time=0.0; rad=0.15;xc=0.5;yc=0.7; % Initial drop size and location % Numerical variables nx=32;ny=32;dt=0.00125;nstep=100; maxit=200;maxError=0.001;beta=1.2; % Zero various arrys u=zeros(nx+1,ny+2); v=zeros(nx+2,ny+1); p=zeros(nx+2,ny+2); ut=zeros(nx+1,ny+2); vt=zeros(nx+2,ny+1); tmp1=zeros(nx+2,ny+2); uu=zeros(nx+1,ny+1); vv=zeros(nx+1,ny+1); tmp2=zeros(nx+2,ny+2); % Set the grid dx=Lx/nx;dy=Ly/ny; for i=1:nx+2; x(i)=dx*(i-1.5);end; for j=1:ny+2; y(j)=dy*(j-1.5);end; % Set density in the domain and the drop r=zeros(nx+2,ny+2)+rho1; for i=2:nx+1,for j=2:ny+1; if ( (x(i)-xc)^2+(y(j)-yc)^2 < rad^2), r(i,j)=rho2; end, end,end %================== SETUP THE FRONT =================== Nf=100; xf=zeros(1,Nf+2);yf=zeros(1,Nf+2); uf=zeros(1,Nf+2);vf=zeros(1,Nf+2); tx=zeros(1,Nf+2);ty=zeros(1,Nf+2); for l=1:Nf+2, xf(l)=xc-rad*sin(2.0*pi*(l-1)/(Nf)); yf(l)=yc+rad*cos(2.0*pi*(l-1)/(Nf));end %================== START TIME LOOP====================================== for is=1:nstep fx=zeros(nx+2,ny+2);fy=zeros(nx+2,ny+2); % Set fx & fy to zero %---------------------------------------------------------% tangential velocity at boundaries u(1:nx+1,1)=2*usouth-u(1:nx+1,2);u(1:nx+1,ny+2)=2*unorth-u(1:nx+1,ny+1); v(1,1:ny+1)=2*vwest-v(2,1:ny+1);v(nx+2,1:ny+1)=2*veast-v(nx+1,1:ny+1); for i=2:nx,for j=2:ny+1 % TEMPORARY u-velocity ut(i,j)=u(i,j)+dt*(-0.25*(((u(i+1,j)+u(i,j))^2-(u(i,j)+ ... u(i-1,j))^2)/dx+((u(i,j+1)+u(i,j))*(v(i+1,j)+ ... v(i,j))-(u(i,j)+u(i,j-1))*(v(i+1,j-1)+v(i,j-1)))/dy)+ ... (m0/(0.5*(r(i+1,j)+r(i,j))) )*( ... (u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2+ ... (u(i,j+1)-2*u(i,j)+u(i,j-1))/dy^2 )+ gx); end,end for i=2:nx+1,for j=2:ny % TEMPORARY v-velocity vt(i,j)=v(i,j)+dt*(-0.25*(((u(i,j+1)+u(i,j))*(v(i+1,j)+ ... v(i,j))-(u(i-1,j+1)+u(i-1,j))*(v(i,j)+v(i-1,j)))/dx+ ... ((v(i,j+1)+v(i,j))^2-(v(i,j)+v(i,j-1))^2)/dy)+ ... (m0/(0.5*(r(i,j+1)+r(i,j))) )*( ... (v(i+1,j)-2*v(i,j)+v(i-1,j))/dx^2+ ... (v(i,j+1)-2*v(i,j)+v(i,j-1))/dy^2 )+ gy); end,end %======================================================================== % Compute source term and the coefficient for p(i,j) rt=r; lrg=1000; rt(1:nx+2,1)=lrg;rt(1:nx+2,ny+2)=lrg; rt(1,1:ny+2)=lrg;rt(nx+2,1:ny+2)=lrg; for i=2:nx+1,for j=2:ny+1 tmp1(i,j)= (0.5/dt)*( (ut(i,j)-ut(i-1,j))/dx+(vt(i,j)-vt(i,j-1))/dy ); tmp2(i,j)=1.0/( (1./dx)*( 1./(dx*(rt(i+1,j)+rt(i,j)))+... 1./(dx*(rt(i-1,j)+rt(i,j))) )+... (1./dy)*(1./(dy*(rt(i,j+1)+rt(i,j)))+... 1./(dy*(rt(i,j-1)+rt(i,j))) ) ); end,end for it=1:maxit % SOLVE FOR PRESSURE oldArray=p; for i=2:nx+1,for j=2:ny+1 p(i,j)=(1.0-beta)*p(i,j)+beta* tmp2(i,j)*(... (1./dx)*( p(i+1,j)/(dx*(rt(i+1,j)+rt(i,j)))+... p(i-1,j)/(dx*(rt(i-1,j)+rt(i,j))) )+... (1./dy)*( p(i,j+1)/(dy*(rt(i,j+1)+rt(i,j)))+... p(i,j-1)/(dy*(rt(i,j-1)+rt(i,j))) ) - tmp1(i,j)); end,end if max(max(abs(oldArray-p))) <maxError, break,end end for i=2:nx,for j=2:ny+1 % CORRECT THE u-velocity u(i,j)=ut(i,j)-dt*(2.0/dx)*(p(i+1,j)-p(i,j))/(r(i+1,j)+r(i,j)); end,end for i=2:nx+1,for j=2:ny % CORRECT THE v-velocity v(i,j)=vt(i,j)-dt*(2.0/dy)*(p(i,j+1)-p(i,j))/(r(i,j+1)+r(i,j)); end,end %================== ADVECT FRONT ===================== for l=2:Nf+1 ip=floor(xf(l)/dx)+1; jp=floor((yf(l)+0.5*dy)/dy)+1; ax=xf(l)/dx-ip;ay=(yf(l)+0.5*dy)/dy-jp; uf(l)=(1.0-ax)*(1.0-ay)*u(ip,jp)+ax*(1.0-ay)*u(ip+1,jp)+... (1.0-ax)*ay*u(ip,jp+1)+ax*ay*u(ip+1,jp+1); ip=floor((xf(l)+0.5*dx)/dx)+1; jp=floor(yf(l)/dy)+1; ax=(xf(l)+0.5*dx)/dx-ip+1;ay=yf(l)/dy-jp+1; vf(l)=(1.0-ax)*(1.0-ay)*v(ip,jp)+ax*(1.0-ay)*v(ip+1,jp)+... (1.0-ax)*ay*v(ip,jp+1)+ax*ay*v(ip+1,jp+1); end for i=2:Nf+1, xf(i)=xf(i)+dt*uf(i); yf(i)=yf(i)+dt*vf(i);end %MOVE THE FRONT xf(1)=xf(Nf+1);yf(1)=yf(Nf+1);xf(Nf+2)=xf(2);yf(Nf+2)=yf(2); %------------ Add points to the front ------------ xfold=xf;yfold=yf; j=1; for l=2:Nf+1 ds=sqrt( ((xfold(l)-xf(j))/dx)^2 + ((yfold(l)-yf(j))/dy)^2); if (ds > 0.5) j=j+1;xf(j)=0.5*(xfold(l)+xf(j-1));yf(j)=0.5*(yfold(l)+yf(j-1)); j=j+1;xf(j)=xfold(l);yf(j)=yfold(l); elseif (ds < 0.25) % DO NOTHING! else j=j+1;xf(j)=xfold(l);yf(j)=yfold(l); end end Nf=j-1; xf(1)=xf(Nf+1);yf(1)=yf(Nf+1);xf(Nf+2)=xf(2);yf(Nf+2)=yf(2); %------------ distribute gradient -------------- fx=zeros(nx+2,ny+2);fy=zeros(nx+2,ny+2); % Set fx & fy to zero for l=2:Nf+1 nfx=-0.5*(yf(l+1)-yf(l-1))*(rho2-rho1); nfy=0.5*(xf(l+1)-xf(l-1))*(rho2-rho1); % Normal vector ip=floor(xf(l)/dx)+1; jp=floor((yf(l)+0.5*dy)/dy)+1; ax=xf(l)/dx-ip+1; ay=(yf(l)+0.5*dy)/dy-jp+1; fx(ip,jp) =fx(ip,jp)+(1.0-ax)*(1.0-ay)*nfx/dx/dy; fx(ip+1,jp) =fx(ip+1,jp)+ax*(1.0-ay)*nfx/dx/dy; fx(ip,jp+1) =fx(ip,jp+1)+(1.0-ax)*ay*nfx/dx/dy; fx(ip+1,jp+1)=fx(ip+1,jp+1)+ax*ay*nfx/dx/dy; ip=floor((xf(l)+0.5*dx)/dx)+1; jp=floor(yf(l)/dy)+1; ax=(xf(l)+0.5*dx)/dx-ip+1; ay=yf(l)/dy-jp+1; fy(ip,jp) =fy(ip,jp)+(1.0-ax)*(1.0-ay)*nfy/dx/dy; fy(ip+1,jp) =fy(ip+1,jp)+ax*(1.0-ay)*nfy/dx/dy; fy(ip,jp+1) =fy(ip,jp+1)+(1.0-ax)*ay*nfy/dx/dy; fy(ip+1,jp+1)=fy(ip+1,jp+1)+ax*ay*nfy/dx/dy; end %------------ construct the density -------------- for iter=1:maxit oldArray=r; for i=2:nx+1,for j=2:ny+1 r(i,j)=0.25*(r(i+1,j)+r(i-1,j)+r(i,j+1)+r(i,j-1)+... dx*fx(i-1,j)-dx*fx(i,j)+... dy*fy(i,j-1)-dy*fy(i,j)); end,end if max(max(abs(oldArray-r))) <maxError, break,end %======================================================================== time=time+dt % plot the results uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,2:ny+2)+u(1:nx+1,1:ny+1)); vv(1:nx+1,1:ny+1)=0.5*(v(2:nx+2,1:ny+1)+v(1:nx+1,1:ny+1)); for i=1:nx+1,xh(i)=dx*(i-1);end; for j=1:ny+1,yh(j)=dy*(j-1);end hold off,contourf(x,y,flipud(rot90(r))),axis equal,axis([0 Lx 0 Ly]); hold on;quiver(xh,yh,flipud(rot90(uu)),flipud(rot90(vv)),'r'); plot(xf(1:Nf+2),yf(1:Nf+2),'k','linewidth',5);pause(0.01) end

[Jan]

Please format your code properly. Follow the "? Help" link to learn the basics about this forum. Thanks.

"The code has error" does not allow to understand the problems and "modify the coordinate of xc, yc" is not a clear description of what you exactly do. Answering to these to details would demand for a lot of wild guessing.

[Ut]

Thanks you for your help, the coordinate at here it mean i will change yc=0, or xc=0,ect, the code has the error. This code here already working.

Thanks

[Jan]

@Ut: Again: Please format your code. It is not readable yet. Concentrate on the relevant part of the code and post the error message also.

Answer 1

[Jan]

The standard method to solve such problems is to investigate the core of the problem at first. For this job two things are important:

1. Raed the error message carefully. Matlab's messages are the best I've ever seen and usually the point in the right direction for a solution already.

2. Use the debugger: Either set a breakpoint into the code and step through the program line by line. Or let Matlab stop, when the error occurs by typing this in the command line:

dbstop if error

Then start the program again and when it stops, check the sizes and values of the locally used variables.

If these two steps did not help to solve the problem, the forum can be useful. Show us the error message, explain as much detail as needed and explain what you want to achieve.

[Ut]

Thank so much for your help, it would be possible, could you help me for this code.

[Jan]

@Ut: I still do not understand the problem and I still cannot read your code.