Unit Conversion

zarbaya · December 20, 2012, 12:28am

Hello TIMM,
I have this code in fortran, and I like that you look at it, my physical problem is I work in 3D, I have a flow in a cube 1m * 1m * 1m, I inject pressure in inlet and outlet (1.1 and 1 unit LB) the problem is that with a mesh 50 * 50 * 50 calculation works, but if I increase the mesh for example 100 * 100 * 100 calculation becomes unstable.
thank you for giving me your comments about my code (the steps of calculation, how can I find the value viteese in physical units, why if we increase the mesh becomes unstable calculates!?)

MODULE geomConst
integer, parameter:: fluid = 1, solide = 0, inlet = 10, outlet = 11
END MODULE geomConst

! ========================================================
! Les constantes de lattice boltzmann
! ========================================================
MODULE D3Q19Const

double precision,parameter:: t(19) = (/1.0d0/18.0d0, 1.0d0/18.0d0, 1.0d0/18.0d0,1.0d0/18.0d0, 1.0d0/18.0d0&
                                       &,1.0d0/18.0d0, 1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0, 1.0d0/36.0d0&
	                           &,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0, 1.0d0/36.0d0&
		                   &, 1.0d0/36.0d0, 1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/3.0d0/)

! D3Q19 Directions
integer, parameter:: cx(19) = (/1,-1,0,0,0,0,1,-1,-1,1,1,1,-1,-1,0,0,0,0,0/)
integer, parameter:: cy(19) = (/0,0,1,-1,0,0,1,1,-1,-1,0,0,0,0,1,1,-1,-1,0/)
integer, parameter:: cz(19) = (/0,0,0,0,1,-1,0,0,0,0,1,-1,-1,1,1,-1,-1,1,0/)
!

integer, parameter:: opposite(19) = (/2,1,4,3,6,5,9,10,7,8,13,14,11,12,17,18,15,16,19/)

END MODULE D3Q19Const

! ========================================================
! Les constantes de la simulation
! ========================================================
MODULE simParam
integer, parameter:: nx = 100
integer, parameter:: ny = 100
integer, parameter:: nz = 100
integer, parameter:: sphX = (nx+1)/2
integer, parameter:: sphY = (ny+1)/2
integer, parameter:: sphZ = (nz+1)/2
integer, parameter:: sphR = (ny+1)/10
integer, parameter:: tMax = 10000
double precision, parameter:: tho = 0.75d0

END MODULE simParam

! ========================================================
! main
! ========================================================

PROGRAM LBM
USE geomConst
USE simParam
USE D3Q19Const
implicit none
integer :: i, j,k,x, y, z, x0,y0,z0
double precision:: fTmp(nx,ny,nz,19)
double precision:: dx,dy,dz,dt,hx,hy,hz,xx, yy, zz
double precision::sd,s,er,tt
double precision:: t1, t2, temps
double precision, dimension(:,:,:,:), allocatable::f, fEq
double precision, dimension(:,:,:), allocatable:: rho, u,ux,uy,uz,uux,uuy,uuz
integer, dimension(:,:,:), allocatable:: geom
double precision:: cu, ccu
integer:: tStep

allocate(f(nx,ny,nz,19))  
allocate(fEq(nx,ny,nz,19))  
allocate(ux(nx,ny,nz))
allocate(uy(nx,ny,nz))
allocate(uz(nx,ny,nz))
allocate(u(nx,ny,nz))
allocate(uux(nx,ny,nz))
allocate(uuy(nx,ny,nz))
allocate(uuz(nx,ny,nz))
allocate(rho(nx,ny,nz))
allocate(geom(nx,ny,nz))

dx=1.d0/nx
dy=1.0d0/ny
dz=1.0d0/nz
dt=dx**2 * (tho-0.5d0)/(3*nu)
ccu=dx/dt

fTmp=0.0d0

geom = fluid
geom(:,:,1) = inlet
geom(:,:,nz) = outlet
geom(1,:, = solide
geom(nx,:, = solide
geom(:,1, = solide
geom(:,ny, = solide

 uux(:,:,:) = 0.0d0
 uuy(:,:,:) = 0.0d0
 uuz(:,:,:) = 0.0d0
 ux(:,:,:) = 0.0d0
 uy(:,:,:) = 0.0d0
 uz(:,:,:) = 0.0d0

rho  = 1.0d0
u = ux(:,:,:) * ux(:,:,:) + uy(:,:,:) * uy(:,:,:) + uz(:,:,:) * uz(:,:,:)

do i = 1, 19
do z = 1, nz
do y = 1, ny
do x=1, nx
cu = 3*( ux(x,y,z) * cx(i) + uy(x,y,z) * cy(i)+ uz(x,y,z) * cz(i))
fEq(x,y,z,i) = t(i) * rho(x,y,z) * (1.0d0 + cu + 0.5d0 * cu * cu - 1.5d0 * u(x,y,z))
end do
end do
end do
end do

f = fEq

temps = 0.0d0

do tStep = 1,tMax

       CALL CPU_TIME(t1)

! ------------------------------------
! Streaming
! ------------------------------------

f(:,:,:,1) = f((/nx, (i, i= 1, nx-1)/),:,:,1)

! -------------------------------------
! -------------------------------------
!

f(:,:,:,2) = f((/(i, i= 2, nx), 1/),:,:,2)

! -------------------------------------
f(:,:,:,3) = f(:,(/ny, (i, i= 1, ny-1)/),:,3)

!-------------------------------------
!

f(:,:,:,4) = f(:,(/(i, i= 2,ny), 1/),:,4)

!---------------------------------------------------
f(:,:,:,5) = f(:, :, (/nz, (i, i= 1, nz-1)/), 5)

!-------------------------------------

f(: ,:,:,6) = f(:,:,(/(i, i= 2, nz), 1/),6)

! -------------------------------------

f(:,:,:,7) = f((/nx, (i, i= 1, nx-1)/),(/ny, (j, j= 1, ny-1)/), :, 7)

!----------------------------------------
f(:,:,:,8)= f((/(i, i= 2, nx), 1/),(/ny, (j, j= 1, ny-1)/),:,8)

!------------

f(:,:,:,9)=f((/(i, i= 2, nx), 1/),(/(j, j= 2, ny), 1/),:,9)


f(:,:,:,10)=f((/nx, (i, i= 1, nx-1)/),(/(j, j= 2, ny), 1/),:,10)

!-------------
f(:,:,:,11)=f((/nx, (i, i= 1, nx-1)/),:,(/nz, (j, j= 1, nz-1)/), 11)

f(:,:,:,12)=f((/nx, (i, i= 1, nx-1)/),:,(/(j, j= 2, nz), 1/), 12)

!-----------
!------------

f(:,:,:,13)=f((/(i, i= 2, nx), 1/),:,(/(j, j= 2, nz), 1/), 13)

!-------------

f(:,:,:,14)=f((/(i, i= 2, nx), 1/),:,(/nz, (j, j= 1, nz-1)/), 14)

f(:,:,:,15)=f(:, (/ny, (i, i= 1,ny-1)/), (/nz, (j, j= 1, nz-1)/), 15)

!------------

f(:,:,:,16)=f(:,(/ny, (i, i= 1,ny-1)/),(/(j, j= 2, nz), 1/),16)

!-----------

f(:,:,:,17)=f(:,(/(i, i= 2, ny), 1/),(/(j, j= 2, nz), 1/),17)
!----------

f(:,:,:,18)=f(:,(/(i, i= 2, ny), 1/),(/nz, (j, j= 1, nz-1)/),18)

!do i=1,19
!write(,) f(2,2,1,17)
!write(,) f(1,2,1,2)
!pause;
!end do

!--------------------------
! Valeur Macro: ux,uy,uz,u
!-------------------------

do z = 1, nz
do y = 1, ny
do x=1, nx
if(geom(x,y,z)/=solide) then

 rho(x,y,z)=  f(x,y,z,1) + f(x,y,z,2) + f(x,y,z,3) + f(x,y,z,4) + f(x,y,z,5) + f(x,y,z,6) &
+ f(x,y,z,7) + f(x,y,z,8)+f(x,y,z,9)+f(x,y,z,10)+f(x,y,z,11)+f(x,y,z,12)+f(x,y,z,13)+f(x,y,z,14)+f(x,y,z,15)&
+f(x,y,z,16)+f(x,y,z,17)+f(x,y,z,18)+f(x,y,z,19) 

 
 ux(x,y,z)  = (f(x,y,z,1) - f(x,y,z,2) + f(x,y,z,7) - f(x,y,z,8) - f(x,y,z,9) + f(x,y,z,10)+&
		 f(x,y,z,11)+f(x,y,z,12)- f(x,y,z,13)- f(x,y,z,14)) / rho(x,y,z)
		
 uy(x,y,z)  = (f(x,y,z,3) - f(x,y,z,4) + f(x,y,z,7) + f(x,y,z,8) - f(x,y,z,9) - f(x,y,z,10)+&
		 f(x,y,z,15)+f(x,y,z,16)- f(x,y,z,17)- f(x,y,z,18)) / rho(x,y,z)

 uz(x,y,z)  = (f(x,y,z,5) - f(x,y,z,6) + f(x,y,z,11) - f(x,y,z,12) - f(x,y,z,13) + f(x,y,z,14)+&
		 f(x,y,z,15)-f(x,y,z,16)- f(x,y,z,17)+ f(x,y,z,18)) / rho(x,y,z)

! u(x,y,z) = sqrt(ux(x,y,z) * ux(x,y,z) + uy(x,y,z) * uy(x,y,z)+ uz(x,y,z) * uz(x,y,z))
 
   
           end if
       end do  
     end do
    end do

!---------------------
! Inlet - outlet
!---------------------

do z = 1, nz
do y = 1, ny
do x=1, nx
if (geom(x,y,z) == inlet) then
ux(x,y,z) = 0.0d0
uy(x,y,z) = 0.0d0
rho(x,y,z)= 1.10d0

            uz(x,y,z)  = 1-(f(x,y,z,4)+f(x,y,z,1)+f(x,y,z,2)+f(x,y,z,3)+f(x,y,z,10)+f(x,y,z,7)+f(x,y,z,8)+ &
                 f(x,y,z,9) +f(x,y,z,19)+2*(f(x,y,z,6)+f(x,y,z,13)+f(x,y,z,12)+f(x,y,z,17)+f(x,y,z,16)))/rho(x,y,z)

!rho(x,y,z) = (f(x,y,z,4)+f(x,y,z,1)+f(x,y,z,2)+f(x,y,z,3)+f(x,y,z,10)+f(x,y,z,7)+f(x,y,z,8)+ &
! f(x,y,z,9) +f(x,y,z,19)+2*(f(x,y,z,6)+f(x,y,z,13)+f(x,y,z,12)+f(x,y,z,17)+f(x,y,z,16)))/(1-uz(x,y,z))

            f(x,y,z,5)  = f(x,y,z,6)  + 1/3.0d0 * rho(x,y,z) * uz(x,y,z) 
            f(x,y,z,15) = f(x,y,z,17) - 1/4.0d0 * (f(x,y,z,3)-f(x,y,z,4)) + 1/6.0d0 * rho(x,y,z) * uz(x,y,z)
            f(x,y,z,18) = f(x,y,z,16) + 1/4.0d0 * (f(x,y,z,3)-f(x,y,z,4)) + 1/6.0d0 * rho(x,y,z) * uz(x,y,z)
            f(x,y,z,11) = f(x,y,z,13) - 1/4.0d0 * (f(x,y,z,1)-f(x,y,z,2)) + 1/6.0d0 * rho(x,y,z) * uz(x,y,z)
            f(x,y,z,14) = f(x,y,z,12) + 1/4.0d0 * (f(x,y,z,1)-f(x,y,z,2)) + 1/6.0d0 * rho(x,y,z) * uz(x,y,z)
 

        else if (geom(x,y,z) == outlet) then

            ux(x,y,z)  = 0.0d0
            uy(x,y,z)  = 0.0d0
            rho(x,y,z) = 1.0d0

            uz(x,y,z)  = -1+(f(x,y,z,4)+f(x,y,z,1)+f(x,y,z,2)+f(x,y,z,3)+f(x,y,z,10)+f(x,y,z,7)+f(x,y,z,8)+ &
                          f(x,y,z,9) +f(x,y,z,19) + 2 * (f(x,y,z,5)+f(x,y,z,11)+f(x,y,z,14)+f(x,y,z,15)+f(x,y,z,18)))/rho(x,y,z)

             f(x,y,z,6)  = f(x,y,z,5)  -  1/3.0d0 * rho(x,y,z) * uz(x,y,z)
             f(x,y,z,17) = f(x,y,z,15) +  1/4.0d0 * (f(x,y,z,3)-f(x,y,z,4)) - 1/6.0d0 * rho(x,y,z)*uz(x,y,z)
             f(x,y,z,16) = f(x,y,z,18) -  1/4.0d0 * (f(x,y,z,3)-f(x,y,z,4)) - 1/6.0d0 * rho(x,y,z)*uz(x,y,z)
             f(x,y,z,13) = f(x,y,z,11) +  1/4.0d0 * (f(x,y,z,1)-f(x,y,z,2)) - 1/6.0d0 * rho(x,y,z)*uz(x,y,z)
             f(x,y,z,12) = f(x,y,z,14) -  1/4.0d0 * (f(x,y,z,1)-f(x,y,z,2)) - 1/6.0d0 * rho(x,y,z)*uz(x,y,z)  
        end if
     end do
    end do
end do

! write(,) ‘vz’, uz(2,2,1)
!pause;

!--------------------
! Colision
!--------------------

do i = 1, 19
do z = 1, nz
do y = 1, ny
do x=1, nx
if(geom(x,y,z)/=solide) then
cu = 3*( ux(x,y,z) * cx(i) + uy(x,y,z) * cy(i)+ uz(x,y,z) * cz(i))
fEq(x,y,z,i) = t(i) * rho(x,y,z) * (1.0d0 + cu + 0.5d0 * cu * cu - &
1.5d0 * (ux(x,y,z) * ux(x,y,z) + uy(x,y,z) * uy(x,y,z)+ uz(x,y,z) *uz(x,y,z)))

            f(x,y,z,i) = f(x,y,z,i) - 1/tho * (f(x,y,z,i) - fEq(x,y,z,i))
           end if
        end do
    end do
end do

end do

!----------------
!Cl : Bounce-back
!----------------

do z = 1, nz
do y = 1, ny
do x=1, nx
if (geom(x,y,z) == solide) then

        tt         = f(x,y,z,1) 
        f(x,y,z,1) = f(x,y,z,2)   
        f(x,y,z,2) = tt

        tt         = f(x,y,z,3) 
        f(x,y,z,3) = f(x,y,z,4)   
        f(x,y,z,4) = tt
        
         tt         = f(x,y,z,5) 
        f(x,y,z,5) = f(x,y,z,6)   
        f(x,y,z,6) = tt

         tt         = f(x,y,z,7) 
        f(x,y,z,7) = f(x,y,z,9)   
        f(x,y,z,9) = tt

         tt         = f(x,y,z,8) 
        f(x,y,z,8) = f(x,y,z,10)   
        f(x,y,z,10) = tt

         tt         = f(x,y,z,11) 
        f(x,y,z,11) = f(x,y,z,13)   
        f(x,y,z,13) = tt

         tt         = f(x,y,z,12) 
        f(x,y,z,12) = f(x,y,z,14)   
        f(x,y,z,14) = tt

         tt         = f(x,y,z,15) 
        f(x,y,z,15) = f(x,y,z,17)   
        f(x,y,z,17) = tt

         tt         = f(x,y,z,16) 
        f(x,y,z,16) = f(x,y,z,18)   
        f(x,y,z,18) = tt
           
        end if 
      end do
    end do
end do

do z = 1, nz
do y = 1, ny
do x=1, nx

u(x,y,z) = sqrt(ux(x,y,z) * ux(x,y,z) + uy(x,y,z) * uy(x,y,z)+ uz(x,y,z) * uz(x,y,z))

end do
end do
end do

!-------------------------
! teste de Convergence
!-------------------------

sd=0.0d0
s =0.0d0

uux=ux-uux
uuy=uy-uuy
uuz=uz-uuz

do z = 1, nz
do y = 1, ny
do x=1, nx
sd=sd+sqrt(uux(x,y,z)**2 + uuy(x,y,z)**2 + uuz(x,y,z)**2)
s=s+sqrt(ux(x,y,z) * ux(x,y,z) + uy(x,y,z) * uy(x,y,z)+ uz(x,y,z) * uz(x,y,z))
end do
end do
end do
write(,) sd/s
if(sd/s < 1.0d-4) then
write(,) ’ nbre des iterations =’, tStep
write(,) 'temps de calcul= ', temps
exit
else
uux=ux
uuy=uy
uuz=uz
end if
CALL CPU_TIME(t2)
temps = temps + (t2-t1)

end do

open(unit=17,file=‘suuu.txt’,status=‘unknown’)
do z=1, nz
do y=1, ny
do x=1, nx
write(17,) x, y, z, sqrt(ux(x,y,z) * ux(x,y,z) + uy(x,y,z) * uy(x,y,z) + uz(x,y,z) * uz(x,y,z))
if (x==nx) write(17,)
end do
end do
end do

close(unit=17)

open(unit=18,file=‘50uzz1.txt’,status=‘unknown’)
do k=1, nz
zz=k*dz

write(18,*) zz, ccu*uz(k,50,50)

end do

close(unit=18)

open(unit=19,file=‘50uzz2.txt’,status=‘unknown’)
do k=1, nz
zz=k*dz

write(19,*) zz, ccu*uz(50,k,50)

end do

close(unit=19)

open(unit=21,file=‘geo.txt’,status=‘unknown’)
do k=1, nz
do j = 1, ny
do i = 1, nx
xx=idx
yy=jdy
zz=kdz
write(21,) xx,yy,zz, geom(i,j,k)
end do
end do
end do
close(unit=21)

deallocate(f)
deallocate(fEq)
deallocate(ux)
deallocate(uy)
deallocate(uz)
deallocate(uux)
deallocate(uuy)
deallocate(uuz)
deallocate(u)
deallocate(rho)
deallocate(geom)

END PROGRAM LBM

Timm · January 2, 2013, 8:10pm

Hello,

I am sorry, but I won’t be able to take a look at your code. You should not expect other people to debug your code.
Regarding the unit conversions: Please have a look at this document. You should go through the examples and try to understand how the unit conversion works. I am pretty sure that the instability is caused by a flawed conversion of units.

Timm