flow3.f90

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

program main

!*****************************************************************************80
!
!! MAIN is the main program for FLOW3.
!
!  Discussion:
!
!    FLOW3 controls a package that solves a fluid flow problem.
!
!    MAXNX =    21,    31,    41,    61,     81,     121,      161,
!    MAXNY =     7,    10,    13,    19,     25,      37,       49,
!    H =       1/4,   1/6,   1/8,  1/12,   1/16,    1/24,     1/32,
!           0.25, 0.166, 0.125, 0.083, 0.0625, 0.04166,  0.03125,
!
!    The assignment of MAXROW should really read (maxrow = 28*min(nx,ny)).
!
!  Usage:
!
!    flow3 < input_file > output_file
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    17 January 2007
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Max Gunzburger,
!    Finite Element Methods for Viscous Incompressible Flows,
!    A Guide to Theory, Practice, and Algorithms,
!    Academic Press, 1989,
!    ISBN: 0-12-307350-2,
!    LC: TA357.G86.
!
  implicit none

  integer, parameter :: maxnx = 41
  integer, parameter :: maxny = 13

  integer, parameter :: liv = 60
  integer, parameter :: maxdim = 3
  integer, parameter :: maxparb = 5
  integer, parameter :: maxparf = 5
  integer, parameter :: npe = 6

  integer, parameter :: maxrow = 28 * maxny
  integer, parameter :: maxelm = 2 * ( maxnx - 1 ) * ( maxny - 1 )
  integer, parameter :: maxeqn = 2*(2*maxnx-1)*(2*maxny-1)+maxnx*maxny
  integer, parameter :: maxnp = (2*maxnx-1)*(2*maxny-1)
  integer, parameter :: maxpar = maxparb + maxparf + 1

  integer, parameter :: lv = 78+maxpar*(maxpar+21)/2

  real ( kind = 8 ) a(maxrow,maxeqn)
  real ( kind = 8 ) area(maxelm,3)
  real ( kind = 8 ) base(maxpar,maxdim)
  real ( kind = 8 ) cost
  real ( kind = 8 ) costar
  real ( kind = 8 ) costb
  real ( kind = 8 ) costp
  real ( kind = 8 ) costu
  real ( kind = 8 ) costv
  real ( kind = 8 ) dir(maxpar,maxdim)
  real ( kind = 8 ) disjac
  real ( kind = 8 ) dopt(maxpar)
  real ( kind = 8 ) dpara3(maxpar)
  real ( kind = 8 ) dparsn(maxpar)
  real ( kind = 8 ) dparfd(maxpar)
  real ( kind = 8 ) dparfdc(maxpar)
  real ( kind = 8 ) dpdyn(maxnp)
  real ( kind = 8 ) dudyn(maxnp)
  real ( kind = 8 ) dvdyn(maxnp)
  real ( kind = 8 ) dydpn(maxnp,maxparb)
  character ( len = 2 ) eqn(maxeqn)
  real ( kind = 8 ) etan(6)
  real ( kind = 8 ) etaq(3)
  character ( len = 80 ) plot_file
  character ( len = 80 ) march_file
  real ( kind = 8 ) flarea
  real ( kind = 8 ) g(maxeqn)
  real ( kind = 8 ) gdif(maxeqn,maxpar)
  real ( kind = 8 ) gdifc(maxeqn,maxpar)
  real ( kind = 8 ) gold(maxeqn)
  real ( kind = 8 ) gopt(maxpar)
  real ( kind = 8 ) gradf(maxeqn,maxpar)
  real ( kind = 8 ) gtar(maxeqn)
  integer i
  integer ibc
  integer ibump
  integer idfd
  integer ids
  integer ierror
  integer ifds
  integer igrad
  integer igrid
  integer igunit
  integer ijac
  integer indx(maxnp,3)
  integer iopt(maxpar)
  integer ip
  integer ipivot(maxeqn)
  integer iplot
  integer ipred
  integer iprdat
  integer ishapb
  integer ishapbt
  integer ishapf
  integer ishapft
  integer ismooth
  integer isotri(maxelm)
  integer istep1
  integer istep2
  integer itar
  integer itemp
  integer itunit
  integer itype
  integer iuval
  integer ival
  integer ivopt(liv)
  integer iwrite
  integer jjac
  integer jstep1
  integer jstep2
  logical lmat
  logical lval
  integer maxnew
  integer maxstp
!
!  TEMPORARY
!
  integer nabor(19,maxnp)
  integer ndim
  integer nelem
  integer neqn
  integer nlband
  integer node(maxelm,npe)
  integer nopt
  integer np
  integer npar
  integer nparb
  integer nparbt
  integer nparf
  integer nprof(2*maxny-1)
  integer nrow
  integer nstep3
  integer numel(maxnp)
  integer numstp
  integer nx
  integer ny
  real ( kind = 8 ) para(maxpar)
  real ( kind = 8 ) para1(maxpar)
  real ( kind = 8 ) para2(maxpar)
  real ( kind = 8 ) para3(maxpar)
  real ( kind = 8 ) parjac(maxpar)
  real ( kind = 8 ) parnew(maxpar)
  real ( kind = 8 ) parold(maxpar)
  real ( kind = 8 ) partar(maxpar)
  real ( kind = 8 ) phi(maxelm,3,6,10)
  real ( kind = 8 ) res(maxeqn)
  real ( kind = 8 ) sens(maxeqn,maxpar)
  real ( kind = 8 ) splbmp(4,maxparb+2,0:maxparb)
  real ( kind = 8 ) splflo(4,maxparf+2,0:maxparf)
  real ( kind = 8 ) stpmax
  character ( len = 20 ) syseqn
  real ( kind = 8 ) taubmp(maxparb+2)
  real ( kind = 8 ) tauflo(maxparf+2)
  character ( len = 80 ) title
  real ( kind = 8 ) tolnew
  real ( kind = 8 ) tolopt
  real ( kind = 8 ) u1
  real ( kind = 8 ) u1max
  real ( kind = 8 ) u2
  real ( kind = 8 ) u2max
  real ( kind = 8 ) vopt(lv)
  real ( kind = 8 ) wateb
  real ( kind = 8 ) wateb1
  real ( kind = 8 ) wateb2
  real ( kind = 8 ) watep
  real ( kind = 8 ) wateu
  real ( kind = 8 ) watev
  real ( kind = 8 ) wquad(3)
  real ( kind = 8 ) xbl
  real ( kind = 8 ) xbleft
  real ( kind = 8 ) xbltar
  real ( kind = 8 ) xbord(maxnx)
  real ( kind = 8 ) xbr
  real ( kind = 8 ) xbrite
  real ( kind = 8 ) xbrtar
  real ( kind = 8 ) xc(maxnp)
  real ( kind = 8 ) xopt(maxpar)
  real ( kind = 8 ) xquad(maxelm,3)
  real ( kind = 8 ) xprof
  real ( kind = 8 ) xsin(6)
  real ( kind = 8 ) xsiq(3)
  real ( kind = 8 ) y2max
  real ( kind = 8 ) ybl
  real ( kind = 8 ) ybleft
  real ( kind = 8 ) ybltar
  real ( kind = 8 ) ybord(maxny)
  real ( kind = 8 ) ybr
  real ( kind = 8 ) ybrite
  real ( kind = 8 ) ybrtar
  real ( kind = 8 ) yc(maxnp)
  real ( kind = 8 ) yquad(maxelm,3)
  real ( kind = 8 ) yval

  call timestamp ( )

  iprdat = 0
  ival = 0
!
!  Print the program name, date, and computer name.
!
  write ( *, * ) ' '
  write ( *, * ) 'FLOW3'
  write ( *, * ) '  FORTRAN90 version'
  write ( *, * ) '  This is the version of 01 January 2001.'
  write ( *, * ) '  The maximum problem size is ',maxnx,' by ',maxny
!
!  Initialize the variables.
!
  call init ( area,cost,costar,costb,costp,costu,costv,disjac,dopt,dpara3, &
    dparsn,dparfd,dparfdc,dpdyn,dudyn,dvdyn,dydpn,eqn,etan,plot_file, &
    march_file,g,gdif,gold,gopt,gradf,gtar,ibc,ibump,idfd,ids,ierror, &
    ifds,igrad,igrid,igunit,ijac,indx,iopt,iplot,ipred,ishapb,ishapbt, &
    ishapf,ishapft,ismooth,isotri,istep1,istep2,itar,itunit,itype,ivopt, &
    iwrite,jjac,jstep1,jstep2,liv,lv,maxelm,maxeqn,maxnew,maxnp,maxnx,maxny, &
    maxpar,maxparb,maxstp,node,nopt,npar,nparb,nparf,npe,nstep3,nx,ny,para1, &
    para2,para3,parjac,partar,sens,stpmax,syseqn,tolnew,tolopt,vopt,wateb, &
    wateb1,wateb2,watep,wateu,watev,xbleft,xbltar,xbord,xbrite,xbrtar,xprof, &
    xsin,ybleft,ybltar,ybord,ybrite,ybrtar)
!
!  Read the user input.
!
  call input ( disjac,plot_file,march_file,ibc,ibump,idfd,ids,ierror,ifds,igrad, &
    igrid,ijac,iopt,iplot,ipred,ishapb,ishapbt,ishapf,ishapft,ismooth,istep1, &
    istep2,itar,itype,iwrite,jjac,jstep1,jstep2,maxnew,maxnx,maxny,maxpar, &
    maxstp,npar,nparb,nparf,nstep3,nx,ny,para1,para2,para3,partar,stpmax, &
    syseqn,tolnew,tolopt,wateb,wateb1,wateb2,watep,wateu,watev,xbleft, &
    xbltar,xbord,xbrite,xbrtar,xprof,ybleft,ybltar,ybord,ybrite,ybrtar)

  if ( ierror == 0 ) then
    call imemry ( 'get', 'Points_Opt', ival )
    ival = 1
    call imemry ( 'inc', 'Restarts', ival )
    write ( *, * ) ' '
    write ( *, * ) 'Flow - GO signal from user input.'
  else
    write ( *, * ) ' '
    write ( *, * ) 'Flow - STOP signal from user.'
    call pr_work
    call after ( plot_file, march_file, igunit, itunit )
  end if
!
!  Set the number of optimization variables.
!
  nopt = sum ( iopt(1:npar) )
  nelem = 2 * ( nx - 1 ) * ( ny - 1 )
  np = ( 2 * nx - 1 ) * ( 2 * ny - 1 )

  write ( *, * ) ' '
  write ( *, * ) 'FLOW3 - Note.'
  write ( *, * ) '  Number of elements  = ', nelem
  write ( *, * ) '  Number of nodes =   ', np
!
!  Check the data.
!
  call chkdat ( ibump,idfd,ids,ifds,igrad,ijac,iopt,ipred,itype,jjac,maxpar, &
    maxparb,maxparf,nopt,npar,nparb,nparf,nstep3,para1,partar,xbleft,xbrite)
!
!  Print the information that the user can change via the input
!  file.  But only do this once.  The output from this routine
!  is extensive, and tedious to see repeatedly.
!
  if ( iprdat == 0 ) then

    iprdat = iprdat+1

    call pr_dat ( disjac,plot_file,march_file,ibc,ibump,idfd,ids,ifds,igrad, &
      igrid,ijac,iopt,iplot,ipred,ishapb,ishapbt,ishapf,ishapft,ismooth, &
      istep1,istep2,itar,itype,iwrite,jjac,jstep1,jstep2,maxnew,maxpar, &
      maxstp,nopt,npar,nparb,nparf,nstep3,nx,ny,para1,para2,para3,partar, &
      stpmax,syseqn,tolnew,tolopt,wateb,wateb1,wateb2,watep,wateu,watev, &
      xbleft,xbltar,xbord,xbrite,xbrtar,xprof,ybleft,ybltar,ybord,ybrite,ybrtar)

  end if
!
!  Open the plot file.
!
  call plot_file_open ( plot_file, igunit, iplot )
!
!  Open the marching file.
!
  call march_file_open ( march_file, itunit )
!
!  Solve the optimization problem,
!  which involves repeated evaluation of the functional
!  at various flow solutions (U,V,P)+(FLO,BUM,NU_INV).
!
!  Set data specific to target calculations.
!
  lval = .true.
  call lmemry('set','target',lval)

  xbl = xbltar
  xbr = xbrtar
  ybl = ybltar
  ybr = ybrtar
  parnew(1:npar) = partar(1:npar)
!
!  Set up the elements and nodes.
!
  call node_set ( eqn, ibump, indx, isotri, maxeqn, nelem, neqn, node, np, npe, &
    nx, ny, xbl, xbr )
!
!  TEMPORARY
!
!  Set the neighbor array.
!
  call setnab ( nabor, np, nx, ny )
!
!  Find points on velocity profile sampling line.
!
  itemp = nint ( ( 2.0D+00 * real ( nx - 1 ) * xprof ) / 10.0D+00 )

  do i = 1,2*ny-1
    ip = itemp*(2*ny-1)+i
    nprof(i) = ip
  end do
!
!  Get matrix bandwidth.
!
  call setban ( indx, maxrow, nelem, neqn, nlband, node, np, npe, nrow )
!
!  Demand that SETXY set internal node coordinates, even if IGRID = 2.
!
  lval = .true.
  call lmemry ( 'set', 'need_xy', lval )
!
!  Demand that SETBAS set the basis functions.
!
  lval = .true.
  call lmemry ( 'set', 'need_phi', lval )

  if ( itar == 0 ) then

    write ( *, * ) ' '
    write ( *, * ) 'FLOW3'
    write ( *, * ) '  Computing the target solution, GTAR.'

    para(1:npar) = partar(1:npar)

    call flosol ( a,area,disjac,eqn,etan,etaq,flarea,g,ierror,igrid,ijac, &
      indx,ipivot,ishapbt,ishapft,isotri,iwrite,jjac,maxnew,nelem,neqn,nlband, &
      node,np,npar,nparb,nparf,npe,nrow,nx,ny,para,parjac,phi,res,splbmp, &
      splflo, &
      syseqn,taubmp,tauflo,tolnew,wquad,xbl,xbord,xbr,xc,xquad,xsin,xsiq, &
      ybl,ybord,ybr,yc,yquad)

    if ( ierror /= 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'FLOW3 - Fatal error!'
      write ( *, * ) '  Could not compute the target solution.'
      stop
    end if

    gtar(1:neqn) = g(1:neqn)
!
!  Compute the finite coefficient differences dG/dPara for the target solution.
!
    call getgrd ( a,area,cost,disjac,dpara3,eqn,etan,etaq,flarea,g, &
      gdif,gtar,idfd,ierror,igrid,ijac,indx,iopt,ipivot,ishapbt,ishapft, &
      isotri,iwrite,jjac,maxnew,nelem,neqn,nlband,node,np,npar,nparb,nparf, &
      npe,nprof,nrow,nx,ny,para,parjac,phi,res,splbmp,splflo,syseqn,taubmp, &
      tauflo,tolnew,wateb,watep,wateu,watev,wquad,xbl,xbord,xbr,xc,xquad, &
      xsin,xsiq,ybl,ybord,ybr,yc,yquad)

    if ( ierror /= 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'FLOW3 - Fatal error!'
      write ( *, * ) '  GETGRD returns IERROR = ',ierror
      stop
    end if

    numstp = 0
!
!  Compute the cost COST associated with the solution GTAR.
!
    call get_cost ( cost, costb, costp, costu, costv, g, gtar, indx, &
      ishapbt, neqn, np, nparb, nprof, ny, splbmp, taubmp, wateb, watep, &
      wateu, watev, xbl, xbr, ybl, ybr, yc )

    if ( iwrite == 1 ) then
      title = 'Cost associated with target solution.'
      call pr_cost1 ( cost, title )
    else if ( iwrite >= 2 ) then
      title = 'Cost associated with target solution.'
      call pr_cost2 ( cost, costb, costp, costu, costv, title, wateb, watep, &
        wateu, watev )
    end if
!
!  Print stuff out.
!
    if ( ifds /= 0 ) then
      call cost_gradient ( dparfd,g,gtar,indx,ishapbt,neqn,np,npar,nparb, &
        nparf,nprof,ny, &
        gdif,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc)
    end if

    if ( ifds /= 0 ) then
      call cost_gradient ( dparfdc,g,gtar,indx,ishapbt,neqn,np,npar,nparb, &
        nparf,nprof, &
        ny,gdifc,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc)
    end if

    if ( ids /= 0 ) then
      call cost_gradient ( dparsn,g,gtar,indx,ishapbt,neqn,np,npar,nparb, &
        nparf,nprof,ny, &
        sens,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc)
    end if

    title = 'Computed target solution'

    call pr_solution ( dpara3, dparfd, dparfdc, dparsn, g, gdif, gdifc, gtar, &
      idfd, ids, ifds, indx, iwrite, neqn, np, npar, nparb, nparf, nprof, &
      numstp, ny, parnew, sens, title, yc )

    if ( iwrite >= 4 ) then
      call xy_print ( maxnp, np, xc, yc )
    end if

    costar = cost

  else if ( itar == 1 ) then

    write ( *, * ) ' '
    write ( *, * ) 'FLOW3'
    write ( *, * ) '  Calling GETTAR2 to get target data.'

    nparbt = 0
    ishapbt = 0

    call xy_set ( igrid, ishapbt, np, nparbt, nx, ny, splbmp, taubmp, &
      xbltar, xbord, xbrtar, xc, ybltar, ybord, ybrtar, yc )

    g(1:neqn) = 0.0D+00

    u1max = 9.0D+00 / 4.0D+00
    y2max = (6.0D+00-2.0D+00*sqrt(3.0D+00))/4.0D+00
    u2max = (3.0D+00-y2max)*(1.5D+00-y2max)*y2max
!
!  WARNING!  The only reason I can get away with referencing
!  INDX in the main program (where it has first dimension MAXNP)
!  is because I reference column 1.  Very lucky!
!
    do i = 1,2*ny-1
      iuval = indx(nprof(i),1)
      yval = real ( i - 1, kind = 8 )*3.0D+00/ real ( 2 * ny - 2, kind = 8 )
      yc(nprof(i)) = yval
      u1 = (3.0D+00-yval)*yval/u1max
      u2 = (3.0D+00-yval)*(1.5D+00-yval)*yval/u2max
      g(iuval) = partar(1)*u1+partar(2)*u2
    end do
!
!  Compute the cost, relative to zero.
!
    gtar(1:neqn) = 0.0D+00

    call disc_cost ( costp, costu, costv, g, gtar, indx, neqn, np, nprof, ny, yc )

    cost = costu
    costar = cost

    call pr_profile ( g, indx, neqn, np, nprof, ny, yc )

    gtar(1:neqn) = g(1:neqn)

  end if
!
!  Write target information to plot file.
!
  if ( iplot /= 0 ) then

    call plot_file_write ( eqn, g, gdif, igunit, indx, isotri, iwrite, nelem, &
      neqn, node, np, npar, npe, nprof, nx, ny, partar, sens, xc, xprof, yc )

  end if
!
!  Turn off target calculation flag.
!
  lval = .false.
  call lmemry ( 'set', 'target', lval )

  call pr_work
!
!  Zero stuff out.
!
  ival = 0
  call imemry('zero','nothing',ival)

  parjac(1:npar) = 0.0D+00

  lmat = .false.