channel.f90

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

      write(*,*)i,r(i)
    end do
  end if

  return
end
function idamax ( n, dx, incx )

!*****************************************************************************80
!
!! IDAMAX finds the index of element having max. absolute value.
!
!     jack dongarra, linpack, 3/11/78.
!     modified 3/93 to return if incx  <=  0.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
  integer idamax
  real ( kind = 8 ) dmax
  real ( kind = 8 ) dx(*)
  integer i,incx,ix,n

  idamax = 0
  if (  n < 1 .or. incx <= 0 ) return
  idamax = 1
  if ( n == 1)return
  if ( incx == 1)go to 20
!
!        code for increment not equal to 1
!
  ix = 1
  dmax = dabs(dx(1))
  ix = ix + incx
  do i = 2,n
     if ( dabs(dx(ix)) <= dmax) go to 5
     idamax = i
     dmax = dabs(dx(ix))
    5    ix = ix + incx
  end do

  return
!
!        code for increment equal to 1
!
   20 continue

  dmax = abs ( dx(1) )

  do i = 2,n
    if ( dmax < abs ( dx(i) ) ) then
      idamax = i
      dmax = abs ( dx(i) )
    end if
  end do

  return
end
function igetl ( k, iline, my )

!*****************************************************************************80
!
!! IGETL gets the local unknown number along the profile line.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer my

  integer igetl
  integer iline(my)
  integer j
  integer k

  do j = 1, my
    if (iline(j) == k) then
      igetl = j
      return
    end if
  end do

  write ( *, * ) ' '
  write (*,*) 'IGETL - fatal error!'
  write (*,*) '  Unable to get local unknown number for '
  write (*,*) '  Global variable number ',k
  igetl = 0
  stop
end
subroutine linsys ( a, area, f, g, indx, insc, ipivot, maxrow, nelemn, &
  neqn, nlband, nnodes, node, np, nquad, nrow, para1, para2, phi, psi, &
  reynld, yc )

!*****************************************************************************80
!
!! LINSYS sets up and solves the linear system.
!
!  Discussion:
!
!    The G array contains the previous solution.
!
!    The F array contains the right hand side initially and then the
!    current solution.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 February 2006
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer maxrow
  integer nelemn
  integer neqn
  integer nnodes
  integer np
  integer nquad
  integer nrow

  real ( kind = 8 ) a(maxrow,neqn)
  real ( kind = 8 ) ar
  real ( kind = 8 ) area(nelemn)
  real ( kind = 8 ) bb
  real ( kind = 8 ) bbb
  real ( kind = 8 ) bbbl
  real ( kind = 8 ) bbl
  real ( kind = 8 ) bbx
  real ( kind = 8 ) bby
  real ( kind = 8 ) bx
  real ( kind = 8 ) by
  real ( kind = 8 ) f(neqn)
  real ( kind = 8 ) g(neqn)
  integer i
  integer idamax
  integer ihor
  integer indx(np,2)
  integer info
  integer insc(np)
  integer ioff
  integer ip
  integer ipivot(neqn)
  integer ipp
  integer iprs
  integer iq
  integer iqq
  integer iquad
  integer it
  integer iuse
  integer iver
  integer j
  integer job
  integer jp
  integer ju
  integer jv
  integer nlband
  integer node(nelemn,nnodes)
  real ( kind = 8 ) para1
  real ( kind = 8 ) para2
  real ( kind = 8 ) phi(nelemn,nquad,nnodes,3)
  real ( kind = 8 ) psi(nelemn,nquad,nnodes)
  real ( kind = 8 ) reynld
  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)

  ioff = nlband + nlband + 1
  visc = 1.0D+00 / reynld
  f(1:neqn) = 0.0D+00
  a(1:nrow,1:neqn) = 0.0D+00
!
!  For each element,
!
  do it = 1, nelemn

    ar = area(it) / 3.0D+00
!
!  and for each quadrature point in the element,
!
    do iquad = 1, nquad
!
!  Evaluate velocities at quadrature point
!
      call uval (g,indx,iquad,it,nelemn,neqn,nnodes,node, &
        np,nquad,para1,phi,un,unx,uny,yc)
!
!  For each basis function,
!
      do 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)
        ihor = indx(ip,1)
        iver = indx(ip,2)
        iprs = insc(ip)

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

        if ( 0 < iver ) then
          f(iver) = f(iver)+ar*bb*(un(1)*unx(2)+un(2)*uny(2))
        end if
!
!  For another basis function,
!
        do 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)
!
!  Horizontal velocity variable
!
          if ( 0 < ju ) then

            if ( 0 < ihor ) then
              iuse = ihor-ju+ioff
              a(iuse,ju) = a(iuse,ju)+ar*(visc*(by*bby+bx*bbx) &
                + bb*(bbb*unx(1)+bbx*un(1)+bby*un(2)))
            end if

            if ( 0 < iver ) then
              iuse = iver-ju+ioff
              a(iuse,ju) = a(iuse,ju)+ar*bb*bbb*unx(2)
            end if

            if ( 0 < iprs ) then
              iuse = iprs-ju+ioff
              a(iuse,ju) = a(iuse,ju)+ar*bbx*bbl
            end if

          else if ( ju < 0 ) then

            uu = ubdry(yc(ipp),para2)
            if ( 0 < ihor ) then 
              f(ihor) = f(ihor)-ar*uu*(visc*(by*bby+bx*bbx) &
                + bb*(bbb*unx(1)+bbx*un(1)+bby*un(2)))
            end if

            if ( 0 < iver ) then
              f(iver) = f(iver)-ar*uu*bb*bbb*unx(2)
            end if

            if ( 0 < iprs ) then
              f(iprs) = f(iprs)-ar*uu*bbx*bbl
            end if

          end if
!
!  Vertical velocity variable
!
          if ( 0 < jv ) then

            if ( 0 < ihor ) then
              iuse = ihor-jv+ioff
              a(iuse,jv) = a(iuse,jv)+ar*bb*bbb*uny(1)
            end if

            if ( 0 < iver ) then
              iuse = iver-jv+ioff
              a(iuse,jv) = a(iuse,jv)+ar*(visc*(by*bby+bx*bbx) &
                +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1)))
            end if

            if ( 0 < iprs ) then
              iuse = iprs-jv+ioff
              a(iuse,jv) = a(iuse,jv)+ar*bby*bbl
            end if

          end if
!
!  Pressure variable
!
          if ( 0 < jp ) then

            if ( 0 < ihor ) then
              iuse = ihor-jp+ioff
              a(iuse,jp) = a(iuse,jp)-ar*bx*bbbl
            end if

            if ( 0 < iver ) then
              iuse = iver-jp+ioff
              a(iuse,jp) = a(iuse,jp)-ar*by*bbbl
            end if

          end if

        end do
      end do
    end do
  end do
!
!  To avoid singularity of the pressure system, the last pressure
!  is simply assigned a value of 0.
!
  f(neqn) = 0.0D+00
  do j = neqn-nlband, neqn-1
    i = neqn-j+ioff
    a(i,j) = 0.0D+00
  end do
  a(ioff,neqn) = 1.0D+00
!
!  Factor the matrix
!
  call dgbfa ( a, maxrow, neqn, nlband, nlband, ipivot, info )

  if ( info /= 0 ) then
    write (*,*) ' '
    write (*,*) 'LINSYS - fatal error!'
    write (*,*) 'DGBFA returns INFO = ',info
    stop
  end if
!
!  Solve the linear system
!
  job = 0
  call dgbsl ( a, maxrow, neqn, nlband, nlband, ipivot, f, job ) 

  return
end
subroutine nstoke (a,area,f,g,indx,insc,ipivot,iwrite,maxnew,maxrow, &
  nelemn,neqn,nlband,nnodes,node,np,nquad,nrow,numnew,para,phi,psi, &
  reynld,tolnew,yc)

!*****************************************************************************80
!
!! NSTOKE solves the Navier Stokes equation using Taylor-Hood elements.
!
!  The G array contains the previous iterate.
!
!  The F array contains the right hand side initially and then the
!  current iterate.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer maxrow
  integer nelemn
  integer neqn
  integer nnodes
  integer np
  integer nquad
  integer nrow

  real ( kind = 8 ) a(maxrow,neqn)
  real ( kind = 8 ) area(nelemn)
  double precision diff
  real ( kind = 8 ) f(neqn)
  real ( kind = 8 ) g(neqn)
  integer idamax
  integer indx(np,2)
  integer insc(np)
  integer ipivot(neqn)
  integer iter
  integer iwrite
  integer maxnew
  integer nlband
  integer node(nelemn,nnodes)
  integer numnew
  real ( kind = 8 ) para
  real ( kind = 8 ) phi(nelemn,nquad,nnodes,3)
  real ( kind = 8 ) psi(nelemn,nquad,nnodes)
  real ( kind = 8 ) reynld
  real ( kind = 8 ) tolnew
  real ( kind = 8 ) yc(np)

  do iter = 1, maxnew

    numnew = numnew + 1

    call linsys (a,area,f,g,indx,insc,ipivot, &
      maxrow,nelemn,neqn,nlband,nnodes,node, &
      np,nquad,nrow,para,para,phi,psi,reynld,yc)
!
!  Check for convergence
!
    g(1:neqn) = g(1:neqn) - f(1:neqn)
    diff = abs ( g(idamax(neqn,g,1)) )

    if ( 1 <= iwrite ) then
      write(*,*)'NSTOKE iteration ',iter,' Mnorm = ',diff
    end if

    g(1:neqn) = f(1:neqn)

    if ( diff <= tolnew ) then
      write (*,*) 'Navier Stokes iteration converged in ', &
        iter,' iterations.'
      return
    end if

  end do

  write (*,*) 'Navier Stokes solution did not converge!'

  return
end
subroutine pval (g,insc,long,mx,my,nelemn,neqn,nnodes,node,np,press)

!*****************************************************************************80
!
!! PVAL computes a table of pressures.
!
!  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

  real ( kind = 8 ) g(neqn)
  integer i
  integer insc(np)
  integer ip
  integer iq
  integer it
  integer ivar
  integer j
  logical long
  integer mx
  integer my
  integer node(nelemn,nnodes)
  real ( kind = 8 ) press(mx,my)

  press(1:mx,1:my) = 0.0D+00
!
!  Read the pressures where they are computed.  
!  These are "(odd, odd)" points.
!
  do it = 1, nelemn
    do iq = 1, 3

      ip = node(it,iq)
      ivar = insc(ip)

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

      if ( 0 < ivar ) then
        press(i,j) = g(ivar)
      else
        press(i,j) = 0.0D+00
      end if

    end do
  end do
!
!  Interpolate the pressures at points (even, odd) and (odd, even).
!
  do i = 2,mx-1,2
    do j = 1,my,2
      press(i,j) = 0.5D+00*(press(i-1,j)+press(i+1,j))
    end do
  end do

  do j = 2,my-1,2
    do i = 1,mx,2
      press(i,j) = 0.5D+00*(press(i,j-1)+press(i,j+1))
    end do
  end do
!
!  Interpolate the pressures at points (even,even).
!
  do j = 2,my-1,2
    do i = 2,mx-1,2
      press(i,j) = 0.5D+00*(press(i-1,j-1)+press(i+1,j+1))
    end do
  end do

  return
end
subroutine qbf (x,y,it,in,bb,bx,by,nelemn,nnodes,node,np,xc,yc)

!*****************************************************************************80
!
!! QBF evaluates a quadratic basis function in a triangle.
!
!  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

  real ( kind = 8 ) bb
  real ( kind = 8 ) bx
  real ( kind = 8 ) by
  real ( kind = 8 ) c
  real ( kind = 8 ) d
  integer i1
  integer i2
  integer i3
  integer in
  integer in1
  integer in2
  integer in3
  integer inn
  integer it