program newburger include 'newbin' c c newburger applies newton's method to solve burger's equations c for u(x,y) and v(x,y) c exbur - exact solution of burgers equations c resbur evaluates the residuals c jacbur evaluates the jacobian c fact factorises the jacobian in l.u c solve solves the linear system for dt c call pread call pwrite call grid call startcon it = 0 itp = 0 c 0. residuals first time call resbur do while( (it .lt. itmx) .and. (rms .gt. eps) + .and. (rms .lt. 1.0e+03) ) if (it .ge. itp) itp = itp + ipn it = it+1 c 1. jacobian call jacbur c 2. factorise jacobian call fact c 3. solve for the correction, dt call solve c 4. increment u and v call newuv c 5. residuals again call resbur end do c 6. final output call fwrite stop end c ------------ : ------------ subroutine pread include 'newbin' c read program parameters open(11,file='newb.dat') open(16,file='newb.out') open(17,file='newb.sto') read(11,1) itmx,ird,ipn read(11,2) eps,re,om,dt read(11,3) (a(j),j=1,na), al 1 format(i5) 2 format(e10.3) 3 format(f8.2) c return end c ------------ : ------------ subroutine pwrite include 'newbin' c write program parameters open(14,file='newb.bin',form='unformatted') write(14) nx,ny write(16,4) naj,itmx,ird,ipn write(16,5) dt,eps,re,om write(16,6) al,(a(j),j=1,na) 4 format(1x,'newtons method for n =',i3,4x,' itmx=', i3, + ' ird= ',i2,' ipn= ',i2) 5 format(1x,'dt =',e12.5,' eps=',e10.3,' re=',e10.3,' om=',e10.3) 6 format(1x,' al=',f5.2,' a =', 5f8.2,//) return end c ------------ : ------------ subroutine grid include 'newbin' open(16,file='newb.out') xmin = -1. xmax = 1. ymin = 0. ymax = pi/6./al dx = (xmax-xmin)/(nx-1) dy = (ymax-ymin)/(ny-1) do 1 j = 1,nx 1 x(j) = xmin + (j-1)*dx do 2 k = 1,ny 2 y(k) = ymin + (k-1)*dy write(14) x,y return end c ------------ : ------------ subroutine startcon include 'newbin' integer j, k call exbur rms=2*eps c ird = 2 -> init sol from newb.sto if (ird .eq. 2) then do 7 k = 1,ny read(17,101) (u(k,j),j=1,nx) read(17,101) (v(k,j),j=1,nx) 7 continue else do 10 j = 1,nx do 10 k = 1,ny u(k,j) = ue(k,j) v(k,j) = ve(k,j) 10 continue do 11 j = 2,nx-1 do 11 k = 2,ny-1 u(k,j) = 0.0 v(k,j) = 0.0 11 continue endif do 12 k = 1,ny 12 write(16,103)(ue(k,j),j=1,nx) do 14 k = 1,ny 14 write(16,105)(ve(k,j),j=1,nx) 101 format(10f8.5) 103 format(' ue=',7f10.4) 105 format(' ve=',7f10.4) return end c ------------ : ------------ subroutine exbur include 'newbin' c exact solution of burger's equations using the c cole hopf transformation real ph(nx,ny),xd(nx),yd(ny),dex(nx),ddx(nx) x0 = 1.0 do 1 j = 1,nx 1 xd(j) = al*(x(j) - x0) do 2 j = 1,nx dex(j) = exp(xd(j)) + exp(-xd(j)) ddx(j) = exp(xd(j)) - exp(-xd(j)) 2 continue do 3 k = 1,ny 3 yd(k) = al*y(k) do 4 k = 1,ny sy = sin(yd(k)) cy = cos(yd(k)) do 4 j = 1,nx ph(k,j) = a(1) + a(2)*x(j) + a(3)*y(k) + + a(4)*x(j)*y(k) + a(5)*dex(j)*cy phx = a(2) + a(4)*y(k) + a(5)*al*ddx(j)*cy phy = a(3) + a(4)*x(j) - a(5)*al*dex(j)*sy ue(k,j) = - 2.*phx/ph(k,j)/re ve(k,j) = - 2.*phy/ph(k,j)/re 4 continue return end c ------------ : ------------ subroutine resbur include 'newbin' c residuals of steady burger's equations real*8 cx,cy,ccx,ccy,dum,dsqrt c cx = 0.5/dx cy = 0.5/dy ccx = 1./re/dx/dx ccy = 1./re/dy/dy nct = 1 do 1 j = 2,nxp jm = j-1 jp = j+1 do 1 k = 2,nyp km = k-1 kp = k+1 dum = cx*u(k,j)*( u(k,jp)-u(k,jm) ) + + cy*v(k,j)*( u(kp,j)-u(km,j) ) r(nct) = dum - ccx*( u(k,jm)-2.*u(k,j)+u(k,jp) ) + - ccy*( u(km,j)-2.*u(k,j)+u(kp,j) ) dum = cx*u(k,j)*( v(k,jp)-v(k,jm) ) + + cy*v(k,j)*( v(kp,j)-v(km,j) ) r(nct+1) = dum - ccx*( v(k,jm)-2.*v(k,j)+v(k,jp) ) + - ccy*( v(km,j)-2.*v(k,j)+v(kp,j) ) nct = nct+2 1 continue c write(16,43)(r(j),j=1,n) sum = 0. do 17 i = 1,naj rd(i) = r(i) 17 sum = sum + r(i)*r(i) 18 rms = dsqrt(sum/naj) if (it .ge. itp) write(16,44) it,rms,(r(j),j=1,3) 43 format(' r=',8f8.4) 44 format(' it=',i3,' rms=',d11.4,' r=',3d11.4) return end c ------------ : ------------ subroutine jacbur include 'newbin' c jacobian of steady burger's equations real dti real*8 cx,cy,ccx,ccy dti = 1./dt n = (nx-2)*(ny-2)*2 do 1 j = 1,n do 1 k = 1,n aj(k,j) = 0.0 1 continue cx = 0.5/dx cy = 0.5/dy ccx = 1./re/dx/dx ccy = 1./re/dy/dy do 7 j = 2,nxp jm = j-1 jp = j+1 do 6 k = 2,nyp km = k-1 kp = k+1 mb = 2*(k-2) + 2*(ny-2)*(j-2) lu = mb+1 lv = mb+2 ml = mb - 2*(ny-2) if (ml .lt. 0) goto 3 aj(lu,ml+1) = -cx*u(k,j) - ccx aj(lv,ml+2) = -cx*u(k,j) - ccx 3 aj(lu,mb+1) = 2.*(ccx+ccy) + cx*( u(k,jp)-u(k,jm) ) + dti aj(lu,mb+2) = cy*( u(kp,j)-u(km,j) ) aj(lv,mb+1) = cx*( v(k,jp)-v(k,jm) ) aj(lv,mb+2) = 2.*(ccx+ccy) + cy*( v(kp,j)-v(km,j) ) + dti mr = mb+2*(ny-2) if (mr .gt. n-2) goto 4 aj(lu,mr+1) = cx*u(k,j) - ccx aj(lv,mr+2) = cx*u(k,j) - ccx 4 mk = mb-2 if (mk .lt. 0) goto 5 aj(lu,mk+1) = - cy*v(k,j) - ccy aj(lv,mk+2) = - cy*v(k,j) - ccy 5 mt = mb+2 if (mt .gt. n-2) goto 6 aj(lu,mt+1) = cy*v(k,j) - ccy aj(lu,mt+2) = cy*v(k,j) - ccy 6 continue 7 continue c do 12 j = 1,n c 12 write(16,13)(aj(k,j),k=1,n) c 13 format(' a=',18f7.1) return end c ------------ : ------------ subroutine fact include 'newbin' c factorises aj into permuted l.u so that perm*a=l*u c jpvt(k) gives index of kth pivot row c sets jpvt(n) = -1 if zero pivot occurs c integer i,j,k,l real*8 s c gauss elimination with partial pivoting do 5 k = 1,n-1 c select pivot l = k do 1 i = k+1,n 1 if ( abs(aj(i,k)) .gt. abs(aj(l,k)) ) l = i jpvt(k) = l s = aj(l,k) aj(l,k) = aj(k,k) aj(k,k) = s c check for zero pivot if (abs(s) .lt. 1.0e-15) goto 6 c calculate multipliers do 2 i = k+1,n 2 aj(i,k) = -aj(i,k)/s c interchange and eliminate by columns do 4 j = k+1,n s = aj(l,j) aj(l,j) = aj(k,j) aj(k,j) = s if (abs(s) .lt. 1.0e-15) goto 4 do 3 i = k+1,n 3 aj(i,j) = aj(i,j) + aj(i,k)*s 4 continue 5 continue return 6 jpvt(n) = -1 if (jpvt(n) .eq. -1) then write(16,20) it stop endif 20 format(1x,'zero pivot detected, it=',i4) return end c ------------ : ------------ subroutine solve include 'newbin' c solve solves the linear system aj*x=r c assumes a is factorised into l.u form (by fact) c returns solution in x integer i,k,l,km,ka real*8 s c forward elimination do 2 k = 1, n-1 l = jpvt(k) s = rd(l) rd(l) = rd(k) rd(k) = s do 1 i = k+1,n 1 rd(i) = rd(i) + aj(i,k)*s 2 continue c back substitution do 4 ka = 1,n-1 km = n - ka k = km + 1 rd(k) = rd(k)/aj(k,k) s = - rd(k) do 3 i = 1,km 3 rd(i) = rd(i) + aj(i,k)*s 4 continue rd(1) = rd(1)/aj(1,1) c 42 write(16,43)(rd(j),j=1,n) c 43 format(' rd=',8f10.6) return end c ------------ : ------------ subroutine newuv include 'newbin' do 22 j = 2,nxp do 22 k = 2,nyp mb = 2*(k-2) + 2*(ny-2)*(j-2) du = -rd(mb+1) dv = -rd(mb+2) u(k,j) = u(k,j) + du*om v(k,j) = v(k,j) + dv*om 22 continue c write(16,*) 'it:',it,' rms:',rms c do 42 k = 1,ny c 42 write(16,43)(u(k,j),j=1,nx) c do 44 k = 1,ny c 44 write(16,45)(v(k,j),j=1,nx) 43 format(' u=',7f10.4) c 45 format(' v=',7f10.4) return end c ------------ : ------------ subroutine fwrite include 'newbin' if (rms .le. eps) write(16,101) it,rms if (rms .gt. eps) write(16,103) itmx,rms if (rms .le. eps) write(6,101) it,rms if (rms .gt. eps) write(6,103) itmx,rms do 25 k = 1,ny 25 write(16,105)(u(k,j),j=1,nx) do 27 k = 1,ny 27 write(16,107)(v(k,j),j=1,nx) if (ird .ne. 0) then do 29 k = 1,ny write(17,109)(u(k,j),j=1,nx) write(17,109)(v(k,j),j=1,nx) 29 continue endif write(14) u,v close(11) close(14) close(16) close(17) 101 format(//,'Convergence after ',i3, + ' iterations. RMS residual is', d12.5) 103 format(//,'No convergence after ',i3, + ' iterations. RMS residual is', d12.5) 105 format(' u=',7f10.4) 107 format(' v=',7f10.4) 109 format(10f8.5) return end