channel.f90

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

program main

!*****************************************************************************80
!
!! MAIN is the main program for CHANNEL.
!
!  Discussion:
!
!    CHANNEL solves the channel flow problem using the finite element method.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    20 January 2007
!
!  Author:
!
!    John Burkardt
!
!  Journal:
!
!    4 September 1992
!    Added LONG to GDUMP output.
!
!    2 September 1992:
!    Set g to zero before secant iteration.
!    Made sure to update g BEFORE returning from NSTOKE.
!    New code converges in 3 secant steps, 13 Newton steps,
!    and 48 seconds.
!
!    27 August 1992: Corrected a mistake in XYDUMP that only occurred
!    for NX <= NY.
!
!    26 August 1992: broke out SETLIN, SETQUD and SETBAS for 
!    similarity with BUMP.
!
!    22 August 1992: Moved bandwidth calculations to SETBAN, mainly for
!    the benefit of the BUMP computation.
!
!    20 August 1992: Moved linear system setup and solve out of NSTOKE,
!    merged it with solvlin, and called it "LINSYS".
!
!    19 August 1992.  After I changed the real ( kind = 8 ) of the matrix
!    A in NSTOKE from (MAXROW,NEQN) to (NROW,NEQN) I found performance
!    was degraded.  The time went up from about 47 seconds to 57 seconds!
!    (I made other, also theoretically harmless changes at the same time).
!    Changing back to MAXROW brought the time down to 50 seconds.  What's
!    going on?  This is definitely a performance problem in LAPACK!
!
!    18 August 1992, the residual calculation has been healed, somehow!
!    That means that I can, if I wish, think about a jacobian.
!
!    13 August 1992, the way that the last pressure was set to 1 was not
!    working properly.  I don't know why, but I changed the code so that
!    the last pressure was set to zero, and the contour plots smoothed out,
!    and the secant iteration converged in three steps, rather than four.
!    So I'm also dropping the entire "PRESET" routine, which forced the
!    average pressure to zero.  What's the point?  Perhaps that's important
!    when you're modifying the grid, though.  So I won't actually delete 
!    the code.
!
!    11 August 1992, the original scheme for the last pressure equation
!    does NOT force the average pressure to be zero.  If you don't
!    believe me, repeat the calculation of PMEAN after the adjustment.
!    I've fixed that.
! 
!
!
!  Control for NAVIER-STOKES on channel using primitive variables.
!
!  This program solves a fluid flow problem on the unit square.
!
!  The fluid flow problem is formulated in terms of primitive variables 
!  - u, v, and p.
!
!  This code tries to match a downstream profile by altering
!  one parameter, the value of the inflow parameter at the mid-height
!  of the channel.
!
!  Piecewise linear functions on triangles approximate
!  the pressure and quadratics on triangles approximate the velocity.
!  This is the "Taylor-Hood" finite element basis.
!
!
!
!  Variables
!
!  A      real ( kind = 8 ) A(MAXROW,MAXEQN), the banded matrix 
!         used in NSTOKE and SOLVLIN.  It is stored according to
!         LINPACK/LAPACK general band storage mode.
!
!  AREA   real ( kind = 8 ) AREA(NELEMN), contains the area of each
!         element.
!
!  F      real ( kind = 8 ) F(MAXEQN).  After the call to NSTOKE,
!         F contains the solution to the current Navier Stokes problem.
!
!  FILEG  CHARACTER*30 FILEG, the filename of the graphics data file
!         to be dumped for the DISPLAY program.
!
!  FILEU  CHARACTER*30 FILEU, the filename of the graphics UV data file
!         to be dumped for the PLOT3D program.
!
!  FILEX  CHARACTER*30 FILEX, the filename of the graphics XY data file
!         to be dumped for the PLOT3D program.
!
!  G      real ( kind = 8 ) G(MAXEQN).  After the call to SOLVLIN,
!         G contains the sensitivities.
!
!  INDX   INTEGER INDX(NP,2), records, for each node, whether
!         any velocity unknowns are associated with the node.
!
!         INDX(I,1) records information for horizontal velocities.
!         If it is 0, then no unknown is associated with the
!         node, and a zero horizontal velocity is assumed.
!         If it is -1, then no unknown is associated with the
!         node, but a velocity is specified via the UBDRY routine
!         If it is positive, then the value is the index in the
!         F and G arrays of the unknown coefficient.
!
!         INDX(I,2) records information for vertical velocities.
!         If it is 0, then no unknown is associated with the
!         node, and a zero vertical velocity is assumed.
!         If it is positive, then the value is the index in the
!         F and G arrays of the unknown coefficient.
!
!  INSC   INTEGER INSC(NP), records, for each node, whether an
!         unknown pressure is associated with the node.
!
!         If INSC(I) is zero, then no unknown is associated with the
!         node, and a pressure of 0 is assumed.  
!         If INSC(I) is positive, then the value is the index in the
!         F and G arrays of the unknown coefficient.
!
!  IOUNIT INTEGER IOUNIT, the FORTRAN unit number for the file
!         into which graphics data will be written.
!
!  IPIVOT INTEGER IPIVOT(MAXEQN), pivot information used by the linear
!         solver.
!
!  ISOTRI INTEGER ISOTRI(NELEMN), records whether or not a given 
!         element is isometric or not.
!
!         0, element is not isometric.
!         1, element is isometric.
!
!  IVUNIT Input, INTEGER IVUNIT, the output unit to which the data is to
!         be written.  The data file is unformated.
!
!  IWRITE INTEGER IWRITE, controls the amount of output produced by the
!         program.
!
!         0, minimal output.
!         1, normal output, plus graphics data files created by
!            GDUMP, UVDUMP and XYDUMP.
!         2, copious output.
!
!  IXUNIT Input, INTEGER IXUNIT, the output unit to which the data is to
!         be written.  The data file is unformated.
!
!  LONG   LOGICAL LONG, 
!
!         .TRUE. if region is "long and thin", and
!         .FALSE. if region is "tall and skinny".
!
!         This determines how we number nodes, elements, and variables.
!
!  MAXEQN INTEGER MAXEQN, the maximum number of equations allowed.
!
!  MAXNEW INTEGER MAXNEW, the maximum number of Newton steps per iteration.
!
!  MAXROW INTEGER MAXROW, the maximum row real ( kind = 8 ) of the coefficient
!         matrix that is allowed.
!
!  MAXSEC INTEGER MAXSEC, the maximum number of secant steps allowed.
!
!  MX     INTEGER MX, MX = 2*NX-1, the total number of grid points on
!         the horizontal side of the region.
!
!  MY     INTEGER MY, MY = 2*NY-1, the total number of grid points on
!         the vertical side of the region.
!
!  NELEMN INTEGER NELEMN, the number of elements used.   
!
!  NEQN   INTEGER NEQN, the number of equations or functions for
!         the full system.
!
!  NLBAND INTEGER NLBAND, the number of diagonals below the main diagonal 
!         of the matrix A which are nonzero. 
!
!  NNODES INTEGER NNODES, the number of nodes per element, 6.
!
!  NODE   INTEGER NODE(NELEMN,6), records the global node numbers of the
!         6 nodes that make up each element.
!
!  NODEX0 INTEGER NODEX0, the lowest numbered node in the column of
!         nodes where the profile is measured.
!
!  NP     INTEGER NP, the number of nodes.
!
!  NQUAD  INTEGER NQUAD, the number of quadrature points, currently
!         set to 3.
!
!  NROW   INTEGER NROW, the used row real ( kind = 8 ) of the coefficient
!         matrix.
!
!  NUMNEW INTEGER NUMNEW, total number of Newton iterations taken in
!         the NSTOKE routine during the entire run.
!
!  NUMSEC INTEGER NUMSEC, the number of secant steps taken.
!
!  NX     INTEGER NX, the number of "main" grid points on the horizontal 
!         side of the region.
!
!  NY     INTEGER NY, the number of "main" grid points on the vertical 
!         side of the region.
!
!  PHI    REAL PHI(NELEMN,NQUAD,NNODES,3).  Each entry of PHI contains
!         the value of a quadratic basis function or its derivative, 
!         evaluated at a quadrature point.
!
!         In particular, PHI(I,J,K,1) is the value of the quadratic basis
!         function associated with local node K in element I, evaluatated
!         at quadrature point J.
!
!         PHI(I,J,K,2) is the X derivative of that same basis function,
!         PHI(I,J,K,3) is the Y derivative of that same basis function.
!
!  PSI    REAL PSI(NELEMN,NQUAD,NNODES).  Each entry of PSI contains
!         the value of a linear basis function evaluated at a 
!         quadrature point.
!
!         PSI(I,J,K) is the value of the linear basis function associated
!         with local node K in element I, evaluated at quadrature point J.
!
!  RES    REAL RES(MAXEQN), contains the residuals.
!
!  REYNLD REAL REYNLD, the value of the Reynolds number.  In the
!         program's system of units, viscosity = 1 / REYNLD.
!
!  RJPNEW REAL RJPNEW, the derivative with respect to the
!         parameter A of the functional J.
!
!  TARRAY REAL TARRAY(2), an array needed in order to store 
!         results of a call to the UNIX CPU timing routine ETIME.
!
!  TOLNEW REAL TOLNEW, the convergence tolerance for the Newton
!         iteration in NSTOKE.
!
!  TOLSEC REAL TOLSEC, the convergence tolerance for the secant
!         iteration in the main program.
!
!  XC     REAL XC(NP), XC(I) is the X coordinate of node I.
!
!  XLNGTH REAL XLNGTH, the length of the region.
!
!  XM     REAL XM(NELEMN,NQUAD), XM(IT,I) is the X coordinate of the 
!         I-th quadrature point in element IT.
!
!  YC     REAL YC(NP), YC(I) is the Y coordinate of node I.
!
!  YLNGTH REAL YLNGTH, the height of the region.
!
!  YM     REAL YM(NELEMN,NQUAD).  YM(IT,I) is the Y coordinate of the
!         I-th quadrature point in element IT.
!
!
!  The Navier Stokes equations:
!
!  The primitive variable formulation involves horizontal velocity U,
!  vertical velocity V, and pressure P.  The equations are:
!
!    U dUdx + V dUdy + dPdx - mu*(ddU/dxdx + ddU/dydy) = F1
!
!    U dVdx + V dVdy + dPdy - mu*(ddV/dxdx + ddV/dydy) = F2
!
!    dUdx + dVdy = 0
!
!
!
!  Finite element form of Navier Stokes equations:
!
!  When reformulated into finite element form, with PHI(i) being the I-th
!  common basis function for U and V, and PSI(i) the I-th basis function 
!  for P, these equations become:
!
!  Integral (U dUdx PHI(i) + V dUdy PHI(I) - P dPHI(i)/dx 
!            + mu (dUdx dPHI(i)/dx + dUdy dPHI(i)/dy) ) = 
!  Integral (F1 * PHI(i))
!
!  Integral (U dVdx PHI(i) + V dVdy PHI(I) - P dPHI(i)/dy 
!            + mu (dVdx dPHI(i)/dx + dVdy dPHI(i)/dy) ) = 
!  Integral (F2 * PHI(i))
!
!  Integral (dUdx PSI(i) + dVdx PSI(i)) = 0
!
  implicit none

  integer, parameter :: nx = 21
  integer, parameter :: ny = 7

  integer, parameter :: maxrow = 27*ny
  integer, parameter :: nelemn = 2*(nx-1)*(ny-1)
  integer, parameter :: mx = 2*nx-1
  integer, parameter :: my = 2*ny-1
  integer, parameter :: np = mx*my
  integer, parameter :: maxeqn = 2*mx*my+nx*ny
  integer, parameter :: nnodes = 6
  integer, parameter :: nquad = 3

  real ( kind = 8 ) a(maxrow,maxeqn)
  real ( kind = 8 ) a2
  real ( kind = 8 ) abound
  real ( kind = 8 ) anew
  real ( kind = 8 ) aold
  real ( kind = 8 ) area(nelemn)
  real ( kind = 8 ) dcda(my)
  real ( kind = 8 ) f(maxeqn)
  character fileg*30
  character fileu*30
  character filex*30
  real ( kind = 8 ) g(maxeqn)
  real ( kind = 8 ) gr(my,my)
  integer i
  integer iline(my)
  integer indx(np,2)
  integer insc(np)
  integer iounit
  integer ipivot(maxeqn)
  integer isotri(nelemn)
  integer iter
  integer ivunit
  integer iwrite
  integer ixunit
  integer j
  logical long
  integer maxnew
  integer maxsec
  integer nband
  integer neqn
  integer nlband
  integer node(nelemn,nnodes)
  integer nodex0
  integer npara
  integer nrow
  integer numnew
  integer numsec
  real ( kind = 8 ) para
  real ( kind = 8 ) phi(nelemn,nquad,nnodes,3)
  real ( kind = 8 ) psi(nelemn,nquad,nnodes)
  real ( kind = 8 ) r(my)
  real ( kind = 8 ) res(maxeqn)
  real ( kind = 8 ) reynld
  real ( kind = 8 ) rjpnew
  real ( kind = 8 ) rjpold
  real ( kind = 8 ) rtemp(my)
  real tarray(2)
  real ( kind = 8 ) temp
  real ( kind = 8 ) test
  real ( kind = 8 ) tolnew
  real ( kind = 8 ) tolsec
  real ( kind = 8 ) ui(my)
  real ( kind = 8 ) unew(my)
  real ( kind = 8 ) xc(np)
  real ( kind = 8 ) xlngth
  real ( kind = 8 ) xm(nelemn,nquad)
  real ( kind = 8 ) yc(np)
  real ( kind = 8 ) ylngth
  real ( kind = 8 ) ym(nelemn,nquad)

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'CHANNEL'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Channel flow control problem'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Last modified:'
  write ( *, '(a)' ) '    4 September 1992.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Flow control problem:'
  write ( *, '(a)' ) '    Inflow controlled by one parameter.'
  write ( *, '(a)' ) '    Velocities measured along vertical line.'
  write ( *, '(a)' ) '    Try to match specified velocity profile.'
!
!  Set input data
!
  fileg = 'display.txt'
  fileu = 'uv.txt'
  filex = 'xy.txt'
  iounit = 2
  ivunit = 4
  iwrite = 2
  ixunit = 3
  maxnew = 10
  maxsec = 8
  npara = 1
  numnew = 0
  numsec = 0
  reynld = 1.0D+00
  rjpnew = 0.0D+00
  tolnew = 1.0D-04
  tolsec = 1.0D-06
  xlngth = 10.0D+00
  ylngth = 3.0D+00

  write (*,*) ' '
  write (*,*) 'NX = ', nx
  write (*,*) 'NY = ', ny
  write (*,*) 'Number of elements = ', nelemn
  write (*,*) 'Reynolds number = ', reynld
  write (*,*) 'Secant tolerance = ', tolsec
  write (*,*) 'Newton tolerance = ', tolnew
  write (*,*) ' '
!
!  SETGRD constructs grid, numbers unknowns, calculates areas,
!  and points for midpoint quadrature rule.
!
  call setgrd (indx, insc, isotri, iwrite, long, maxeqn, mx, my, &
    nelemn, neqn, nnodes, node, np, nx, ny )
!
!  Compute the bandwidth
!
  call setban(indx,insc,maxrow,nband,nelemn,nlband,nnodes, &
    node,np,nrow)
!
!  Record variable numbers along profile sampling line.
!
  call setlin(iline,indx,iwrite,long,mx,my,nodex0,np, &
    nx,ny,xlngth)
!
!  Set the coordinates of grid points.
!
  call setxy(iwrite,long,mx,my,np,nx,ny,xc,xlngth,yc,ylngth)
!
!  Set quadrature points
!
  call setqud(area,nelemn,nnodes,node,np,nquad,xc,xm,yc,ym)
!
!  Evaluate basis functions at quadrature points
!
  call setbas(nelemn,nnodes,node,np,nquad,phi,psi,xc,xm,yc,ym)
!
!  NSTOKE now solves the Navier Stokes problem for an inflow 
!  parameter of 1.0.
!
  para = 1.0D+00
  write (*,*) ' '
  write (*,*) 'Solve Navier Stokes problem with parameter = ',para
  write (*,*) 'for profile at x = ', xc(nodex0)
  g(1:neqn) = 1.0D+00

  call nstoke (a,area,f,g,indx,insc,ipivot,iwrite, &
    maxnew,maxrow,nelemn,neqn,nlband,nnodes,node, &
    np,nquad,nrow,numnew,para,phi,psi,reynld,tolnew,yc)
!
!  RESID computes the residual at the given solution
!
  if ( 1 <= iwrite ) then
    call resid (area,f,indx,insc,iwrite,nelemn,neqn, &
      nnodes,node,np,nquad,para,phi,psi,res,reynld,yc)
  end if
!
!  GETG computes the internal velocity profile at X = XC(NODEX0), which will
!  be used to measure the goodness-of-fit of the later solutions.
!
  call getg ( f, iline, my, neqn, ui )

  if ( 1 <= iwrite ) then
    write (*,*) ' '
    write (*,*) 'U profile:'
    write (*,*) ' '
    write (*,'(5g14.6)') ui(1:my)
  end if
!
!  GRAM generates the Gram matrix GR and the vector 
!  R = line integral of ui*phi
!
  call gram (gr,iline,indx,iwrite,my,nelemn,nnodes,node, &
    nodex0,np,para,r,ui,xc,yc)
!
!  GDUMP dumps information for graphics display by DISPLAY.
!
  if ( .false. ) then
    write (*,*) 'Writing graphics data to file '//fileg
    call delete(fileg)
    open (unit = iounit,file=fileg,form='formatted',status='new', &
      err = 50)
    rjpnew = 0.0D+00
    call gdump (f,indx,insc,iounit,isotri,long,nelemn,neqn, &
      nnodes,node,np,npara,nx,ny,para,reynld,rjpnew,xc,yc)
  end if
!
!  Write the XY data to a file.
!
  if ( .false. ) then
    call delete(filex)
    open(unit = ixunit,file=filex,form='formatted',status='new')
    call xy_plot3d (ixunit,long,np,nx,ny,xc,yc)
    close(unit = ixunit)
  else
    call delete ( filex )
    open ( unit = ixunit, file = filex, form = 'formatted', &
      status = 'new')
    call xy_table ( ixunit, np, xc, yc )
    close ( unit = ixunit )
  end if
!
!  Write the velocity data to a file.
!
  if ( .false. ) then
    call delete(fileu)
    open(unit = ivunit,file=fileu,form='formatted',status='new')
    call uv_plot3d (f,indx,insc,ivunit,long,mx,my, &
      nelemn,neqn,nnodes,node,np,para,a,reynld,yc)
    close(unit = ivunit)
  else
    call delete(fileu)
    open(unit = ivunit,file=fileu,form='formatted',status='new')
    call uv_table ( f, indx, ivunit, neqn, np, para, yc )
    close(unit = ivunit)
  end if
!
!  Destroy information about true solution
!
  f(1:neqn) = 0.0D+00
  g(1:neqn) = 0.0D+00
!
!  Secant iteration loop
!
  aold = 0.0D+00
  rjpold = 0.0D+00
  anew = 0.1D+00

  do iter = 1, maxsec

    numsec = numsec+1
    write (*,*) ' '
    write (*,*) 'Secant iteration ',iter
!
!  Solve for unew at new value of parameter anew
!
    write (*,*) ' '
    write (*,*) 'Solving Navier Stokes problem for parameter = ',anew
!
!  Use solution F at previous value of parameter for starting point.
!
    call dcopy(neqn,f,1,g,1)
    para = anew

    call nstoke (a,area,f,g,indx,insc,ipivot,iwrite, &
      maxnew,maxrow,nelemn,neqn,nlband,nnodes,node, &
      np,nquad,nrow,numnew,para,phi,psi,reynld,tolnew,yc)
!
!  Get velocity profile
!
    call getg ( f, iline, my, neqn, unew )

    if ( 1 <= iwrite ) then
      write (*,*) ' '
      write (*,*) 'Velocity profile:'
      write (*,*) ' '
      write (*,'(5g14.6)') unew(1:my)