flow3.f90

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

  real ( kind = 8 ) h1
  real ( kind = 8 ) h2
  integer indx(np,3)
  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 ) u0
  real ( kind = 8 ) u1
  real ( kind = 8 ) u2
  real ( kind = 8 ) ubc
  real ( kind = 8 ) v0
  real ( kind = 8 ) v1
  real ( kind = 8 ) v2
  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.
!
  h1 = yc(ip+1)-yc(ip)
  h2 = yc(ip+2)-yc(ip)

  u0 = g(indx(ip,1))
  u1 = g(indx(ip+1,1))
  u2 = g(indx(ip+2,1))

  dudy = (h1**2*u2 - h2**2*u1 + (h2**2-h1**2)*u0)/(h1*h2*(h1-h2))

  v0 = g(indx(ip,2))
  v1 = g(indx(ip+1,2))
  v2 = g(indx(ip+2,2))

  dvdy = (h1**2*v2 - h2**2*v1 + (h2**2-h1**2)*v0)/(h1*h2*(h1-h2))
!
!  Set the boundary conditions.
!
  ubc = - dudy * shape
  vbc = - dvdy * shape

  return
end
subroutine bump_cost ( costb,ishapb,nparb,splbmp,taubmp,xbl,xbr,ybl,ybr)

!*****************************************************************************80
!
!! BUMP_COST evaluates the cost of the bump control.
!
!  Discussion:
!
!    The bump connects the points (XBL,YBL) and (XBR,YBR).
!
!    Compute its "cost" by comparing its slope to the slope of the
!    straight line that connects those two points.
!
!      COSTB = Integral (XBL <= X <= XBR) (Bump'(X) - Line'(X))**2 dX
!
!    Here, Bump(X) represents the function describing the shape
!    of the bump, and Line(X) represents the straight line which
!    simply joins the two endpoints, (XBL,YBL) and (XBR,YBR).
!
!    This integral is approximated by numerical integration.
!
!    The interval between XBL and XBR is divided into NPARB+1
!    intervals, over each of which the bump's height is described
!    by a cubic.
!
!    For each such interval, pick NQUAD1 quadrature points,
!    evaluate the derivative of the bump function there, and
!    subtract the slope of the straight line.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, real ( kind = 8 ) COSTB, the integral of the difference of the
!    derivatives of the straight line joining the two straight line
!    line segments of the bottom, and the bump that is
!    actually drawn there.
!
!    This measures the cost of bump control.
!
  implicit none

  integer nparb
  integer, parameter :: nquad1 = 5

  real ( kind = 8 ) costb
  real ( kind = 8 ) cprime
  integer i
  integer ishapb
  integer j
  integer jderiv
  real ( kind = 8 ) slope
  real ( kind = 8 ) splbmp(4,nparb+2,0:nparb)
  real ( kind = 8 ) taubmp(nparb+2)
  real ( kind = 8 ) wquad1(nquad1)
  real ( kind = 8 ) xbl
  real ( kind = 8 ) xbr
  real ( kind = 8 ) xleft
  real ( kind = 8 ) xquad1(nquad1)
  real ( kind = 8 ) xrite
  real ( kind = 8 ) xx
  real ( kind = 8 ) ybl
  real ( kind = 8 ) ybr
  real ( kind = 8 ) yvec(nparb+2)

  costb = 0.0D+00

  if ( nparb == 0 ) then
    return
  end if

  if ( xbl >= xbr ) then
    return
  end if
!
!  Get the Gauss weights and abscissas.
!
  call gquad1 ( nquad1, wquad1, xquad1 )
!
!  Get the slope of the line joining the endpoints of the bump.
!
  slope = (ybr-ybl)/(xbr-xbl)

  do i = 1,nparb+1

    xleft = ( real ( nparb - i + 2, kind = 8 ) * xbl &
            + real (         i - 1, kind = 8 ) * xbr ) &
            / real ( nparb +     1, kind = 8 )

    xrite = ( real ( nparb - i + 1, kind = 8 ) * xbl &
            + real (         i,     kind = 8 ) * xbr ) &
            / real ( nparb     + 1, kind = 8)

    do j = 1,nquad1

      xx = 0.5D+00*((1.0D+00+xquad1(j))*xrite+(1.0D+00-xquad1(j))*xleft)

      if ( ishapb == 1 ) then
        yvec(1:nparb+2) = splbmp(1,1:nparb+2,0)
        call pldx(nparb+2,xx,taubmp,cprime,yvec)
      else if ( ishapb == 2 ) then
        yvec(1:nparb+2) = splbmp(1,1:nparb+2,0)
        call pqdx(nparb+2,xx,taubmp,cprime,yvec)
      else if ( ishapb == 3 ) then
        jderiv = 1
        call ppvalu(taubmp,splbmp,nparb+1,4,xx,jderiv,cprime)
      end if

      costb = costb + 0.5D+00 * wquad1(j) * ( xrite - xleft ) &
        * ( cprime - slope )**2

    end do

  end do

  return
end
subroutine bump_der ( dpara, ishapb, npar, nparb, nparf, splbmp, taubmp, &
  wateb, xbl, xbr, ybl, ybr )

!*****************************************************************************80
!
!! BUMP_DER differentiates the bump control cost with respect to the parameters.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    02 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!
  implicit none

  integer npar
  integer nparb
  integer, parameter :: nquad1 = 5

  real ( kind = 8 ) cprime
  real ( kind = 8 ) dpara(npar)
  integer i
  integer ishapb
  integer j
  integer jderiv
  integer k
  integer nparf
  real ( kind = 8 ) slope
  real ( kind = 8 ) splbmp(4,nparb+2,0:nparb)
  real ( kind = 8 ) taubmp(nparb+2)
  real ( kind = 8 ) temp
  real ( kind = 8 ) wateb
  real ( kind = 8 ) wquad1(nquad1)
  real ( kind = 8 ) xbl
  real ( kind = 8 ) xbr
  real ( kind = 8 ) xleft
  real ( kind = 8 ) xquad1(nquad1)
  real ( kind = 8 ) xrite
  real ( kind = 8 ) xx
  real ( kind = 8 ) ybl
  real ( kind = 8 ) ybr
  real ( kind = 8 ) yvec(nparb+2)

  if ( nparb == 0 ) then
    return
  end if

  if ( xbr  <=  xbl )then
    return
  end if
!
!  Get the Gauss weights and abscissas.
!
  call gquad1(nquad1,wquad1,xquad1)
!
!  Get the slope of the straight line connecting the endpoints.
!
  slope = (ybr-ybl)/(xbr-xbl)

  do i = 1,nparb+1

    xleft = ( real ( nparb - i + 2, kind = 8 ) * xbl &
            + real (         i - 1, kind = 8 ) * xbr ) &
            / real ( nparb     + 1, kind = 8 )

    xrite = ( real ( nparb - i + 1, kind = 8 ) * xbl &
            + real (         i,     kind = 8 ) * xbr ) &
            / real ( nparb     + 1, kind = 8 )

    do j = 1,nquad1

      xx = 0.5D+00*((1.0D+00+xquad1(j))*xrite+(1.0D+00-xquad1(j))*xleft)

      if ( ishapb == 1 ) then
        yvec(1:nparb+2) = splbmp(1,1:nparb+2,0)
        call pldx(nparb+2,xx,taubmp,cprime,yvec)
      else if ( ishapb == 2 ) then
        yvec(1:nparb+2) = splbmp(1,1:nparb+2,0)
        call pqdx(nparb+2,xx,taubmp,cprime,yvec)
      else if ( ishapb == 3 ) then
        jderiv = 1
        call ppvalu(taubmp,splbmp,nparb+1,4,xx,jderiv,cprime)
      end if

      do k = 1,nparb

        if ( ishapb == 1 ) then
          call pldx1(k+1,nparb+2,xx,taubmp,temp)
        else if ( ishapb == 2 ) then
          call pqdx1(k+1,nparb+2,xx,taubmp,temp)
        else if ( ishapb == 3 ) then
          jderiv = 1
          call ppvalu(taubmp,splbmp(1,1,k),nparb+1,4,xx,jderiv,temp)
        end if

        dpara(nparf+k) = dpara(nparf+k)+wateb*wquad1(j)*(xrite-xleft)* &
          (cprime-slope)*temp

      end do

    end do

  end do

  return
end
subroutine bump_fx ( area,dudyn,dvdyn,eqn,g,ibc,indx,ipar,ishapb,nelem,neqn, &
  node,np,npar,nparb,nparf,npe,para,phi,res,sens,splbmp,taubmp,xc,yc)

!*****************************************************************************80
!
!! BUMP_FX measures the residual in the bump sensitivity equations.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
  implicit none

  integer nelem
  integer neqn
  integer np
  integer npar
  integer nparf
  integer npe

  real ( kind = 8 ) ar
  real ( kind = 8 ) area(nelem,3)
  real ( kind = 8 ) dpdx
  real ( kind = 8 ) dpdy
  real ( kind = 8 ) dppdx
  real ( kind = 8 ) dppdy
  real ( kind = 8 ) dudx
  real ( kind = 8 ) dudy
  real ( kind = 8 ) dudyn(np)
  real ( kind = 8 ) dupdx
  real ( kind = 8 ) dupdy
  real ( kind = 8 ) dvdx
  real ( kind = 8 ) dvdy
  real ( kind = 8 ) dvdyn(np)
  real ( kind = 8 ) dvpdx
  real ( kind = 8 ) dvpdy
  real ( kind = 8 ) dwidx
  real ( kind = 8 ) dwidy
  character ( len = 2 ) eqn(neqn)
  real ( kind = 8 ) g(neqn)
  integer i
  integer ibc
  integer ielem
  integer ihor
  integer indx(np,3)
  integer ip
  integer ipar
  integer iparb
  integer iprs
  integer iq
  integer iquad
  integer ishapb
  integer iver
  integer node(nelem,npe)
  integer nparb
  real ( kind = 8 ) p
  real ( kind = 8 ) para(npar)
  real ( kind = 8 ) phi(nelem,3,6,10)
  real ( kind = 8 ) pp
  real ( kind = 8 ) qi
  real ( kind = 8 ) res(neqn)
  real ( kind = 8 ) nu_inv
  real ( kind = 8 ) sens(neqn,npar)
  real ( kind = 8 ) shape
  real ( kind = 8 ) splbmp(4,nparb+2,0:nparb)
  real ( kind = 8 ) taubmp(nparb+2)
  real ( kind = 8 ) u
  real ( kind = 8 ) ubc
  real ( kind = 8 ) up
  real ( kind = 8 ) v
  real ( kind = 8 ) vbc
  real ( kind = 8 ) vp
  real ( kind = 8 ) wi
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) yc(np)

  nu_inv = para(nparf+nparb+1)
  iparb = ipar-nparf

  res(1:neqn) = 0.0D+00
!
!  Approximate the integral by summing over all elements.
!
  do ielem = 1, nelem
!
!  Evaluate the integrand at the quadrature points.
!
    do iquad = 1,3

      ar = area(ielem,iquad)
!
!  For the given quadrature point, evaluate P, U, V.
!
      call uvalq ( dpdx,dpdy,dudx,dudy,dvdx,dvdy,g,ielem,indx,iquad,nelem,neqn, &
        node,np,npe,p,phi,u,v)

      call upvalq ( dppdx,dppdy,dupdx,dupdy,dvpdx,dvpdy,sens(1,ipar),ielem, &
        indx,iquad,nelem,neqn,node,np,npe,phi,pp,up,vp)
!
!  The only equations that have a contribution from this element
!  are those associated with basis functions for the element.
!  These, in turn, are associated with the nodes of the element.
!
!  So now we consider each node in the element.
!
      do iq = 1,6

        ip = node(ielem,iq)

        wi = phi(ielem,iquad,iq,1)
        dwidx = phi(ielem,iquad,iq,2)
        dwidy = phi(ielem,iquad,iq,3)
        qi = phi(ielem,iquad,iq,4)

        ihor = indx(ip,1)
        if ( eqn(ihor) == 'U' ) then
          res(ihor) = res(ihor)+ar*(dupdx*dwidx+dupdy*dwidy &
            +nu_inv*(up*dudx+u*dupdx+vp*dudy+v*dupdy+dppdx)*wi)
        else 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
          res(ihor) = sens(ihor,ipar)-ubc
        else if ( eqn(ihor) == 'UI' ) then
          res(ihor) = sens(ihor,ipar)
        else if ( eqn(ihor) == 'UW' ) then
          res(ihor) = sens(ihor,ipar)
        end if

        iver = indx(ip,2)
        if ( eqn(iver) == 'V' ) then
          res(iver) = res(iver)+ar*(dvpdx*dwidx+dvpdy*dwidy &
            +nu_inv*(up*dvdx+u*dvpdx+vp*dvdy+v*dvpdy+dppdy)*wi)
        else 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
          res(iver) = sens(iver,ipar)-vbc
        else if ( eqn(iver) == 'VI' ) then
          res(iver) = sens(iver,ipar)
        else if ( eqn(iver) == 'VW' ) then
          res(iver) = sens(iver,ipar)
        end if

        iprs = indx(ip,3)
        if ( iprs>0 ) then
          if ( nu_inv == 0.0D+00 ) then
            res(iprs) = sens(iprs,ipar)
          else if ( eqn(iprs) == 'P' ) then
            res(iprs) = res(iprs)+ar*(dupdx+dvpdy)*qi
          else if ( eqn(iprs) == 'PB' ) then
            res(iprs) = sens(iprs,ipar)
          end if
        end if

      end do
    end do
  end do

  return
end
subroutine bump_sen ( dudyn,dvdyn,eqn,f,g,ibc,indx,ipar,ishapb,neqn,np,nparb, &
  nparf,splbmp,taubmp,xc,yc)

!*****************************************************************************80
!
!! BUMP_SEN sets up the right hand side F associated with the
!  sensitivities of a given flow solution (U,V,P) with respect to the
!  IPAR-th bump parameter.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    01 January 2001
!
!  Author:
!
!    John Burkardt