flow3.f90

FLOW采用有限单元法fortran90编写的求解不可压缩流体的稳态流速和压力场的程序

!
!  Parameters:
!
  implicit none

  integer neqn
  integer np
  integer nparf

  real ( kind = 8 ) dudyn(np)
  real ( kind = 8 ) dvdyn(np)
  character ( len = 2 ) eqn(neqn)
  real ( kind = 8 ) f(neqn)
  real ( kind = 8 ) g(neqn)
  integer i
  integer ibc
  integer ihor
  integer indx(np,3)
  integer ip
  integer ipar
  integer iparb
  integer ishapb
  integer iver
  integer nparb
  real ( kind = 8 ) shape
  real ( kind = 8 ) splbmp(4,nparb+2,0:nparb)
  real ( kind = 8 ) taubmp(nparb+2)
  real ( kind = 8 ) ubc
  real ( kind = 8 ) vbc
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) yc(np)

  iparb = ipar-nparf

  f(1:neqn) = 0.0D+00

  do ip = 1,np

    ihor = indx(ip,1)
    if ( eqn(ihor) == 'UB' ) then
      if ( ibc == 0 ) then
        call bump_bc(dudyn,dvdyn,ip,iparb,ishapb,np,nparb,shape,splbmp,taubmp, &
          ubc,vbc,xc)
      else if ( ibc == 1 ) then
        call bump_bc1(g,indx,ip,iparb,ishapb,neqn,np,nparb,shape,splbmp,taubmp, &
          ubc,vbc,xc,yc)
      else if ( ibc == 2 ) then
        call bump_bc2(g,indx,ip,iparb,ishapb,neqn,np,nparb,shape,splbmp,taubmp, &
          ubc,vbc,xc,yc)
      else if ( ibc == 3 ) then
        call bump_bc(dudyn,dvdyn,ip,iparb,ishapb,np,nparb,shape,splbmp,taubmp, &
          ubc,vbc,xc)
      end if
      f(ihor) = ubc
    end if

    iver = indx(ip,2)
    if ( eqn(iver) == 'VB' ) then
      if ( ibc == 0 ) then
        call bump_bc(dudyn,dvdyn,ip,iparb,ishapb,np,nparb,shape,splbmp,taubmp, &
          ubc,vbc,xc)
      else if ( ibc == 1 ) then
        call bump_bc1(g,indx,ip,iparb,ishapb,neqn,np,nparb, &
          shape,splbmp,taubmp,ubc,vbc,xc,yc)
      else if ( ibc == 2 ) then
        call bump_bc2(g,indx,ip,iparb,ishapb,neqn,np,nparb, &
          shape,splbmp,taubmp,ubc,vbc,xc,yc)
      else if ( ibc == 3 ) then
        call bump_bc(dudyn,dvdyn,ip,iparb,ishapb,np,nparb, &
          shape,splbmp,taubmp,ubc,vbc,xc)
      end if
      f(iver) = vbc
    end if

  end do

  return
end
subroutine bump_spl ( ishapb, npar, nparb, nparf, par, splbmp, taubmp, xbl, &
  xbr, ybl, ybr )

!*****************************************************************************80
!
!! BUMP_SPL sets up or updates the spline data that describes the bump.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    22 December 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ISHAPB.
!    1, the bump is modeled by C0 linear splines.
!    2, the bump is modeled by C0 quadratic splines.
!    3, the bump is modeled by C1 cubic splines.
!
  implicit none

  integer npar
  integer nparb
  integer nparf

  integer i
  integer ibcbeg
  integer ibcend
  integer ipar
  integer ishapb
  real ( kind = 8 ) par(npar)
  real ( kind = 8 ) splbmp(4,nparb+2,0:nparb)
  real ( kind = 8 ) taubmp(nparb+2)
  real ( kind = 8 ) xbl
  real ( kind = 8 ) xbr
  real ( kind = 8 ) ybl
  real ( kind = 8 ) ybr

  if ( nparb <= 0 ) then
    return
  end if
!
!  Set up the bump arrays, including:
!
!    TAUBMP, containing the abscissas, which never change,
!    SPLBMP(1,I,0), the location of the bump at abscissa I,
!    SPLBMP(1,I,IPAR), the I-th coefficient of the partial derivative
!    of the bump shape with respect to the IPAR-th bump parameter.
!
  do i = 1, nparb+2
    taubmp(i) = ( real ( nparb + 2 - i ) * xbl + real ( i - 1 ) * xbr ) &
      / real ( nparb + 1 )
  end do

  do i = 1, nparb+2

    if ( i == 1 ) then
      splbmp(1,i,0) = ybl
    else if ( 2  <= i .and. i <= nparb+1 ) then
      splbmp(1,i,0) = par(nparf+i-1)
    else if ( i == nparb+2 ) then
      splbmp(1,i,0) = ybr
    end if

  end do

  if ( ishapb == 3 ) then

    do i = 1, nparb+2

      do ipar = 1, nparb

        if ( ipar+1 /=  i ) then
          splbmp(1,i,ipar) = 0.0D+00
        else
          splbmp(1,i,ipar) = 1.0D+00
        end if

      end do

    end do

    ibcbeg = 0
    ibcend = 0
    do i = 0, nparb
      call cubspl ( taubmp, splbmp(1,1,i), nparb+2, ibcbeg, ibcend )
    end do

  end if

  return
end
subroutine ch_cap ( c )

!*****************************************************************************80
!
!! CH_CAP capitalizes a single character.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    19 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none

  character c
  integer itemp

  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
subroutine chkdat(ibump,idfd,ids,ifds,igrad,ijac,iopt,ipred,itype,jjac, &
  maxpar,maxparb,maxparf,nopt,npar,nparb,nparf,nstep3,para1,partar,xbleft, &
  xbrite)

!*****************************************************************************80
!
!! CHKDAT performs some simple checks on the input data.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IGRAD, the cost gradient approximation option.
!    0, No cost gradient approximation is made.
!    1, Chain rule on discretized sensitivities.
!    2, Chain rule on finite coefficient differences.
!    3, Chain rule on corrected finite coefficient differences.
!    4, Direct finite cost differences.
!
  implicit none

  integer maxpar

  integer i
  integer ibump
  integer idfd
  integer ids
  integer ifds
  integer igrad
  integer ijac
  integer iopt(maxpar)
  integer ipred
  integer isum
  integer itype
  integer jjac
  integer maxparb
  integer maxparf
  integer nopt
  integer npar
  integer nparb
  integer nparf
  integer nstep3
  real ( kind = 8 ) para1(maxpar)
  real ( kind = 8 ) partar(maxpar)
  real ( kind = 8 ) psum
  real ( kind = 8 ) xbleft
  real ( kind = 8 ) xbrite
!
!  IBUMP
!
  if ( ibump < 0 .or. ibump > 3 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IBUMP is out of bounds!'
    write ( *, * ) '  Current value of IBUMP = ',ibump
    write ( *, * ) '  Legal values must be between 0 and 3.'
    stop
  end if
!
!  IGRAD
!
  if ( igrad < 0 .or. igrad > 4 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IGRAD must be between 0 and 4,'
    write ( *, * ) '  but your value is ',igrad
    stop
  end if

  if ( igrad == 0 ) then
    if ( itype == 3 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  A cost gradient approximation MUST be made'
      write ( *, * ) '  when ITYPE = ',itype
      write ( *, * ) '  but you have chosen IGRAD = ',igrad
      stop
    end if
  else if ( igrad == 1 ) then
    if ( ids == 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  Your cost gradient choice IGRAD = ',igrad
      write ( *, * ) '  requires a nonzero value of IDS,'
      write ( *, * ) '  but your value is IDS = ',ids
      stop
    end if
  else if ( igrad == 2 ) then
    if ( ifds == 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  Your cost gradient choice IGRAD = ',igrad
      write ( *, * ) '  requires a nonzero value of IFDS,'
      write ( *, * ) '  but your value is IFDS = ',ifds
      stop
    end if
  else if ( igrad == 3 ) then
    if ( ifds == 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  Your cost gradient choice IGRAD = ',igrad
      write ( *, * ) '  requires a nonzero value of IFDS,'
      write ( *, * ) '  but your value is IFDS = ',ifds
      stop
    end if
  else if ( igrad == 4 ) then
    if ( idfd == 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  Your cost gradient choice IGRAD = ',igrad
      write ( *, * ) '  requires a nonzero value of IDFD,'
      write ( *, * ) '  but your value is IDFD = ',idfd
      stop
    end if
  end if
!
!  IJAC
!
  if ( ijac<1 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IJAC is out of bounds!'
    write ( *, * ) '  Current value of IJAC = ',ijac
    write ( *, * ) '  Legal values must be 1 or greater.'
    stop
  end if
!
!  IOPT
!
  do i = 1,npar
    if ( iopt(i)/= 0 .and. iopt(i)/=1 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  IOPT(*) must be 0 or 1, but'
      write ( *, * ) '  IOPT(',I,') = ',iopt(i)
      stop
    end if
  end do
!
!  IPRED
!
  if ( ipred<0 .or. ipred>4 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IPRED must be between 0 and 4,'
    write ( *, * ) '  but your value is ',ipred
    stop
  end if

  if ( ipred == 1 .and. ids==0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IPRED = ',ipred
    write ( *, * ) '  but IDS = ',ids
    stop
  else if ( ipred == 2 .and. ifds==0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IPRED = ',ipred
    write ( *, * ) '  but IFDS = ',ifds
    stop
  else if ( ipred == 3 .and. ifds==0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IPRED = ',ipred
    write ( *, * ) '  but IFDS = ',ifds
    stop
  else if ( ipred == 4 .and. ids==0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  IPRED = ',ipred
    write ( *, * ) '  but IDS = ',ids
    stop
  end if
!
!  ITYPE
!
  if ( itype<1 .or. itype>7 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  Illegal value of ITYPE = ',itype
    write ( *, * ) '  Legal values are between 1 and 7.'
    write ( *, * ) '  ChkDat forces a STOP!'
    stop
  end if
!
!  JJAC
!
  if ( jjac<0 .or. jjac>3 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  Illegal value of JJAC = ',jjac
    write ( *, * ) '  Legal values are between 0 and 3.'
    write ( *, * ) '  ChkDat forces a STOP!'
    stop
  end if
!
!  NOPT
!
  if ( nopt <= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  The number of free parameters, NOPT = ',nopt
    write ( *, * ) '  but this value must be positive!'
    stop
  end if
!
!  NPAR
!
  if ( npar>maxpar ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  NPAR is out of bounds!'
    write ( *, * ) '  Current value of NPAR = ',npar
    write ( *, * ) '  Maximum legal value, MAXPARA = ',maxpar
    stop
  end if
!
!  NPARB
!
  if ( nparb<0 .or. nparb>maxparb ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  NPARB is out of bounds!'
    write ( *, * ) '  Input value of NPARB = ',nparb
    write ( *, * ) '  Maximum legal value, MAXPARB = ',maxparb
    stop
  end if
!
!  NPARF
!
  if ( nparf <= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  NPARF is out of bounds!'
    write ( *, * ) '  The input value of NPARF is ',nparf
    write ( *, * ) '  But NPARF must be at least 1.'
    stop
  end if

  if ( nparf>maxparf ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Fatal error!'
    write ( *, * ) '  NPARF is too big!'
    write ( *, * ) '  Input value of NPARF = ',nparf
    write ( *, * ) '  Maximum legal value, MAXPARF = ',maxparf
    stop
  end if
!
!  NSTEP3
!
  if ( itype == 4 ) then
    if ( nstep3<1 ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) ' Nonpositive value for NSTEP3 = ',nstep3
      stop
    else if ( nstep3 == 1 ) then
      write ( *, * ) 'CHKDAT - Warning!'
      write ( *, * ) '  NSTEP3 = 1 is an unusual value!'
    end if
  end if
!
!  PARA1(1:NPARF)
!
  if ( itype == 3 ) then

    psum = sum ( abs ( para1(1:nparf) ) )

    if ( psum == 0.0D+00 ) then

      isum = 0
      do i = 1,nparf
        isum = isum+abs(iopt(i))
      end do

      if ( isum == 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'CHKDAT - Fatal error!'
        write ( *, * ) '  You are optimizing, ITYPE = ',itype
        write ( *, * ) '  But all inflows are zero,'
        write ( *, * ) '  and no inflow may vary.'
        stop
      end if

    end if
  end if
!
!  PARA1(NPARF+NPARB+1), the value of NU_INV.