program fivol2 include 'fivin' integer n c c fivol applies the finite volume method to Laplaces Equation c in cartesian coordinates on a polar grid. c the disrcretized equation is solved by SOR c call pread call pwrite call grid call startcon call bdcon call matrixcoef n = 0 do while( (n .lt. niter) .and. (rms .gt. eps) ) call sor n = n+1 end do if (rms .gt. eps) write(16,15) niter,rms if (rms .le. eps) write(16,16) n,rms 15 format(/,/,' Convergence NOT achieved in',i5,' steps',5x, + ' rms=',e12.5) 16 format(/,/,' Convergence achieved in',i5,' steps',5x, + ' rms=',e12.5) call check write(14) n,rms,om close(14) c stop end c ------------ : ------------ subroutine pread include 'fivin' c input data: c parameters explained in input file sim1.dat open(11,file='fivol.dat') open(16,file='fivol.out') read(11,1) niter read(11,2) rw,rx,ry,rz,theb,then,eps,om 1 format(i8) 2 format(e10.3) close(11) return end c ------------ : ------------ subroutine pwrite include 'fivin' open(14,file='fivol.bin',form='unformatted') write(14) nr,ntheta write(16,3) write(16,4) nr,ntheta,niter,eps,om write(16,5) rw,rx,ry,rz,theb,then 3 format(' Laplace Equation by Finite Volume method', //) 4 format(' nr=', i2,' ntheta=', i2,' niter=', i5, 5x, + ' eps=', e10.3, ' om=', f5.3) 5 format(' rw=', f5.3,' rx=', f5.3,' ry=', f5.3, + ' rz=', f5.3, 5x,' theb=', f5.1,' then=', f5.1,//) return end c ------------ : ------------ subroutine grid include 'fivin' integer j, k real thetm1, drwx, drzy, thrad, cosk, sink thetm1 = ntheta-1. drwx = (rx-rw)/(nr-1.) drzy = (ry-rz)/thetm1 dtheta = (then-theb)/thetm1 c c set x, y, exact and initial phi c do 7 k = 1,ntheta thrad = (theb + (k-1.)*dtheta)*pi/180. cosk = cos(thrad) sink = sin(thrad) dr(k) = drwx + (drzy-drwx)*(k-1.)/thetm1 rwz = rw + (rz-rw)*(k-1.)/thetm1 do 6 j = 1, nr r(j,k) = rwz + (j-1.)*dr(k) x(j,k) = r(j,k)*cosk y(j,k) = r(j,k)*sink 6 continue c write(6,98) k, thrad, (x(j,k),j=1,nr) c write(6,99) (y(j,k),j=1,nr) 7 continue write(14) x,y c 98 format(' k, theta, x: ',i2,7f7.3) c 99 format(' y: ',6f7.3) return end c ------------ : ------------ subroutine startcon include 'fivin' integer j, k real thrad, sink rms=2*eps do 7 k = 1,ntheta thrad = (theb + (k-1.)*dtheta)*pi/180. sink = sin(thrad) do 6 j = 1, nr phix(j,k) = sink/r(j,k) phi(j,k) = phix(j,k) 6 continue c write(6,97) (phix(j,k),j=1,nr) 7 continue c 97 format(' phix: ',6f7.3) return end c ------------ : ------------ subroutine bdcon include 'fivin' integer j, k c c set boundary values of phi c c write(16,16) do 8 j = 1,nr phi(j,1) = 0. phi(j,ntheta) = phix(j,ntheta) 8 continue do 9 k = 1,ntheta phi(1,k) = phix(1,k) phi(nr,k) = phix(nr,k) c write(16,17) k c write(16,19) (phi(j,k),j=1,nr) 9 continue c 16 format(/,' Initial Configuration:') c 17 format(/,' k=', i2) c 19 format(' phi=', 10f8.4) return end c ------------ : ------------ subroutine matrixcoef include 'fivin' integer j, k, jm, jp, km, kp c c set grid related parameters c do 11 k = 2, ntheta-1 km = k-1 kp = k+1 do 10 j = 2,nr-1 jm = j-1 jp = j+1 xa = 0.25*( x(j,k)+x(jm,k)+x(jm,km)+x(j,km) ) ya = 0.25*( y(j,k)+y(jm,k)+y(jm,km)+y(j,km) ) xb = 0.25*( x(j,k)+x(j,km)+x(jp,km)+x(jp,k) ) yb = 0.25*( y(j,k)+y(j,km)+y(jp,km)+y(jp,k) ) xc = 0.25*( x(j,k)+x(jp,k)+x(jp,kp)+x(j,kp) ) yc = 0.25*( y(j,k)+y(jp,k)+y(jp,kp)+y(j,kp) ) xd = 0.25*( x(j,k)+x(j,kp)+x(jm,kp)+x(jm,k) ) yd = 0.25*( y(j,k)+y(j,kp)+y(jm,kp)+y(jm,k) ) c c side AB dxa = xb - xa dya = yb - ya dxk = x(j,k) - x(j,k-1) dyk = y(j,k) - y(j,k-1) sab = abs( dxa*dyk-dxk*dya ) qab(j,k) = (dxa*dxa+dya*dya)/sab pab(j,k) = (dxa*dxk+dya*dyk)/sab c c side BC dxb = xc - xb dyb = yc - yb dxj = x(j,k) - x(j+1,k) dyj = y(j,k) - y(j+1,k) sbc = abs( dyj*dxb-dxj*dyb ) qbc(j,k) = (dxb*dxb+dyb*dyb)/sbc pbc(j,k) = (dxb*dxj+dyb*dyj)/sbc c c side CD dxc = xd - xc dyc = yd - yc dxk = x(j,k) - x(j,k+1) dyk = y(j,k) - y(j,k+1) scd = abs( dxc*dyk-dxk*dyc ) qcd(j,k) = (dxc*dxc+dyc*dyc)/scd pcd(j,k) = (dxc*dxk+dyc*dyk)/scd c c side DA dxd = xa - xd dyd = ya - yd dxj = x(j,k) - x(j-1,k) dyj = y(j,k) - y(j-1,k) sda = abs( dxj*dyd-dyj*dxd ) qda(j,k) = (dxd*dxd+dyd*dyd)/sda pda(j,k) = (dxd*dxj+dyd*dyj)/sda 10 continue 11 continue return end c ------------ : ------------ subroutine sor include 'fivin' integer j, k, jm, jp, km, kp real sum,phd,dif c iteration using sor sum = 0. do 13 k = 2,ntheta-1 km = k-1 kp = k+1 do 12 j = 2,nr-1 jm = j-1 jp = j+1 phd = 0.25*( pcd(j,k)-pda(j,k) ) *phi(jm,kp) phd = phd + (qcd(j,k)+0.25*( pbc(j,k)-pda(j,k) ))*phi(j,kp) phd = phd + 0.25*( pbc(j,k)-pcd(j,k) ) *phi(jp,kp) phd = phd + (qda(j,k)+0.25*( pcd(j,k)-pab(j,k) ))*phi(jm,k) phd = phd + (qbc(j,k)+0.25*( pab(j,k)-pcd(j,k) ))*phi(jp,k) phd = phd + 0.25*( pda(j,k)-pab(j,k) ) *phi(jm,km) phd = phd + (qab(j,k)+0.25*( pda(j,k)-pbc(j,k) ))*phi(j,km) phd = phd + 0.25*( pab(j,k)-pbc(j,k) ) *phi(jp,km) phd = phd / ( qab(j,k)+qbc(j,k)+qcd(j,k)+qda(j,k) ) dif = phd - phi(j,k) sum = sum + dif*dif phi(j,k) = phi(j,k) + om*dif 12 continue 13 continue rms = sqrt(sum/(nr-2.)/(ntheta-2.)) return end c ------------ : ------------ subroutine check include 'fivin' integer j, k real sum,dif c compare with exact solution write(16,16) sum = 0. do 21 k = 1,ntheta write(16,17) k do 18 j = 1,nr dif = phi(j,k) - phix(j,k) sum = sum + dif*dif 18 continue write(16,19) (phi(j,k),j=1,nr) write(16,20) (phix(j,k),j=1,nr) 21 continue rms = sqrt(sum/(nr-2.)/(ntheta-2.)) write(16,22) rms write(14) phix,phi 16 format(/,' Final Configuration and Exact Solution:') 17 format(/,' k=', i2) 19 format(' phi=', 10f8.4) 20 format(' phix=', 10f8.4) 22 format(/,' Error compared to exact solution: rms=',e12.5) return end c ------------ : ------------