flow3.f90

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

!
  if ( para1(nparf+nparb+1)<0.0D+00 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Warning!'
    write ( *, * ) '  NU_INV entry of PARA1 is negative.'
    write ( *, * ) '  This is being changed to one.'
    para1(nparf+nparb+1) = 1.0D+00
  end if
!
!  PARTAR(NPARF+NPARB+1), the value of NU_INV.
!
  if ( partar(nparf+nparb+1)<0.0D+00 ) then
    write ( *, * ) ' '
    write ( *, * ) 'CHKDAT - Warning!'
    write ( *, * ) '  NU_INV entry of PARTAR is negative.'
    write ( *, * ) '  This is being changed to one.'
    partar(nparf+nparb+1) = 1.0D+00
  end if
!
!  XBLEFT, XBRITE
!
  if ( nparb>0 ) then
    if ( xbleft >= xbrite ) then
      write ( *, * ) ' '
      write ( *, * ) 'CHKDAT - Fatal error!'
      write ( *, * ) '  XBLEFT >= XBRITE.'
      write ( *, * ) '  XBLEFT = ',xbleft
      write ( *, * ) '  XBRITE = ',xbrite
      stop
    end if
  end if

  return
end
subroutine chkopt(cost,costar,dparfd,dparfdc,dparsn,g,gtar,ids,ifds,neqn, &
  npar,para,partar)

!*****************************************************************************80
!
!! CHKOPT is called at the end of an optimization to check the results.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
  implicit none

  integer neqn
  integer npar

  real ( kind = 8 ) cdist
  real ( kind = 8 ) cost
  real ( kind = 8 ) costar
  real ( kind = 8 ) cpfd
  real ( kind = 8 ) cpfdc
  real ( kind = 8 ) cpsn
  real ( kind = 8 ) dcfd
  real ( kind = 8 ) dcfdc
  real ( kind = 8 ) dcsn
  real ( kind = 8 ) dparfd(npar)
  real ( kind = 8 ) dparfdc(npar)
  real ( kind = 8 ) dparsn(npar)
  real ( kind = 8 ) g(neqn)
  real ( kind = 8 ) gdist
  real ( kind = 8 ) gtar(neqn)
  integer i
  integer ids
  integer ifds
  real ( kind = 8 ) para(npar)
  real ( kind = 8 ) partar(npar)
  real ( kind = 8 ) pdist

  cdist = cost-costar

  gdist = 0.0D+00
  do i = 1,neqn
    gdist = gdist+(gtar(i)-g(i))**2
  end do
  gdist = sqrt(gdist)

  pdist = 0.0D+00
  do i = 1,npar
    pdist = pdist+(partar(i)-para(i))**2
  end do
  pdist = sqrt(pdist)

  write ( *, * ) ' '
  write ( *, * ) 'ChkOpt:'
  write ( *, * ) '  L2 Distance from target solution =   ',gdist
  write ( *, * ) '  L2 Distance from target parameters = ',pdist
  write ( *, * ) '  Distance from target cost =          ',cdist
  write ( *, * ) ' '
  write ( *, * ) '  Estimated cost change if we moved to target:'
  write ( *, * ) ' '

  if ( ids/= 0 ) then

    dcsn = 0.0D+00
    do i = 1,npar
      dcsn = dcsn+dparsn(i)*(partar(i)-para(i))
    end do
    write ( *, * ) '  Sensitivities:               ',dcsn
    cpsn = 0.0D+00
    do i = 1,npar
      cpsn = cpsn+dparsn(i)**2
    end do
    cpsn = sqrt(cpsn)
    write ( *, * ) '  L2 norm of disc. sens. cost gradients:  ',cpsn
  end if

  if ( ifds/= 0 ) then

    dcfd = 0.0D+00
    do i = 1,npar
      dcfd = dcfd+dparfd(i)*(partar(i)-para(i))
    end do
    write ( *, * ) '  Finite differences:          ',dcfd

    cpfd = 0.0D+00
    do i = 1,npar
      cpfd = cpfd+dparfd(i)**2
    end do
    cpfd = sqrt(cpfd)
    write ( *, * ) '  L2 norm of fd cost gradients:           ',cpfd

  end if

  dcfdc = 0.0D+00
  do i = 1,npar
    dcfdc = dcfdc+dparfdc(i)*(partar(i)-para(i))
  end do
  write ( *, * ) '  Corrected finite differences:',dcfdc

  cpfdc = 0.0D+00
  do i = 1,npar
    cpfdc = cpfdc+dparfdc(i)**2
  end do
  cpfdc = sqrt(cpfdc)
  write ( *, * ) '  L2 norm of corrected fd cost gradients: ',cpfdc

  return
end
subroutine chrctd(string,dval,ierror,lchar)

!*****************************************************************************80
!
!! CHRCTD accepts a string of characters, and tries to extract a
!  real real number from the initial part of the
!  string.
!
!  CHRCTD will read as many characters as possible until it reaches
!  the end of the string, or encounters a character which cannot be
!  part of the number.
!
!  Legal input is:
!
!     1 blanks,
!     2 '+' or '-' sign,
!     3 integer part,
!     4 decimal point,
!     5 fraction part,
!     6 'E' or 'e' or 'D' or 'd', exponent marker,
!     7 exponent sign,
!     8 exponent integer part,
!     9 exponent decimal point,
!    10 exponent fraction part,
!    11 blanks,
!    12 final comma,
!
!  with most quantities optional.
!
!  Examples:
!
!    STRING            DVAL
!
!    '1'               1.0D+00
!    '     1   '       1.0D+00
!    '1A'              1.0D+00
!    '12,34,56'        12.0D+00
!    '  34 7'          34.0D+00
!    '-1E2ABCD'        -100.0D+00
!    '-1X2ABCD'        -1.0D+00
!    ' 2E-1'           0.2
!    '23.45'           23.45
!    '-4.2E+2'         -420.0D+00
!    '17d2'            1700.0D+00
!    '-14e-2'         -0.14
!    'e2'              100.0D+00
!    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!  STRING Input, CHARACTER*(*) STRING, the string containing the
!         data to be read.  Reading will begin at position 1 and
!         terminate at the end of the string, or when no more
!         characters can be read to form a legal real.  Blanks,
!         commas, or other nonnumeric data will, in particular,
!         cause the conversion to halt.
!
!  DVAL   Output, real ( kind = 8 ) DVAL, the value that was read
!         from the string.
!
!  IERROR Output, integer IERROR, error flag.
!
!         0, no errors occurred.
!
!         1, 2, 6 or 7, the input number was garbled.  The
!         value of IERROR is the last type of input successfully
!         read.  For instance, 1 means initial blanks, 2 means
!         a plus or minus sign, and so on.
!
!  LCHAR  Output, integer LCHAR, the number of characters read from
!         STRING to form the number, including any terminating
!         characters such as a trailing comma or blanks.
!
  implicit none

  character chrtmp
  real ( kind = 8 ) dval
  integer ierror
  integer ihave
  integer isgn
  integer iterm
  integer jbot
  integer jsgn
  integer jtop
  integer lchar
  logical s_eqi
  integer nchar
  integer ndig
  real ( kind = 8 ) rbot
  real ( kind = 8 ) rexp
  real ( kind = 8 ) rtop
  character ( len = * ) string

  nchar = len(string)

  ierror = 0
  dval = 0.0D+00
  lchar = -1
  isgn = 1
  rtop = 0.0D+00
  rbot = 1.0D+00
  jsgn = 1
  jtop = 0
  jbot = 1
  ihave = 1
  iterm = 0

10    continue

  lchar = lchar+1
  chrtmp = string(lchar+1:lchar+1)
!
!  Blank character.
!
  if ( chrtmp == ' ' ) then
    if ( ihave == 2 .or. ihave==6.or.ihave==7 ) then
      iterm = 1
    else if ( ihave>1 ) then
      ihave = 11
    end if
!
!  Comma
!
  else if ( chrtmp == ',' ) then
    if ( ihave/= 1 ) then
      iterm = 1
      ihave = 12
      lchar = lchar+1
    end if
!
!  Minus sign.
!
  else if ( chrtmp == '-' ) then
    if ( ihave == 1 ) then
      ihave = 2
      isgn = -1
    else if ( ihave == 6 ) then
      ihave = 7
      jsgn = -1
    else
      iterm = 1
    end if
!
!  Plus sign.
!
  else if ( chrtmp == '+' ) then
    if ( ihave == 1 ) then
      ihave = 2
    else if ( ihave == 6 ) then
      ihave = 7
    else
      iterm = 1
    end if
!
!  Decimal point.
!
  else if ( chrtmp == '.' ) then
    if ( ihave<4 ) then
      ihave = 4
    else if ( ihave >= 6 .and. ihave<=8 ) then
      ihave = 9
    else
      iterm = 1
    end if
!
!  Exponent marker.
!
  else if ( s_eqi(chrtmp,'e') .or. s_eqi(chrtmp,'d')  ) then
    if ( ihave<6 ) then
      ihave = 6
    else
      iterm = 1
    end if
!
!  Digit.
!
  else if ( ihave<11 .and.  lge(chrtmp,'0').and.lle(chrtmp,'9')  ) then

    if ( ihave <= 2 ) then
      ihave = 3
    else if ( ihave == 4 ) then
      ihave = 5
    else if ( ihave == 6 .or. ihave==7 ) then
      ihave = 8
    else if ( ihave == 9 ) then
      ihave = 10
    end if

    read(chrtmp,'(i1)')ndig

    if ( ihave == 3 ) then
      rtop = 10*rtop+ndig
    else if ( ihave == 5 ) then
      rtop = 10*rtop+ndig
      rbot = 10*rbot
    else if ( ihave == 8 ) then
      jtop = 10*jtop+ndig
    else if ( ihave == 10 ) then
      jtop = 10*jtop+ndig
      jbot = 10*jbot
    end if
!
!  Anything else is regarded as a terminator.
!
  else
    iterm = 1
  end if
!
!  If we haven't seen a terminator, and we haven't examined the
!  entire string, go get the next character.
!
  if ( iterm/= 1 .and. lchar+1<nchar)go to 10
!
!  If we haven't seen a terminator, and we have examined the
!  entire string, then we're done, and LCHAR is equal to NCHAR.
!
  if ( iterm/= 1 .and. lchar+1 == nchar)lchar=nchar
!
!  Number seems to have terminated.  Have we got a legal number?
!  Not if we terminated in states 1, 2, 6 or 7!
!
  if ( ihave == 1 .or. ihave==2.or.ihave==6.or.ihave==7 ) then
    ierror = ihave
    write ( *, * ) ' '
    write ( *, * ) 'CHRCTD - Fatal error!'
    write ( *, * ) '  Illegal or nonnumeric input!'
    return
  end if
!
!  Number seems OK.  Form it.
!
  if ( jtop == 0 ) then
    rexp = 1.0D+00
  else
    if ( jbot == 1 ) then
      rexp = 10.0D+00**(jsgn*jtop)
    else
      rexp = real ( jsgn * jtop, kind = 8 )
      rexp = rexp / real ( jbot, kind = 8 )
      rexp = 10.0D+00**rexp
    end if
  end if

  dval = real ( isgn, kind = 8 ) * rexp * rtop / rbot

  return
end
subroutine chrcti(string,intval,ierror,lchar)

!*****************************************************************************80
!
!! CHRCTI accepts a STRING of characters and reads an integer
!  from STRING into INTVAL.  The STRING must begin with an integer
!  but that may be followed by other information.
!
!  CHRCTI will read as many characters as possible until it reaches
!  the end of the STRING, or encounters a character which cannot be
!  part of the number.
!
!  Legal input is
!
!    blanks,
!    initial sign,
!    integer part,
!    blanks,
!    final comma,
!
!  with most quantities optional.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2009
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!  STRING Input, CHARACTER*(*) STRING, the string containing the
!         data to be read.  Reading will begin at position 1 and
!         terminate at the end of the string, or when no more
!         characters can be read to form a legal integer.  Blanks,
!         commas, or other nonnumeric data will, in particular,
!         cause the conversion to halt.
!
!         Sample results:
!
!         STRING            INTVAL
!
!         '1'               1
!         '     1   '       1
!         '1A'              1
!         '12,34,56'        12
!         '  34 7'          34
!         '-1E2ABCD'        -100
!         '-1X2ABCD'        -1
!         ' 2E-1'           0
!         '23.45'           23
!
!  INTVAL Output, integer INTVAL, the integer read from the string.
!
!  IERROR Output, integer IERROR, error flag.
!         0 if no errors,
!         Value of IHAVE when error occurred otherwise.
!
!  LCHAR  Output, integer LCHAR, number of characters read from
!         STRING to form the number.
!
  implicit none

  character chrtmp
  integer ierror
  integer ihave
  integer intval
  integer isgn
  integer iterm
  integer itop
  integer lchar
  integer nchar
  integer ndig
  character ( len = * ) string

  nchar = len(string)

  ierror = 0
  intval = 0
  lchar = -1
  isgn = 1
  itop = 0
  ihave = 1
  iterm = 0

10    continue

  lchar = lchar+1
  chrtmp = string(lchar+1:lchar+1)

  if ( chrtmp == ' ' ) then
    if ( ihave == 2 ) then
      iterm = 1
    else if ( ihave == 3 ) then
      ihave