PROGRAM flow IMPLICIT NONE INTEGER :: i INTEGER :: j INTEGER,PARAMETER :: Nx = 70 INTEGER,PARAMETER :: Ny = 20 INTEGER,PARAMETER :: ix = 25!Position of the obstacle INTEGER,PARAMETER :: iT = 9!x dimension of the obstacle INTEGER,PARAMETER :: jW = 5!y dimension of the obstacle INTEGER,PARAMETER :: M = 200!No. of iterations after which the convergence criterion is checked INTEGER :: im,jm!iteration counter REAL*16 :: psi(0:Nx+1,0:Ny+1) REAL*16 :: zeta(0:Nx+1,0:Ny+1) REAL*16 :: psiold(0:Nx+1,0:Ny+1) REAL*16 :: R = 0.0001q0!lattice Reynolds number REAL*16 :: error = 10.q0!error for the convergence criterion REAL*16,PARAMETER :: precisao = 1.q-8!requested precision for the convergence criterion REAL*16,PARAMETER :: w = 0.3q0!relaxation parameter LOGICAL :: readfromfile OPEN(100,file='flow.dat',status='unknown') write(100,*)'Program flow.f95' write(100,*)' ' write(100,*)'Solves the stationary flow past a rectangular obstacle via the iterative method on a dicretized lattice' write(100,*)'Lattice size Nx =',Nx write(100,*)'Lattice size Ny =',Ny write(100,*)'Obstacle width 2W =',2*jW write(100,*)'Obstacle depth T =',iT write(100,*)'Obstacle distance from the upstream =',ix write(100,*)'Lattice Reynolds number = v0 / viscosity =',R write(100,*)'Physical Reynolds number = 2W v0 / viscosity =',R/real(2*jW) write(100,*)' ' write(100,*)'No. of iterations after which interactions are checked =',M write(100,*)'Relaxation parameter =',w !!INITIAL CONDITION readfromfile = .TRUE. ! readfromfile = .FALSE. if(readfromfile)then OPEN(152,file='stream-field.dat',status='unknown') OPEN(153,file='vortex-field.dat',status='unknown') do i=0,Nx+1 do j=0,Ny+1 read(152,*)im,im,psi(i,j) read(153,*)im,im,zeta(i,j) enddo enddo CLOSE(152) CLOSE(153) else do j=1,Ny psi(:,j) = real(j,16)!initial condition enddo zeta(:,:) = 0.q0!initial condition and boundary condition endif psi(:,0) = 0.q0!Boundary AB and DE zeta(:,0) = 0.q0!Boundary AB and DE psi(:,Ny+1) = 2.q0 + psi(:,Ny-1)!Boundary GH zeta(:,Ny) = 0.q0!Boundary GH psi(0,:) = psi(2,:)!Boundary AH psi(Nx+1,:) = psi(Nx-1,:)!Boundary FG zeta(Nx+1,:) = zeta(Nx-1,:)!Boundary FG do j=0,jW psi(ix,j) = 0.q0 !Boundary BC zeta(ix,j) = -2.q0 * psi(ix-1,j) !Boundary BC psi(ix+iT,j) = 0.q0 !Boundary DE zeta(ix+iT,j) = -2.q0 * psi(ix+1,j) !Boundary DE enddo do i=ix+1,ix+iT-1 psi(i,jW) = 0.q0 !Boundary CD zeta(i,jW) = -2.q0 * psi(i,jW+1) !Boundary CD enddo !!ITERATE THE METHOD R = R / 16.q0!divide by 16 for latter convenience psiold(:,:) = psi(:,:) jm = 0 do while(error.gt.precisao.and.jm*M.lt.60000) do im=1,M!perform the iterations CALL updatepsi(Nx,Ny,ix,iT,jW,w,psi,zeta) CALL updateBCDE(Nx,Ny,ix,iT,jW,psi,zeta) CALL updatezeta(Nx,Ny,ix,iT,jW,w,R,psi,zeta) enddo CALL checkerror(Nx,Ny,ix,iT,jW,psi,psiold,error) jm = jm + 1 enddo !!EXPORT DATA CALL exportdata(Nx,Ny,ix,iT,jW,psi,zeta) write(100,*)' ' write(100,*)'Precision required for stop =',precisao write(100,*)'error =',error write(100,*)'No. of iterations required =',M*jm write(100,*)'Done.' CLOSE(100) write(6,*)'Done.' STOP END PROGRAM flow !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!END OF PROGRAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE updatepsi(Nx,Ny,ix,iT,jW,w,psi,zeta) IMPLICIT NONE INTEGER,INTENT(IN) :: Nx INTEGER,INTENT(IN) :: Ny INTEGER,INTENT(IN) :: ix INTEGER,INTENT(IN) :: iT INTEGER,INTENT(IN) :: jW REAL*16,INTENT(IN) :: w REAL*16,INTENT(INOUT) :: psi(0:Nx+1,0:Ny+1) REAL*16,INTENT(IN) :: zeta(0:Nx+1,0:Ny+1) INTEGER :: i!x-counter INTEGER :: j!y-counter REAL*16 :: aux !!!!!!UPDATE psi ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!! do i=1,ix-1 do j=1,Ny aux = psi(i+1,j) + psi(i-1,j) + psi(i,j+1) + psi(i,j-1) + zeta(i,j) aux = aux / 4.q0 psi(i,j) = psi(i,j)*(1.q0 - w) + w*aux!apply the update enddo enddo !!!!!!UPDATE psi ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!! do i=ix,ix+iT do j=jW+1,Ny aux = psi(i+1,j) + psi(i-1,j) + psi(i,j+1) + psi(i,j-1) + zeta(i,j) aux = aux / 4.q0 psi(i,j) = psi(i,j)*(1.q0 - w) + w*aux!apply the update enddo enddo !!!!!!UPDATE psi ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!!! do i=ix+iT+1,Nx do j=1,Ny aux = psi(i+1,j) + psi(i-1,j) + psi(i,j+1) + psi(i,j-1) + zeta(i,j) aux = aux / 4.q0 psi(i,j) = psi(i,j)*(1.q0 - w) + w*aux!apply the update enddo enddo !!!!!!BOUNDARY CONDITIONS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! psi(0,:) = psi(2,:)!update boundary AH psi(:,Ny+1) = 2.q0 + psi(:,Ny-1)!update boundary GH psi(Nx+1,:) = psi(Nx-1,:)!update boundary FG RETURN END SUBROUTINE updatepsi !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE updatezeta(Nx,Ny,ix,iT,jW,w,R,psi,zeta) IMPLICIT NONE INTEGER,INTENT(IN) :: Nx INTEGER,INTENT(IN) :: Ny INTEGER,INTENT(IN) :: ix INTEGER,INTENT(IN) :: iT INTEGER,INTENT(IN) :: jW REAL*16,INTENT(IN) :: w REAL*16,INTENT(IN) :: R REAL*16,INTENT(IN) :: psi(0:Nx+1,0:Ny+1) REAL*16,INTENT(INOUT) :: zeta(0:Nx+1,0:Ny+1) INTEGER :: i!x-counter INTEGER :: j!y-counter REAL*16 :: aux !!!!!!UPDATE zeta ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!! do i=1,ix-1 do j=1,Ny-1 aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1)) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update enddo enddo !!!!!!UPDATE zeta ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!! !here, we need to be careful because of the ambiguity of the boundary condition at the corners C and D do i=ix+1,ix+iT-1!the corners C and D do not enter in this sum do j=jW+1,Ny-1 aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1)) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update enddo enddo i=ix do j=jW+2,Ny-1!the corners C and D do not enter in this sum aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1)) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update enddo i=ix+iT do j=jW+2,Ny-1!the corners C and D do not enter in this sum aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1)) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update enddo i = ix j = jW+1!the updtate involves corner C aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) + 2.q0*psi(i,j) )!here, we used the value of zeta(corner C) = -2 psi(i,j+1) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update to the site just besides corner C i = ix+iT j = jW+1!the updtate involves corner D aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) + 2.q0*psi(i,j) )!here, we used the value of zeta(corner D) = -2 psi(i,j+1) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update to the site just besides corner D !!!!!!UPDATE zeta ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!! do i=ix+iT+1,Nx do j=1,Ny-1 aux = (psi(i+1,j) - psi(i-1,j)) * (zeta(i,j+1) - zeta(i,j-1)) aux = aux - (psi(i,j+1) - psi(i,j-1)) * (zeta(i+1,j) - zeta(i-1,j)) aux = R*aux aux = aux + ( zeta(i+1,j) + zeta(i-1,j) + zeta(i,j+1) + zeta(i,j-1) ) / 4.q0 zeta(i,j) = zeta(i,j)*(1.q0 - w) + w*aux!apply the update enddo enddo !!!!!!BOUNDARY CONDITIONS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! zeta(Nx+1,:) = zeta(Nx-1,:)!update boundary FG RETURN END SUBROUTINE updatezeta !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE updateBCDE(Nx,Ny,ix,iT,jW,psi,zeta) IMPLICIT NONE INTEGER,INTENT(IN) :: Nx INTEGER,INTENT(IN) :: Ny INTEGER,INTENT(IN) :: ix INTEGER,INTENT(IN) :: iT INTEGER,INTENT(IN) :: jW REAL*16,INTENT(IN) :: psi(0:Nx+1,0:Ny+1) REAL*16,INTENT(INOUT) :: zeta(0:Nx+1,0:Ny+1) INTEGER :: i!x-counter INTEGER :: j!y-counter !!!!!!UPDATE zeta ON BOUNDARIES BC AND DE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do j=1,jW zeta(ix,j) = -2.q0 * psi(ix-1,j) zeta(ix+iT,j) = -2.q0 * psi(ix+1,j) enddo !!!!!!UPDATE zeta ON BOUNDARY CD!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do i=ix+1,ix+iT-1 zeta(i,jW) = -2.q0 * psi(i,jW+1) enddo RETURN END SUBROUTINE updateBCDE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE exportdata(Nx,Ny,ix,iT,jW,psi,zeta) IMPLICIT NONE INTEGER,INTENT(IN) :: Nx INTEGER,INTENT(IN) :: Ny INTEGER,INTENT(IN) :: ix INTEGER,INTENT(IN) :: iT INTEGER,INTENT(IN) :: jW REAL*16,INTENT(IN) :: psi(0:Nx+1,0:Ny+1) REAL*16,INTENT(IN) :: zeta(0:Nx+1,0:Ny+1) REAL*16 :: vx,vy,dpts,aux INTEGER :: i!x-counter INTEGER :: j!y-counter INTEGER :: k,npts npts = 20 dpts = real(npts,16) OPEN(150,file='velocity-field.dat',status='unknown') OPEN(151,file='velocity-field2.dat',status='unknown') !!!!!!EXPORT DATA ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!! do i=1,ix-1 do j=1,Ny vx = ( psi(i,j+1) - psi(i,j-1) ) / 2.q0 vy = ( psi(i-1,j) - psi(i+1,j) ) / 2.q0 write(151,*)i,j,vx,vy do k=0,npts aux = real(k,16)/dpts write(150,*)real(i,16)+vx*aux,real(j,16)+vy*aux enddo do k=0,npts/2 aux = real(k,16)/dpts*2.q0 write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)-0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)+0.2q0*vx*aux write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)+0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)-0.2q0*vx*aux enddo enddo enddo !!!!!!EXPORT DATA ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!! do i=ix,ix+iT do j=jW+1,Ny vx = ( psi(i,j+1) - psi(i,j-1) ) / 2.q0 vy = ( psi(i-1,j) - psi(i+1,j) ) / 2.q0 write(151,*)i,j,vx,vy do k=0,npts aux = real(k,16)/dpts write(150,*)real(i,16)+vx*aux,real(j,16)+vy*aux enddo do k=0,npts/2 aux = real(k,16)/dpts*2.q0 write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)-0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)+0.2q0*vx*aux write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)+0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)-0.2q0*vx*aux enddo enddo enddo !!!!!!EXPORT DATA ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!! do i=ix+iT+1,Nx do j=1,Ny vx = ( psi(i,j+1) - psi(i,j-1) ) / 2.q0 vy = ( psi(i-1,j) - psi(i+1,j) ) / 2.q0 write(151,*)i,j,vx,vy do k=0,npts aux = real(k,16)/dpts write(150,*)real(i,16)+vx*aux,real(j,16)+vy*aux enddo do k=0,npts/2 aux = real(k,16)/dpts*2.q0 write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)-0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)+0.2q0*vx*aux write(150,*)real(i,16)+vx*(1.q0-0.2q0*aux)+0.2q0*vy*aux,real(j,16)+vy*(1.q0-0.2q0*aux)-0.2q0*vx*aux enddo enddo enddo CLOSE(150) CLOSE(151) OPEN(152,file='stream-field.dat',status='unknown') OPEN(153,file='vortex-field.dat',status='unknown') do i=0,Nx+1 do j=0,Ny+1 write(152,*)i,j,psi(i,j) write(153,*)i,j,zeta(i,j) enddo enddo CLOSE(152) CLOSE(153) RETURN END SUBROUTINE exportdata !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE checkerror(Nx,Ny,ix,iT,jW,psi,psiold,error) IMPLICIT NONE INTEGER,INTENT(IN) :: Nx INTEGER,INTENT(IN) :: Ny INTEGER,INTENT(IN) :: ix INTEGER,INTENT(IN) :: iT INTEGER,INTENT(IN) :: jW REAL*16,INTENT(IN) :: psi(0:Nx+1,0:Ny+1) REAL*16,INTENT(INOUT) :: psiold(0:Nx+1,0:Ny+1) REAL*16,INTENT(OUT) :: error REAL*16 :: vx,vy,dpts,aux INTEGER :: i!x-counter INTEGER :: j!y-counter INTEGER :: k,npts error = 0.q0 !!!!!!EXPORT DATA ON REGION 1: BEFORE THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!! do i=1,ix-1 do j=1,Ny error = error + abs( psi(i,j) - psiold(i,j) ) enddo enddo !!!!!!EXPORT DATA ON REGION 2: BESIDES THE OBSTACLE!!!!!!!!!!!!!!!!!!!!! do i=ix,ix+iT do j=jW+1,Ny error = error + abs( psi(i,j) - psiold(i,j) ) enddo enddo !!!!!!EXPORT DATA ON REGION 3: AFTER THE OBSTACLE!!!!!!!!!!!!!!!!!!!!!!! do i=ix+iT+1,Nx do j=1,Ny error = error + abs( psi(i,j) - psiold(i,j) ) enddo enddo error = error / real(Nx*Ny-it*jw,16) psiold(:,:) = psi(:,:) RETURN END SUBROUTINE checkerror !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!