matrix3.f

二维热能方程的求解

!=======================================================================
      subroutine matrix3(imode,nx,ny,level,ax,bx,ay,by,a,ierr)
!=======================================================================
      implicit none
      integer imode,nx,ny,level,ierr
      real*8 ax,bx,ay,by
      real*8 a(3,*)

! for imode=1, a(3,*)=a(3,nx,ny)=(b1 Ux)x-0.5cU
! for imode=2, a(3,*)=a(3,ny,nx)=(b2 Uy)y-0.5cU
    
      include 'commons.h'
      integer ix,iy,ntmp,node,info,iout
      real*8 hx,hy,px,py,twohx,twohy,rhx2,rhy2,d1,d2,diag
      real*8 half,ftnc
      character*2 matrix

      if(imode.eq.1) then
         matrix='A1'
      elseif(imode.eq.2) then
         matrix='A2'
      else
         print'("Error from MATRIX3: imode=",i2)',imode
         ierr=1
         return
      endif

      if(level.ge.1) then
         print'("MATRIX3: Matrix=",a2,"; imode=",i1)',matrix,imode
      endif

      hx=(bx-ax)/dble(nx-1)
      hy=(by-ay)/dble(ny-1)
      twohx=2.d0*hx
      twohy=2.d0*hy
      rhx2=1.d0/(hx*hx)
      rhy2=1.d0/(hy*hy)
      half=0.5d0

!--- MAKE A -----------------
 
      if(imode.eq.1) then
         do iy=1,ny
            py=hy*dble(iy-1)+ay
            ntmp=(iy-1)*nx
         do ix=1,nx
            node=ntmp+ix
            px=hx*dble(ix-1)+ax
            call setpar3(nx,ny,ix,iy,imode,px,py,hx,hy,twohx,twohy,
     &                   d1,d2,diag)
            a(1,node)= rhx2*d1
            a(2,node)=-rhx2*diag-half*ftnc(px,py)
            a(3,node)= rhx2*d2
         end do
         end do
      else
         do ix=1,nx
            px=hx*dble(ix-1)+ax
            ntmp=(ix-1)*ny
         do iy=1,ny
            py=hy*dble(iy-1)+ay
            node=ntmp+iy
            call setpar3(nx,ny,ix,iy,imode,px,py,hx,hy,twohx,twohy,
     &                   d1,d2,diag)
            a(1,node)= rhy2*d1
            a(2,node)=-rhy2*diag-half*ftnc(px,py)
            a(3,node)= rhy2*d2
         end do
         end do
      endif

      if(level.ge.4) then
         info=10+2*imode-1
         iout=10+2*imode
         if(imode.eq.1) then
            do iy=1,ny
               ntmp=(iy-1)*nx
            do ix=1,nx
               node=ntmp+ix
               write(iout,'(3(f10.5))') (a(ntmp,node),ntmp=1,3)
            enddo
            enddo
         else
            do ix=1,nx
               ntmp=(ix-1)*ny
            do iy=1,ny
               node=ntmp+iy
               write(iout,'(3(f10.5))') (a(ntmp,node),ntmp=1,3)
            enddo
            enddo
         endif
         write(info,*) "nx=",nx
         write(info,*) "ny=",ny
         write(info,*) "imode=",imode
         write(info,*) "Matrix=",matrix
         print*,matrix,": info,   saved in unit=",info
         print*,matrix,": matrix, saved in unit=",iout
      endif

      return
      end

c ======================================================================
      subroutine setpar3(nx,ny,ix,iy,imode,px,py,hx,hy,twohx,twohy,
     &                   d1,d2,diag)
c ======================================================================
      implicit none
      integer nx,ny,ix,iy,imode
      real*8  px,py,hx,hy,twohx,twohy,d1,d2,diag

      real*8  a0,a1,a2,bdry,b1,b2,beta

      bdry =0.d0

      if(imode.eq.1) then

         a0=b1(px,py)
         a1=b1(px-hx,py)
         a2=b1(px+hx,py)
         d1=0.5d0*(a0+a1)
         d2=0.5d0*(a0+a2)
         if (ix.eq.1) then
            d2=d1+d2
            d1=0.0d0
            bdry=beta(px,py,1)*twohx
         else if (ix.eq.nx) then
            d1=d1+d2
            d2=0.0d0
            bdry=beta(px,py,2)*twohx
         end if

      else

         a0=b2(px,py)
         a1=b2(px,py-hy)
         a2=b2(px,py+hy)
         d1=0.5d0*(a0+a1)
         d2=0.5d0*(a0+a2)
         if (iy.eq.1) then
            d2=d1+d2
            d1=0.0d0
            bdry=beta(px,py,3)*twohy
         else if (iy.eq.ny) then
            d1=d1+d2
            d2=0.0d0
            bdry=beta(px,py,4)*twohy
         end if

      endif

      diag=d1+d2+bdry

      return
      end