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 ! ! 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. 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. ! ! 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 ! ! 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. ! ! 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. ! ! 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. ! ! 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) 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. ! ! 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. ! ! 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. ! ! 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. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none integer neqn integer np integer nparf real ( kind = 8 ) dudyn(np) real ( kind = 8 ) dvdyn(np) character ( len = 2 ) eqn(neqn) real ( kind = 8 ) f(neqn) real ( kind = 8 ) g(neqn) integer i integer ibc integer ihor integer indx(np,3) integer ip integer ipar integer iparb integer ishapb integer iver 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) iparb = ipar-nparf f(1:neqn) = 0.0D+00 do ip = 1,np ihor = indx(ip,1) 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 f(ihor) = ubc end if iver = indx(ip,2) 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 f(iver) = vbc end if end do return end subroutine bump_spl ( ishapb, npar, nparb, nparf, par, splbmp, taubmp, xbl, & xbr, ybl, ybr ) !*****************************************************************************80 ! !! BUMP_SPL sets up or updates the spline data that describes the bump. ! ! Modified: ! ! 22 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ISHAPB. ! 1, the bump is modeled by C0 linear splines. ! 2, the bump is modeled by C0 quadratic splines. ! 3, the bump is modeled by C1 cubic splines. ! implicit none integer npar integer nparb integer nparf integer i integer ibcbeg integer ibcend integer ipar integer ishapb real ( kind = 8 ) par(npar) real ( kind = 8 ) splbmp(4,nparb+2,0:nparb) real ( kind = 8 ) taubmp(nparb+2) real ( kind = 8 ) xbl real ( kind = 8 ) xbr real ( kind = 8 ) ybl real ( kind = 8 ) ybr if ( nparb <= 0 ) then return end if ! ! Set up the bump arrays, including: ! ! TAUBMP, containing the abscissas, which never change, ! SPLBMP(1,I,0), the location of the bump at abscissa I, ! SPLBMP(1,I,IPAR), the I-th coefficient of the partial derivative ! of the bump shape with respect to the IPAR-th bump parameter. ! do i = 1, nparb+2 taubmp(i) = ( real ( nparb + 2 - i ) * xbl + real ( i - 1 ) * xbr ) & / real ( nparb + 1 ) end do do i = 1, nparb+2 if ( i == 1 ) then splbmp(1,i,0) = ybl else if ( 2 <= i .and. i <= nparb+1 ) then splbmp(1,i,0) = par(nparf+i-1) else if ( i == nparb+2 ) then splbmp(1,i,0) = ybr end if end do if ( ishapb == 3 ) then do i = 1, nparb+2 do ipar = 1, nparb if ( ipar+1 /= i ) then splbmp(1,i,ipar) = 0.0D+00 else splbmp(1,i,ipar) = 1.0D+00 end if end do end do ibcbeg = 0 ibcend = 0 do i = 0, nparb call cubspl ( taubmp, splbmp(1,1,i), nparb+2, ibcbeg, ibcend ) end do end if return end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end subroutine chkdat(ibump,idfd,ids,ifds,igrad,ijac,iopt,ipred,itype,jjac, & maxpar,maxparb,maxparf,nopt,npar,nparb,nparf,nstep3,para1,partar,xbleft, & xbrite) !*****************************************************************************80 ! !! CHKDAT performs some simple checks on the input data. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IGRAD, the cost gradient approximation option. ! 0, No cost gradient approximation is made. ! 1, Chain rule on discretized sensitivities. ! 2, Chain rule on finite coefficient differences. ! 3, Chain rule on corrected finite coefficient differences. ! 4, Direct finite cost differences. ! implicit none integer maxpar integer i integer ibump integer idfd integer ids integer ifds integer igrad integer ijac integer iopt(maxpar) integer ipred integer isum integer itype integer jjac integer maxparb integer maxparf integer nopt integer npar integer nparb integer nparf integer nstep3 real ( kind = 8 ) para1(maxpar) real ( kind = 8 ) partar(maxpar) real ( kind = 8 ) psum real ( kind = 8 ) xbleft real ( kind = 8 ) xbrite ! ! IBUMP ! if ( ibump < 0 .or. ibump > 3 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IBUMP is out of bounds!' write ( *, * ) ' Current value of IBUMP = ',ibump write ( *, * ) ' Legal values must be between 0 and 3.' stop end if ! ! IGRAD ! if ( igrad < 0 .or. igrad > 4 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IGRAD must be between 0 and 4,' write ( *, * ) ' but your value is ',igrad stop end if if ( igrad == 0 ) then if ( itype == 3 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' A cost gradient approximation MUST be made' write ( *, * ) ' when ITYPE = ',itype write ( *, * ) ' but you have chosen IGRAD = ',igrad stop end if else if ( igrad == 1 ) then if ( ids == 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Your cost gradient choice IGRAD = ',igrad write ( *, * ) ' requires a nonzero value of IDS,' write ( *, * ) ' but your value is IDS = ',ids stop end if else if ( igrad == 2 ) then if ( ifds == 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Your cost gradient choice IGRAD = ',igrad write ( *, * ) ' requires a nonzero value of IFDS,' write ( *, * ) ' but your value is IFDS = ',ifds stop end if else if ( igrad == 3 ) then if ( ifds == 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Your cost gradient choice IGRAD = ',igrad write ( *, * ) ' requires a nonzero value of IFDS,' write ( *, * ) ' but your value is IFDS = ',ifds stop end if else if ( igrad == 4 ) then if ( idfd == 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Your cost gradient choice IGRAD = ',igrad write ( *, * ) ' requires a nonzero value of IDFD,' write ( *, * ) ' but your value is IDFD = ',idfd stop end if end if ! ! IJAC ! if ( ijac<1 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IJAC is out of bounds!' write ( *, * ) ' Current value of IJAC = ',ijac write ( *, * ) ' Legal values must be 1 or greater.' stop end if ! ! IOPT ! do i = 1,npar if ( iopt(i)/= 0 .and. iopt(i)/=1 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IOPT(*) must be 0 or 1, but' write ( *, * ) ' IOPT(',I,') = ',iopt(i) stop end if end do ! ! IPRED ! if ( ipred<0 .or. ipred>4 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IPRED must be between 0 and 4,' write ( *, * ) ' but your value is ',ipred stop end if if ( ipred == 1 .and. ids==0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IPRED = ',ipred write ( *, * ) ' but IDS = ',ids stop else if ( ipred == 2 .and. ifds==0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IPRED = ',ipred write ( *, * ) ' but IFDS = ',ifds stop else if ( ipred == 3 .and. ifds==0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IPRED = ',ipred write ( *, * ) ' but IFDS = ',ifds stop else if ( ipred == 4 .and. ids==0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' IPRED = ',ipred write ( *, * ) ' but IDS = ',ids stop end if ! ! ITYPE ! if ( itype<1 .or. itype>7 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Illegal value of ITYPE = ',itype write ( *, * ) ' Legal values are between 1 and 7.' write ( *, * ) ' ChkDat forces a STOP!' stop end if ! ! JJAC ! if ( jjac<0 .or. jjac>3 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Illegal value of JJAC = ',jjac write ( *, * ) ' Legal values are between 0 and 3.' write ( *, * ) ' ChkDat forces a STOP!' stop end if ! ! NOPT ! if ( nopt <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' The number of free parameters, NOPT = ',nopt write ( *, * ) ' but this value must be positive!' stop end if ! ! NPAR ! if ( npar>maxpar ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' NPAR is out of bounds!' write ( *, * ) ' Current value of NPAR = ',npar write ( *, * ) ' Maximum legal value, MAXPARA = ',maxpar stop end if ! ! NPARB ! if ( nparb<0 .or. nparb>maxparb ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' NPARB is out of bounds!' write ( *, * ) ' Input value of NPARB = ',nparb write ( *, * ) ' Maximum legal value, MAXPARB = ',maxparb stop end if ! ! NPARF ! if ( nparf <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' NPARF is out of bounds!' write ( *, * ) ' The input value of NPARF is ',nparf write ( *, * ) ' But NPARF must be at least 1.' stop end if if ( nparf>maxparf ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' NPARF is too big!' write ( *, * ) ' Input value of NPARF = ',nparf write ( *, * ) ' Maximum legal value, MAXPARF = ',maxparf stop end if ! ! NSTEP3 ! if ( itype == 4 ) then if ( nstep3<1 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' Nonpositive value for NSTEP3 = ',nstep3 stop else if ( nstep3 == 1 ) then write ( *, * ) 'CHKDAT - Warning!' write ( *, * ) ' NSTEP3 = 1 is an unusual value!' end if end if ! ! PARA1(1:NPARF) ! if ( itype == 3 ) then psum = sum ( abs ( para1(1:nparf) ) ) if ( psum == 0.0D+00 ) then isum = 0 do i = 1,nparf isum = isum+abs(iopt(i)) end do if ( isum == 0 ) then write ( *, * ) ' ' write ( *, * ) 'CHKDAT - Fatal error!' write ( *, * ) ' You are optimizing, ITYPE = ',itype write ( *, * ) ' But all inflows are zero,' write ( *, * ) ' and no inflow may vary.' stop end if end if end if ! ! PARA1(NPARF+NPARB+1), the value of NU_INV. ! 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. ! ! 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) ! ! 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+11)go to 110 ! ! Not-a-knot and N >= 3, and either N>3 or also not-a-knot ! at left end point. ! if ( n/= 3 .or. ibcbeg/=0 ) then g = c(3,n-1)+c(3,n) c(2,n) = ((c(3,n)+2.0D+00*g)*c(4,n)*c(3,n-1)+c(3,n)**2*(c(1,n-1)-c(1,n-2))/c(3,n-1))/g g = -g/c(4,n-1) c(4,n) = c(3,n-1) c(4,n) = g*c(3,n-1)+c(4,n) c(2,n) = (g*c(2,n-1)+c(2,n))/c(4,n) go to 160 end if ! ! Either (N = 3 and not-a-knot also at left) or (N=2 and not not-a- ! knot at left end point). ! 100 continue c(2,n) = 2.0D+00*c(4,n) c(4,n) = 1.0D+00 g = -1.0D+00/c(4,n-1) c(4,n) = g*c(3,n-1)+c(4,n) c(2,n) = (g*c(2,n-1)+c(2,n))/c(4,n) go to 160 ! ! Second derivative prescribed at right endpoint. ! 110 continue c(2,n) = 3.0D+00*c(4,n)+c(3,n)/2.0D+00*c(2,n) c(4,n) = 2.0D+00 g = -1.0D+00/c(4,n-1) c(4,n) = g*c(3,n-1)+c(4,n) c(2,n) = (g*c(2,n-1)+c(2,n))/c(4,n) go to 160 120 continue if ( ibcend/= 1 ) then if ( ibcend>1 ) then c(2,n) = 3.0D+00*c(4,n)+c(3,n)/2.0D+00*c(2,n) c(4,n) = 2.0D+00 g = -1.0D+00/c(4,n-1) c(4,n) = g*c(3,n-1)+c(4,n) c(2,n) = (g*c(2,n-1)+c(2,n))/c(4,n) go to 160 end if if ( ibcbeg>0)go to 100 ! ! Not-a-knot at right endpoint and at left endpoint and N = 2. ! c(2,n) = c(4,n) end if ! ! Carry out the back substitution ! 160 continue do j = n-1,1,-1 c(2,j) = (c(2,j)-c(3,j)*c(2,j+1))/c(4,j) end do ! ! Generate cubic coefficients in each interval, that is, the ! derivatives at its left endpoint, from value and slope at its ! endpoints. ! do i = 2,n dtau = c(3,i) divdf1 = (c(1,i)-c(1,i-1))/dtau divdf3 = c(2,i-1)+c(2,i)-2.0D+00*divdf1 c(3,i-1) = 2.0D+00*(divdf1-c(2,i-1)-divdf3)/dtau c(4,i-1) = (divdf3/dtau)*(6.0D+00/dtau) end do return end subroutine disc_cost ( costp, costu, costv, g, gtar, indx, neqn, np, nprof, & ny, yc ) !*****************************************************************************80 ! !! DISC_COST computes the discrepancy cost integrals. ! ! Discussion: ! ! The discrepancy cost integrals measure the discrepancies between the ! target and computed pressure, horizontal and vertical velocities along ! the profile line. ! ! This integration scheme assumes that the profile line, and ! the element sides that define it, are straight. Otherwise, ! the integration scheme used is not correct. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none integer, parameter :: nquad1 = 5 integer neqn integer np integer ny real ( kind = 8 ) bval real ( kind = 8 ) costp real ( kind = 8 ) costu real ( kind = 8 ) costv real ( kind = 8 ) g(neqn) real ( kind = 8 ) gtar(neqn) integer i integer ii integer indx(np,3) integer j integer k integer npol integer nprof(2*ny-1) real ( kind = 8 ) pcof(2) real ( kind = 8 ) pval real ( kind = 8 ) ucof(3) real ( kind = 8 ) uval real ( kind = 8 ) vcof(3) real ( kind = 8 ) vval real ( kind = 8 ) wquad1(nquad1) real ( kind = 8 ) xquad1(nquad1) real ( kind = 8 ) yc(np) real ( kind = 8 ) yhi real ( kind = 8 ) ylo real ( kind = 8 ) ypol(3) real ( kind = 8 ) yval ! ! Get the weights and abscissas to approximate a line integral. ! call gquad1 ( nquad1, wquad1, xquad1 ) ! ! Compute the integral of the difference squared between the ! current velocity and the target values. ! costu = 0.0D+00 costv = 0.0D+00 ! ! The line along which we integrate is broken into NY-1 ! subintervals, over each of which, U and V are represented ! by quadratic functions. ! do i = 1, ny-1 ! ! Get the values of U and V at the beginning, middle, and ! end of the subinterval. Use these to compute the quadratic ! representation of U and V for any point on the subinterval. ! ylo = yc(nprof(2*i-1)) yhi = yc(nprof(2*i+1)) npol = 3 do k = 1, npol ii = 2 * i - 2 + k ypol(k) = yc(nprof(ii)) j = indx(nprof(ii),1) ucof(k) = g(j) - gtar(j) j = indx(nprof(ii),2) vcof(k) = g(j) - gtar(j) end do ! ! Evaluate the discrepancy at each quadrature point. ! do j = 1, nquad1 call rint_to_rint ( -1.0D+00, +1.0D+00, xquad1(j), ylo, yhi, yval ) uval = 0.0D+00 vval = 0.0D+00 do k = 1, npol call lbase ( k, npol, bval, ypol, yval ) uval = uval + bval * ucof(k) vval = vval + bval * vcof(k) end do costu = costu + 0.5D+00 * wquad1(j) * ( yhi - ylo ) * uval**2 costv = costv + 0.5D+00 * wquad1(j) * ( yhi - ylo ) * vval**2 end do end do ! ! Compute the square root of the integral of the difference ! squared between the current pressure and the target values. ! costp = 0.0D+00 do i = 1, ny-1 ylo = yc(nprof(2*i-1)) yhi = yc(nprof(2*i+1)) npol = 2 do k = 1, npol ii = 2*i-3+2*k ypol(k) = yc(nprof(ii)) j = indx(nprof(ii),3) if ( j <= 0 ) then pcof(k) = 0.0D+00 else pcof(k) = g(j)-gtar(j) end if end do do j = 1, nquad1 call rint_to_rint ( -1.0D+00, +1.0D+00, xquad1(j), ylo, yhi, yval ) pval = 0.0D+00 do k = 1, npol call lbase ( k, npol, bval, ypol, yval ) pval = pval + bval * pcof(k) end do costp = costp + 0.5D+00 * wquad1(j) * ( yhi - ylo ) * pval**2 end do end do return end subroutine disc_der ( dpara, g, gtar, indx, neqn, np, npar, nprof, ny, sens, & watep, wateu, watev, yc ) !*****************************************************************************80 ! !! DISC_DER computes the derivative of the discrepancy cost integrals ! ! Discussion: ! ! The discrepancy integrals for pressure, horizontal and vertical ! velocities are to be differentiated with respect to the ! various parameters. ! ! In practice, the dummy parameter DPARA will actually be DPARA1 or ! DPARA2, and the dummy input parameter SENS will actually be SENS or ! GDIF, depending on which estimate of the cost derivatives is ! desired. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real ( kind = 8 ) DPARA(NPAR), contains the derivatives of ! the cost function with respect to the various parameters. In particular, ! DPARA(I) = D cost / D parameter(I). ! ! Input, real ( kind = 8 ) G(NEQN). ! G is the current solution vector, in which are stored ! pressures and velocities. ! ! Input, real ( kind = 8 ) GTAR(NEQN), ! the target solution vector. ! ! Input, integer INDX(NP,3). ! INDX contains, for each node I, the index of U, V and P at ! that node, or 0 or a negative value. ! If K=INDX(I,J) is positive, then the value of the degree ! of freedom is stored in the solution vector entry G(K). ! If INDX(I,J) is positive, then that means that a degree of ! freedom for variable J (U, V or P) is associated with node ! I, and an equation will be generated to determine its value. ! If INDX(I,J) is zero, then that means the the value of variable ! J (U, V or P) has been specified at node I. No equation is ! generated to determine its value. ! ! Input, integer NEQN, the number of equations, and functions. ! ! Input, integer NP, the number of nodes. NP=(2*NX-1)*(2*NY-1). ! ! Input, integer NPAR. ! The number of parameters. NPAR = NPARF + NPARB + 1. ! The parameters control the shape of the inflow, ! the shape of the bump obstacle, and the strength of the ! flow. ! ! Input, integer NPROF(2*MAXNY-1), the indices of the nodes along ! the profile line. ! ! Input, integer NY. ! NY controls the spacing of nodes and elements in ! the Y direction. There are 2*NY-1 nodes along various ! lines in the Y direction. ! Roughly speaking, NY is the number of elements along ! a line in the Y direction. ! ! Input, real ( kind = 8 ) SENS(MAXEQN,NPAR), the sensitivities. ! SENS(I,J) contains the sensitivity of the I-th unknown ! with respect to the J-th parameter. ! ! Input, real ( kind = 8 ) WATEP, WATEU, WATEV. ! These are weights used in computing the overall cost ! function based on the costs of the flow discrepancy. ! ! Input, real ( kind = 8 ) YC(NP), the Y coordinates of the nodes ! implicit none integer mpar integer nquad1 parameter (mpar = 10) parameter (nquad1 = 5) integer neqn integer np integer npar integer ny real ( kind = 8 ) bval real ( kind = 8 ) dpara(npar) real ( kind = 8 ) dpcof(mpar,mpar) real ( kind = 8 ) dpval(mpar) real ( kind = 8 ) ducof(mpar,mpar) real ( kind = 8 ) duval(mpar) real ( kind = 8 ) dvcof(mpar,mpar) real ( kind = 8 ) dvval(mpar) real ( kind = 8 ) g(neqn) real ( kind = 8 ) gtar(neqn) integer i integer ii integer indx(np,3) integer j integer k integer l integer npol integer nprof(2*ny-1) real ( kind = 8 ) pcof(2) real ( kind = 8 ) pval real ( kind = 8 ) sens(neqn,npar) real ( kind = 8 ) ucof(3) real ( kind = 8 ) uval real ( kind = 8 ) vcof(3) real ( kind = 8 ) vval real ( kind = 8 ) watep real ( kind = 8 ) wateu real ( kind = 8 ) watev real ( kind = 8 ) wquad1(nquad1) real ( kind = 8 ) xquad1(nquad1) real ( kind = 8 ) yc(np) real ( kind = 8 ) yhi real ( kind = 8 ) ylo real ( kind = 8 ) ypol(3) real ( kind = 8 ) yval if ( npar > mpar ) then write ( *, * ) ' ' write ( *, * ) 'DISC_DER - Fatal error!' write ( *, * ) ' The number of parameters is too high.' write ( *, * ) ' This routine can handle NPAR = ',mpar write ( *, * ) ' Your problem has NPAR = ',npar stop end if ! ! Get the Gauss weights and abscissas to approximate a line integral. ! call gquad1 ( nquad1, wquad1, xquad1 ) ! ! The line along which we integrate is broken into NY-1 ! subintervals, over each of which, U and V are represented ! by quadratic functions. ! do i = 1, ny-1 ! ! Get the values of U and V at the beginning, middle, and ! end of the subinterval. Use these to compute the quadratic ! representation of U and V for any point on the subinterval. ! ylo = yc(nprof(2*i-1)) yhi = yc(nprof(2*i+1)) npol = 3 do k = 1,npol ii = 2*i-2+k ypol(k) = yc(nprof(ii)) j = indx(nprof(ii),1) ucof(k) = g(j)-gtar(j) do l = 1,npar ducof(k,l) = sens(j,l) end do j = indx(nprof(ii),2) vcof(k) = g(j)-gtar(j) do l = 1,npar dvcof(k,l) = sens(j,l) end do end do ! ! Evaluate the discrepancy at each quadrature point. ! do j = 1,nquad1 call rint_to_rint ( -1.0D+00, +1.0D+00, xquad1(j), ylo, yhi, yval ) uval = 0.0D+00 vval = 0.0D+00 do l = 1,npar duval(l) = 0.0D+00 dvval(l) = 0.0D+00 end do do k = 1,npol call lbase(k,npol,bval,ypol,yval) uval = uval+bval*ucof(k) vval = vval+bval*vcof(k) do l = 1,npar duval(l) = duval(l)+bval*ducof(k,l) dvval(l) = dvval(l)+bval*dvcof(k,l) end do end do do l = 1,npar dpara(l) = dpara(l)+0.5D+00*wquad1(j)*(yhi-ylo) & *2.0D+00*(wateu*uval*duval(l)+watev*vval*dvval(l)) end do end do end do ! ! Compute the square root of the integral of the difference ! squared between the current pressure and the target values. ! do i = 1,ny-1 ylo = yc(nprof(2*i-1)) yhi = yc(nprof(2*i+1)) npol = 2 do k = 1,npol ii = 2*i-3+2*k ypol(k) = yc(nprof(ii)) j = indx(nprof(ii),3) if ( j <= 0 ) then pcof(k) = 0.0D+00 else pcof(k) = g(j)-gtar(j) end if do l = 1,npar if ( j <= 0 ) then dpcof(k,l) = 0.0D+00 else dpcof(k,l) = sens(j,l) end if end do end do do j = 1,nquad1 call rint_to_rint ( -1.0D+00, +1.0D+00, xquad1(j), ylo, yhi, yval ) pval = 0.0D+00 do l = 1,npar dpval(l) = 0.0D+00 end do do k = 1,npol call lbase(k,npol,bval,ypol,yval) pval = pval+bval*pcof(k) do l = 1,npar dpval(l) = dpval(l)+bval*dpcof(k,l) end do end do do l = 1,npar dpara(l) = dpara(l) + 0.5D+00 * wquad1(j) * ( yhi - ylo ) & * 2.0D+00 * watep * pval * dpval(l) end do end do end do return end subroutine dmemry(action,name,dval) !*****************************************************************************80 ! !! DMEMRY allows the user to define the name of a real ! variable, set it, increment it, or get the value. ! ! ACTION Input, Character*(*) ACTION, desired action. ! ! 'Init', reset all values to zero, wipe out all names. ! 'Name', add a variable of the given name. ! 'Inc', increment variable NAME by DVAL. ! 'Set', set variable NAME to DVAL. ! 'Get', return value of NAME in DVAL. ! 'Zero', reset all values to zero. ! ! NAME Input, Character*(*) NAME, the name of the variable. ! ! DVAL Input/output, real ( kind = 8 ) DVAL. ! ! For the 'Inc' and 'Set' commands, DVAL must contain the ! increment or set value. ! ! For the 'Get' command, DVAL will contain the value of the ! named variable on output. ! implicit none integer maxnam parameter (maxnam = 100) character ( len = * ) action real ( kind = 8 ) dval real ( kind = 8 ) dvals(maxnam) integer i integer lenc logical s_eqi character ( len = * ) name character ( len = 20 ) names(maxnam) integer numnam save dvals save names save numnam ! ! Initialize everything. ! if ( s_eqi(action,'init') ) then numnam = 0 do i = 1,maxnam dvals(i) = 0.0D+00 names(i) = ' ' end do ! ! Name something. ! else if ( s_eqi(action,'name') ) then do i = 1,numnam if ( s_eqi(name,names(i)) ) then write ( *, * ) ' ' write ( *, * ) 'DMemry - Warning!' write(*,'('' There is ALREADY a variable '',a)')name return end if end do if ( numnam0 ) then if ( eqn(iprs) == 'P' ) then iuse = iprs-jhor+2*nlband+1 term = ar*dwjdx*qi a(iuse,jhor) = a(iuse,jhor)+term end if end if ! ! Contributions of the JVER vertical velocity variable to the ! U, V and P equations. ! if ( eqn(ihor) == 'U' ) then if ( s_eqi(syseqn,'NavierStokes') ) then iuse = ihor-jver+2*nlband+1 term = ar*nu_inv*wj*dudy*wi a(iuse,jver) = a(iuse,jver)+term end if end if iuse = iver-jver+2*nlband+1 if ( eqn(iver) == 'V' ) then if ( s_eqi(syseqn,'NavierStokes') ) then term = ar*(dwjdx*dwidx+dwjdy*dwidy & +nu_inv*(uold*dwjdx+wj*dvdy+vold*dwjdy)*wi) else if ( s_eqi(syseqn,'Stokes') ) then term = ar*(dwjdx*dwidx+dwjdy*dwidy) end if a(iuse,jver) = a(iuse,jver)+term end if if ( iprs>0 ) then if ( eqn(iprs) == 'P' ) then iuse = iprs-jver+2*nlband+1 term = ar*dwjdy*qi a(iuse,jver) = a(iuse,jver)+term end if end if ! ! Contributions of the JPRS pressure to the U and V equations. ! if ( jprs>0 ) then if ( eqn(ihor) == 'U' ) then iuse = ihor-jprs+2*nlband+1 term = ar*nu_inv*dqjdx*wi a(iuse,jprs) = a(iuse,jprs)+term end if if ( eqn(iver) == 'V' ) then iuse = iver-jprs+2*nlband+1 term = ar*nu_inv*dqjdy*wi a(iuse,jprs) = a(iuse,jprs)+term end if end if end do end do end do end do ! ! Set up the equations that enforce boundary conditions. ! do ip = 1,np ihor = indx(ip,1) iver = indx(ip,2) iprs = indx(ip,3) if ( eqn(ihor) == 'UB' .or. eqn(ihor)=='UI'.or.eqn(ihor)=='UW' ) then a(2*nlband+1,ihor) = 1.0D+00 end if if ( eqn(iver) == 'VB' .or. eqn(iver)=='VI'.or.eqn(iver)=='VW' ) then a(2*nlband+1,iver) = 1.0D+00 end if if ( iprs>0 ) then if ( eqn(iprs) == 'PB' ) then a(2*nlband+1,iprs) = 1.0D+00 end if end if end do return end subroutine fprnab ( anab,area,eqn,g,indx,nelem,neqn,node,np,npar,nparb,nparf, & npe,ny,para,phi,syseqn) !*****************************************************************************80 ! !! FPRNAB computes the jacobian and stores it in a neighbor matrix. ! ! Discussion: ! ! This routine computes the jacobian of the Navier Stokes or Stokes ! residual functions and stores it in the neighbor matrix ANAB. ! ! The differentiated Navier Stokes functions have the form: ! ! ! d U-Eqn/d U-Coef: ! ! Integral ! ! dWj/dx * dWi/dx + dWj/dy * dWi/dy ! + nu_inv * (Wj*dUold/dx + Uold*dWj/dx+ Vold*dWj/dy) * Wi dx dy ! ! d U-Eqn/d V-Coef: ! ! Integral ! ! nu_inv * Wj*dUold/dy * Wi dx dy ! ! d U-Eqn/d P-Coef: ! ! Integral ! ! nu_inv * dQj/dx * Wi dx dy ! ! d V-Eqn/d U-Coef: ! ! Integral ! ! nu_inv * Wj*dVold/dx * Wi dx dy ! ! d V-Eqn/d V-Coef: ! ! Integral ! ! dWj/dx * dWi/dx + dWj/dy * dWi/dy ! + nu_inv * (Uold*dWj/dx + Wj*dVold/dy + Vold*dWj/dy) * Wi dx dy ! ! d V-Eqn/d P-Coef: ! ! Integral ! ! nu_inv * dQj/dy * Wi dx dy ! ! d P-Eqn/d U-Coef: ! ! Integral ! ! dWj/dx * Qi dx dy ! ! Integral ! ! dWj/dy * Qi dx dy ! ! ! The differentiated Stokes functions have the form: ! ! ! d U-Eqn/d U-Coef: ! ! Integral ! ! dWj/dx * dWi/dx + dWj/dy * dWi/dy dx dy ! ! d U-Eqn/d P-Coef: ! ! Integral ! ! nu_inv * dQj/dx * Wi dx dy ! ! d V-Eqn/d V-Coef: ! ! Integral ! ! dWj/dx * dWi/dx + dWj/dy * dWi/dy dx dy ! ! d V-Eqn/d P-Coef: ! ! Integral ! ! nu_inv * dQj/dy * Wi dx dy ! ! d P-Eqn/d U-Coef: ! ! Integral ! ! dWj/dx * Qi dx dy ! ! Integral ! ! dWj/dy * Qi dx dy ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) ANAB(3,3,19,NP), contains the value of D F(I)/D X(J) ! for each of the NEQN residual functions F(I) with respect to each ! of the unknown coefficients X(J), stored in a neighbor format. ! implicit none integer nelem integer neqn integer np integer npar integer nparf integer npe real ( kind = 8 ) anab(3,3,19,np) real ( kind = 8 ) ar real ( kind = 8 ) area(nelem,3) real ( kind = 8 ) dpdx real ( kind = 8 ) dpdy real ( kind = 8 ) dqjdx real ( kind = 8 ) dqjdy real ( kind = 8 ) dudx real ( kind = 8 ) dudy real ( kind = 8 ) dvdx real ( kind = 8 ) dvdy real ( kind = 8 ) dwidx real ( kind = 8 ) dwidy real ( kind = 8 ) dwjdx real ( kind = 8 ) dwjdy character ( len = 2 ) eqn(neqn) real ( kind = 8 ) g(neqn) integer i integer ielem integer ihor integer indx(np,3) integer ip integer iprs integer iq integer iquad integer iuse integer ival integer iver integer j integer jcol integer jp integer jq integer jrow integer k integer l logical s_eqi integer my integer node(nelem,npe) integer nparb integer ny real ( kind = 8 ) pold real ( kind = 8 ) para(npar) real ( kind = 8 ) phi(nelem,3,6,10) real ( kind = 8 ) qi real ( kind = 8 ) nu_inv character ( len = * ) syseqn real ( kind = 8 ) term real ( kind = 8 ) uold real ( kind = 8 ) vold real ( kind = 8 ) wi real ( kind = 8 ) wj external s_eqi ival = 1 call imemry('inc','Fprime_calls',ival) nu_inv = para(nparf+nparb+1) anab(1:3,1:3,1:19,1:np) = 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 and V. ! call uvalq ( dpdx,dpdy,dudx,dudy,dvdx,dvdy,g,ielem,indx, & iquad,nelem,neqn,node,np,npe,pold,phi,uold,vold) ! ! 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) iver = indx(ip,2) iprs = indx(ip,3) ! ! Now compute the derivatives of the functions associated ! with U, V and P, with respect to the coefficients associated ! with basis vectors at each node of the element. ! do jq = 1,6 jp = node(ielem,jq) jcol = ((jp-1)/(2*ny-1))+1 jrow = jp-(jcol-1)*(2*ny-1) wj = phi(ielem,iquad,jq,1) dwjdx = phi(ielem,iquad,jq,2) dwjdy = phi(ielem,iquad,jq,3) dqjdx = phi(ielem,iquad,jq,5) dqjdy = phi(ielem,iquad,jq,6) call nodnab(ip,jp,my,iuse) ! ! Contributions of the JHOR variable to various equations. ! if ( eqn(ihor) == 'U' ) then if ( s_eqi(syseqn,'NavierStokes') ) then term = ar*(dwjdx*dwidx+dwjdy*dwidy+ & nu_inv*(wj*dudx+uold*dwjdx+vold*dwjdy)*wi) else if ( s_eqi(syseqn,'Stokes') ) then term = ar*(dwjdx*dwidx+dwjdy*dwidy) end if anab(1,1,iuse,ip) = anab(1,1,iuse,ip)+term end if if ( eqn(iver) == 'V' ) then if ( s_eqi(syseqn,'NavierStokes') ) then term = ar*(nu_inv*wj*dvdx*wi) anab(2,1,iuse,ip) = anab(2,1,iuse,ip)+term end if end if if ( iprs>0 ) then if ( eqn(iprs) == 'P' ) then term = ar*dwjdx*qi anab(3,1,iuse,ip) = anab(3,1,iuse,ip)+term end if end if ! ! Contributions of the JVER variable to various equations. ! if ( eqn(ihor) == 'U' ) then if ( s_eqi(syseqn,'NavierStokes') ) then term = ar*nu_inv*wj*dudy*wi anab(1,2,iuse,ip) = anab(1,2,iuse,ip)+term end if end if if ( eqn(iver) == 'V' ) then if ( s_eqi(syseqn,'NavierStokes') ) then term = ar*(dwjdx*dwidx+dwjdy*dwidy & +nu_inv*(uold*dwjdx+wj*dvdy+vold*dwjdy)*wi) else if ( s_eqi(syseqn,'Stokes') ) then term = ar*(dwjdx*dwidx+dwjdy*dwidy) end if anab(2,2,iuse,ip) = anab(2,2,iuse,ip)+term end if if ( iprs>0 ) then if ( eqn(iprs) == 'P' ) then term = ar*dwjdy*qi anab(3,2,iuse,ip) = anab(3,2,iuse,ip)+term end if end if ! ! Contributions of the JPRS variable to various equations. ! if ( mod(jrow,2) == 1 .and. mod(jcol,2)==1 ) then if ( eqn(ihor) == 'U' ) then term = ar*nu_inv*dqjdx*wi anab(1,3,iuse,ip) = anab(1,3,iuse,ip)+term end if if ( eqn(iver) == 'V' ) then term = ar*nu_inv*dqjdy*wi anab(2,3,iuse,ip) = anab(2,3,iuse,ip)+term end if end if end do end do end do end do ! ! Set up the equations that enforce boundary conditions. ! do ip = 1,np ihor = indx(ip,1) iver = indx(ip,2) iprs = indx(ip,3) if ( eqn(ihor) == 'UB' .or. eqn(ihor)=='UI'.or.eqn(ihor)=='UW' ) then anab(1,1,10,ip) = 1.0D+00 end if if ( eqn(iver) == 'VB' .or. eqn(iver)=='VI'.or.eqn(iver)=='VW' ) then anab(2,2,10,ip) = 1.0D+00 end if if ( iprs>0 ) then if ( eqn(iprs) == 'PB' ) then anab(3,3,10,ip) = 1.0D+00 end if end if end do return end subroutine fx ( area,eqn,g,indx,ishapf,nelem,neqn,node,np,npar,nparb,nparf, & npe,para,phi,res,splflo,syseqn,tauflo,yc) !*****************************************************************************80 ! !! FX computes the residual of the Navier Stokes or Stokes equations. ! ! Discussion: ! ! The Navier Stokes equations have the form: ! ! Integral ! ! dU/dx * dW/dx + dU/dy * dW/dy ! + nu_inv * (U*dU/dx + V*dU/dy + dP/dx) * W dx dy = 0 ! ! Integral ! ! dV/dx * dW/dx + dV/dy * dW/dy ! + nu_inv * (U*dV/dx + V*dV/dy + dP/dy) * W dx dy = 0 ! ! Integral ! ! (dU/dx + dV/dy) * Q dx dy = 0 ! ! The Stokes equations have the form: ! ! Integral ! ! dU/dx * dW/dx + dU/dy * dW/dy + nu_inv * dP/dx * W dx dy = 0 ! ! Integral ! ! dV/dx * dW/dx + dV/dy * dW/dy + nu_inv * dP/dy * W dx dy = 0 ! ! Integral ! ! (dU/dx + dV/dy) * Q dx dy = 0 ! ! Here W is a basis function for U and V, and Q is a basis ! function for P. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) RES(NEQN), contains the value ! of the residual. ! 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 ) dudx real ( kind = 8 ) dudy real ( kind = 8 ) dvdx real ( kind = 8 ) dvdy real ( kind = 8 ) dwidx real ( kind = 8 ) dwidy character ( len = 2 ) eqn(neqn) real ( kind = 8 ) g(neqn) integer i integer ielem integer ihor integer indx(np,3) integer ip integer iprs integer iq integer iquad integer ishapf integer ival integer iver logical s_eqi 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 ) qi real ( kind = 8 ) res(neqn) real ( kind = 8 ) nu_inv real ( kind = 8 ) splflo(4,nparf+2,0:nparf) character ( len = 20 ) syseqn real ( kind = 8 ) tauflo(nparf+2) real ( kind = 8 ) u real ( kind = 8 ) ubc real ( kind = 8 ) v real ( kind = 8 ) vbc real ( kind = 8 ) wi real ( kind = 8 ) yc(np) external s_eqi ival = 1 call imemry('inc','Fx_calls',ival) nu_inv = para(nparf+nparb+1) res(1:neqn) = 0.0D+00 ! ! Consider an element. ! do ielem = 1,nelem ! ! Evaluate the integrand at the quadrature points. ! do iquad = 1,3 ar = area(ielem,iquad) ! ! Evaluate P, U and V. ! call uvalq ( dpdx,dpdy,dudx,dudy,dvdx,dvdy,g,ielem,indx, & iquad,nelem,neqn,node,np,npe,p,phi,u,v) ! ! Look at nearby basis functions. ! 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) ! ! The horizontal velocity equations. ! ihor = indx(ip,1) if ( eqn(ihor) == 'U' ) then if ( s_eqi(syseqn,'NavierStokes') ) then res(ihor) = res(ihor)+ar*(dudx*dwidx + dudy*dwidy & +nu_inv*(u*dudx+v*dudy+dpdx)*wi ) else if ( s_eqi(syseqn,'Stokes') ) then res(ihor) = res(ihor)+ar*(dudx*dwidx + dudy*dwidy & +nu_inv*dpdx*wi ) end if else if ( eqn(ihor) == 'UB' ) then res(ihor) = g(ihor) else if ( eqn(ihor) == 'UI' ) then call flo_spl_val ( ishapf, nparf, splflo, tauflo, ubc, vbc, yc(ip) ) res(ihor) = g(ihor)-ubc else if ( eqn(ihor) == 'UW' ) then res(ihor) = g(ihor) end if ! ! The vertical velocity equations. ! iver = indx(ip,2) if ( eqn(iver) == 'V' ) then if ( s_eqi(syseqn,'NavierStokes') ) then res(iver) = res(iver)+ar*(dvdx*dwidx + dvdy*dwidy & +nu_inv*(u*dvdx+v*dvdy+dpdy)*wi ) else if ( s_eqi(syseqn,'Stokes') ) then res(iver) = res(iver)+ar*(dvdx*dwidx + dvdy*dwidy & +nu_inv*dpdy*wi ) end if else if ( eqn(iver) == 'VB' ) then res(iver) = g(iver) else if ( eqn(iver) == 'VI' ) then call flo_spl_val ( ishapf, nparf, splflo, tauflo, ubc, vbc, yc(ip) ) res(iver) = g(iver)-vbc else if ( eqn(iver) == 'VW' ) then res(iver) = g(iver) end if ! ! The pressure equations. ! iprs = indx(ip,3) if ( iprs>0 ) then if ( eqn(iprs) == 'P' ) then res(iprs) = res(iprs)+ar*(dudx+dvdy)*qi else if ( eqn(iprs) == 'PB' ) then res(iprs) = g(iprs) end if end if end do end do end do return end subroutine get_cost ( cost, costb, costp, costu, costv, g, gtar, indx, & ishapb, neqn, np, nparb, nprof, ny, splbmp, taubmp, wateb, watep, wateu, & watev, xbl, xbr, ybl, ybr, yc ) !*****************************************************************************80 ! !! GET_COST determines the cost of a solution G given a target solution GTAR. ! ! Modified: ! ! 31 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none integer neqn integer np integer nparb integer ny real ( kind = 8 ) cost real ( kind = 8 ) costb real ( kind = 8 ) costp real ( kind = 8 ) costu real ( kind = 8 ) costv real ( kind = 8 ) g(neqn) real ( kind = 8 ) gtar(neqn) integer indx(np,3) integer ishapb integer ival integer nprof(2*ny-1) real ( kind = 8 ) splbmp(4,nparb+2,0:nparb) real ( kind = 8 ) taubmp(nparb+2) real ( kind = 8 ) wateb real ( kind = 8 ) watep real ( kind = 8 ) wateu real ( kind = 8 ) watev real ( kind = 8 ) xbl real ( kind = 8 ) xbr real ( kind = 8 ) ybl real ( kind = 8 ) ybr real ( kind = 8 ) yc(np) ival = 1 call imemry ( 'inc', 'GET_COST_calls', ival ) call bump_cost ( costb, ishapb, nparb, splbmp, taubmp, xbl, xbr, ybl, ybr ) call disc_cost ( costp, costu, costv, g, gtar, indx, neqn, np, nprof, ny, yc ) cost = wateb * costb + watep * costp + wateu * costu + watev * costv return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5 and 6). ! ! Otherwise, IUNIT is an integer between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! implicit none integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine getdu ( dpdyn,dudyn,dvdyn,etan,g,indx,isotri,nelem,neqn,node,np, & npe,numel,xc,xsin,yc) !*****************************************************************************80 ! !! GETDU estimates spatial derivatives of state variables at nodes. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none integer nelem integer neqn integer np integer npe real ( kind = 8 ) det real ( kind = 8 ) detadx real ( kind = 8 ) detady real ( kind = 8 ) dpdx real ( kind = 8 ) dpdy real ( kind = 8 ) dpdyn(np) real ( kind = 8 ) dudx real ( kind = 8 ) dudy real ( kind = 8 ) dudyn(np) real ( kind = 8 ) dvdx real ( kind = 8 ) dvdy real ( kind = 8 ) dvdyn(np) real ( kind = 8 ) dxsidx real ( kind = 8 ) dxsidy real ( kind = 8 ) eta real ( kind = 8 ) etan(6) real ( kind = 8 ) g(neqn) integer ielem integer indx(np,3) integer ip integer iq integer isotri(nelem) integer node(nelem,npe) integer numel(np) real ( kind = 8 ) p real ( kind = 8 ) u real ( kind = 8 ) v real ( kind = 8 ) xc(np) real ( kind = 8 ) xq real ( kind = 8 ) xsi real ( kind = 8 ) xsin(6) real ( kind = 8 ) yc(np) real ( kind = 8 ) yq ! ! Estimate dUdY, dVdY, dPdY at each node. ! ! We have a problem, because these quantities are not continuously ! defined in all directions at nodes, which lie on the border between ! two or more elements. ! ! In order to assign some reasonable value to these quantities, ! we look at each node, count the number of elements in which it ! lies, evaluate the quantity within each element, and average the ! result. ! numel(1:np) = 0 dpdyn(1:np) = 0.0D+00 dudyn(1:np) = 0.0D+00 dvdyn(1:np) = 0.0D+00 do ielem = 1,nelem do iq = 1,6 ip = node(ielem,iq) numel(ip) = numel(ip)+1 xq = xc(ip) yq = yc(ip) eta = etan(iq) xsi = xsin(iq) if ( isotri(ielem) == 2 ) then call trans ( det,detadx,detady,dxsidx,dxsidy,eta,ielem,nelem,node,np, & npe,xc,xsi,yc) end if call uval ( detadx,detady,dpdx,dpdy,dudx,dudy,dvdx,dvdy, & dxsidx,dxsidy,eta,g,ielem,indx,isotri,nelem,neqn, & node,np,npe,p,u,v,xc,xq,xsi,yc,yq) dpdyn(ip) = dpdyn(ip) + dpdy dudyn(ip) = dudyn(ip) + dudy dvdyn(ip) = dvdyn(ip) + dvdy end do end do ! ! Take the average value of the quantities over all the ! different elements along whose boundaries they are defined. ! do ip = 1, np dpdyn(ip) = dpdyn(ip) / real ( numel(ip), kind = 8 ) dudyn(ip) = dudyn(ip) / real ( numel(ip), kind = 8 ) dvdyn(ip) = dvdyn(ip) / real ( numel(ip), kind = 8 ) end do return end subroutine getdu4(dudyn,np,nx,ny) !*****************************************************************************80 ! !! GETDU4 uses the Zienkiewicz-Zhou technique to attempt to improve the ! accuracy of the computed values of dPdY, dUdY and dVdY at corner ! nodes. ! ! It differs from GETDU2 in that it also modifies the values of ! dUdY at the midside nodes. ! ! Here is a picture of the element patch around an interior node "C": ! ! .----.----NN---NNE--NNEE ! | /| /| ! | / | / | ! | / | / | ! | / | / | ! . NW N NE NEE ! | / | / | ! | / | / | ! | / | / | ! |/ |/ | ! WW---W----C----E----EE ! | /| /| ! | / | / | ! | / | / | ! | / | / | ! SWW SW S SE . ! | / | / | ! | / | / | ! | / | / | ! |/ |/ | ! SSWW-SSW--SS---.----. ! ! The midside nodes SWW, SSW, SW, W, NW, S, N, SE, E, NE, NNE, NEE ! are all quadrature points, at which a higher rater of convergence ! is expected in the solution. The nodes labeled SSWW, WW, SS, C, ! NN, EE and NNEE are corner nodes, whose values of dUdY we wish ! to improve. On this step, we will try to compute improved values ! of the node "C", and the side nodes. The step will be repeated ! so that every pressure node gets to be the "C" node of a patch ! (except for corner nodes). Nodes labeled "." are not included in ! the patch, and are shown only for completeness. ! ! The value of the quantity dUdY at each of the midside nodes is ! taken as data, to be fitted by a least squares polynomial based ! on the six basis functions (1, x, y, x*x, x*y, y*y). For the ! case where the node C is in the interior of the region, this ! makes 12 equations in 6 unknowns. ! ! The above diagram is assigned "logical" coordinates that ! range from 0 <= x,y <= 4, with SWW having coordinates (0,1) ! and C having coordinates (2,2). (These are really a sort of ! XSI, ETA coordinate system, rather than a "physical" X, Y ! system.) ! ! Using standard discrete least squares methods, we write our ! original set of equations as ! ! AC x = b ! ! where x is the set of 6 unknown polynomial coefficients, ! row i of matrix AC is the value of (1, x, y, x*x, x*y, y*y) ! at the i-th of our 12 quadrature points. ! ! To get a solvable system, we multiply by the transpose of AC, ! getting the square system known as the normal equations: ! ! ACT AC x = ACT b ! ! The matrix ACT AC only has to be set up and factored once, and ! we can use the factors over and over for all the patches ! associated with interior nodes. (If we really like this method, ! we should stop using the normal equations and use instead ! a QR routine like the LINPACK SQRDC/SQRSL). ! ! Once the coefficients of the polynomial are found, the least ! squares polynomial may be evaluated at any point in the ! "patch" of six elements. In this routine, we choose to ! evaluate the polynomial only at the point C, and use this ! value in place of the value computed by the finite element ! method. The values at the midside nodes are unchanged. ! ! Special treatment must be made for nodes which lie on the ! north, east, west or south sides, or in the corners. ! ! In the case of side nodes, we simply adapt the process to ! account for 7 data values, using the six basis functions. ! ! In the case of a corner node, we simply apply the relevant ! value computed for the patch along the north or south boundary ! nearest to the node. ! ! In particular, separate matrices AN, AE, AW and AS must be set ! up to represent the least squares linear systems that occur for ! nodes which lie along the north, east, west or south boundaries. ! ! ! Since the midside nodes can occur in more than one element, ! we ADD the improved estimate to a running total, and average ! later. ! ! Note: this routine explicitly assumes that the nodes are ordered ! from top to bottom, then from left to right, as in this example: ! ! 5 10 15 ! 4 9 14 ! 3 8 13 ! 2 7 12 ! 1 6 11 ! ! It also assumes that EVERY node has an associated value of DUDYN, ! and that the pressure nodes, at which we will smooth the value of ! DUDYN, are exactly those nodes in an odd row and odd column. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, REAL DUDYN(NP). ! On input, DUDYN contains the original value of DUDY at every node. ! On output, DUDYN contains the smoothed values of dUdY. ! ! Input, INTEGER NP, the number of nodes. ! ! Input, INTEGER NX, the number of columns of pressure nodes ! in the X direction. The actual number of columns of nodes ! is 2*NX-1. ! ! Input, INTEGER NY, the number of rows of pressure nodes ! in the Y direction. The actual number of rows of nodes ! is 2*NY-1. ! implicit none integer, parameter :: ncoef = 6 integer, parameter :: npmax = 7889 integer np real ( kind = 8 ) ac(12,ncoef) real ( kind = 8 ) ae(7,ncoef) real ( kind = 8 ) an(7,ncoef) real ( kind = 8 ) as(7,ncoef) real ( kind = 8 ) aw(7,ncoef) real ( kind = 8 ) atac(ncoef,ncoef) real ( kind = 8 ) atae(ncoef,ncoef) real ( kind = 8 ) atan(ncoef,ncoef) real ( kind = 8 ) atas(ncoef,ncoef) real ( kind = 8 ) ataw(ncoef,ncoef) real ( kind = 8 ) b(ncoef) real ( kind = 8 ) dat(12) real ( kind = 8 ) dudyn(np) real ( kind = 8 ) dudyn2(npmax) integer i integer ic integer icol integer ie integer iee integer in integer ine integer inee integer info integer inn integer inne integer inw integer ipivc(ncoef) integer ipive(ncoef) integer ipivn(ncoef) integer ipivs(ncoef) integer ipivw(ncoef) integer irow integer is integer ise integer iss integer issw integer isw integer isww integer iw integer iww integer j integer k integer ndat integer nodat(12) integer nrhs integer nrep(npmax) integer nx integer ny real ( kind = 8 ) temp real ( kind = 8 ) x real ( kind = 8 ) xdatc(12) real ( kind = 8 ) xdatn(7) real ( kind = 8 ) xdate(7) real ( kind = 8 ) xdatw(7) real ( kind = 8 ) xdats(7) real ( kind = 8 ) y real ( kind = 8 ) ydatc(12) real ( kind = 8 ) ydatn(7) real ( kind = 8 ) ydate(7) real ( kind = 8 ) ydatw(7) real ( kind = 8 ) ydats(7) if ( np>npmax ) then write ( *, * ) ' ' write ( *, * ) 'GETDU4 - Fatal error!' write ( *, * ) ' NP must be less than NPMAX, but' write ( *, * ) ' NP = ',np write ( *, * ) ' NPMAX = ',npmax write ( *, * ) ' Change NP or NPMAX and try again!' stop end if do i = 1,np nrep(i) = 0 dudyn2(i) = 0.0D+00 end do ! ! Here is a picture of the element patch around a node on the ! northern boundary: ! ! ! WW---W----C----E----EE ! | /| /| ! | / | / | ! | / | / | ! | / | / | ! SWW SW S SE . ! | / | / | ! | / | / | ! | / | / | ! |/ |/ | ! SSWW-SSW--SS---.----. ! ! NORTH ! ndat = 7 xdatn(1) = 0.0D+00 ydatn(1) = 1.0D+00 xdatn(2) = 1.0D+00 ydatn(2) = 0.0D+00 xdatn(3) = 1.0D+00 ydatn(3) = 1.0D+00 xdatn(4) = 1.0D+00 ydatn(4) = 2.0D+00 xdatn(5) = 2.0D+00 ydatn(5) = 1.0D+00 xdatn(6) = 3.0D+00 ydatn(6) = 1.0D+00 xdatn(7) = 3.0D+00 ydatn(7) = 2.0D+00 ! ! Compute matrix of basis polynomials (1,x,y,x**2,xy,y**2) evaluated ! at the data points. ! do i = 1,ndat an(i,1) = 1.0D+00 an(i,2) = xdatn(i) an(i,3) = ydatn(i) an(i,4) = xdatn(i)**2 an(i,5) = xdatn(i)*ydatn(i) an(i,6) = ydatn(i)**2 end do ! ! Compute the normal equations matrix. ! do i = 1,ncoef do j = 1,ncoef atan(i,j) = 0.0D+00 do k = 1,ndat atan(i,j) = atan(i,j)+an(k,i)*an(k,j) end do end do end do ! ! Factor the normal equations matrix. ! call sgetrf(ncoef,ncoef,atan,ncoef,ipivn,info) if ( info/= 0 ) then write ( *, * ) ' ' write ( *, * ) 'GETDU2 - Fatal error!' write ( *, * ) ' SGETRF returned INFO = ',info write ( *, * ) ' while factoring matrix ATAN.' stop end if ! ! EAST ! ! Here is a picture of the element patch around a node ! on the eastern boundary: ! ! .----.----NN ! | /| ! | / | ! | / | ! | / | ! . NW N ! | / | ! | / | ! | / | ! |/ | ! WW---W----C ! | /| ! | / | ! | / | ! | / | ! SWW SW S ! | / | ! | / | ! | / | ! |/ | ! SSWW-SSW--SS ! ndat = 7 xdate(1) = 0.0D+00 ydate(1) = 1.0D+00 xdate(2) = 1.0D+00 ydate(2) = 0.0D+00 xdate(3) = 1.0D+00 ydate(3) = 1.0D+00 xdate(4) = 1.0D+00 yd