channel.f90

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

  end do

  ipvt(n) = n

  if ( abd(m,n) == 0.0D+00 ) then
    info = n
  end if

  return
end
subroutine dgbsl ( abd, lda, n, ml, mu, ipvt, b, job )

!*****************************************************************************80
!
!! DGBSL solves a banded system factored by DGBFA.
!
!  Discussion:
!
!    SGBSL can solve either a * x = b  or  trans(a) * x = b.
!
!  Parameters:
!
!     on entry
!
!        abd     real ( kind = 8 )(lda, n)
!                the output from dgbco or dgbfa.
!
!        lda     integer
!                the leading dimension of the array  abd .
!
!        n       integer
!                the order of the original matrix.
!
!        ml      integer
!                number of diagonals below the main diagonal.
!
!        mu      integer
!                number of diagonals above the main diagonal.
!
!        ipvt    integer(n)
!                the pivot vector from dgbco or dgbfa.
!
!        b       real ( kind = 8 )(n)
!                the right hand side vector.
!
!        job     integer
!                = 0         to solve  a*x = b ,
!                = nonzero   to solve  trans(a)*x = b , where
!                            trans(a)  is the transpose.
!
!     on return
!
!        b       the solution vector  x .
!
!     error condition
!
!        a division by zero will occur if the input factor contains a
!        zero on the diagonal.  technically this indicates singularity
!        but it is often caused by improper arguments or improper
!        setting of lda .  it will not occur if the subroutines are
!        called correctly and if dgbco has set 0.0 < RCOND
!        or dgbfa has set info == 0 .
!
!     to compute  inverse(a) * c  where  c  is a matrix
!     with  p  columns
!           call dgbco ( abd, lda, n, ml, mu, ipvt, rcond, z )
!           if (rcond is too small) go to ...
!           do j = 1, p
!              call dgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 )
!           end do
!
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
  integer lda
  integer n

  real ( kind = 8 ) abd(lda,n)
  real ( kind = 8 ) b(n)
  integer ipvt(n)
  integer job
  integer k
  integer l
  integer la
  integer lb
  integer lm
  integer m
  integer ml
  integer mu
  real ( kind = 8 ) ddot
  real ( kind = 8 ) t

  m = mu + ml + 1
!
!  JOB = 0, Solve  a * x = b.
!
!  First solve l*y = b.
!
  if ( job == 0 ) then

     if ( 0 < ml ) then

        do k = 1, n-1
           lm = min ( ml, n-k )
           l = ipvt(k)
           t = b(l)
           if ( l /= k ) then
              b(l) = b(k)
              b(k) = t
           end if
           call daxpy ( lm, t, abd(m+1,k), 1, b(k+1), 1 )
        end do

     end if
!
!  Now solve u*x = y.
!
     do k = n, 1, -1
        b(k) = b(k) / abd(m,k)
        lm = min ( k, m ) - 1
        la = m - lm
        lb = k - lm
        t = -b(k)
        call daxpy ( lm, t, abd(la,k), 1, b(lb), 1 )
     end do
!
!  JOB nonzero, solve  trans(a) * x = b.
!
!  First solve  trans(u)*y = b.
!
  else

     do k = 1, n
        lm = min ( k, m ) - 1
        la = m - lm
        lb = k - lm
        t = ddot ( lm, abd(la,k), 1, b(lb), 1 )
        b(k) = ( b(k) - t ) / abd(m,k)
     end do
!
!  Now solve trans(l)*x = y
!
     if ( 0 < ml ) then

        do k = n-1, 1, -1
           lm = min ( ml, n-k )
           b(k) = b(k) + ddot ( lm, abd(m+1,k), 1, b(k+1), 1 )
           l = ipvt(k)
           if ( l /= k ) then
              t = b(l)
              b(l) = b(k)
              b(k) = t
           end if
        end do

      end if

  end if

  return
end
subroutine dscal ( n, da, dx, incx )

!*****************************************************************************80
!
!! DSCAL scales a vector by a constant.
!
!     uses unrolled loops for increment equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 3/93 to return if incx  <=  0.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
  real ( kind = 8 ) da,dx(*)
  integer i,incx,m,n,nincx
!
  if (  n <= 0 .or. incx <= 0 )return
  if ( incx == 1)go to 20
!
!        code for increment not equal to 1
!
  nincx = n*incx
  do i = 1,nincx,incx
    dx(i) = da*dx(i)
  end do
  return
!
!        code for increment equal to 1
!
!
!        clean-up loop
!
   20 continue

  m = mod(n,5)
  if (  m  ==  0 ) go to 40
  dx(1:m) = da*dx(1:m)
  if (  n  <  5 ) return

40 continue

  do i = m+1,n,5
    dx(i) = da*dx(i)
    dx(i + 1) = da*dx(i + 1)
    dx(i + 2) = da*dx(i + 2)
    dx(i + 3) = da*dx(i + 3)
    dx(i + 4) = da*dx(i + 4)
  end do

  return
end
subroutine gdump (f,indx,insc,iounit,isotri,long,nelemn,neqn, &
  nnodes,node,np,npara,nx,ny,para,reynld,rjpnew,xc,yc)

!*****************************************************************************80
!
!! GDUMP writes information to a file.
!
!  Discussion:
!
!    The information can be used to create
!    graphics images.  In order to keep things simple, exactly one
!    value, real or integer, is written per record.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, INTEGER NPARA, the number of parameters.  Fixed at 1
!    for now.
!
!    Input, real ( kind = 8 ) PARA(MAXPAR), the parameters.
!
  integer nelemn
  integer neqn
  integer nnodes
  integer np

  real ( kind = 8 ) f(neqn)
  real ( kind = 8 ) fval
  integer i
  integer indx(np,2)
  integer insc(np)
  integer iounit
  integer, save :: iset = 0
  integer isotri(nelemn)
  integer j
  logical long
  integer node(nelemn,nnodes)
  integer npara
  integer nx
  integer ny
  real ( kind = 8 ) para
  real ( kind = 8 ) reynld
  real ( kind = 8 ) rjpnew
  real ( kind = 8 ) ubdry
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) yc(np)

  iset = iset+1

  write(iounit,*)long
  write (iounit,*) nelemn
  write (iounit,*) np
  write (iounit,*) npara
  write (iounit,*) nx
  write (iounit,*) ny
!
!  Pressures
!
  do i = 1, np
    j = insc(i)
    if (j <= 0) then
      fval = 0.0D+00
    else
      fval = f(j)
    end if
    write (iounit,*) fval
  end do
!
!  Horizontal velocities, U
!
  do i = 1, np
    j = indx(i,1)
    if (j == 0) then
      fval = 0.0D+00
    else if (j < 0) then
      fval = ubdry(yc(i),para)
    else
      fval = f(j)
    end if
    write (iounit,*) fval
  end do
!
!  Vertical velocities, V
!
  do i = 1, np
    j = indx(i,2)
    if (j <= 0) then
      fval = 0.0D+00
    else
      fval = f(j)
    end if
    write (iounit,*) fval
  end do

  do i = 1, np
    write (iounit,*) indx(i,1)
    write (iounit,*) indx(i,2)
  end do

  do i = 1, np
    write (iounit,*) insc(i)
  end do

  do i = 1, nelemn
    write (iounit,*) isotri(i)
  end do

  do i = 1, nelemn
    do j = 1, 6
      write (iounit,*) node(i,j)
    end do
  end do

  write (iounit,*) para
  write (iounit,*) reynld
  write (iounit,*) rjpnew

  do i = 1, np
    write (iounit,*) xc(i)
  end do

  do i = 1, np
    write (iounit,*) yc(i)
  end do

  write (*,*) 'GDUMP wrote data set ',iset,' to file.'

  return
end
subroutine getg ( f, iline, my, neqn, u )

!*****************************************************************************80
!
!! GETG outputs field values along the profile line X = XZERO.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer neqn
  integer my

  real ( kind = 8 ) f(neqn)
  integer iline(my)
  integer j
  integer k
  real ( kind = 8 ) u(my)

  do j = 1, my
    k = iline(j)
    if ( 0 < k ) then
      u(j) = f(k)
    else
      u(j) = 0.0D+00
    end if
  end do

  return
end
subroutine gram ( gr, iline, indx, iwrite, my, nelemn, nnodes, node, &
  nodex0, np, para, r, ui, xc, yc )

!*****************************************************************************80
!
!! GRAM computes the Gram matrix, GR(I,J) = INTEGRAL PHI(I)*PHI(J).
!
!  and the vector R(I) = INTEGRAL UI*PHI(I).
!
!  The integrals are computed along the line where the profile is
!  specified.  The three point Gauss quadrature rule is used for the
!  line integral.
!
!  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

  real ( kind = 8 ) ar
  real ( kind = 8 ) bb
  real ( kind = 8 ) bbb
  real ( kind = 8 ) bbx
  real ( kind = 8 ) bby
  real ( kind = 8 ) bma2
  real ( kind = 8 ) bx
  real ( kind = 8 ) by
  real ( kind = 8 ) gr(my,my)
  integer i
  integer igetl
  integer ii
  integer iline(my)
  integer indx(np,2)
  integer ip
  integer ipp
  integer iq
  integer iqq
  integer iquad
  integer it
  integer iun
  integer iwrite
  integer j
  integer jj
  integer k
  integer kk
  integer node(nelemn,nnodes)
  integer nodex0
  real ( kind = 8 ) para
  real ( kind = 8 ) r(my)
  real ( kind = 8 ) ubc
  real ( kind = 8 ) ubdry
  real ( kind = 8 ) ui(my)
  real ( kind = 8 ) uiqdpt
  real ( kind = 8 ) wt(3)
  real ( kind = 8 ) x
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) xzero
  real ( kind = 8 ) y
  real ( kind = 8 ) yq(3)
  real ( kind = 8 ) yc(np)
!
!  Input values for 3 point Gauss quadrature
!
  wt(1) = 5.0D+00 / 9.0D+00
  wt(2) = 8.0D+00 / 9.0D+00
  wt(3) = wt(1)
  yq(1) = -0.7745966692D+00
  yq(2) = 0.0D+00
  yq(3) = -yq(1)
!
!  zero arrays
!
  r(1:my) = 0.0D+00
  gr(1:my,1:my) = 0.0D+00
!
!  Compute line integral by looping over intervals along line
!  using three point Gauss quadrature
!
  xzero = xc(nodex0)
  do 70 it = 1, nelemn
!
!  Check to see if we are in a triangle with a side along line
!  x = xzero.  If not, skip out
!
    k = node(it,1)
    kk = node(it,2)

    if ( 1.0D-04 < abs(xc(k)-xzero) ) then
      cycle
    end if

    if ( 1.0D-04 < abs(xc(kk)-xzero) ) then
      cycle
    end if

    do 60 iquad = 1, 3
      bma2 = (yc(kk)-yc(k))/2.0
      ar = bma2*wt(iquad)
      x = xzero
      y = yc(k)+bma2*(yq(iquad)+1.0D+00 )
!
!  Compute u internal at quadrature points
!
      uiqdpt = 0
      do 30 iq = 1, nnodes
        if ( 4 < iq ) go to 30
        if (iq == 3) go to 30
        call qbf (x,y,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc)
        ip = node(it,iq)
        iun = indx(ip,1)
        if ( 0 < iun ) then
          ii = igetl(iun,iline,my)
          uiqdpt = uiqdpt+bb*ui(ii)
        else if (iun < 0) then
          ubc = ubdry(yc(ip),para)
          uiqdpt = uiqdpt+bb*ubc
        end if
   30     continue
!
!  Only loop over nodes lying on line x = xzero
!
      do 50 iq = 1, nnodes
        if ( iq == 1.or.iq == 2.or.iq == 4 ) then
          ip = node(it,iq)
          call qbf (x,y,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc)
          i = indx(ip,1)
          if (i <= 0) go to 50
          ii = igetl(i,iline,my)
          r(ii) = r(ii)+bb*uiqdpt*ar

          do iqq = 1, nnodes
            if ( iqq == 1.or.iqq == 2.or.iqq == 4 ) then
              ipp = node(it,iqq)
              call qbf (x,y,it,iqq,bbb,bbx,bby,nelemn,nnodes, &
                node,np,xc,yc)
              j = indx(ipp,1)
              if (j /= 0) then
                jj = igetl(j,iline,my)
                gr(ii,jj) = gr(ii,jj)+bb*bbb*ar
              end if
            end if
          end do

        end if
   50     continue
   60   continue
   70 continue

  if ( 2 <= iwrite ) then
    write(*,*)' '
    write(*,*)'Gram matrix:'
    write(*,*)' '
    do i = 1,my
      do j = 1,my
        write(*,*)i,j,gr(i,j)
      end do
    end do
    write(*,*)' '
    write(*,*)'R vector:'
    write(*,*)' '
    do i = 1,my