Welcome! Log In Create A New Profile

Advanced

Error in code

Posted by AtulBhagatAtulBhagat  
Error in code
January 26, 2017 07:18AM
Dear all ,
I am new to this forum, and want to get help to correct the code I have written,
This is very simple lid driven cavity code using D2Q9 model.
The Zou and He boundary condition is used at moving lid. The compilation of code is fine but after running the
code I am getting nan values, which is due to reduction in density at moving lid,
but I have coded according to Zou and He paper but still getting error.
Can anybody please help me to site the error in it.
Thanks in advance



Edited 1 time(s). Last edit at 01/26/2017 09:44AM by AtulBhagat.
Re: Error in code
January 26, 2017 09:44AM
Bug found!
code corrected.
Re: Error in code
February 03, 2017 07:59AM
Dear Atulbhagat,

You are probably dividing by 0 somewhere in the code.
But if you can show the code we might be able to fix it.

Regards
Githin
Re: Error in code
March 12, 2017 01:48PM
i am working on a project fdlbm(finite difference latiice boltezmam )
and i need a code for lid cavity
and i need it very much
best regards
mahmoud
Re: Error in code
March 12, 2017 04:34PM
Hi mahmoud,
This is the code I have written rather say converted from fortran to c. This code you can get in A.A.Mohammad Book.
Hope this helps
Language: C++
#include <cstdio> #include <cstdlib> #include <cstring> #include <cmath>   #include <time.h>   const int nx = 100, ny = 100, q = 9; double csq, w0 = 4.0/9.0, w1 = 1.0/9.0, w2 = 1.0/36.0; double ***f, ***feq, **rho, **u, **v, **strf; double dx, dt, alpha, omega, rhon, u0 = 0.1, v0 = 0.0; double R, U, V, u2, v2, uv, s2, uv2, uvm, uvm2; FILE *data2d;   // Memory allocation void memoryallocation() { feq = new double **[q]; f = new double **[q];   for (int i = 0; i < q; i++) { feq[i] = new double *[nx]; f[i] = new double *[nx];   for (int j = 0; j < nx; j++) { feq[i][j] = new double [ny]; memset (feq[i][j], 0.0, (ny)*sizeof(double));   f[i][j] = new double [ny]; memset (f[i][j], 0.0, (ny)*sizeof(double)); } }   rho = new double *[nx]; u = new double *[nx]; v = new double *[nx]; strf = new double *[nx];   for (int i = 0; i < nx; i++) { rho[i] = new double [ny]; memset (rho[i], 0.0, (ny)*sizeof(double)); u[i] = new double [ny]; memset (u[i], 0.0, (ny)*sizeof(double)); v[i] = new double [ny]; memset (v[i], 0.0, (ny)*sizeof(double)); strf[i] = new double [ny]; memset (strf[i], 0.0, (ny)*sizeof(double)); } }   // Deletion of Memory void deletememory() { for (int i = 0; i < q; i++) { for (int j = 0; j < nx; j++) { delete[] feq[i][j]; delete[] f[i][j]; } delete[] feq[i]; delete[] f[i]; }   for (int i = 0; i < nx; i++) { delete[] rho[i]; delete[] u[i]; delete[] v[i]; delete[] strf[i]; } }   // Equilibrium state void equilibrium() { for(int i = 0; i < nx; i++) { for(int j = 0; j < ny; j++) { U = u[i][j]; V = v[i][j]; R = rho[i][j]; u2 = U*U; v2 = V*V; s2 = -1.5*(u2+v2); uv = U+V; uv2 = (U+V)*(U+V); uvm = U-V; uvm2 = (U-V)*(U-V); feq[0][i][j] = w0*R*(1. + s2); feq[1][i][j] = w1*R*(1. + 3.*U + 4.5*u2 + s2); feq[2][i][j] = w1*R*(1. + 3.*V + 4.5*v2 + s2); feq[3][i][j] = w1*R*(1. - 3.*U + 4.5*u2 + s2); feq[4][i][j] = w1*R*(1. - 3.*V + 4.5*v2 + s2); feq[5][i][j] = w2*R*(1. + 3.*uv + 4.5*uv2 + s2); feq[6][i][j] = w2*R*(1. - 3.*uvm + 4.5*uvm2 + s2); feq[7][i][j] = w2*R*(1. - 3.*uv + 4.5*uv2 + s2); feq[8][i][j] = w2*R*(1. + 3.*uvm + 4.5*uvm2 + s2); } } }   // Initial conditions void initialstate() { for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) { rho[i][j] = 5.0; strf[i][j] = 0.0; } }   equilibrium();   for(int k = 0; k < q; k++) { for(int i = 0; i < nx; i++) { for(int j = 0; j < ny; j++) { f[k][i][j] = feq[k][i][j]; } } } }   // Macro-variable calculations void macrovar() { for(int i = 0; i < nx; i++) { for(int j = 0; j < ny; j++) { rho[i][j] = 0.0; for(int k = 0; k < q; k++) { rho[i][j] += f[k][i][j]; } } } for(int i = 0; i < nx; i++) { for(int j = 0; j < ny; j++) { u[i][j] = (f[1][i][j] - f[3][i][j] + f[5][i][j] - f[6][i][j] - f[7][i][j] + f[8][i][j])/rho[i][j]; v[i][j] = (f[5][i][j] + f[2][i][j] + f[6][i][j] -f[7][i][j] -f[4][i][j] - f[8][i][j])/rho[i][j]; } } } // collision void collision() { for(int k = 0; k < q; k++) { for(int i = 0; i < nx; i++) { for(int j = 0; j < ny; j++) { f[k][i][j] = f[k][i][j]*(1.0-omega) + omega*feq[k][i][j]; } } } }   //streaming void streaming() { for(int i = (nx-1); i > 0; i--) { for(int j = (ny-1); j > 0; j--) { f[1][i][j] = f[1][i-1][j]; f[5][i][j] = f[5][i-1][j-1]; } }   for(int i = 0; i < (nx-1); i++) { for(int j = 0; j < (ny-1); j++) { f[3][i][j] = f[3][i+1][j]; f[7][i][j] = f[7][i+1][j+1]; } }   for(int j = (ny-1); j > 0; j--) { for(int i = 0; i < (nx-1); i++) { f[2][i][j] = f[2][i][j-1]; f[6][i][j] = f[6][i+1][j-1]; } }   for(int i = (nx-1); i > 0; i--) { for(int j = 0; j < (ny-1); j++) { f[4][i][j] = f[4][i][j+1]; f[8][i][j] = f[8][i-1][j+1]; } } }   // boundary condition void boundarycondition() { for(int j = 0; j < ny; j++) { // left layer i.e. west boundary   f[1][0][j] = f[3][0][j]; f[5][0][j] = f[7][0][j]; f[8][0][j] = f[6][0][j];   // right layer i.e. east boundary   f[3][nx-1][j] = f[1][nx-1][j]; f[6][nx-1][j] = f[8][nx-1][j]; f[7][nx-1][j] = f[5][nx-1][j]; }   // bottom layer i.e. south boundary for(int i = 0; i < nx; i++) { f[2][i][0] = f[4][i][0]; f[5][i][0] = f[7][i][0]; f[6][i][0] = f[8][i][0]; }   // top layer i.e. north boundary (a moving plate/lid) for(int i = 1; i < (nx-1); i++) { rhon = (1.0/(1.0+v0))*(f[0][i][ny-1]+f[1][i][ny-1]+f[3][i][ny-1]+2.0*(f[2][i][ny-1]+f[6][i][ny-1]+f[5][i][ny-1])); f[4][i][ny-1] = f[2][i][ny-1] - ((2.0/3.0)*rhon*v0); f[7][i][ny-1] = f[5][i][ny-1]-(rhon*v0/6.0)-((1.0/2.0)*rhon*u0); f[8][i][ny-1] = f[6][i][ny-1]-(rhon*v0/6.0)+((1.0/2.0)*rhon*u0); } }   int main() { // Time calculations clock_t t; t = clock();   // Variable declaration dt = 1.0; dx = 1.0; alpha = 0.01;   // Parameter calculation csq = dx*dx/(dt*dt); omega = 1.0/(3.*alpha/(dt*csq)+0.5);   memoryallocation(); initialstate();   int nstep = 200000;   // main loop of LBM for (int m = 0; m < nstep; m++) { equilibrium(); collision(); streaming(); boundarycondition(); macrovar(); if(m % 100 == 0) { printf("nsteps = %d , U(nx/2,ny/2) = %e\n", m, u[nx/2][ny/2]); } }   data2d = fopen("lid_mid_plane_output.txt", "w"); for (int i = 0; i < nx; i++) { fprintf(data2d, "%i\t%f\t%f\t%f\n", i, rho[i][ny/2], u[i][ny/2], v[i][ny/2]); //X,y,rho,u,v; } fclose(data2d);     for (int i = 1; i < nx; i++) { double rhoav = 0.5*(rho[i-1][0]+rho[i][0]); if (i != 0) strf[i][0] = strf[i-1][0] - rhoav*0.5*(v[i-1][0]+v[i][0]);   for (int j = 1; j < ny; j++) { double rhom = 0.5*(rho[i][j]+rho[i][j-1]); strf[i][j] = strf[i][j-1] + rhom*0.5*(u[i][j-1]+u[i][j]); } } data2d = fopen("stremfun_output.txt", "w"); for(int j = 0; j < ny; j++) { for (int i = 0; i < nx; i++) { fprintf(data2d, "%f\t", strf[i][j]); } fprintf(data2d, "\n"); } fclose(data2d);   deletememory();   t = clock() - t; double time_taken = ((double)t)/CLOCKS_PER_SEC; printf ("Time of execution = %f\n", time_taken);   return 0; }
Re: Error in code
May 30, 2017 01:22PM
hi AtulBhagat
sorry for replying late,thanks for the code
but the code i need is finite difference lbm code
which have not streaming in that method
once again very very thanks for the code
if you have a code for linearly heating natural convection cavity
i run the heating natural convection cavity at A.A.Mohammad Book. and it work fine
and when i modify it to run linearly heating natural convection cavity it did not work fine
and i do not why
i changed the boundary and i think i changed write
6.6.1 Example: Natural Convection in a Differentially Heated Cavity page 98 and it run with me good
and here is my code
! this code for linealy heated cavity with all wall at static with linealy left heated wall and adiabatic at north
! and heated at bottom and low temp at right



[code="fortran"parameter (n=100,m=100)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
real g(0:8,0:n,0:m), geq(0:8,0:n,0:m),th(0:n,0:m)
integer i
open(2,file='uvfield')
open(3,file='uvely')
open(4,file='tmx')

cx(:)=(/0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,-1.0,1.0/)
cy(:)=(/0.0,0.0,1.0,0.0,-1.0,1.0,1.0,-1.0,-1.0/)
w(:)=(/4./9.,1./9.,1./9.,1./9.,1./9.,1./36.,1./36.,1./36.,1./36./)
uo=0.0
sumvelo=0.0
rhoo=6.00
dx=1.0
dy=dx
dt=1.0
tw=1.0
th=0.0
ra=1.0e3
pr=0.7
visco=0.02
alpha=visco/pr
pr=visco/alpha
gbeta=ra*visco*alpha/(float(m*m*m))
Re=uo*m/alpha
print *, "Re=", Re
omega=1.0/(3.*visco+0.5)
omegat=1.0/(3.*alpha+0.5)
mstep=100000
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do i=0,n
u(i,m)=uo
v(i,m)=0.0
end do
! main loop
1 do kk=1,mstep
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m,th,gbeta)
call streaming(f,n,m)
call bounceb(f,n,m)
call rhouv(f,rho,u,v,cx,cy,n,m)
! ——————————–
! collestion for scalar
call collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
! streaming for scalar
call streaming(g,n,m)
call gbound(g,tw,w,n,m)
do j=0,m
do i=0,n
sumt=0.0
do k=0,8
sumt=sumt+g(k,i,j)
end do
th(i,j)=sumt
end do
end do
!print *, th(n/2,m/2),v(5,m/2),rho(0,m/2),u(n/2,m/2),v(n/2,m/2),rho(n,m/2)
END DO
! end of the main loop
call result(u,v,rho,th,uo,n,m,ra)
stop
end
! end of the main program
subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m,th,gbeta)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
real th(0:n,0:m)
tref=0.50
DO i=0,n
DO j=0,m
t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
force=3.*w(k)*gbeta*(th(i,j)-tref)*cy(k)*rho(i,j)
if(i.eq.0.or.i.eq.n) force =0.0
if(j.eq.0.or.j.eq.m) force =0.0
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)+force
END DO
END DO
END DO
return
end
subroutine collt(u,v,g,geq,th,omegat,w,cx,cy,n,m)
real g(0:8,0:n,0:m),geq(0:8,0:n,0:m),th(0:n,0:m)
real w(0:8),cx(0:8),cy(0:8)
real u(0:n,0:m),v(0:n,0:m)
do i=0,n
do j=0,m
do k=0,8
geq(k,i,j)=th(i,j)*w(k)*(1.0+3.0*(u(i,j)*cx(k)+v(i,j)*cy(k)))
g(k,i,j)=omegat*geq(k,i,j)+(1.0-omegat)*g(k,i,j)
end do
end do
end do
return
end
subroutine streaming(f,n,m)
real f(0:8,0:n,0:m)
! streaming
DO j=0,m
DO i=n,1,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=0,n-1 !LEFT TO RIGHT
f(3,i,j)=f(3,i+1,j)
END DO
END DO
DO j=m,1,-1 !TOP TO BOTTOM
DO i=0,n
f(2,i,j)=f(2,i,j-1)
END DO
DO i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO i=0,n-1
f(6,i,j)=f(6,i+1,j-1)
END DO
END DO
DO j=0,m-1 !BOTTOM TO TOP
DO i=0,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,1,-1
f(8,i,j)=f(8,i-1,j+1)
END DO
END DO
return
end
subroutine bounceb(f,n,m)
real f(0:8,0:n,0:m)
do j=0,m
!west boundary
f(1,0,j)=f(3,0,j)
f(5,0,j)=f(7,0,j)
f(8,0,j)=f(6,0,j)
!east boundary
f(3,n,j)=f(1,n,j)
f(7,n,j)=f(5,n,j)
f(6,n,j)=f(8,n,j)
end do
do i=0,n
!south boundary
f(2,i,0)=f(4,i,0)
f(5,i,0)=f(7,i,0)
f(6,i,0)=f(8,i,0)
!north boundary
f(4,i,m)=f(2,i,m)
f(8,i,m)=f(6,i,m)
f(7,i,m)=f(5,i,m)
end do
return
end
subroutine gbound(g,tw,w,n,m)
real g(0:8,0:n,0:m)
real w(0:8),tw
! Boundary conditions
! West boundary condition, T=1.
do j=0,m
g(1,0,j)=tw*(1-(j/m))*(w(1)+w(3))-g(3,0,j)
g(5,0,j)=tw*(1-(j/m))*(w(5)+w(7))-g(7,0,j)
g(8,0,j)=tw*(1-(j/m))*(w(8)+w(6))-g(6,0,j)
!g(1,0,j)=tw*((j/m))*(w(1)+w(3))-g(3,0,j)
!g(5,0,j)=tw*((j/m))*(w(5)+w(7))-g(7,0,j)
!g(8,0,j)=tw*((j/m))*(w(8)+w(6))-g(6,0,j)

! g(1,0,j)=tw*(w(1)+w(3))-g(3,0,j)
!g(5,0,j)=tw*(w(5)+w(7))-g(7,0,j)
!g(8,0,j)=tw*(w(8)+w(6))-g(6,0,j)
end do
! East boundary condition, T=0.
do j=0,m
g(6,n,j)=-g(8,n,j)
g(3,n,j)=-g(1,n,j)
g(7,n,j)=-g(5,n,j)
end do
! Top boundary conditions, Adiabatic
do i=0,n
!g(8,i,m)=g(8,i,m-1)
!g(7,i,m)=g(7,i,m-1)
!g(6,i,m)=g(6,i,m-1)
!g(5,i,m)=g(5,i,m-1)
!g(4,i,m)=g(4,i,m-1)
!g(3,i,m)=g(3,i,m-1)
!g(2,i,m)=g(2,i,m-1)
!g(1,i,m)=g(1,i,m-1)
!g(0,i,m)=g(0,i,m-1)
g(4,i,m)=g(2,i,m)
g(8,i,m)=g(6,i,m)
g(7,i,m)=g(5,i,m)

end do
!Bottom boundary conditions, Adiabatic
do i=0,n
!g(1,i,0)=g(1,i,1)
!g(2,i,0)=g(2,i,1)
!g(3,i,0)=g(3,i,1)
!g(4,i,0)=g(4,i,1)
!g(5,i,0)=g(5,i,1)
!g(6,i,0)=g(6,i,1)
!g(7,i,0)=g(7,i,1)
!g(8,i,0)=g(8,i,1)
!g(0,i,0)=g(0,i,1)
g(2,i,0)=tw*(w(2)+w(4))-g(4,i,0)
g(5,i,0)=tw*(w(5)+w(7))-g(7,i,0)
g(6,i,0)=tw*(w(8)+w(6))-g(8,i,0)



end do
return
end
subroutine tcalcu(g,th,n,m)
real g(0:8,0:n,0:m),th(0:n,0:m)
do j=0,m
do i=0,n
ssumt=0.0
do k=0,8
ssumt=ssumt+g(k,i,j)
end do
th(i,j)=ssumt
end do
end do
return
end
subroutine rhouv(f,rho,u,v,cx,cy,n,m)
real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m),cx(0:8),cy(0:8)
do j=0,m
do i=0,n
ssum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
end do
end do
DO i=0,n
DO j=0,m
usum=0.0
vsum=0.0
DO k=0,8
usum=usum+f(k,i,j)*cx(k)
vsum=vsum+f(k,i,j)*cy(k)
END DO
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
END DO
END DO
return
end
subroutine result(u,v,rho,th,uo,n,m,ra)
real u(0:n,0:m),v(0:n,0:m),th(0:n,0:m)
real strf(0:n,0:m),rho(0:n,0:m)
open(5,file='streamt')
open(6,file='nuav')
! streamfunction calculations
strf(0,0)=0.
do i=0,n
rhoav=0.5*(rho(i-1,0)+rho(i,0))
if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
do j=1,m
rhom=0.5*(rho(i,j)+rho(i,j-1))
strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
end do
end do
! ———————————–
write(2,*)"VARIABLES =X, Y, U, V"
write(2,*)"ZONE ","I=",n+1,"J=",m+1,",","F=BLOCK"
do j=0,m
write(2,*)(i/float(m),i=0,n)
end do
do j=0,m
write(2,*)(j/float(m),i=0,n)
end do
do j=0,m
write(2,*)(u(i,j),i=0,n)
end do
do j=0,m
write(2,*)(v(i,j),i=0,n)
end do
do j=0,m
!write(3,*)j/float(m),u(5,j)/uo,u(n/2,j)/uo,u(n-10,j)/uo
end do
do i=0,n
write(4,*) i/float(n),th(i,m/2)
end do
write(5,*)"VARIABLES =X, Y, S, T"
write(5,*)"ZONE ","I=",n+1,"J=",m+1,",","F=BLOCK"
do j=0,m
write(5,*)(i/float(m),i=0,n)
end do
do j=0,m
write(5,*)(j/float(m),i=0,n)
end do
do j=0,m
write(5,*)(strf(i,j),i=0,n)
end do
do j=0,m
write(5,*)(th(i,j),i=0,n)
end do
! Nusselt number Calculation
snul=0.0
snur=0.0
do j=0,m
rnul=(th(0,j)-th(1,j))*float(n)
rnur=(th(n-1,j)-th(n,j))*float(n)
snul=snul+rnul
snur=snur+rnur
write(5,*)j/float(m),rnul,rnur
end do
avnl=snul/float(m)
avnr=snur/float(m)
write(6,*)ra,avnl,avnr

OPEN(21,FILE='DATx.m')
OPEN(22,FILE='the.m')
OPEN(23,FILE='DATu.m')
OPEN(24,FILE='DATv.m')
OPEN(25,FILE='steam.m')
REWIND(21)
REWIND(22)
REWIND(23)
REWIND(24)
REWIND(25)
WRITE(21,*)'x = [...'
WRITE(22,*)'the = [...'
WRITE(23,*)'u = [...'
WRITE(24,*)'v = [...'
WRITE(25,*)'stream = [...'
DO j=0,m
WRITE(21,*)(i,i=0,n)

write(22,100)(th(i,j),i=0,n)

WRITE(23,100)(u(i,j),i=0,n)
WRITE(24,100)(v(i,j),i=0,n)
WRITE(25,100)(strf(i,j),i=0,n)
ENDDO
WRITE(21,*)'] ;'
WRITE(22,*)'] ;'
WRITE(23,*)'] ;'
WRITE(24,*)'] ;'
WRITE(25,*)'] ;'

100 FORMAT(1X,256(E12.6,2X))






return
end]
please correct it if you can
[/code]
Sorry, you do not have permission to post/reply in this forum.