channel.f90

主要用来计算水槽中中粘性不可压缩流动的程序

  integer j1
  integer j2
  integer j3
  integer node(nelemn,nnodes)
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) x
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) y
  real ( kind = 8 ) yc(np)

  if (in <= 3) then
    in1 = in
    in2 = mod(in,3)+1
    in3 = mod(in+1,3)+1
    i1 = node(it,in1)
    i2 = node(it,in2)
    i3 = node(it,in3)
    d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1))
    t = 1.0+((yc(i2)-yc(i3))*(x-xc(i1))+(xc(i3)-xc(i2))*(y-yc(i1)))/d
    bb = t*(2.0D+00*t-1.0D+00)
    bx = (yc(i2)-yc(i3))*(4.0D+00*t-1.0D+00)/d
    by = (xc(i3)-xc(i2))*(4.0D+00*t-1.0D+00)/d
  else
    inn = in-3
    in1 = inn
    in2 = mod(inn,3)+1
    in3 = mod(inn+1,3)+1
    i1 = node(it,in1)
    i2 = node(it,in2)
    i3 = node(it,in3)
    j1 = i2
    j2 = i3
    j3 = i1
    d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1))
    c = (xc(j2)-xc(j1))*(yc(j3)-yc(j1))-(xc(j3)-xc(j1))*(yc(j2)-yc(j1))
    t = 1.0D+00+((yc(i2)-yc(i3))*(x-xc(i1))+(xc(i3)-xc(i2))*(y-yc(i1)))/d
    s = 1.0D+00+((yc(j2)-yc(j3))*(x-xc(j1))+(xc(j3)-xc(j2))*(y-yc(j1)))/c
    bb = 4.0D+00*s*t
    bx = 4.0D+00*(t*(yc(j2)-yc(j3))/c+s*(yc(i2)-yc(i3))/d)
    by = 4.0D+00*(t*(xc(j3)-xc(j2))/c+s*(xc(i3)-xc(i2))/d)
  end if

  return
end
subroutine resid (area,g,indx,insc,iwrite,nelemn,neqn,nnodes, &
  node,np,nquad,para,phi,psi,res,reynld,yc)

!*****************************************************************************80
!
!! RESID computes the residual.
!
!  Discussion:
!
!    The G array contains the current iterate.
!
!    The RES array will contain the value of the residual.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer nelemn
  integer neqn
  integer nnodes
  integer np
  integer nquad

  real ( kind = 8 ) aijpu
  real ( kind = 8 ) aijpv
  real ( kind = 8 ) aijup
  real ( kind = 8 ) aijuu
  real ( kind = 8 ) aijuv
  real ( kind = 8 ) aijvp
  real ( kind = 8 ) aijvu
  real ( kind = 8 ) aijvv
  real ( kind = 8 ) ar
  real ( kind = 8 ) area(nelemn)
  real ( kind = 8 ) bb
  real ( kind = 8 ) bbb
  real ( kind = 8 ) bbbl
  real ( kind = 8 ) bbbx
  real ( kind = 8 ) bbby
  real ( kind = 8 ) bbl
  real ( kind = 8 ) bbx
  real ( kind = 8 ) bby
  real ( kind = 8 ) bx
  real ( kind = 8 ) by
  real ( kind = 8 ) g(neqn)
  integer i
  integer ibad
  integer ihor
  integer imax
  integer indx(np,2)
  integer insc(np)
  integer ip
  integer ipp
  integer iprs
  integer iq
  integer iqq
  integer iquad
  integer it
  integer iver
  integer iwrite
  integer j
  integer jp
  integer ju
  integer jv
  integer node(nelemn,nnodes)
  real ( kind = 8 ) para
  real ( kind = 8 ) phi(nelemn,nquad,nnodes,3)
  real ( kind = 8 ) psi(nelemn,nquad,nnodes)
  real ( kind = 8 ) res(neqn)
  real ( kind = 8 ) reynld
  real ( kind = 8 ) rmax
  real ( kind = 8 ) test
  real ( kind = 8 ) ubdry
  real ( kind = 8 ) un(2)
  real ( kind = 8 ) unx(2)
  real ( kind = 8 ) uny(2)
  real ( kind = 8 ) uu
  real ( kind = 8 ) visc
  real ( kind = 8 ) yc(np)

  visc = 1.0D+00 / reynld

  res(1:neqn) = 0.0D+00
!
!  For each element,
!
  do 90 it = 1, nelemn

    ar = area(it) / 3.0D+00
!
!  and for each quadrature point in the element,
!
    do 80 iquad = 1, nquad
!
!  Evaluate velocities at quadrature point
!
      call uval (g,indx,iquad,it,nelemn,neqn,nnodes,node, &
        np,nquad,para,phi,un,unx,uny,yc)
!
!  For each basis function,
!
      do 70 iq = 1, nnodes
        ip = node(it,iq)
        bb = phi(it,iquad,iq,1)
        bx = phi(it,iquad,iq,2)
        by = phi(it,iquad,iq,3)
        bbl = psi(it,iquad,iq)
        iprs = insc(ip)
        ihor = indx(ip,1)
        iver = indx(ip,2)

        if ( 0 < ihor ) then
          res(ihor) = res(ihor)+(un(1)*unx(1)+un(2)*uny(1))*bb*ar
        end if

        if ( 0 < iver ) then
          res(iver) = res(iver)+(un(1)*unx(2)+un(2)*uny(2))*bb*ar
        end if
!
!  For another basis function,
!
        do 50 iqq = 1, nnodes
          ipp = node(it,iqq)
          bbb = phi(it,iquad,iqq,1)
          bbx = phi(it,iquad,iqq,2)
          bby = phi(it,iquad,iqq,3)
          bbbl = psi(it,iquad,iqq)
          ju = indx(ipp,1)
          jv = indx(ipp,2)
          jp = insc(ipp)

          if ( 0 < ju ) then
            if ( 0 < ihor ) then
              aijuu = visc*(by*bby+bx*bbx) &
                + bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))
              res(ihor) = res(ihor)+aijuu*ar*g(ju)
            end if
            if ( 0 < iver ) then
              aijvu = bb*bbb*unx(2)
              res(iver) = res(iver)+aijvu*ar*g(ju)
            end if
            if ( 0 < iprs ) then
              aijpu = bbx*bbl
              res(iprs) = res(iprs)+aijpu*ar*g(ju)
            end if
          else if ( ju < 0 ) then
            uu = ubdry(yc(ipp),para)
            if ( 0 < ihor ) then
              aijuu = visc*(by*bby+bx*bbx) &
                + bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))
              res(ihor) = res(ihor)+ar*aijuu*uu
            end if
            if ( 0 < iver ) then
              aijvu = bb*bbb*unx(2)
              res(iver) = res(iver)+ar*aijvu*uu
            end if
            if ( 0 < iprs ) then
              aijpu = bbx*bbl
              res(iprs) = res(iprs)+ar*aijpu*uu
            end if
          end if

          if ( 0 < jv ) then
            if ( 0 < ihor ) then
              aijuv = bb*bbb*uny(1)
              res(ihor) = res(ihor)+aijuv*ar*g(jv)
            end if
            if ( 0 < iver ) then
              aijvv = visc*(by*bby+bx*bbx) &
                +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))
              res(iver) = res(iver)+aijvv*ar*g(jv)
            end if
            if ( 0  < iprs ) then
              aijpv = bby*bbl
              res(iprs) = res(iprs)+aijpv*ar*g(jv)
            end if
          end if

          if ( 0 < jp ) then
            if ( 0 < ihor ) then
              aijup = -bx*bbbl
              res(ihor) = res(ihor)+aijup*ar*g(jp)
            end if
            if ( 0 < iver ) then
              aijvp = -by*bbbl
              res(iver) = res(iver)+aijvp*ar*g(jp)
            end if
          end if

   50       continue
   70     continue
   80   continue
   90 continue

  res(neqn) = g(neqn)

  rmax = 0.0D+00
  imax = 0
  ibad = 0

  do i = 1,neqn

    test = abs(res(i))

    if ( rmax < test ) then
      rmax = test
      imax = i
    end if

    if ( 1.0D-03 < test ) then
      ibad = ibad+1
    end if

  end do

  if ( 1 <= iwrite ) then
    write(*,*)' '
    write(*,*)'RESIDUAL INFORMATION:'
    write(*,*)' '
    write(*,*)'Worst residual is number ',IMAX
    write(*,*)'of magnitude ',RMAX
    write(*,*)' '
    write(*,*)'Number of "bad" residuals is ',IBAD,' out of ',NEQN 
    write(*,*)' '
  end if

  if ( 2 <= iwrite ) then
    write(*,*)'Raw residuals:'
    write(*,*)' '
    i = 0

    do j = 1,np

      if ( 0 < indx(j,1) ) then
        i = i+1
        if ( abs(res(i)) <= 1.0D-03 ) then
          write(*,'(1x,a1,2i5,g14.6)')'U',i,j,res(i)
        else
          write(*,'(a1,a1,2i5,g14.6)')'*','U',i,j,res(i)
        end if
      end if

      if ( 0 < indx(j,2) ) then
        i = i+1
        if ( abs(res(i)) <= 1.0D-03 ) then
          write(*,'(1x,a1,2i5,g14.6)')'V',i,j,res(i)
        else
          write(*,'(a1,a1,2i5,g14.6)')'*','V',i,j,res(i)
        end if
      end if

      if ( 0 < insc(j) ) then
        i = i+1
        if ( abs(res(i)) <= 1.0D-03 ) then
          write(*,'(1x,a1,2i5,g14.6)')'P',i,j,res(i)
        else
          write(*,'(a1,a1,2i5,g14.6)')'*','P',i,j,res(i)
        end if
      end if

    end do

  end if

  return
end
subroutine setban(indx,insc,maxrow,nband,nelemn,nlband,nnodes, &
  node,np,nrow)

!*****************************************************************************80
!
!! SETBAN computes the half band width.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer nelemn
  integer nnodes
  integer np

  integer i
  integer indx(np,2)
  integer insc(np)
  integer ip
  integer ipp
  integer iq
  integer iqq
  integer it
  integer iuk
  integer iukk
  integer j
  integer maxrow
  integer nband
  integer nlband
  integer node(nelemn,nnodes)
  integer nrow

  nlband = 0

  do it = 1, nelemn

    do iq = 1, nnodes

      ip = node(it,iq)

      do iuk = 1, 3

        if (iuk == 3) then
          i = insc(ip)
        else
          i = indx(ip,iuk)
        end if

        if ( 0 < i ) then
          do iqq = 1, nnodes
            ipp = node(it,iqq)
            do iukk = 1, 3
              if (iukk == 3) then
                j = insc(ipp)
              else
                j = indx(ipp,iukk)
              end if
              nlband = max(nlband,j-i)
            end do
          end do
        end if

      end do

    end do
  end do

  nband = nlband+nlband+1
  nrow = nlband+nlband+nlband+1

  write (*,*) 'Lower bandwidth = ',nlband
  write (*,*) 'Total bandwidth = ',nband
  write (*,*) 'NROW  = ',nrow
  if ( maxrow < nrow ) then
    write(*,*)'SETBAN - NROW is too large!'
    write(*,*)'The maximum allowed is ',maxrow
    stop
  end if

  return
end
subroutine setbas(nelemn,nnodes,node,np,nquad,phi,psi,xc,xm,yc,ym)

!*****************************************************************************80
!
!! SETBAS computes the basis functions at each integration point.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer nelemn
  integer nnodes
  integer np
  integer nquad

  real ( kind = 8 ) bb
  real ( kind = 8 ) bsp
  real ( kind = 8 ) bx
  real ( kind = 8 ) by
  integer iq
  integer it
  integer j
  integer node(nelemn,nnodes)
  real ( kind = 8 ) phi(nelemn,nquad,nnodes,3)
  real ( kind = 8 ) psi(nelemn,nquad,nnodes)
  real ( kind = 8 ) x
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) xm(nelemn,nquad)
  real ( kind = 8 ) y
  real ( kind = 8 ) yc(np)
  real ( kind = 8 ) ym(nelemn,nquad)

  do it = 1,nelemn
    do j = 1,nquad
      x = xm(it,j)
      y = ym(it,j)
      do iq = 1,6
        psi(it,j,iq) = bsp(x,y,it,iq,1,nelemn,nnodes,node,np,xc,yc)
        call qbf(x,y,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc)
        phi(it,j,iq,1) = bb
        phi(it,j,iq,2) = bx
        phi(it,j,iq,3) = by
      end do
    end do
  end do
  return
end
subroutine setgrd (indx,insc,isotri,iwrite,long,maxeqn,mx,my, &
  nelemn,neqn,nnodes,node,np,nx,ny)

!*****************************************************************************80
!
!! SETGRD sets up the grid for the problem..
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer my
  integer nelemn
  integer nnodes
  integer np

  integer i
  integer ic
  integer icnt
  integer ielemn
  integer indx(np,2)
  integer insc(np)
  integer ip
  integer ip1
  integer ip2
  integer isotri(nelemn)
  integer it
  integer iwrite
  integer jc
  integer jcnt
  logical long
  integer maxeqn
  integer mx
  integer neqn
  integer node(nelemn,nnodes)
  integer nx
  integer ny
!
!  Determine whether region is long or skinny.  This will determine
!  how we number the nodes and elements.
!
  if ( ny < nx ) then
    long = .true.
    write(*,*)'Using vertical ordering.'
  else
    long = .false.
    write(*,*)'Using horizontal ordering.'
  end if
!
!  Set parameters for Taylor Hood element
!
  write (*,*) ' '
  write (*,*) 'SETGRD: Taylor Hood element'
!
!  Construct grid coordinates, elements, and ordering of unknowns
!
  neqn = 0
  ielemn = 0

  do ip = 1,np

    if ( long ) then
      ic = ((ip-1)/my)+1
      jc = mod((ip-1),my)+1
    else
      ic = mod((ip-1),mx)+1
      jc = ((ip-1)/mx)+1
    end if

    icnt = mod(ic,2)
    jcnt = mod(jc,2)
!
!  If both the row count and the column count are odd,
!  and we're not in the last row or top column, 
!  then we can define two new triangular elements based at the node.
!
!  For horizontal ordering, 
!  given the following arrangement of nodes, for instance:
! 
!    21 22 23 24 25
!    16 17 18 19 20
!    11 12 13 14 15
!    06 07 08 09 10
!    01 02 03 04 05
!
!  when we arrive at node 13, we will define