flow3.f90

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

  call lmemry('set','have_fp',lmat)
!
!  Demand that the internal nodes coordinates be reset by SETXY,
!  even if IGRID = 2, for the next calculation.
!
  lval = .true.
  call lmemry('set','need_xy',lval)

  lval = .true.
  call lmemry('set','need_phi',lval)
!
!  Set data specific to the optimization.
!
  xbl = xbleft
  xbr = xbrite
  ybl = ybleft
  ybr = ybrite
!
!  Set the elements and nodes.
!
  if ( xbleft/= xbltar .or. xbrite /= xbrtar .or. ybleft /= ybltar .or. &
    ybrite /= ybrtar ) then

    call node_set ( eqn,ibump,indx,isotri,maxeqn,nelem,neqn,node,np,npe, &
      nx,ny,xbl,xbr)

  end if
!
!  Is this a march?
!
  if ( itype == 1 .or. itype == 2 .or. itype == 4 ) then

    if ( itype == 1 ) then
      ndim = 1
    else if ( itype == 2 ) then
      ndim = 2
    else if ( itype == 4 ) then
      ndim = 3
    end if

    call march ( a,area,base,dir,disjac,dpara3,dparsn,dparfd,dparfdc,dpdyn, &
      dudyn,dvdyn,dydpn,eqn,etan,etaq,g,gdif,gdifc,gold, &
      gradf,gtar,ibc,idfd,ids,ifds,igrid,igunit,ijac,indx,iopt,ipivot,iplot, &
      ipred,ishapb,ishapf,ismooth,isotri,istep1,istep2,itunit,itype,iwrite, &
      jjac,jstep1,jstep2,maxnew,maxstp,ndim,nelem,neqn,nlband,node,np,npar, &
      nparb,nparf,npe,nprof,nrow,nstep3,numel,nx,ny,para,para1,para2,para3, &
      parjac, &
      parnew,parold,phi,res,sens,splbmp,splflo,stpmax,syseqn,taubmp,tauflo, &
      tolnew,wateb,wateb1,wateb2,watep,wateu,watev,wquad,xbl,xbord,xbr,xc, &
      xprof,xquad,xsin,xsiq,ybl,ybord,ybr,yc,yquad)
!
!  ...or a sensitivity optimization?
!
  else if ( itype == 3 ) then

    call qsolve ( a,area,costar,disjac,dopt,dpara3,dparfd,dparfdc,dparsn,dpdyn, &
      dudyn,dvdyn,dydpn,eqn,etan,etaq,g,gdif,gdifc,gold,gopt, &
      gradf,gtar,ibc,idfd,ids,ifds,igrad,igrid,igunit,ijac,indx,iopt,ipivot, &
      iplot,ipred,ishapb,ishapf,ismooth,isotri,itunit,itype,ivopt,iwrite, &
      jjac,liv,lv,maxnew,maxstp,nelem,neqn,nlband,node,nopt,np,npar,nparb, &
      nparf,npe,nprof,nrow,numel,nx,ny,para,para1,parjac,parnew,partar, &
      phi,res,sens,splbmp,splflo,stpmax,syseqn,taubmp,tauflo,tolnew,tolopt, &
      vopt,wateb,watep,wateu,watev,wquad,xbl,xbord,xbr,xc,xopt,xprof,xquad, &
      xsin,xsiq,ybl,ybord,ybr,yc,yquad)
!
!  ...or a one point computation?
!
  else if ( itype == 5 ) then

    call rsolve ( a,area,disjac,dpara3,dparfd,dparfdc,dparsn,dpdyn,dudyn, &
      dvdyn,dydpn,eqn,etan,etaq,g,gdif,gdifc,gradf,gtar, &
      ibc,idfd,ids,ifds,igrid,igunit,ijac,indx,iopt,ipivot,iplot,ipred, &
      ishapb,ishapf,ismooth,isotri,itype,iwrite,jjac,maxnew,nelem,neqn, &
      nlband,node,np,npar,nparb,nparf,npe,nprof,nrow,numel,nx,ny,para,para1, &
      parjac,phi,res,sens,splbmp,splflo,stpmax,syseqn,taubmp, &
      tauflo,tolnew,wateb,watep,wateu,watev,wquad,xbl,xbord,xbr,xc,xprof, &
      xquad,xsin,xsiq,ybl,ybord,ybr,yc,yquad)
!
!  ...or an optimization using function values only?
!
  else if ( itype == 7 ) then

    call osolve ( a,area,disjac,dopt,dpara3,dparfd,dparfdc,dparsn,dpdyn,dudyn, &
      dvdyn,dydpn,eqn,etan,etaq,g,gdif,gdifc,gold,gradf, &
      gtar,ibc,idfd,ids,ifds,igrid,igunit,ijac,indx,iopt,ipivot,iplot,ipred, &
      ishapb,ishapf,ismooth,isotri,itunit,itype,ivopt,iwrite,jjac,liv,lv, &
      maxnew,maxstp,nelem,neqn,nlband,node,nopt,np,npar,nparb,nparf,npe, &
      nprof, &
      nrow,numel,nx,ny,para,para1,parjac,parnew,parold,phi,res,sens,splbmp, &
      splflo,stpmax,syseqn,taubmp,tauflo,tolnew,tolopt,vopt,wateb,watep, &
      wateu,watev,wquad,xbl,xbord,xbr,xc,xopt,xprof,xquad,xsin,xsiq,ybl, &
      ybord,ybr,yc,yquad)

  else

    write ( *, * ) ' '
    write ( *, * ) 'FLOW3 - Fatal error!'
    write ( *, * ) '  Unknown value of ITYPE = ', itype

  end if

  call pr_work

  call after ( plot_file, march_file, igunit, itunit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FLOW3'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine after ( plot_file, march_file, igunit, itunit )

!*****************************************************************************80
!
!! AFTER shuts down files and stops.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
  implicit none

  character ( len = * ) plot_file
  character ( len = * ) march_file
  integer igunit
  integer itunit
!
!  Close the marching file.
!
  if ( itunit /= 0 ) then
    close ( unit = itunit )
    write ( *, * ) ' '
    write ( *, * ) 'AFTER - Closing the marching file "' // &
      trim ( march_file ) // '".'
  end if
!
!  Close the graphics file.
!
  if ( igunit /= 0 ) then
    close ( unit = igunit )
    write ( *, * ) ' '
    write ( *, * ) 'AFTER - Closing the graphics file "' // &
      trim ( plot_file ) // '".'
  end if

  return
end
subroutine bsp ( q, dqdx, dqdy, ielem, iq, nelem, node, np, npe, xc, xq, yc, &
  yq )

!*****************************************************************************80
!
!! BSP evaluates the basis functions associated with pressure.
!
!  Discussion:
!
!    Here is a picture of a typical finite element associated with pressure:
!
!        2
!       /|
!      / |
!     /  |
!    1---3
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = 8 ) Q, the value of the IQ-th basis
!    function at the point with global coordinates (XQ,YQ).
!
!    Output, real ( kind = 8 ) DQDX, DQDY, the X and Y
!    derivatives of the IQ-th basis function at the point
!    with global coordinates (XQ,YQ).
!
!    Input, integer IELEM, the global element number about which
!    we are inquiring.
!
!    Input, integer IQ, the index of the desired basis
!    function.  This is also the node of the reference
!    triangle which is associated with the basis function.
!
!    Basis function IQ is 1 at node IQ, and zero at the
!    other two nodes.
!
!    Input, integer NELEM, the number of elements.
!
!    Input, integer NODE(MAXELM,6).  NODE(I,J) is
!    the global node number of the J-th node in the I-th element.
!
!    Input, integer NP, the number of nodes.
!
!    Input, real ( kind = 8 ) XC(NP), the global X coordinates
!    of the element nodes.
!
!    Input, real ( kind = 8 ) XQ, the global X coordinate of
!    the point in which we are interested.
!
!    Input, real ( kind = 8 ) YC(NP), the global Y coordinates
!    of the element nodes.
!
!    Input, real ( kind = 8 ) YQ, the global Y coordinate of
!    the point in which we are interested.
!
  implicit none

  integer nelem
  integer np
  integer npe

  real ( kind = 8 ) q
  real ( kind = 8 ) dqdx
  real ( kind = 8 ) dqdy
  real ( kind = 8 ) d
  integer i1
  integer i2
  integer i3
  integer ielem
  integer iq
  integer iq1
  integer iq2
  integer iq3
  integer node(nelem,npe)
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) xq
  real ( kind = 8 ) yc(np)
  real ( kind = 8 ) yq

  if ( iq < 1 .or. npe < iq ) then
    write ( *, * ) ' '
    write ( *, * ) 'BSP - Fatal error!'
    write ( *, * ) '  The requested basis function is IQ = ',iq
    write ( *, * ) '  but only values from 1 to 6 are legal.'
    stop
  else if ( iq >= 4 .and. iq <= npe ) then
    q = 0.0D+00
    dqdx = 0.0D+00
    dqdy = 0.0D+00
    return
  end if

  iq1 = iq
  iq2 = mod ( iq, 3 ) + 1
  iq3 = mod ( iq+1, 3 ) + 1

  i1 = node(ielem,iq1)
  i2 = node(ielem,iq2)
  i3 = node(ielem,iq3)

  d = ( xc(i2) - xc(i1) ) * ( yc(i3) - yc(i1) ) &
    - ( xc(i3) - xc(i1) ) * ( yc(i2) - yc(i1) )

  dqdx = ( yc(i2) - yc(i3) ) / d
  dqdy = ( xc(i3) - xc(i2) ) / d

  q = 1.0D+00 + dqdx * ( xq - xc(i1) ) + dqdy * ( yq - yc(i1) )

  return
end
subroutine bump_bc ( dudyn,dvdyn,ip,iparb,ishapb,np,nparb,shape,splbmp,taubmp, &
  ubc,vbc,xc)

!*****************************************************************************80
!
!! BUMP_BC computes the boundary conditions for velocity sensitivities.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 December 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer IP, the global node number.
!
!    Input, integer IPARB, the bump parameter with respect to
!    which the sensitivities are desired.
!
!    Output, real ( kind = 8 ) UBC, the boundary condition for the 
!    horizontal velocity sensitivity.
!
!    Output, real ( kind = 8 ) VBC, the boundary condition for the vertical 
!    velocity sensitivity.
!
  implicit none

  integer np

  real ( kind = 8 ) dudyn(np)
  real ( kind = 8 ) dvdyn(np)
  integer ip
  integer iparb
  integer ishapb
  integer jderiv
  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)
!
!  Determine the value of the "basis" shape function associated
!  with parameter IPARB, evaluated at node IP.
!
  if ( ishapb == 1 ) then
    call plval1 ( iparb+1, nparb+2, xc(ip), taubmp, shape )
  else if ( ishapb == 2 ) then
    call pqval1 ( iparb+1, nparb+2, xc(ip), taubmp, shape )
  else if ( ishapb == 3 ) then
    jderiv = 0
    call ppvalu ( taubmp, splbmp(1,1,iparb), nparb+1, 4, xc(ip), jderiv, shape )
  end if

  ubc = - dudyn(ip) * shape
  vbc = - dvdyn(ip) * shape

  return
end
subroutine bump_bc1 ( g,indx,ip,iparb,ishapb,neqn,np,nparb,shape,splbmp,taubmp, &
  ubc,vbc,xc,yc)

!*****************************************************************************80
!
!! BUMP_BC1 computes the value of the boundary conditions for the
!  horizontal and vertical velocity sensitivities with respect
!  to a given shape parameter using a two point finite difference
!  formula.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    28 December 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) G(MAXEQN), the finite element coefficients.
!
!    Input, integer IP, the global node number.
!
!    Input, integer IPARB, the bump parameter with respect to
!    which the sensitivities are desired.
!
!    Output, real ( kind = 8 ) UBC, the boundary condition for the
!    horizontal velocity sensitivity.
!
!    Output, real ( kind = 8 ) VBC, the boundary condition for the vertical 
!    velocity sensitivity.
!
  implicit none

  integer neqn
  integer np

  real ( kind = 8 ) dudy
  real ( kind = 8 ) dvdy
  real ( kind = 8 ) g(neqn)
  integer ihor
  integer ihorn
  integer indx(np,3)
  integer ip
  integer iparb
  integer ishapb
  integer iver
  integer ivern
  integer jderiv
  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)
!
!  Determine the value of the "basis" shape function associated
!  with parameter IPARB, evaluated at node IP.
!
  if ( ishapb == 1 ) then
    call plval1 ( iparb+1, nparb+2, xc(ip), taubmp, shape )
  else if ( ishapb == 2 ) then
    call pqval1 ( iparb+1, nparb+2, xc(ip), taubmp, shape )
  else if ( ishapb == 3 ) then
    jderiv = 0
    call ppvalu ( taubmp, splbmp(1,1,iparb), nparb+1, 4, xc(ip), jderiv, shape )
  end if
!
!  Estimate dUdY and dVdY at the node.
!
  ihor = indx(ip,1)
  ihorn = indx(ip+1,1)
  dudy = (g(ihorn)-g(ihor))/(yc(ip+1)-yc(ip))

  iver = indx(ip,2)
  ivern = indx(ip+1,2)
  dvdy = (g(ivern)-g(iver))/(yc(ip+1)-yc(ip))
!
!  Set the boundary conditions.
!
  ubc = - dudy * shape
  vbc = - dvdy * shape

  return
end
subroutine bump_bc2 ( g, indx, ip, iparb, ishapb, neqn, np, nparb, shape, &
  splbmp, taubmp, ubc, vbc, xc, yc )

!*****************************************************************************80
!
!! BUMP_BC2 evaluates the velocity sensitivity boundary conditions.
!
!  Discussion:
!
!    The routine evaluates the boundary conditions for the
!    horizontal and vertical velocity sensitivities with respect
!    to a given shape parameter using a three point finite difference
!    formula.
!
!    We derive that formula from the following Taylor series:
!
!      u0 = u0
!      u1 = u0 + h1 u0' + h1**2 u0"/2 + O(h1**3)
!      u2 = u0 + h2 u0' + h2**2 u0"/2 + O(h2**3)
!
!    arriving at:
!
!      u0' = (h1**2 u2 - h2**2 u1 + (h2**2-h1**2) u0) / (h1*h2*(h1-h2))
!        + O(max(h1,h2)**2)
!
!    Note that these Taylor series are really only valid when all three
!    nodes lie in one element, so that U is a smooth polynomial.  This
!    is not true for midside nodes, though to correct this would be
!    painful.
!
!    When all three nodes do lie in one element, and if u is at
!    most a quadratic function, then the formula should be exact.
!    This is the case for the velocities U and V at a node.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) G(MAXEQN), the finite element coefficients.
!
!    Input, integer IP, the global node number.
!
!    Input, integer IPARB, the bump parameter with respect to
!    which the sensitivities are desired.
!
!    Output, real ( kind = 8 ) UBC, the boundary condition for the 
!    horizontal velocity sensitivity.
!
!    Output, real ( kind = 8 ) VBC, the boundary condition for the 
!    vertical velocity sensitivity.
!
  implicit none

  integer neqn
  integer np

  real ( kind = 8 ) dudy
  real ( kind = 8 ) dvdy
  real ( kind = 8 ) g(neqn)