channel.f90

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

    end if
!
!  Solve linear system for du/da
!
    para = anew
    abound = 1.0D+00
    call linsys (a,area,g,f,indx,insc,ipivot, &
      maxrow,nelemn,neqn,nlband,nnodes,node, &
      np,nquad,nrow,para,abound,phi,psi,reynld,yc)
!
!  Output in DCDA
!
    call getg ( g, iline, my, neqn, dcda )

    if ( 2 <= iwrite ) then
      write (*,*) ' '
      write (*,*) 'Sensitivities:'
      write (*,*) ' '
      write (*,'(5g14.6)') dcda(1:my)
    end if
!
!  Evaluate J prime at current value of parameter where J is 
!  functional to be minimized.
!
!  JPRIME = 2.0 * DCDA(I) * (GR(I,J)*UNEW(J)-R(I))
!
    rjpnew = 0.0D+00
    do i = 1, my
      temp = -r(i)
      do j = 1, my
        temp = temp + gr(i,j) * unew(j)
      end do
      rjpnew = rjpnew + 2.0D+00 * dcda(i) * temp
    end do

    write (*,*) ' '
    write (*,*) 'Parameter  = ',anew,' J prime = ',rjpnew
!
!  Dump information for graphics
!
    if ( .false. ) then
      para = anew
      call gdump (f,indx,insc,iounit,isotri,long,nelemn,neqn, &
        nnodes,node,np,npara,nx,ny,para,reynld,rjpnew,xc,yc)
    end if
!
!  Update the estimate of the parameter using the secant step
!
    if (iter == 1) then
      a2 = 0.5D+00
    else
      a2 = aold-rjpold*(anew-aold)/(rjpnew-rjpold)
    end if

    aold = anew
    anew = a2
    rjpold = rjpnew
    test = abs(anew-aold)/abs(anew)

    write (*,*) 'New value of parameter = ',anew
    write(*,*)'Convergence test = ',test

    if (abs(anew-aold) < abs(anew)*tolsec) then
      write (*,*) 'Secant iteration converged.'
      go to 40
    end if

  end do

  write (*,*) 'Secant iteration failed to converge.'
   40 continue
!
!  Produce total CPU time used.
!
!     tarray(2) = 0.0D+00
!     tarray(1) = second()
!
  call etime (tarray)
  tarray(1) = tarray(1)+tarray(2)
  tarray(2) = tarray(1)/60.0D+00

  write (*,*) ' '
  write (*,*) 'Total execution time = ', tarray(1), ' seconds.'
  write (*,*) '                     = ', tarray(2), ' minutes.'
  write (*,*)'Number of secant steps = ', NUMSEC
  write (*,*)'Number of Newton steps = ', NUMNEW
!
!  Close graphics file
!
  if ( .false. ) then
    close ( unit = iounit )
  end if
  stop
!
!  Error opening graphics file
!
   50 continue
  write (*,*) 'CHANNEL could not open the graphics file!'
  stop
end
function bsp ( xq, yq, it, iq, id, nelemn, nnodes, node, np, xc, yc )

!*****************************************************************************80
!
!! BSP evaluates the linear basis functions associated with pressure.
!
!  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 ) bsp
  real ( kind = 8 ) d
  integer i1
  integer i2
  integer i3
  integer id
  integer iq
  integer iq1
  integer iq2
  integer iq3
  integer it
  integer node(nelemn,nnodes)
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) xq
  real ( kind = 8 ) yc(np)
  real ( kind = 8 ) yq

  iq1 = iq
  iq2 = mod(iq,3)+1
  iq3 = mod(iq+1,3)+1
  i1 = node(it,iq1)
  i2 = node(it,iq2)
  i3 = node(it,iq3)
  d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1))

  if (id == 1) then
 
    bsp = 1.0+((yc(i2)-yc(i3))*(xq-xc(i1))+(xc(i3)-xc(i2))*(yq-yc(i1)))/d
 
  else if (id == 2) then
 
    bsp = (yc(i2)-yc(i3))/d
 
  else if (id == 3) then
 
    bsp = (xc(i3)-xc(i2))/d
 
  else
 
    write (*,*) 'BSP - fatal error!'
    write (*,*) 'unknown value of id = ',id
    stop

  end if

  return
end
subroutine daxpy ( n, da, dx, incx, dy, incy )

!*****************************************************************************80
!
!! DAXPY computes constant times a vector plus a vector.
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    Jack Dongarra
!
  real ( kind = 8 ) dx(*),dy(*),da
  integer i,incx,incy,ix,iy,m,n
!
  if ( n <= 0)return
  if (da  ==  0.0d+00 ) return
  if ( incx == 1.and.incy == 1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
  ix = 1
  iy = 1
  if ( incx < 0)ix = (-n+1)*incx + 1
  if ( incy < 0)iy = (-n+1)*incy + 1
  do i = 1,n
    dy(iy) = dy(iy) + da*dx(ix)
    ix = ix + incx
    iy = iy + incy
  end do
  return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,4)
  if (  m  ==  0 ) go to 40
  do 30 i = 1,m
    dy(i) = dy(i) + da*dx(i)
   30 continue
  if (  n  <  4 ) return
   40 continue

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

  return
end
subroutine dcopy ( n, dx, incx, dy, incy )

!*****************************************************************************80
!
!! DCOPY copies a vector, x, to a vector, y.
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    Jack Dongarra
!
  real ( kind = 8 ) dx(*),dy(*)
  integer i,incx,incy,ix,iy,m,n

  if ( n <= 0)return
  if ( incx == 1.and.incy == 1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
  ix = 1
  iy = 1
  if ( incx < 0)ix = (-n+1)*incx + 1
  if ( incy < 0)iy = (-n+1)*incy + 1
  do i = 1,n
    dy(iy) = dx(ix)
    ix = ix + incx
    iy = iy + incy
  end do
  return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,7)
  if (  m  ==  0 ) go to 40

  dy(1:m) = dx(1:m)

  if (  n  <  7 ) return
   40 continue
  do i = m+1, n ,7
    dy(i) = dx(i)
    dy(i + 1) = dx(i + 1)
    dy(i + 2) = dx(i + 2)
    dy(i + 3) = dx(i + 3)
    dy(i + 4) = dx(i + 4)
    dy(i + 5) = dx(i + 5)
    dy(i + 6) = dx(i + 6)
  end do

  return
end
function ddot ( n, dx, incx, dy, incy )

!*****************************************************************************80
!
!! DDOT forms the dot product of two vectors.
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    Jack Dongarra
!
  real ( kind = 8 ) ddot
  real ( kind = 8 ) dx(*),dy(*),dtemp
  integer i,incx,incy,ix,iy,m,n

  ddot = 0.0d+00
  dtemp = 0.0d+00
  if ( n <= 0)return
  if ( incx == 1.and.incy == 1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
  ix = 1
  iy = 1
  if ( incx < 0)ix = (-n+1)*incx + 1
  if ( incy < 0)iy = (-n+1)*incy + 1
  do i = 1,n
    dtemp = dtemp + dx(ix)*dy(iy)
    ix = ix + incx
    iy = iy + incy
  end do

  ddot = dtemp
  return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,5)
  if (  m  ==  0 ) go to 40
  do i = 1,m
    dtemp = dtemp + dx(i)*dy(i)
  end do

  if (  n  <  5 ) go to 60
   40 continue
  do i = m+1, n, 5
    dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + &
        dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
  end do
   60 ddot = dtemp
  return
end
subroutine delete (filnam)

!*****************************************************************************80
!
!! DELETE deletes a file.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
  character filnam*(*)
  integer ios

  open (unit = 99,file=filnam,status='old', iostat = ios )

  if ( ios /= 0 ) then
    return
  end if

  close (unit = 99,status='delete', iostat = ios)

  return
end
subroutine dgbfa ( abd, lda, n, ml, mu, ipvt, info )

!*****************************************************************************80
!
!! DGBFA factors a band matrix by elimination.
!
!  Discussion:
!
!     dgbfa is usually called by dgbco, but it can be called
!     directly with a saving in time if  rcond  is not needed.
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    Cleve Moler
!
!  Parameters:
!
!     on entry
!
!        abd     real ( kind = 8 )(lda, n)
!                contains the matrix in band storage.  the columns
!                of the matrix are stored in the columns of  abd  and
!                the diagonals of the matrix are stored in rows
!                ml+1 through 2*ml+mu+1 of  abd .
!                see the comments below for details.
!
!        lda     integer
!                the leading dimension of the array  abd .
!                2*ml + mu + 1 <= LDA.
!
!        n       integer
!                the order of the original matrix.
!
!        ml      integer
!                number of diagonals below the main diagonal.
!                0 <= ml < n .
!
!        mu      integer
!                number of diagonals above the main diagonal.
!                0 <= mu < n .
!                more efficient if  ml <= mu .
!     on return
!
!        abd     an upper triangular matrix in band storage and
!                the multipliers which were used to obtain it.
!                the factorization can be written  a = l*u  where
!                l  is a product of permutation and unit lower
!                triangular matrices and  u  is upper triangular.
!
!        ipvt    integer(n)
!                an integer vector of pivot indices.
!
!        info    integer
!                = 0  normal value.
!                = k  if  u(k,k) == 0.0D+00 .  this is not an error
!                     condition for this subroutine, but it does
!                     indicate that dgbsl will divide by zero if
!                     called.  use  rcond  in dgbco for a reliable
!                     indication of singularity.
!
!     band storage
!
!           if  a  is a band matrix, the following program segment
!           will set up the input.
!
!                   ml = (band width below the diagonal)
!                   mu = (band width above the diagonal)
!                   m = ml + mu + 1
!                   do j = 1, n
!                      i1 = max ( 1, j-mu )
!                      i2 = min ( n, j+ml )
!                      do i = i1, i2
!                         k = i - j + m
!                         abd(k,j) = a(i,j)
!                      end do
!                   end do
!
!           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
!           in addition, the first  ml  rows in  abd  are used for
!           elements generated during the triangularization.
!           the total number of rows needed in  abd  is  2*ml+mu+1 .
!           the  ml+mu by ml+mu  upper left triangle and the
!           ml by ml  lower right triangle are not referenced.
!
!     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)
  integer i
  integer i0
  integer info
  integer ipvt(n)
  integer idamax
  integer j
  integer j0
  integer j1
  integer ju
  integer jz
  integer k
  integer l
  integer lm
  integer m
  integer ml
  integer mm
  integer mu
  real ( kind = 8 ) t

  m = ml + mu + 1
  info = 0
!
!  Zero initial fill-in columns.
!
  j0 = mu + 2
  j1 = min ( n, m ) - 1

  do jz = j0, j1
     i0 = m + 1 - jz
     do i = i0, ml
        abd(i,jz) = 0.0D+00
     end do
  end do

  jz = j1
  ju = 0
!
!  Gaussian elimination with partial pivoting.
!
  do k = 1, n-1
!
!  Zero next fill-in column.
!
     jz = jz + 1
     if ( jz <= n ) then
       abd(1:ml,jz) = 0.0D+00
     end if
!
!  Find L = pivot index.
!
     lm = min ( ml, n-k )
     l = idamax ( lm+1, abd(m,k), 1 ) + m - 1
     ipvt(k) = l + k - m
!
!  Zero pivot implies this column already triangularized.
!
     if ( abd(l,k) == 0.0D+00 ) then

       info = k
!
!  Interchange if necessary.
!
     else

        if ( l /= m ) then
           t = abd(l,k)
           abd(l,k) = abd(m,k)
           abd(m,k) = t
        end if
!
!  Compute multipliers.
!
        t = -1.0D+00 / abd(m,k)
        call dscal ( lm, t, abd(m+1,k), 1 )
!
!  Row elimination with column indexing.
!
        ju = min ( max ( ju, mu+ipvt(k) ), n )
        mm = m

        do j = k+1, ju
           l = l - 1
           mm = mm - 1
           t = abd(l,j)
           if ( l /= mm ) then
              abd(l,j) = abd(mm,j)
              abd(mm,j) = t
           end if
           call daxpy ( lm, t, abd(m+1,k), 1, abd(mm+1,j), 1 )
        end do

     end if