c======================================================================= c NOTE: The code below has been "hacked" from code for the solution c of the Laplace equation in 1-, 2- or 3-D, in order to be used c for the solution of the Crank-Nicholson scheme for the diffusion c equation. However, only the 1- and 2-D cases have been been so c hacked. c======================================================================= c----------------------------------------------------------------------- c Residual Tracing routine c----------------------------------------------------------------------- subroutine trace_res(string,l,resnrm) implicit none character*(*) string integer l real*8 resnrm write(0,1000) string, l, resnrm 1000 format(a,' Level: ',i2,' |res.|: ',1pe12.3) return end c----------------------------------------------------------------------- c Does one-step 1d Point-GS relaxation and returns norm of c running ('instantaneous') residuals---checkerboard ordering. c----------------------------------------------------------------------- real*8 function relaxc1(u,rhs,nx,h) implicit none include 'mgcom.inc' real*8 dvnrm2 integer nx real*8 u(nx), rhs(nx) real*8 h real*8 hm2, jacel, rresl integer i, pass logical ltrace parameter ( ltrace = .false. ) relaxc1 = 0.0d0 jacel = 1.0d0 + cdt / h**2 if( ltrace ) then write(0,*) 'relaxc1: |u pre| = ', dvnrm2(u,nx) end if do pass = 1 , 2 do i = pass + 1 , nx - 1 , 2 rresl = u(i) - & 0.5d0 * cdt * (u(i+1) - 2.0d0 * u(i) + u(i-1)) * / (h**2) - & rhs(i) u(i) = u(i) - rresl / jacel relaxc1 = relaxc1 + rresl * rresl end do end do if( ltrace ) then write(0,*) 'relaxc1: |u pst| = ', dvnrm2(u,nx) write(0,*) 'relaxc1: |rhs| = ', dvnrm2(rhs,nx) write(0,*) 'relaxc1: |rres| = ', relaxc1 end if relaxc1 = sqrt(relaxc1 / nx) return end c----------------------------------------------------------------------- c Computes negative of 1d residuals ... c----------------------------------------------------------------------- subroutine cmres1(mres,u,rhs,nx,h) implicit none include 'mgcom.inc' integer nx real*8 mres(nx), u(nx), * rhs(nx) real*8 h integer i real*8 hm2 mres(1) = rhs(1) - u(1) do i = 2 , nx - 1 mres(i) = rhs(i) - & (u(i) - & 0.5d0 * cdt * (u(i+1) - 2.0d0 * u(i) + u(i-1)) & / (h**2)) end do mres(nx) = rhs(nx) - u(nx) return end c----------------------------------------------------------------------- c Does one-step 2d Point-GS relaxation and returns norm of c running ('instantaneous') residuals---checkerboard ordering. c----------------------------------------------------------------------- real*8 function relaxc2(u,rhs,nx,ny,h) implicit none include 'mgcom.inc' real*8 dmnrm2 integer nx, ny real*8 u(nx,ny), rhs(nx,ny) real*8 h real*8 hm2, mjelm1, rresl integer i, j, * isw, jsw, pass logical ltrace parameter ( ltrace = .false. ) relaxc2 = 0.0d0 hm2 = 1.0d0 / (h * h) mjelm1 = -1.0d0 / (1.0d0 + 2.0d0 * cdt / h**2) if( ltrace ) then write(0,*) 'relaxc2: |u pre| = ', dmnrm2(u,nx,ny) end if jsw = 1 do pass = 1 , 2 isw = jsw do j = 2 , ny - 1 do i = isw + 1 , nx - 1, 2 rresl = u(i,j) - 0.5d0 * cdt * & ( u(i+1,j) - 2.0d0 * u(i,j) + u(i-1,j) + & u(i,j+1) - 2.0d0 * u(i,j) + u(i,j-1) & ) / h**2 * - rhs(i,j) u(i,j) = u(i,j) + mjelm1 * rresl relaxc2 = relaxc2 + rresl * rresl end do isw = 3 - isw end do jsw = 3 - jsw end do relaxc2 = sqrt(relaxc2 / (nx * ny)) if( ltrace ) then write(0,*) 'relaxc2: |u pst| = ', dmnrm2(u,nx,ny) write(0,*) 'relaxc2: |rhs| = ', dmnrm2(rhs,nx,ny) write(0,*) 'relaxc2: |rres| = ', relaxc2 end if return end c----------------------------------------------------------------------- c Computes negative of 2d residuals ... c----------------------------------------------------------------------- subroutine cmres2(mres,u,rhs,nx,ny,h) implicit none include 'mgcom.inc' integer nx, ny real*8 mres(nx,ny), u(nx,ny), * rhs(nx,ny) real*8 h integer i, j real*8 hm2 hm2 = 1.0d0 / (h * h) do j = 2 , ny - 1 do i = 2 , nx - 1 mres(i,j) = rhs(i,j) - & (u(i,j) - 0.5d0 * cdt * & ( u(i+1,j) - 2.0d0 * u(i,j) + u(i-1,j) + & u(i,j+1) - 2.0d0 * u(i,j) + u(i,j-1) & ) / h**2 & ) end do end do call dmlbs(mres,0.0d0,nx,ny) return end c----------------------------------------------------------------------- c Does one-step 3d Point-GS relaxation and returns norm of c running ('instantaneous') residuals---checkerboard ordering. c----------------------------------------------------------------------- real*8 function relaxc3(u,rhs,nx,ny,nz,h) implicit none integer nx, ny, nz real*8 u(nx,ny,nz), rhs(nx,ny,nz) real*8 h real*8 hm2, mjelm1, rresl integer i, j, k, * isw, ksw, pass relaxc3 = 0.0d0 hm2 = 1.0d0 / (h * h) mjelm1 = h * h / 6.0d0 ksw = 2 do pass = 1 , 2 isw = ksw do k = 2 , nz - 1 do j = 2 , ny - 1 do i = isw + 1 , nx - 1, 2 rresl = hm2 * (u(i+1,j,k) + u(i-1,j,k) + * u(i,j+1,k) + u(i,j-1,k) + * u(i,j,k+1) + u(i,j,k-1) - * 6.0d0 * u(i,j,k)) - * rhs(i,j,k) u(i,j,k) = u(i,j,k) + mjelm1 * rresl relaxc3 = relaxc3 + rresl * rresl end do isw = 3 - isw end do end do ksw = 3 - ksw end do relaxc3 = sqrt(relaxc3 / (nx * ny * nz)) return end c----------------------------------------------------------------------- c Computes negative of 3d residuals ... c----------------------------------------------------------------------- subroutine cmres3(mres,u,rhs,nx,ny,nz,h) implicit none integer nx, ny, nz real*8 mres(nx,ny,nz), u(nx,ny,nz), * rhs(nx,ny,nz) real*8 h integer i, j, k real*8 hm2 hm2 = 1.0d0 / (h * h) do k = 2 , nz - 1 do j = 2 , ny - 1 do i = 2 , nx - 1 mres(i,j,k) = rhs(i,j,k) - * hm2 * (u(i+1,j,k) + u(i-1,j,k) + * u(i,j+1,k) + u(i,j-1,k) + * u(i,j,k+1) + u(i,j,k-1) - * 6.0d0*u(i,j,k)) end do end do end do call d3lbs(mres,0.0d0,nx,ny,nz) return end c----------------------------------------------------------------------- c Computes and dumps true residuals c----------------------------------------------------------------------- subroutine cres1_trace(label,unit,res,u,rhs,nx,h) implicit none character*(*) label integer nx, unit real*8 res(*), u(*), rhs(*) real*8 h integer shape(1) real*8 vh(1) shape(1) = nx vh(1) = h call cres(res,u,rhs,1,shape,vh) write(unit,1000) label 1000 format('cres1_trace: ',a) call dvdump(u,nx,'u',unit) call dvdump(res,nx,'residual',unit) return end c----------------------------------------------------------------------- c Computes true (non-negated) residual. c----------------------------------------------------------------------- subroutine cres(res,u,rhs,rank,shape,vh) implicit none integer rank real*8 res(*), u(*), rhs(*), * vh(rank) integer shape(rank) call cmres(res,u,rhs,rank,shape,vh) if( rank .eq. 1 ) then call dvsm(res,-1.0d0,res,shape(1)) else if( rank .eq. 2 ) then call dvsm(res,-1.0d0,res,shape(1)*shape(2)) else if( rank .eq. 3 ) then call dvsm(res,-1.0d0,res,shape(1)*shape(2)*shape(3)) else write(0,*) 'cmres: rank = ', rank,' not implemented' end if return end c----------------------------------------------------------------------- c Front end to cmres1, cmres2, cmres3 ... c----------------------------------------------------------------------- subroutine cmres(mres,u,rhs,rank,shape,vh) implicit none integer rank real*8 mres(*), u(*), rhs(*), * vh(rank) integer shape(rank) if( rank .eq. 1 ) then call cmres1(mres,u,rhs,shape(1), vh(1)) else if( rank .eq. 2 ) then call cmres2(mres,u,rhs,shape(1),shape(2), vh(1)) else if( rank .eq. 3 ) then call cmres3(mres,u,rhs,shape(1),shape(2),shape(3),vh(1)) else write(0,*) 'cmres: rank = ', rank,' not implemented' end if return end