program main !******************************************************************************* ! !! FLOW4 solves a steady viscous 2D Navier-Stokes flow using finite elements. ! ! ! NX = 11, 21, 31, 41, 61, 81, 121, 161, ! NY = 4, 7, 10, 13, 19, 25, 37, 49, ! H = 1, 1/4, 1/6, 1/8, 1/12, 1/16, 1/24, 1/32, ! 0.5, 0.25, 0.166, 0.125, 0.083, 0.0625, 0.04166, 0.03125, ! ! 17 August 1999 ! ! Retrieved FLOW4 from floppy disk. The data had "rotted". There were ! some losses, particulary in subroutine GETDER. Right now, I just ! want to see how the options I used in the Newton iteration, so I'm ! not panicked. ! ! 08 May 1996 ! ! Trying to solve at RE = 100, I'm getting failures. I think I ! could do this, if I turn off the Newton rejection... ! ! 07 May 1996 ! ! Well, I've forgotten where I was on this code, but I need to make ! some runs for optimization over RE and over the Bump. I'm having ! a little trouble getting it started up again. ! ! Getting an index out of bounds error...Tried to run on SUN, ! but missing damn "_CG92_USED_" symbolic name, somehow. ! The only way I could avoid it was to compile everything at once. ! ! 10 November 1995 ! ! Well, I slaved away on the routine to compute the differences ! between U,V,P on a coarse and a fine grid (really so that I ! could, in turn, look at sensitivities) and I'm getting poor ! convergence. I haven't plotted it yet, but it doesn't look good. ! I think I'll can this whole project shortly. ! ! 09 November 1995 ! ! I am worried that my cubic bump is not "consistent" with my ! quadratic elements. So for my comparison runs I'm going to ! look at quadratic bumps only, for a while. ! ! 02 November 1995 ! ! 7 node quadrature routine is in, ready to test. ! ! 01 November 1995 ! ! My studies of the convergence behavior of U don't show H**2 convergence. ! I am going to check this by adding the option of 7 point Gaussian ! quadrature, but this will be a mess. ! ! 26 October 1995 ! ! Got back Cray 615 run: ! ! NX = 11, NY=4, CPU=4 seconds. ! DS-FDC (U,V,P) = (0.177348, 0.132730, 0.124550) ! Try full size problem. ! ! 25 October 1995 ! ! From the SIAM meeting: ! ! "Can decompose any vector field into the gradient of a scalar field ! and something else." -- I could use this fact to compute a ! stream function of vector fields that don't quite satisfy the ! divergence condition, such as Jennifer's least squares flows. ! ! From the SIAM meeting: ! ! "Blah, blah, blah, convert the problem to a quadratic program, ! which only needs the gradients of the cost functional, and solve ! it, so you don't care about consistency between the cost functional ! and the gradients." ! ! From the SIAM meeting: (Matthias) ! ! "SQP method only uses gradients of cost functional." ! ! Lagrangian L(y,g,lambda) = J + (lambda,C(y,g)) ! ! Kuhn Tucker conditions: ! ! Grady J + Cy Lambda = 0 (adjoint) ! Gradg J + Cg Lambda = 0 (gradient) ! C(y,g) = 0 (state) ! ! Use Newton method on this system. ! ! Big goal for San Francisco: ! Slides of O(H) or O(H^2) behavior of (U_ALPHA)^H and (U^H)_ALPHA. ! ! 20 October 1995 ! ! Renamed: ! ! ISHAPB -> IBSCAN ! ISHAPBT -> IBSTAR ! ISHAPF -> IFSCAN ! ISHAPFT -> IFSTAR ! XBLEFT -> XBLCAN ! XBRITE -> XBRCAN ! YBLEFT -> YBLCAN ! YBLEFT -> YBRCAN ! ! 13 October 1995 ! ! I am interested in examining the Zienkiewicz-Zhou method again. ! I suspect that the reason I saw little improvement was that I ! was comparing the finite differences to the discretized sensitivities, ! when I don't know that the finite differences converge to the ! true solution at any particular rate, and I suspect the discretized ! sensitivities converge at an O(h*h) rate. ! ! To verify this, I am going to solve the problem at a sequence of ! values of h, compute dudy in various ways, and look at the behavior ! of dudy as h goes to zero. ! ! If this doesn't work, then I will look at the behavior of U ! as H goes to zero, which should surely show H**2 behavior. ! ! ALSO, I think my serendipity scheme is open to question. I compute ! dUdY on the original grid, then interpolate using serendipity. ! But I don't pass in the X, Y coordinates, so I'm treating dUdY as ! a "basic" variable, when only U should be so treated. So I may ! have to make a new version that can compute dUdY on a serendipity ! element, given U. ! ! 19 August 1995 ! ! Y12M was a big loss. For large problems, the amount of fillin must ! have increased dramatically. I ended up needed as much space as with ! a banded solver. ! ! Then I tried SLAP, and got, as usual with iterative schemes, no ! convergence. ! ! Then I considered GBSOL, and discovered that I had already given up ! on that in June, because GBSOL does not save the LU factors. It is ! a one shot solver only. ! ! Now I am back to a very early project, which is to rewrite SGBFA ! and SGBSL to use a small in core workspace, and a disk file. I have ! actually got it to solve a small problem. Now I have to test it out ! on FLOW. ! ! Switched NODE(I,J) to NODE(J,I) in DISPLAY, FLOW3 and FLOW4. ! ! 14 August 1995 ! ! I am having some success with Y12M. I have gotten it to solve ! a small system, and to solve two systems with the same matrix. ! ! 10 August 1995 ! ! I'm getting an error from NDRV, which looks as though it's ! caused by the fact that the pressure equations don't have ! a pressure coefficient, in other words, the diagonal entry ! of the matrix is zero. ! ! I am going to run a simple test example, but if that bears ! out my suspicion, then either I can't use YALMAT, or I'll have ! to try some sort of fancy manipulation to solve for ! the pressures first. (Matthias did suggest something along ! these lines). ! ! YALMAT routine NDRV fails on the matrix ! ! 1 0 0 ! 0 0 1 ! 0 1 1 ! ! complaining that it has a zero pivot in column 2. ! ! The second alternative is to try ITPACK or NSPCG. ! ! 09 August 1995 ! ! After some exasperation, redimensioned INDX from INDX(MAXNP,3) to ! INDX(3,MAXNP) which makes life easier. ! ! 27 July 1995 ! ! Starting FLOW4.F, archiving FLOW3.F, which had all that optimization ! stuff in it. ! integer, parameter :: maxnpe = 6 integer, parameter :: maxnx = 21 integer, parameter :: maxny = 7 ! ! The assignment of MAXROW should really read (maxrow = 28*min(nx,ny)). ! integer, parameter :: maxdim = 3 integer, parameter :: maxrow = 29 * 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 :: maxparf = 5 integer, parameter :: maxparb = 5 integer, parameter :: maxpar = maxparb+maxparf+1 integer, parameter :: liv = 60 integer, parameter :: lv = 78+maxpar*(maxpar+21)/2 integer, parameter :: maxquad = 7 real ( kind = 8 ) a(maxrow,maxeqn) real ( kind = 8 ) area(maxquad,maxelm) 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 ) 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 ) dpdy2(maxnp) real ( kind = 8 ) dpdy3(maxnp) real ( kind = 8 ) dpdy4(maxnp) real ( kind = 8 ) dudyn(maxnp) real ( kind = 8 ) dudy2(maxnp) real ( kind = 8 ) dudy3(maxnp) real ( kind = 8 ) dudy4(maxnp) real ( kind = 8 ) dudy5(maxnp) real ( kind = 8 ) dvdyn(maxnp) real ( kind = 8 ) dvdy2(maxnp) real ( kind = 8 ) dvdy3(maxnp) real ( kind = 8 ) dvdy4(maxnp) real ( kind = 8 ) dydpn(maxnp,maxparb) real ( kind = 8 ) epsdif character ( len = 2 ) eqn(maxeqn) real ( kind = 8 ) etan(6) real ( kind = 8 ) etaq(maxquad) character ( len = 30 ) fileg character ( len = 30 ) filet real ( kind = 8 ) g(maxeqn) real ( kind = 8 ) g1(maxeqn) real ( kind = 8 ) g2(maxeqn) real ( kind = 8 ) g3(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) real ( kind = 8 ) hx real ( kind = 8 ) hy integer i integer ibc integer icunit integer ifs integer ifscan integer ifstar integer ibs integer ibscan integer ibstar integer ibump integer idfd integer ido integer ids integer ierror integer ifds integer igrad integer igunit integer ijac integer indx(3,maxnp) integer info integer iopt(maxpar) integer ip integer ipivot(maxeqn) integer iplot integer isotri(maxelm) integer istep1 integer istep2 integer itar integer itemp integer itunit integer itype integer iuval integer ivopt(liv) integer iwrite integer j integer jjac integer jstep1 integer jstep2 logical lmat logical lval integer maxnew integer maxstp integer ndim integer nelem integer neqn integer nlband integer node(maxnpe,maxelm) integer nopt integer np integer npar integer nparb integer nparbt integer nparf integer npe integer nprof(2*maxny-1) integer nquad 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 ) partar(maxpar) real ( kind = 8 ) phi(maxquad,6,10,maxelm) 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) character ( len = 20 ) syseqn real ( kind = 8 ) taubmp(maxparb+2) real ( kind = 8 ) tauflo(maxparf+2) character ( len = 30 ) title real ( kind = 8 ) tolnew real ( kind = 8 ) tolopt real ( kind = 8 ) u(maxnp) 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(maxquad) real ( kind = 8 ) xbl real ( kind = 8 ) xblcan real ( kind = 8 ) xbltar real ( kind = 8 ) xbr real ( kind = 8 ) xbrcan real ( kind = 8 ) xbrtar real ( kind = 8 ) xc(maxnp) real ( kind = 8 ) xopt(maxpar) real ( kind = 8 ) xquad(maxquad,maxelm) real ( kind = 8 ) xprof real ( kind = 8 ) xsin(6) real ( kind = 8 ) xsiq(maxquad) real ( kind = 8 ) y2max real ( kind = 8 ) ybl real ( kind = 8 ) yblcan real ( kind = 8 ) ybltar real ( kind = 8 ) ybr real ( kind = 8 ) ybrcan real ( kind = 8 ) ybrtar real ( kind = 8 ) yc(maxnp) real ( kind = 8 ) yquad(maxquad,maxelm) real ( kind = 8 ) yval ! call timestamp ( ) ! ! Say hello. ! call hello ( maxnx, maxny ) ! ! Initialize the variables. ! call init(area,cost,costar,costb,costp,costu,costv,dopt,dpara3,dparsn, & dparfd,dparfdc,dpdyn,dudyn,dvdyn,dydpn,epsdif,eqn,etan,fileg,filet, & g,g1,g2,g3,gdif,gold,gopt,gradf,gtar,ibc,ifscan,ifstar,ibscan,ibstar, & ibump,idfd,ids,ierror,ifds,igrad,igunit,ijac,indx,iopt,iplot,isotri, & istep1,istep2,itar,itunit,itype,ivopt,iwrite,jjac,jstep1,jstep2,liv,lv, & maxelm,maxeqn,maxnew,maxnp,maxpar,maxparb,maxquad,maxstp,node,nopt,npar, & nparb,nparf,npe,nquad,nstep3,nx,ny,para1,para2,para3,partar,sens,syseqn, & tolnew,tolopt,vopt,wateb,wateb1,wateb2,watep,wateu,watev,xblcan,xbltar, & xbrcan,xbrtar,xprof,xsin,yblcan,ybltar,ybrcan,ybrtar) ! ! Read the user input. ! 20 continue call input(epsdif,fileg,filet,ibc,ifscan,ifstar,ibscan,ibstar,ibump,idfd, & ids,ierror,ifds,igrad,ijac,iopt,iplot,istep1,istep2,itar,itype,iwrite, & jjac,jstep1,jstep2,maxnew,maxpar,maxstp,npar,nparb,nparf,nquad,nstep3,nx, & ny,para1,para2,para3,partar,syseqn,tolnew,tolopt,wateb,wateb1,wateb2, & watep,wateu,watev,xblcan,xbltar,xbrcan,xbrtar,xprof,yblcan,ybltar,ybrcan, & ybrtar) if ( ierror == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW4 - GO signal from user input.' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW4 - STOP signal from user.' if ( itunit /= 0 ) then close ( unit = itunit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLOW4 - Note:' write ( *, '(a)' ) ' Closing the marching file ' // trim ( filet ) end if if ( igunit /= 0 ) then close ( unit = igunit ) write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Note:' write ( *, * ) ' Closing the graphics file ' // trim ( fileg ) end if write ( *, * ) ' ' write ( *, * ) 'FLOW4 - End of execution.' write ( *, * ) ' ' stop 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 ( *, * ) 'FLOW4 - Note.' write ( *, * ) ' Number of elements, NELEM = ',nelem write ( *, * ) ' Number of nodes, NP = ',np write ( *, * ) ' Number of nodes per element, NPE = ',npe ! ! Check the data. ! call chkdat(ibump,idfd,ids,ifds,igrad,ijac,iopt,itype,jjac,maxpar,maxparb, & maxparf,nopt,npar,nparb,nparf,nstep3,para1,partar,xblcan,xbrcan) ! ! Print the problem information. ! call prdat(epsdif,fileg,filet,ibc,ifscan,ifstar,ibscan,ibstar,ibump,idfd, & ids,ifds,igrad,ijac,iopt,iplot,istep1,istep2,itar,itype,iwrite,jjac, & jstep1,jstep2,maxnew,maxpar,maxstp,nopt,npar,nparb,nparf,nquad,nstep3,nx, & ny,para1,para2,para3,partar,syseqn,tolnew,tolopt,wateb,wateb1,wateb2, & watep,wateu,watev,xblcan,xbltar,xbrcan,xbrtar,xprof,yblcan,ybltar,ybrcan, & ybrtar) ! ! Open the plot file. ! call pltopn(fileg,igunit,iplot) ! ! Open the marching file. ! call maropn(filet,itunit) ! ! Solve the optimization problem, ! which involves repeated evaluation of the functional ! at various flow solutions (U,V,P)+(FLO,BUM,REY). ! lval = .true. call lmemry('set','target',lval) ifs = ifstar ibs = ibstar xbl = xbltar xbr = xbrtar ybl = ybltar ybr = ybrtar ! ! Set up the elements and nodes. ! write ( *, * ) 'ibump = ',ibump write ( *, * ) 'maxeqn = ',maxeqn write ( *, * ) 'nelem = ',nelem write ( *, * ) 'np = ',np write ( *, * ) 'nx = ',nx write ( *, * ) 'ny = ',ny hx = 10.0D+00 / real (2*(nx-1), kind = 8 ) hy = 3.0D+00 / real (2*(ny-1), kind = 8 ) write ( *, * ) 'HX = ',hx write ( *, * ) 'HY = ',hy call setnod(eqn,ibump,indx,isotri,maxeqn,nelem,neqn,node,np,nx,ny,xbl,xbr) ! ! Find points on velocity profile sampling line. ! itemp = nint ( ( 2.0D+00 * real(nx-1, kind = 8 ) * 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,nlband,node,np,nrow) write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Note:' write ( *, * ) ' Lower bandwidth NLBAND = ',nlband write ( *, * ) ' Total bandwidth = ',2*nlband+1 write ( *, * ) ' Required matrix rows NROW = ',nrow write ( *, * ) ' MAXROW = ',maxrow if ( itar == 0 ) then write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Note:' write ( *, * ) ' Computing the target solution.' para(1:npar) = partar(1:npar) g(1:neqn) = 0.0D+00 ! ! Get the flow solution G. ! call flosol(a,area,eqn,etaq,g,g2,ifs,ibs,ierror,ijac,indx,ipivot,isotri, & iwrite,jjac,maxnew,nelem,neqn,nlband,node,np,npar,nparb,nparf,nquad, & nrow,nx,ny,para,phi,res,splbmp,splflo,syseqn,taubmp,tauflo,tolnew, & wquad,xbl,xbr,xc,xquad,xsiq,ybl,ybr,yc,yquad) if ( ierror/= 0 ) then write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Fatal error!' write ( *, * ) ' FloSol could not compute the target solution.' go to 16 else gtar(1:neqn) = g(1:neqn) end if ! ! Compute the cost function J. ! call getcst(cost,costb,costp,costu,costv,g,gtar,ibs,indx,neqn,np,nparb, & nprof,ny,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) if ( iwrite == 1 ) then write ( *, * ) 'Cost:',cost else if ( iwrite >= 2 ) then title = 'Cost' call prcost2(cost,costb,costp,costu,costv,title,wateb,watep,wateu,watev) end if ! ! Compute the finite difference vectors GDIF. ! call getgrd(a,area,cost,dpara3,epsdif,eqn,etaq,g,g1,g2,gdif,gtar,ifs, & ibs,ierror,ijac,indx,iopt,ipivot,isotri,iwrite,jjac,maxeqn,maxnew, & nelem,neqn,nlband,node,np,npar,nparb,nparf,nprof,nquad,nrow,nx,ny, & para,phi,res,splbmp,splflo,syseqn,taubmp,tauflo,tolnew,wateb,watep, & wateu,watev,wquad,xbl,xbr,xc,xquad,xsiq,ybl,ybr,yc,yquad) if ( ierror /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Fatal error!' write ( *, * ) ' GetGrd returns IERROR = ',ierror stop end if numstp = 0 ! ! Get dPdY, dUdY and dVdY. ! call getdu(dpdyn,dudyn,dvdyn,etan,g,indx,isotri,nelem,neqn,node,np, & numel,xc,xsin,yc) ! ! Get dPdY2, dUdY2, and dVdY2, the ZZ corrections. ! call getdu2(dudyn,dudy2,np,nx,ny) call getdu2(dvdyn,dvdy2,np,nx,ny) call getdu2(dpdyn,dpdy2,np,nx,ny) ! ! Get dPdY3, dUdY3, and dVdY3, the serendipity corrections. ! call getdu3(dudyn,dudy3,np,nx,ny) call getdu3(dvdyn,dvdy3,np,nx,ny) call getdu3(dpdyn,dpdy3,np,nx,ny) ! ! Get dPdY4, dUdY4, and dVdY4, the ZZ modified corrections. ! call getdu4(dudyn,dudy4,np,nx,ny) call getdu4(dvdyn,dvdy4,np,nx,ny) call getdu4(dpdyn,dpdy4,np,nx,ny) ! ! Get dUdY5, the second serendipity corrections. ! do i = 1,np u(i) = g(indx(1,i)) end do call getdu5(dudyn,dudy5,np,nx,ny,u,xc,yc) ! ! Compute the correction term GRADF, and the GDIFC, the corrected ! finite difference estimate of the sensitivity. ! call getfix(dpdyn,dudyn,dvdyn,dydpn,gradf,ibs,indx,iopt,maxeqn,np,npar, & nparb,nparf,splbmp,taubmp,xbl,xbr,xc,yc) do i = 1,neqn do j = 1,npar gdifc(i,j) = gdif(i,j) - gradf(i,j) end do end do ! ! Get the discretized sensitivities. ! lmat = .false. call lmemry('get','have_fp',lmat) if ( .not.lmat ) then call fp(a,area,eqn,g,indx,nelem,neqn,nlband,node,np,npar,nparb,nparf, & nquad,nrow,para,phi,syseqn) call sgbtrf(neqn,neqn,nlband,nlband,a,nrow,ipivot,info) if ( info/= 0 ) then write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Fatal error!' write ( *, * ) ' Jacobian factorization failed.' write ( *, * ) ' SGBTRF returns nonzero INFO = ',info else lmat = .true. call lmemry('set','have_fp',lmat) end if end if if ( lmat ) then call getsen(a,area,dudyn,dvdyn,eqn,g,ibc,ifs,ibs,indx,iopt,ipivot, & maxeqn,nelem,neqn,nlband,node,np,npar,nparb,nparf,nquad,nrow, & phi,sens,splbmp,splflo,taubmp,tauflo,xc,yc) end if gold(1:neqn) = g(1:neqn) if ( iwrite == 1 ) then write ( *, * ) 'Cost: ',cost else if (iwrite >= 2 ) then title = 'Cost' call prcost2(cost,costb,costp,costu,costv,title,wateb,watep,wateu,watev) end if call getder(dparfd,g,gtar,ibs,indx,maxeqn,neqn,np,npar,nparb,nparf, & nprof,ny,gdif,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) call getder(dparfdc,g,gtar,ibs,indx,maxeqn,neqn,np,npar,nparb,nparf, & nprof,ny,gdifc,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) call getder(dparsn,g,gtar,ibs,indx,maxeqn,neqn,np,npar,nparb,nparf,nprof, & ny,sens,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) call prsol(dpara3,dparfd,dparfdc,dparsn,g,g2,gdif,gdifc,idfd,ids,ifds, & indx,iwrite,maxeqn,neqn,np,npar,nparb,nparf,numstp,para,sens) costar = cost else if (itar == 1 ) then write ( *, * ) 'FLOW4 - Note:' write ( *, * ) ' Calling GetTar2 to get target data.' nparbt = 0 ibs = 0 call setxy(ibs,np,nparbt,nx,ny,splbmp,taubmp,xbltar,xbrtar,xc,ybltar, & 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 do i = 1,2*ny-1 iuval = indx(1,nprof(i)) 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 discst(costp,costu,costv,g,gtar,indx,neqn,np,nprof,ny,yc) cost = costu costar = cost gtar(1:neqn) = g(1:neqn) end if ! ! HARDCODE OPTION: ! Write comparison data to file. ! ido = -1 if ( ido == 0 ) then icunit = 17 open ( unit = icunit, file = 'compare.dat', form = 'formatted', & access = 'sequential', status = 'replace' ) call writes(ibs,icunit,nelem,neqn,np,nparb,nquad,nx,ny,xbl,xbr) call writev(area,g,indx,isotri,icunit,nelem,neqn,node,np,nparb,nquad, & splbmp,taubmp,xc,yc) close ( unit = icunit ) write ( *, * ) ' ' write ( *, * ) 'Wrote comparison data to file COMP.DAT.' write ( *, * ) 'Stopping now.' stop end if ! ! Write target information to plot file. ! ! HARDCODE OPTION: ! Plot uncorrected finite difference sensitivities GDIF, ! or corrected finite difference sensitivities GDIFC. ! if ( iplot /= 0 ) then ido = 1 if ( ido == 0 ) then call pltwrt(eqn,g,gdif,igunit,indx,isotri,iwrite,maxeqn,nelem,neqn, & node,np,npar,npe,nprof,nx,ny,para,sens,xc,xprof,yc) else if (ido == 1 ) then call pltwrt(eqn,g,gdifc,igunit,indx,isotri,iwrite,maxeqn,nelem,neqn, & node,np,npar,npe,nprof,nx,ny,para,sens,xc,xprof,yc) end if end if ! ! Turn off target calculation flag. ! lval = .false. call lmemry('set','target',lval) lmat = .false. call lmemry('set','have_fp',lmat) ! ! Set data specific to the optimization. ! ifs = ifscan ibs = ibscan xbl = xblcan xbr = xbrcan ybl = yblcan ybr = ybrcan ! ! Set the elements and nodes. ! if ( xblcan /= xbltar .or. xbrcan /= xbrtar .or. & yblcan /= ybltar .or. ybrcan /= ybrtar ) then call setnod(eqn,ibump,indx,isotri,maxeqn,nelem,neqn,node,np,nx,ny,xbl,xbr) end if ! ! The user only wanted to compute the target point, and nothing more. ! if ( itype == 0 ) then go to 20 ! ! Is this a march? ! else 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,dpara3,dparsn,dparfd,dparfdc,dpdyn,dudyn, & dvdyn,dydpn,epsdif,eqn,etan,etaq,g,g1,g2,g3,gdif,gdifc,gold,gradf, & gtar,ibc,ifs,ibs,idfd,ids,ifds,igunit,ijac,indx,iopt,ipivot,iplot, & isotri,istep1,istep2,itunit,iwrite,jjac,jstep1,jstep2,maxeqn,maxnew, & maxstp,ndim,nelem,neqn,nlband,node,np,npar,nparb,nparf,npe,nprof, & nquad,nrow,nstep3,numel,nx,ny,para,para1,para2,para3,phi,res,sens, & splbmp,splflo,syseqn,taubmp,tauflo,tolnew,wateb,wateb1,wateb2,watep, & wateu,watev,wquad,xbl,xbr,xc,xprof,xquad,xsin,xsiq,ybl,ybr,yc,yquad ) ! ! ...or a sensitivity optimization? ! ...or an optimization using function values only? ! else if (itype ==3.or.itype==7 ) then do i = 1,npar para(i) = para1(i) end do call qsolve(a,area,dopt,dpara3,dparfd,dparfdc,dparsn,dpdyn,dudyn,dvdyn, & dydpn,epsdif,eqn,etan,etaq,g,g1,g2,gdif,gdifc,gold,gopt,gradf,gtar, & ibc,ifs,ibs,idfd,ids,ifds,igrad,igunit,ijac,indx,iopt,ipivot,iplot, & isotri,itunit,itype,ivopt,iwrite,jjac,liv,lv,maxeqn,maxnew,maxstp, & nelem,neqn,nlband,node,nopt,np,npar,nparb,nparf,npe,nprof,nquad,nrow, & numel,nx,ny,para,phi,res,sens,splbmp,splflo,syseqn,taubmp,tauflo, & tolnew,tolopt,vopt,wateb,watep,wateu,watev,wquad,xbl,xbr,xc,xopt,xprof, & xquad,xsin,xsiq,ybl,ybr,yc,yquad) else write ( *, * ) ' ' write ( *, * ) 'FLOW4 - Fatal error!' write ( *, * ) ' Unknown value of ITYPE = ',itype stop end if 16 continue ! ! If an error occurred, then wipe out the current data. ! if ( ierror/= 0 ) then ierror = 0 write ( *, * ) 'FLOW4 - Clearing error flag after failure.' end if ! ! See if user wishes to change variables and continue run. ! go to 20 end subroutine bmpbc(dudyn,dvdyn,ibs,ip,iparb,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc) !******************************************************************************* ! !! BMPBC computes the value of the boundary conditions for the ! horizontal and vertical velocity sensitivities with respect ! to a given shape parameter. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! IP Input, the global node number. ! ! IPARB Input, integer IPARB, the bump parameter with respect to ! which the sensitivities are desired. ! ! UBC Output, the boundary condition for the horizontal velocity ! sensitivity. ! ! VBC Output, the boundary condition for the vertical velocity ! sensitivity. ! ! integer np ! real ( kind = 8 ) dudyn(np) real ( kind = 8 ) dvdyn(np) integer ibs integer ip integer iparb 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 ( ibs == 1 ) then call plval1(iparb+1,nparb+2,xc(ip),taubmp,shape) else if (ibs == 2 ) then call pqval1(iparb+1,nparb+2,xc(ip),taubmp,shape) else if (ibs ==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 bmpbc1(g,ibs,indx,ip,iparb,neqn,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc,yc) ! !******************************************************************************* ! !! BMPBC1 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! G Input, real G(MAXEQN), the finite element ! coefficients. ! ! IP Input, the global node number. ! ! IPARB Input, integer IPARB, the bump parameter with respect to ! which the sensitivities are desired. ! ! UBC Output, the boundary condition for the horizontal velocity ! sensitivity. ! ! VBC Output, the boundary condition for the vertical velocity ! sensitivity. ! ! integer neqn integer np ! real ( kind = 8 ) dudy real ( kind = 8 ) dvdy real ( kind = 8 ) g(neqn) integer ibs integer ihor integer ihorn integer indx(3,np) integer ip integer iparb 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 ( ibs == 1 ) then call plval1(iparb+1,nparb+2,xc(ip),taubmp,shape) else if (ibs == 2 ) then call pqval1(iparb+1,nparb+2,xc(ip),taubmp,shape) else if (ibs ==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(1,ip) ihorn = indx(1,ip+1) dudy = (g(ihorn)-g(ihor))/(yc(ip+1)-yc(ip)) iver = indx(2,ip) ivern = indx(2,ip+1) dvdy = (g(ivern)-g(iver))/(yc(ip+1)-yc(ip)) ! ! Set the boundary conditions. ! ubc = -dudy*shape vbc = -dvdy*shape return end subroutine bmpbc2(g,ibs,indx,ip,iparb,neqn,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc,yc) ! !******************************************************************************* ! !! BMPBC2 computes the value of 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! G Input, real G(MAXEQN), the finite element ! coefficients. ! ! IP Input, the global node number. ! ! IPARB Input, integer IPARB, the bump parameter with respect to ! which the sensitivities are desired. ! ! UBC Output, the boundary condition for the horizontal velocity ! sensitivity. ! ! VBC Output, the boundary condition for the vertical velocity ! sensitivity. ! ! 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 ibs integer indx(3,np) integer ip integer iparb 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 ( ibs == 1 ) then call plval1(iparb+1,nparb+2,xc(ip),taubmp,shape) else if (ibs == 2 ) then call pqval1(iparb+1,nparb+2,xc(ip),taubmp,shape) else if (ibs ==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(1,ip)) u1 = g(indx(1,ip+1)) u2 = g(indx(1,ip+2)) dudy = (h1**2*u2-h2**2*u1+(h2**2-h1**2)*u0)/(h1*h2*(h1-h2)) v0 = g(indx(2,ip)) v1 = g(indx(2,ip+1)) v2 = g(indx(2,ip+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 bmpcst(costb,ibs,nparb,splbmp,taubmp,xbl,xbr,ybl,ybr) ! !******************************************************************************* ! !! BMPCST evaluates the cost of the bump control. ! ! 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! COSTB Output, real COSTB. ! ! COSTB is 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. ! ! integer nparb ! integer, parameter :: nquad1 = 5 ! real ( kind = 8 ) costb real ( kind = 8 ) cprime integer i integer ibs 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)return if ( xbl>= xbr)return ! ! 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+2-i, kind = 8 )*xbl & +real(i-1, kind = 8 )*xbr) & /real(nparb+1, kind = 8 ) xrite = (real(nparb+1-i, kind = 8 )*xbl+real(i, kind = 8 )*xbr)/real(nparb+1, kind = 8 ) do j = 1,nquad1 xx = 0.5*((1.0+xquad1(j))*xrite+(1.0-xquad1(j))*xleft) if ( ibs == 1 ) then yvec(1:nparb+2) = splbmp(1,1:nparb+2,0) call pldx(nparb+2,xx,taubmp,cprime,yvec) else if (ibs == 2 ) then yvec(1:nparb+2) = splbmp(1,1:nparb+2,0) call pqdx(nparb+2,xx,taubmp,cprime,yvec) else if (ibs ==3 ) then jderiv = 1 call ppvalu(taubmp,splbmp,nparb+1,4,xx,jderiv,cprime) end if costb = costb+0.5*wquad1(j)*(xrite-xleft)*(cprime-slope)**2 end do end do return end subroutine bmpder(dpara,ibs,npar,nparb,nparf,splbmp,taubmp,wateb,xbl,xbr, & ybl,ybr) ! !******************************************************************************* ! !! BMPDER evaluates the derivative of the cost of the bump control ! with respect to the parameters. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! integer npar integer nparb ! integer, parameter :: nquad1 = 5 ! real ( kind = 8 ) cprime real ( kind = 8 ) dpara(npar) integer i integer ibs 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)return if ( xbr<= xbl)return ! ! 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+2-i, kind = 8 )*xbl+real(i-1, kind = 8 )*xbr)/real(nparb+1, kind = 8 ) xrite = (real(nparb+1-i, kind = 8 )*xbl+real(i, kind = 8 )*xbr)/real(nparb+1, kind = 8 ) do j = 1,nquad1 xx = 0.5*((1.0+xquad1(j))*xrite+(1.0-xquad1(j))*xleft) if ( ibs == 1 ) then yvec(1:nparb+2) = splbmp(1,1:nparb+2,0) call pldx(nparb+2,xx,taubmp,cprime,yvec) else if (ibs == 2 ) then yvec(1:nparb+2) = splbmp(1,1:nparb+2,0) call pqdx(nparb+2,xx,taubmp,cprime,yvec) else if (ibs ==3 ) then jderiv = 1 call ppvalu(taubmp,splbmp,nparb+1,4,xx,jderiv,cprime) end if do k = 1,nparb if ( ibs == 1 ) then call pldx1(k+1,nparb+2,xx,taubmp,temp) else if (ibs == 2 ) then call pqdx1(k+1,nparb+2,xx,taubmp,temp) else if (ibs ==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 bmpsen(dudyn,dvdyn,eqn,f,g,ibc,ibs,indx,ipar,neqn,np,nparb,nparf, & splbmp,taubmp,xc,yc) ! !******************************************************************************* ! !! BMPSEN 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! 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 ibs integer ihor integer indx(3,np) integer ip integer ipar integer iparb 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 do i = 1,neqn f(i) = 0.0D+00 end do do ip = 1,np ihor = indx(1,ip) if ( eqn(ihor) =='UB' ) then if ( ibc == 0 ) then call bmpbc(dudyn,dvdyn,ibs,ip,iparb,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc) else if (ibc == 1 ) then call bmpbc1(g,ibs,indx,ip,iparb,neqn,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc,yc) else if (ibc == 2 ) then call bmpbc2(g,ibs,indx,ip,iparb,neqn,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc,yc) else if (ibc ==3 ) then call bmpbc(dudyn,dvdyn,ibs,ip,iparb,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc) end if f(ihor) = ubc end if iver = indx(2,ip) if ( eqn(iver) =='VB' ) then if ( ibc == 0 ) then call bmpbc(dudyn,dvdyn,ibs,ip,iparb,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc) else if (ibc == 1 ) then call bmpbc1(g,ibs,indx,ip,iparb,neqn,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc,yc) else if (ibc == 2 ) then call bmpbc2(g,ibs,indx,ip,iparb,neqn,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc,yc) else if (ibc ==3 ) then call bmpbc(dudyn,dvdyn,ibs,ip,iparb,np,nparb,shape,splbmp,taubmp, & ubc,vbc,xc) end if f(iver) = vbc end if end do return end subroutine bmpspl(ibs,npar,nparb,nparf,par,splbmp,taubmp,xbl,xbr,ybl,ybr) ! !******************************************************************************* ! !! BMPSPL sets up or updates the spline data that describes the bump. ! ! It does this for the target parameters and the feasible parameters. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! IBS Input, integer IBS. ! 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. ! ! integer npar integer nparb integer nparf ! integer i integer ibcbeg integer ibcend integer ibs integer ipar 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)return ! ! 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, kind = 8 )*xbl+real(i-1, kind = 8 )*xbr)/real(nparb+1, kind = 8 ) end do splbmp(1,1,0) = ybl do i = 2, nparb+1 splbmp(1,i,0) = par(nparf+i-1) end do splbmp(1,nparb+2,0) = ybr if ( ibs == 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 bsp(q,dqdx,dqdy,ielem,iq,nelem,node,np,xc,xq,yc,yq) ! !******************************************************************************* ! !! BSP computes the value and spatial derivatives of the linear basis ! functions associated with pressure. ! ! ! Here is a picture of a typical finite element associated with ! pressure: ! ! 2 ! /| ! / | ! / | ! 1---3 ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Q Output, real Q, the value of the IQ-th basis ! function at the point with global coordinates (XQ,YQ). ! ! DQDX, ! DQDY Output, real DQDX, DQDY, the X and Y ! derivatives of the IQ-th basis function at the point ! with global coordinates (XQ,YQ). ! ! IELEM Input, integer IELEM, the global element number about which ! we are inquiring. ! ! IQ 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. ! ! NELEM Input, integer NELEM, the number of elements. ! ! NODE Input, integer NODE(6,NELEM). NODE(J,I) is ! the global node number of the J-th node in the I-th ! element. ! ! NP Input, integer NP, the number of nodes. ! ! XC Input, real XC(NP), the global X coordinates ! of the element nodes. ! ! XQ Input, real XQ, the global X coordinate of ! the point in which we are interested. ! ! YC Input, real YC(NP), the global Y coordinates ! of the element nodes. ! ! YQ Input, real YQ, the global Y coordinate of ! the point in which we are interested. ! ! integer nelem integer np ! 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(6,nelem) real ( kind = 8 ) xc(np) real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) yq ! if ( iq<1.or.iq>6 ) 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<=6 ) 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(iq1,ielem) i2 = node(iq2,ielem) i3 = node(iq3,ielem) 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.0 + dqdx*(xq-xc(i1)) + dqdy*(yq-yc(i1)) return end subroutine ch_cap ( c ) ! !******************************************************************************* ! !! CH_CAP capitalizes a single character. ! ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! 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,itype,jjac,maxpar, & maxparb,maxparf,nopt,npar,nparb,nparf,nstep3,para1,partar,xblcan,xbrcan) ! !******************************************************************************* ! !! CHKDAT performs some simple checks on the input data. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! integer maxpar ! integer i integer ibump integer idfd integer ids integer ifds integer igrad integer ijac integer iopt(maxpar) 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 ) sum real ( kind = 8 ) xblcan real ( kind = 8 ) xbrcan ! ! 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 ! ! ITYPE ! if ( itype<0.or.itype>7 ) then write ( *, * ) ' ' write ( *, * ) 'ChkDat - Fatal error!' write ( *, * ) ' Illegal value of ITYPE = ',itype write ( *, * ) ' Legal values are between 0 and 7.' write ( *, * ) ' ChkDat forces a STOP!' stop end if ! ! JJAC ! if ( jjac<0.or.jjac>2 ) then write ( *, * ) ' ' write ( *, * ) 'ChkDat - Fatal error!' write ( *, * ) ' Illegal value of JJAC = ',jjac write ( *, * ) ' Legal values are between 0 and 2.' write ( *, * ) ' ChkDat forces a STOP!' stop end if ! ! NOPT ! if ( nopt<= 0 ) then write ( *, * ) ' ' write ( *, * ) 'ChkDat - Warning!' write ( *, * ) ' The number of free parameters, NOPT = ',nopt write ( *, * ) ' but this value should be positive!' 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) ! sum = 0.0D+00 do i = 1,nparf sum = sum+abs(para1(i)) end do if ( sum == 0.0 ) then write ( *, * ) ' ' write ( *, * ) 'ChkDat - Fatal error!' write ( *, * ) ' All PARA1 inflows are zero.' stop end if ! ! PARA1(NPARF+NPARB+1), the value of NU_INV. ! if ( para1(nparf+nparb+1)<0.0 ) 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(1:NPARF) ! sum = 0.0D+00 do i = 1,nparf sum = sum+abs(partar(i)) end do if ( sum == 0.0 ) then write ( *, * ) ' ' write ( *, * ) 'ChkDat - Fatal error!' write ( *, * ) ' All PARTAR inflows are zero.' stop end if ! ! PARTAR(NPARF+NPARB+1), the value of NU_INV. ! if ( partar(nparf+nparb+1)<0.0 ) 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 ( xblcan>= xbrcan ) then write ( *, * ) ' ' write ( *, * ) 'ChkDat - Fatal error!' write ( *, * ) ' XBLEFT >= XBRITE.' write ( *, * ) ' XBLEFT = ',xblcan write ( *, * ) ' XBRITE = ',xbrcan stop end if end if return end subroutine chrctd(string,dval,ierror,lchar) ! !******************************************************************************* ! !! 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) ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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 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. ! ! 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+1len1)return if ( i>len2)return if ( i+len3-1<= len1 ) then if ( string(i:i+len3-1) ==strng3)return end if strng2(i:i) = string(i:i) go to 10 end subroutine cubspl(tau,c,n,ibcbeg,ibcend) ! !******************************************************************************* ! !! CUBSPL is given data and boundary conditions for a cubic ! spline, and computes information that defines that cubic ! spline. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! TAU Input, real TAU(N), the abscissas or X values of ! the data points. The entries of TAU are assumed to be ! strictly increasing. ! ! N Input, integer N, the number of data points. N is ! assumed to be at least 2. ! ! C Input/output, real C(4,N). ! ! On input, if IBCBEG or IBCBEG is 1 or 2, then C(2,1) ! or C(2,N) should have been set to the desired derivative ! values, as described further under IBCBEG and IBCEND. ! ! On output, C contains the polynomial coefficients of ! the cubic interpolating spline with interior knots ! TAU(2) through TAU(N-1). ! ! In the interval interval (TAU(I), TAU(I+1)), the spline ! F is given by ! ! F(X) = C(1,I)+H*C(2,I)+(1/2)*H*H*C(3,I) ! +(1/6)*H*H*H*C(4,I) ! ! where H = X - TAU(I). The routine PPVALU may be used to ! evaluate F or its derivatives from TAU, C, L = N-1, ! and K = 4. ! ! IBCBEG, ! IBCEND Input, integer IBCBEG, IBCEND, boundary condition ! indicators. ! ! IBCBEG = 0 means no boundary condition at TAU(1) is given. ! In this case, the "not-a-knot condition" is used. That ! is, the jump in the third derivative across TAU(2) is ! forced to zero. Thus the first and the second cubic ! polynomial pieces are made to coincide. ! ! IBCBEG = 1 means that the slope at TAU(1) is to equal the ! input value C(2,1). ! ! IBCBEG = 2 means that the second derivative at TAU(1) is ! to equal C(2,1). ! ! IBCEND = 0, 1, or 2 has analogous meaning concerning the ! boundary condition at TAU(N), with the additional ! information taken from C(2,N). ! ! integer n ! real ( kind = 8 ) c(4,n) real ( kind = 8 ) divdf1 real ( kind = 8 ) divdf3 real ( kind = 8 ) dtau real ( kind = 8 ) g integer i integer ibcbeg integer ibcend integer j integer m real ( kind = 8 ) tau(n) ! ! A tridiagonal linear system for the unknown slopes S(I) of ! F at TAU(I), I = 1,..., N, is generated and then solved by Gauss ! elimination, with S(I) ending up in C(2,I), for all I. ! ! C(3,*) and C(4,*) are used initially for temporary storage. ! ! Store first differences of the TAU sequence in C(3,*). ! ! Store first divided difference of data in C(4,*). ! do m = 2,n c(3,m) = tau(m)-tau(m-1) c(4,m) = (c(1,m)-c(1,m-1))/c(3,m) end do ! ! Construct the first equation from the boundary condition, of ! the form: ! ! C(4,1)*S(1) + C(3,1)*S(2) = C(2,1) ! if ( ibcbeg == 1 ) then c(4,1) = 1.0D+00 c(3,1) = 0.0D+00 go to 60 end if if ( ibcbeg<= 1 ) then ! ! No condition at left end and N = 2. ! if ( n<= 2 ) then c(4,1) = 1.0D+00 c(3,1) = 1.0D+00 c(2,1) = 2.0*c(4,2) go to 120 end if ! ! Not-a-knot condition at left end and N is greater than 2. ! c(4,1) = c(3,3) c(3,1) = c(3,2)+c(3,3) c(2,1) = ((c(3,2)+2.0*c(3,1))*c(4,2)*c(3,3)+c(3,2)**2*c(4,3))/c(3,1) go to 70 end if ! ! Second derivative prescribed at left end. ! c(4,1) = 2.0D+00 c(3,1) = 1.0D+00 c(2,1) = 3.0*c(4,2)-c(3,2)/2.0*c(2,1) 60 continue if ( n == 2)go to 120 ! ! If there are interior knots, generate the corresponding ! equations and carry out the forward pass of Gauss elimination, ! after which the M-th equation reads: ! ! C(4,M)*S(M) + C(3,M)*S(M+1) = C(2,M). ! 70 continue do m = 2,n-1 g = -c(3,m+1)/c(4,m-1) c(2,m) = g*c(2,m-1)+3.0*(c(3,m)*c(4,m+1)+c(3,m+1)*c(4,m)) c(4,m) = g*c(3,m-1)+2.0*(c(3,m)+c(3,m+1)) end do ! ! Construct last equation from the second boundary condition, of ! the form ! ! (-G*C(4,N-1))*S(N-1) + C(4,N)*S(N) = C(2,N) ! ! If slope is prescribed at right end, one can go directly to ! back-substitution, since the C array happens to be set up just ! right for it at this point. ! if ( ibcend == 1)go to 160 if ( ibcend>1)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.0*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.0*c(4,n) c(4,n) = 1.0D+00 g = -1.0/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.0*c(4,n)+c(3,n)/2.0*c(2,n) c(4,n) = 2.0D+00 g = -1.0/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.0*c(4,n)+c(3,n)/2.0*c(2,n) c(4,n) = 2.0D+00 g = -1.0/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.0*divdf1 c(3,i-1) = 2.0*(divdf1-c(2,i-1)-divdf3)/dtau c(4,i-1) = (divdf3/dtau)*(6.0/dtau) end do return end subroutine discst(costp,costu,costv,g,gtar,indx,neqn,np,nprof,ny,yc) ! !******************************************************************************* ! !! DISCST computes the discrepancy integrals along the profile line. ! ! ! Discussion: ! ! There are discrepancy integrals for the pressure, ! horizontal and vertical velocities. ! ! 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! 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(3,np) 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(1,nprof(ii)) ucof(k) = g(j)-gtar(j) j = indx(2,nprof(ii)) 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.5 * wquad1(j)*(yhi-ylo)*uval**2 costv = costv + 0.5 * 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(3,nprof(ii)) 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.5*wquad1(j)*(yhi-ylo)*pval**2 end do end do return end subroutine disder(dpara,g,gtar,indx,maxeqn,neqn,np,npar,nprof,ny,sens, & watep,wateu,watev,yc) ! !******************************************************************************* ! !! DISDER computes the derivative of the discrepancy integrals for the ! pressure, horizontal and vertical velocities, with respect to the ! various parameters. ! ! NOTES: ! ! 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! DPARA Input/output, real DPARA(NPAR). ! ! DPARA contains the derivatives of the cost function with ! respect to the various parameters. In particular, ! DPARA(I) = D cost / D parameter(I). ! ! integer, parameter :: mpar = 10 integer, parameter :: nquad1 = 5 ! integer maxeqn 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(3,np) 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(maxeqn,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 ( *, * ) 'DisDer - 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(1,nprof(ii)) ucof(k) = g(j)-gtar(j) do l = 1,npar ducof(k,l) = sens(j,l) end do j = indx(2,nprof(ii)) 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.5*wquad1(j)*(yhi-ylo) & *2.0*(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(3,nprof(ii)) 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.5*wquad1(j)*(yhi-ylo)*2.0*watep*pval*dpval(l) end do end do end do return end subroutine floduv(ifs,ipar,nparf,splflo,tauflo,ubc,vbc,yy) ! !*************************************************************************** ! !! FLODUV sets the derivative of the parabolic inflow in terms of ! the value of the inflow parameters. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! UBC Output, real UBC. ! The value of the partial derivative of the horizontal ! boundary velocity component with respect to a given ! parameter at the specified point. ! ! VBC Output, real VBC. ! The value of the partial derivative of the vertical ! boundary velocity component with respect to a given ! parameter at the specified point. ! ! integer nparf ! integer ifs integer ipar integer jderiv real ( kind = 8 ) splflo(4,nparf+2,0:nparf) real ( kind = 8 ) tauflo(nparf+2) real ( kind = 8 ) ubc real ( kind = 8 ) vbc real ( kind = 8 ) yy ! ubc = 0.0D+00 vbc = 0.0D+00 ! ! Handle the case where there is only one inflow parameter, and it ! is specified as a constant. ! if ( nparf<= 0.and.ipar== 0 ) then ubc = 4.0*(3.0-yy)*yy/(9.0) return end if if ( ipar<1.or.ipar>nparf)return ! ! Evaluate the basis function which is zero at the other parameters ! and 1 at parameter IPAR. ! if ( ifs == 1 ) then call plval1(ipar+1,nparf+2,yy,tauflo,ubc) else if (ifs == 2 ) then call pqval1(ipar+1,nparf+2,yy,tauflo,ubc) else if (ifs ==3 ) then jderiv = 0 call ppvalu(tauflo,splflo(1,1,ipar),nparf+1,4,yy,jderiv,ubc) else write ( *, * ) ' ' write ( *, * ) 'FloDUV - Fatal error!' write ( *, * ) ' Illegal value of ISHAPF = ',ifs stop end if return end subroutine flosen(eqn,f,ifs,indx,ipar,neqn,np,nparf,splflo,tauflo,yc) ! !******************************************************************************* ! !! FLOSEN 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 inflow parameter. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! F Output, real F(NEQN), the right hand ! side of the sensitivity equations associated with ! the IPAR-th inflow parameter. ! ! integer neqn integer np integer nparf ! character ( len = 2 ) eqn(neqn) real ( kind = 8 ) f(neqn) integer i integer ifs integer ihor integer indx(3,np) integer ip integer ipar integer iver real ( kind = 8 ) splflo(4,nparf+2,0:nparf) real ( kind = 8 ) tauflo(nparf+2) real ( kind = 8 ) ubc real ( kind = 8 ) vbc real ( kind = 8 ) yc(np) ! do i = 1,neqn f(i) = 0.0D+00 end do do ip = 1,np ihor = indx(1,ip) if ( eqn(ihor) =='UI' ) then call floduv(ifs,ipar,nparf,splflo,tauflo,ubc,vbc,yc(ip)) f(ihor) = ubc end if iver = indx(2,ip) if ( eqn(iver) =='VI' ) then call floduv(ifs,ipar,nparf,splflo,tauflo,ubc,vbc,yc(ip)) f(iver) = vbc end if end do return end subroutine flosol(a,area,eqn,etaq,g,g2,ifs,ibs,ierror,ijac,indx,ipivot, & isotri,iwrite,jjac,maxnew,nelem,neqn,nlband,node,np,npar,nparb,nparf, & nquad,nrow,nx,ny,para,phi,res,splbmp,splflo,syseqn,taubmp,tauflo,tolnew, & wquad,xbl,xbr,xc,xquad,xsiq,ybl,ybr,yc,yquad) ! !******************************************************************************* ! !! FLOSOL is given a set of flow parameters in PARA, and an ! approximate solution vector G, and proceeds to set up the ! constraints associated with PARA, and use Newton iteration ! to correct G to a solution that satisfies the constraints ! to within some tolerance. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real G(NEQN), the computed solution vector, in which are stored ! pressures and velocities. ! ! Output, integer IERROR, an error flag. ! 0 means no error, nonzero means an error. ! integer nelem integer neqn integer np integer npar integer nparb integer nparf integer nquad integer nrow integer nx integer ny ! real ( kind = 8 ) a(nrow,neqn) real ( kind = 8 ) area(nquad,nelem) character ( len = 2 ) eqn(neqn) real ( kind = 8 ) etaq(nquad) real ( kind = 8 ) g(neqn) real ( kind = 8 ) g2(neqn) integer ifs integer ibs integer ierror integer ijac integer indx(3,np) integer ipivot(neqn) integer isotri(nelem) integer iwrite integer jjac logical s_eqi integer maxnew integer nlband integer node(6,nelem) real ( kind = 8 ) para(npar) real ( kind = 8 ) phi(nquad,6,10,nelem) real ( kind = 8 ) res(neqn) real ( kind = 8 ) splbmp(4,nparb+2,0:nparb) real ( kind = 8 ) splflo(4,nparf+2,0:nparf) character ( len = 20 ) syseqn real ( kind = 8 ) taubmp(nparb+2) real ( kind = 8 ) tauflo(nparf+2) real ( kind = 8 ) tolnew real ( kind = 8 ) wquad(nquad) real ( kind = 8 ) xbl real ( kind = 8 ) xbr real ( kind = 8 ) xc(np) real ( kind = 8 ) xquad(nquad,nelem) real ( kind = 8 ) xsiq(nquad) real ( kind = 8 ) ybl real ( kind = 8 ) ybr real ( kind = 8 ) yc(np) real ( kind = 8 ) yquad(nquad,nelem) ! external s_eqi ! ! Set the spline coefficients for the bump. ! call bmpspl(ibs,npar,nparb,nparf,para,splbmp,taubmp,xbl,xbr,ybl,ybr) ! ! Set the spline coefficients for the inflow. ! call flospl(ifs,npar,nparf,para,splflo,tauflo) ! ! Set the X and Y coordinates of the nodes that form the grid. ! call setxy(ibs,np,nparb,nx,ny,splbmp,taubmp,xbl,xbr,xc,ybl,ybr,yc) ! ! Set the quadrature points, which move every step if there ! are bump parameters. ! if ( nquad ==3 ) then call setq3(area,etaq,isotri,nelem,node,np,nquad,wquad,xc, & xquad,xsiq,yc,yquad) else if (nquad ==7 ) then call setq7(area,etaq,isotri,nelem,node,np,nquad,wquad,xc, & xquad,xsiq,yc,yquad) end if ! ! Set the value of the basis functions at all quadrature points. ! call setbas(area,etaq,isotri,nelem,node,np,nquad,phi,xc,xquad,xsiq,yc,yquad) ! ! Try to solve the full nonlinear system. ! call newton(a,area,eqn,g,g2,ifs,ierror,ijac,indx,ipivot,iwrite,jjac,maxnew, & nelem,neqn,nlband,node,np,npar,nparb,nparf,nquad,nrow,para,phi,res,splflo, & syseqn,tauflo,tolnew,yc) if ( ierror/= 0.and.s_eqi(syseqn,'NavierStokes') ) then write ( *, * ) ' ' write ( *, * ) 'FloSol - Warning!' write ( *, * ) ' Newton failed on '//syseqn write ( *, * ) ' The parameters at which failure occurred:' write ( *, * ) ' ' call prpar(nparb,nparf,para) syseqn = 'Stokes' write ( *, * ) ' Switching to '//syseqn ! ! Try to solve the linearized system. ! call newton ( a, area, eqn, g, g2, ifs, ierror, ijac, indx, ipivot, & iwrite,jjac,maxnew, & nelem,neqn,nlband,node,np,npar,nparb,nparf,nquad,nrow,para,phi,res, & splflo,syseqn,tauflo,tolnew,yc) syseqn = 'NavierStokes' if ( ierror/= 0 ) then write ( *, * ) ' ' write ( *, * ) 'FloSol - Fatal error!' write ( *, * ) ' Newton failed!' write ( *, * ) ' The parameters at which failure occurred:' write ( *, * ) ' ' call prpar(nparb,nparf,para) ierror = 1 return end if end if return end subroutine flospl(ifs,npar,nparf,par,splflo,tauflo) ! !******************************************************************************* ! !! FLOSPL sets up or updates the spline data that describes the inflow. ! ! It does this for the target parameters and the feasible parameters. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! integer npar integer nparf ! integer i integer ibcbeg integer ibcend integer ifs integer ihi integer ipar real ( kind = 8 ) par(npar) real ( kind = 8 ) splflo(4,nparf+2,0:nparf) real ( kind = 8 ) tauflo(nparf+2) ! ! Set up the abscissas in TAUFLO, which will be evenly spaced. ! If there are no flow parameters, then we generate three evenly ! spaced nodes. ! ihi = nparf+2 do i = 1,ihi tauflo(i) = real(i-1, kind = 8 )*3.0/real(ihi-1, kind = 8 ) end do ! ! Set up the coefficient arrays, including ! SPLFLO(1,I,0), the strength of the inflow at abscissa I, ! SPLFLO(1,I,IPAR), the I-th coefficient of the partial derivative ! of the inflow with respect to the IPAR-th inflow parameter. ! if ( ifs == 1.or.ifs== 2 ) then do i = 1,nparf+2 if ( i == 1 ) then splflo(i,1,0) = 0.0D+00 else if (2<= i.and.i<=nparf+1 ) then splflo(i,1,0) = par(i-1) else if (i ==nparf+2 ) then splflo(i,1,0) = 0.0D+00 end if end do else if (ifs ==3 ) then do i = 1,nparf+2 if ( i == 1 ) then splflo(1,i,0) = 0.0D+00 else if (2<= i.and.i<=nparf+1 ) then splflo(1,i,0) = par(i-1) else if (i ==nparf+2 ) then splflo(1,i,0) = 0.0D+00 end if do ipar = 1,nparf if ( ipar+1/= i ) then splflo(1,i,ipar) = 0.0D+00 else splflo(1,i,ipar) = 1.0D+00 end if end do end do ibcbeg = 0 ibcend = 0 do i = 0,nparf call cubspl(tauflo,splflo(1,1,i),nparf+2,ibcbeg,ibcend) end do end if return end subroutine flouv(ifs,nparf,splflo,tauflo,ubc,vbc,yy) ! !*************************************************************************** ! !! FLOUV computes the specified boundary values of velocity for a ! given position as determined by the value of the parameters. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! integer nparf ! integer ifs integer jderiv real ( kind = 8 ) splflo(4,nparf+2,0:nparf) real ( kind = 8 ) tauflo(nparf+2) real ( kind = 8 ) ubc real ( kind = 8 ) vbc real ( kind = 8 ) yy ! if ( ifs == 1 ) then call plval(nparf+2,yy,tauflo,ubc,splflo) else if (ifs == 2 ) then call pqval(nparf+2,yy,tauflo,ubc,splflo) else if (ifs ==3 ) then jderiv = 0 call ppvalu(tauflo,splflo(1,1,0),nparf+1,4,yy,jderiv,ubc) else write ( *, * ) ' ' write ( *, * ) 'FloUV - Fatal error!' write ( *, * ) ' Illegal value of ISHAPF = ',ifs stop end if vbc = 0.0D+00 return end subroutine fp ( a, area, eqn, g, indx, nelem, neqn, nlband, node, np, & npar, nparb, nparf, nquad, nrow, para, phi, syseqn ) ! !******************************************************************************* ! !! FP sets up a matrix associated with Navier-Stokes or Stokes flow. ! ! ! Discussion: ! ! If the calling routine is trying to solve the Stokes equations, ! (which are linearized Navier-Stokes equations), then ! this routine sets up the appropriate system matrix. ! ! If the calling routine is trying to solve the nonlinear Navier-Stokes ! equations, then this routine computes the jacobian of the ! finite element form of the Navier Stokes equations. ! ! The differentiated finite element 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 ! ! d P-Eqn/d V-Coef: ! ! Integral ! ! dWj/dy * Qi dx dy ! ! ! The Stokes equations 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 ! ! d P-Eqn/d V-Coef: ! ! Integral ! ! dWj/dy * Qi dx dy ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! A Output, real A(NROW,NEQN), 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). ! ! AREA Input, real AREA(3,NELEM). ! ! AREA contains a common factor multiplying the term associated ! with a quadrature point in a given element, namely, ! ! AREA(IQUAD,IELEM) = Ar(IELEM) * WQUAD(IQUAD) ! ! or, if the element is isoperimetric, ! ! AREA(IQUAD,IELEM) = DET * Ar(IELEM) * WQUAD(IQUAD) ! ! Here Ar(IELEM) represents the area of element IELEM. ! ! EQN Input, character ( len = 2 ) EQN(NEQN). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. Note that most boundary ! conditions do not result in an equation. The current values are: ! ! 'U' The horizontal momentum equation. ! 'UB' The condition U = 0 applied at a node on the bump. ! 'UI' The condition U = UInflow(Y,Lambda) at the inflow. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! ! 'V' The vertical momentum equation. ! 'VB' The condition V = 0 applied at a node on the bump. ! 'VI' The condition V = VInflow(Y,Lambda) at the inflow. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! ! 'P' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! G Input, real G(NEQN). ! ! G is the current solution vector, in which are stored ! the finite element coefficients that define the velocity ! and pressure functions, U, V and P. ! ! INDX Input, integer INDX(3,NP). ! ! INDX(I,J) contains, for each node J, 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), ! and an equation will be generated to determine its value. ! ! If INDX(I,J) is not positive, then no equation is ! generated to determine for variable I at node J, either because ! the variable is specified in some other way, or because ! (in the case of pressure), there is no coefficient associated ! with that node. ! ! NELEM Input, integer NELEM, the number of elements. ! ! NEQN Input, integer NEQN, the number of finite element equations used ! to define the horizontal and vertical velocities and the ! pressure. ! ! NLBAND Input, integer NLBAND. ! ! The lower bandwidth of the matrix A. The zero structure of A ! is assumed to be symmetric, and so NLBAND is also the upper ! bandwidth of A. ! ! NODE Input, integer NODE(6,NELEM). ! ! NODE(I,J) contains, for an element J, the global node index of ! the element node whose local number is I. ! ! The local ordering of the nodes is suggested by this diagram: ! ! 2 ! /| ! 4 5 ! / | ! 1-6-3 ! ! NP Input, integer NP, the number of nodes used to define the finite ! element mesh. NP = (2*NX-1)*(2*NY-1). ! ! NPAR 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. ! ! NPARB Input, integer NPARB. ! ! The number of parameters associated with the position and ! shape of the bump. ! ! Note that if NPARB = 0, the bump is replaced by a flat wall. ! ! NPARF Input, integer NPARF. ! ! NPARF is the number of parameters associated with the ! inflow. NPARF must be at least 1. ! ! NROW Input, integer NROW. ! ! The number of rows need to store the matrix A, using the ! LINPACK/LAPACK general banded storage format. NROW must be ! at least 3*NLBAND+1. ! ! PARA Input, real PARA(NPAR). ! ! PARA is the current set of parameter values, including the ! inverse viscosity, the flow parameters, and the bump parameters. ! ! PHI Input, real PHI(NQUAD,6,10,NELEM). ! ! PHI contains the value of a basis function, its derivative, ! or other information, evaluated at a quadrature point. ! ! For a particular element I, quadrature point J, and basis ! function K, we use the following shorthand for the 10 ! entries of PHI: ! ! W, dWdX, dWdY ! Q, dQdX, dQdY ! dXsidX, dXsidY, dEtadX, dEtadY ! ! W is the quadratic basis function associated with velocity, ! Q the linear basis function associated with pressure, ! Xsi and Eta the reference coordinates for the point. ! ! In particular, PHI(J,K,1,I) is the value of the quadratic ! basis function associated with local node K in element I, ! evaluated at quadrature point J. ! ! Note that PHI(J,K,4,I) = PHI(J,K,5,I)=PHI(J,K,6,I)=0 for ! K = 4, 5, or 6, since there are only three linear basis ! functions. ! ! SYSEQN Input, character ( len = 20 ) SYSEQN. ! ! Input to the FX and FP routines, SYSEQN is either 'NavierStokes' ! or 'Stokes', and specifies which state equation is to be set up. ! ! integer nelem integer neqn integer np integer npar integer nparf integer nquad integer nrow ! real ( kind = 8 ) a(nrow,neqn) real ( kind = 8 ) ar real ( kind = 8 ) area(nquad,nelem) 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(3,np) integer ip integer iprs integer iq integer iquad integer iuse integer iver integer j integer jhor integer jp integer jprs integer jq integer jver logical s_eqi integer nlband integer node(6,nelem) integer nparb real ( kind = 8 ) pold real ( kind = 8 ) para(npar) real ( kind = 8 ) phi(nquad,6,10,nelem) 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 ! nu_inv = para(nparf+nparb+1) do i = 1,nrow do j = 1,neqn a(i,j) = 0.0D+00 end do end do ! ! Approximate the integral by summing over all elements. ! do ielem = 1,nelem ! ! Evaluate the integrand at the quadrature points. ! do iquad = 1,nquad ar = area(iquad,ielem) ! ! 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,nquad,pold,phi,uold,vold) ! ! Consider each node in the element. ! do iq = 1,6 ip = node(iq,ielem) wi = phi(iquad,iq,1,ielem) dwidx = phi(iquad,iq,2,ielem) dwidy = phi(iquad,iq,3,ielem) qi = phi(iquad,iq,4,ielem) ihor = indx(1,ip) iver = indx(2,ip) iprs = indx(3,ip) ! ! 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(jq,ielem) wj = phi(iquad,jq,1,ielem) dwjdx = phi(iquad,jq,2,ielem) dwjdy = phi(iquad,jq,3,ielem) dqjdx = phi(iquad,jq,5,ielem) dqjdy = phi(iquad,jq,6,ielem) jhor = indx(1,jp) jver = indx(2,jp) jprs = indx(3,jp) ! ! Contributions of the JHOR horizontal velocity to the U, V, and ! P equations. ! iuse = ihor-jhor+2*nlband+1 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 a(iuse,jhor) = a(iuse,jhor) + term end if if ( eqn(iver) == 'V' ) then if ( s_eqi(syseqn,'NavierStokes') ) then iuse = iver - jhor+2*nlband+1 term = ar * nu_inv * wj * dvdx * wi a(iuse,jhor) = a(iuse,jhor)+term end if end if if ( 0 < iprs ) 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(1,ip) iver = indx(2,ip) iprs = indx(3,ip) 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 fx ( area,eqn,g,ifs,indx,nelem,neqn,node,np,npar,nparb,nparf, & nquad,para,phi,res,splflo,syseqn,tauflo,yc) ! !******************************************************************************* ! !! 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: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real AREA(3,NELEM), contains a common factor multiplying the ! term associated with a quadrature point in a given element, namely, ! AREA(IQUAD,IELEM) = Ar(IELEM) * WQUAD(IQUAD) ! or, if the element is isoperimetric, ! AREA(IQUAD,IELEM) = DET * Ar(IELEM) * WQUAD(IQUAD) ! Here Ar(IELEM) represents the area of element IELEM. ! ! EQN Input, character ( len = 2 ) EQN(NEQN). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. Note that most boundary ! conditions do not result in an equation. The current values are: ! ! 'U' The horizontal momentum equation. ! 'UB' The condition U = 0 applied at a node on the bump. ! 'UI' The condition U = UInflow(Y,Lambda) at the inflow. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! ! 'V' The vertical momentum equation. ! 'VB' The condition V = 0 applied at a node on the bump. ! 'VI' The condition V = VInflow(Y,Lambda) at the inflow. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! ! 'P' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! G Input, real G(NEQN). ! ! G is the current solution vector, in which are stored ! the finite element coefficients that define the velocity ! and pressure functions, U, V and P. ! ! IFS Input, integer IFS. ! 1, the inflow is modeled by C0 linear splines. ! 2, the inflow is modeled by C0 quadratic splines. ! 3, the inflow is modeled by C1 cubic splines. ! ! INDX Input, integer INDX(3,NP). ! ! INDX(I,J) contains, for each node J, 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), ! and an equation will be generated to determine its value. ! ! If INDX(I,J) is not positive, then no equation is ! generated to determine for variable I at node J, either because ! the variable is specified in some other way, or because ! (in the case of pressure), there is no coefficient associated ! with that node. ! ! NELEM Input, integer NELEM, the number of elements. ! ! NEQN Input, integer NEQN, the number of finite element equations used ! to define the horizontal and vertical velocities and the ! pressure. ! ! NODE Input, integer NODE(6,NELEM). ! ! NODE(I,J) contains, for an element J, the global node index of ! the element node whose local number is I. ! ! The local ordering of the nodes is suggested by this diagram: ! ! 2 ! /| ! 4 5 ! / | ! 1-6-3 ! ! NP Input, integer NP, the number of nodes used to define the finite ! element mesh. NP = (2*NX-1)*(2*NY-1). ! ! NPAR 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. ! ! NPARB Input, integer NPARB. ! ! The number of parameters associated with the position and ! shape of the bump. ! ! Note that if NPARB = 0, the bump is replaced by a flat wall. ! ! NPARF Input, integer NPARF. ! ! NPARF is the number of parameters associated with the ! inflow. NPARF must be at least 1. ! ! NQUAD Input, integer NQUAD, the number of quadrature points. ! ! PARA Input, real PARA(NPAR). ! ! PARA is the current set of parameter values, including the ! inverse viscosity, the flow parameters, and the bump parameters. ! ! PHI Input, real PHI(NQUAD,6,10,NELEM). ! ! PHI contains the value of a basis function, its derivative, ! or other information, evaluated at a quadrature point. ! ! For a particular element I, quadrature point J, and basis ! function K, we use the following shorthand for the 10 ! entries of PHI: ! ! W, dWdX, dWdY ! Q, dQdX, dQdY ! dXsidX, dXsidY, dEtadX, dEtadY ! ! W is the quadratic basis function associated with velocity, ! Q the linear basis function associated with pressure, ! Xsi and Eta the reference coordinates for the point. ! ! In particular, PHI(J,K,1,I) is the value of the quadratic ! basis function associated with local node K in element I, ! evaluated at quadrature point J. ! ! Note that PHI(J,K,4,I) = PHI(J,K,5,I)=PHI(J,K,6,I)=0 for ! K = 4, 5, or 6, since there are only three linear basis ! functions. ! ! RES Output, real RES(NEQN), contains the value ! of the residual. ! ! SPLFLO Input, real SPLFLO(4,NPARF+2,0:NPARF). ! ! SPLFLO contains the spline coefficients for the inflow ! in SPLFLO(*,*,0). ! ! SYSEQN Input, character ( len = 20 ) SYSEQN. ! ! Input to the FX and FP routines, SYSEQN is either 'NavierStokes' ! or 'Stokes', and specifies which state equation is to be set up. ! ! TAUFLO Input, real TAUFLO(NPARF+2). ! ! TAUFLO contains the location of the spline abscissas for ! the inflow. There are NPARF+2 of them, because the end ! values of the spline are constrained to have particular ! values. ! ! YC Input, real YC(NP). ! ! The Y coordinates of the nodes. ! ! integer nelem integer neqn integer np integer npar integer nparf integer nquad ! real ( kind = 8 ) ar real ( kind = 8 ) area(nquad,nelem) 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 ifs integer ielem integer ihor integer indx(3,np) integer ip integer iprs integer iq integer iquad integer iver logical s_eqi integer node(6,nelem) integer nparb real ( kind = 8 ) p real ( kind = 8 ) para(npar) real ( kind = 8 ) phi(nquad,6,10,nelem) 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 ! 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,nquad ar = area(iquad,ielem) ! ! Evaluate P, U and V. ! call uvalq(dpdx,dpdy,dudx,dudy,dvdx,dvdy,g,ielem,indx, & iquad,nelem,neqn,node,np,nquad,p,phi,u,v) ! ! Look at nearby basis functions. ! do iq = 1,6 ip = node(iq,ielem) wi = phi(iquad,iq,1,ielem) dwidx = phi(iquad,iq,2,ielem) dwidy = phi(iquad,iq,3,ielem) qi = phi(iquad,iq,4,ielem) ! ! The horizontal velocity equations. ! ihor = indx(1,ip) 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 flouv ( ifs, 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(2,ip) 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 flouv ( ifs, 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(3,ip) if ( 0 < iprs ) 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 getcst(cost,costb,costp,costu,costv,g,gtar,ibs,indx,neqn,np, & nparb,nprof,ny,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) ! !******************************************************************************* ! !! GETCST is given the value of the solution, G, the target ! solution GTAR, and information about the shape of the bump, ! and returns the value of the overall and individual cost ! functions. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! 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 ibs integer indx(3,np) 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) ! call bmpcst(costb,ibs,nparb,splbmp,taubmp,xbl,xbr,ybl,ybr) call discst(costp,costu,costv,g,gtar,indx,neqn,np,nprof,ny,yc) cost = wateb*costb+watep*costp+wateu*costu+watev*costv return end subroutine getder(dpara,g,gtar,ibs,indx,maxeqn,neqn,np,npar,nparb,nparf, & nprof,ny,sens,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) ! !******************************************************************************* ! !! GETDER returns the derivatives of the cost with respect to the parameters. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! DPARA Output, real DPARA(NPAR), the partial ! derivative of the cost with respect to each parameter, ! computed using sensitivities or finite differences. ! integer maxeqn integer neqn integer np integer npar integer nparb integer nparf integer ny ! real ( kind = 8 ) dpara(npar) real ( kind = 8 ) g(neqn) real ( kind = 8 ) gtar(neqn) integer i integer ibs integer indx(3,np) integer nprof(2*ny-1) real ( kind = 8 ) sens(maxeqn,npar) 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) return end subroutine getdu(dpdyn,dudyn,dvdyn,etan,g,indx,isotri,nelem,neqn,node,np, & numel,xc,xsin,yc) ! !******************************************************************************* ! !! GETDU uses the finite element method to estimate the value of ! the spatial derivatives dPdY, dUdY and dVdY at every node. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! integer nelem integer neqn integer np ! 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,6) 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. ! do ip = 1,np numel(ip) = 0 dpdyn(ip) = 0.0D+00 dudyn(ip) = 0.0D+00 dvdyn(ip) = 0.0D+00 end do 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, & 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,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 getdu2 ( dudyn, dudy2, np, nx, ny ) ! !******************************************************************************* ! !! GETDU2 is damaged code. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! integer, parameter :: maxdat = 12 integer, parameter :: ncoef = 100 ! integer np ! real ( kind = 8 ) ac(maxdat,6) real ( kind = 8 ) ae(maxdat,6) real ( kind = 8 ) an(maxdat,6) real ( kind = 8 ) as(maxdat,6) 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 ) aw(maxdat,6) real ( kind = 8 ) b(ncoef) real ( kind = 8 ) dat(maxdat) real ( kind = 8 ) dudyn(np) real ( kind = 8 ) dudy2(np) integer i integer ic integer icol integer ie integer iee integer in integer ine integer inee integer info 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(maxdat) integer nrhs integer nx integer ny real ( kind = 8 ) x real ( kind = 8 ) xdat(maxdat) real ( kind = 8 ) y real ( kind = 8 ) ydat(maxdat) ! ! Copy the data ! do i = 1,np dudy2(i) = dudyn(i) 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 xdat(1) = 0.0D+00 ydat(1) = 1.0D+00 xdat(2) = 1.0D+00 ydat(2) = 0.0D+00 xdat(3) = 1.0D+00 ydat(3) = 1.0D+00 xdat(4) = 1.0D+00 ydat(4) = 2.0D+00 xdat(5) = 2.0D+00 ydat(5) = 1.0D+00 xdat(6) = 3.0D+00 ydat(6) = 1.0D+00 xdat(7) = 3.0D+00 ydat(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) = xdat(i) an(i,3) = ydat(i) an(i,4) = xdat(i)**2 an(i,5) = xdat(i)*ydat(i) an(i,6) = ydat(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 xdat(1) = 0.0D+00 ydat(1) = 1.0D+00 xdat(2) = 1.0D+00 ydat(2) = 0.0D+00 xdat(3) = 1.0D+00 ydat(3) = 1.0D+00 xdat(4) = 1.0D+00 ydat(4) = 2.0D+00 xdat(5) = 1.0D+00 ydat(5) = 3.0D+00 xdat(6) = 2.0D+00 ydat(6) = 1.0D+00 xdat(7) = 2.0D+00 ydat(7) = 3.0D+00 ! ! Compute matrix of basis polynomials (1,x,y,x**2,xy,y**2) evaluated ! at the data points. ! do i = 1,ndat ae(i,1) = 1.0D+00 ae(i,2) = xdat(i) ae(i,3) = ydat(i) ae(i,4) = xdat(i)**2 ae(i,5) = xdat(i)*ydat(i) ae(i,6) = ydat(i)**2 end do ! ! Compute the normal equations matrix. ! do i = 1,ncoef do j = 1,ncoef atae(i,j) = 0.0D+00 do k = 1,ndat atae(i,j) = atae(i,j)+ae(k,i)*ae(k,j) end do end do end do ! ! Factor the normal equations matrix. ! call sgetrf(ncoef,ncoef,atae,ncoef,ipive,info) if ( info/= 0 ) then write ( *, * ) ' ' write ( *, * ) 'GetDU2 - Fatal error!' write ( *, * ) ' SGETRF returned INFO = ',info write ( *, * ) ' while factoring matrix ATAE.' stop end if ! ! WEST ! ! Here is a picture of the element patch around a node on the ! western boundary: ! ! NN---NNE--NNEE ! | /| ! | / | ! | / | ! | / | ! N NE NEE ! | / | ! | / | ! | / | ! |/ | ! C----E----EE ! | /| ! | / | ! | / | ! | / | ! S SE . ! | / | ! | / | ! | / | ! |/ | ! SS---.----. ! ndat = 7 xdat(1) = 2.0D+00 ydat(1) = 1.0D+00 xdat(2) = 2.0D+00 ydat(2) = 3.0D+00 xdat(3) = 3.0D+00 ydat(3) = 1.0D+00 xdat(4) = 3.0D+00 ydat(4) = 2.0D+00 xdat(5) = 3.0D+00 ydat(5) = 3.0D+00 xdat(6) = 3.0D+00 ydat(6) = 4.0D+00 xdat(7) = 4.0D+00 ydat(7) = 3.0D+00 ! ! Compute matrix of basis polynomials (1,x,y,x**2,xy,y**2) evaluated ! at the data points. ! do i = 1,ndat aw(i,1) = 1.0D+00 aw(i,2) = xdat(i) aw(i,3) = ydat(i) aw(i,4) = xdat(i)**2 aw(i,5) = xdat(i)*ydat(i) aw(i,6) = ydat(i)**2 end do ! ! Compute the normal equations matrix. ! do i = 1,ncoef do j = 1,ncoef ataw(i,j) = 0.0D+00 do k = 1,ndat ataw(i,j) = ataw(i,j)+aw(k,i)*aw(k,j) end do end do end do ! ! Factor the normal equations matrix. ! call sgetrf(ncoef,ncoef,ataw,ncoef,ipivw,info) if ( info/= 0 ) then write ( *, * ) ' ' write ( *, * ) 'GetDU2 - Fatal error!' write ( *, * ) ' SGETRF returned INFO = ',info write ( *, * ) ' while factoring matrix ATAW.' stop end if ! ! SOUTH ! ! Here is a picture of the element patch around a node on the ! southern boundary: ! ! .----.----NN---NNE--NNEE ! | /| /| ! | / | / | ! | / | / | ! | / | / | ! . NW N NE NEE ! | / | / | ! | / | / | ! | / | / | ! |/ |/ | ! WW---W----C----E----EE ! ndat = 7 xdat(1) = 1.0D+00 ydat(1) = 2.0D+00 xdat(2) = 1.0D+00 ydat(2) = 3.0D+00 xdat(3) = 2.0D+00 ydat(3) = 3.0D+00 xdat(4) = 3.0D+00 ydat(4) = 2.0D+00 xdat(5) = 3.0D+00 ydat(5) = 3.0D+00 xdat(6) = 3.0D+00 ydat(6) = 4.0D+00 xdat(7) = 4.0D+00 ydat(7) = 3.0D+00 ! ! Compute matrix of basis polynomials (1,x,y,x**2,xy,y**2) evaluated ! at the data points. ! do i = 1,ndat as(i,1) = 1.0D+00 as(i,2) = xdat(i) as(i,3) = ydat(i) as(i,4) = xdat(i)**2 as(i,5) = xdat(i)*ydat(i) as(i,6) = ydat(i)**2 end do ! ! Compute the normal equations matrix. ! do i = 1,ncoef do j = 1,ncoef atas(i,j) = 0.0D+00 do k = 1,ndat atas(i,j) = atas(i,j)+as(k,i)*as(k,j) end do end do end do ! ! Factor the normal equations matrix. ! call sgetrf(ncoef,ncoef,atas,ncoef,ipivs,info) if ( info /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'GetDU2 - Fatal error!' write ( *, * ) ' SGETRF returned INFO = ',info write ( *, * ) ' while factoring matrix ATAS.' stop end if ! ! CENTRAL. ! ndat = 12 xdat(1) = 0.0D+00 ydat(1) = 1.0D+00 xdat(2) = 1.0D+00 ydat(2) = 0.0D+00 xdat(3) = 1.0D+00 ydat(3) = 1.0D+00 xdat(4) = 1.0D+00 ydat(4) = 2.0D+00 xdat(5) = 1.0D+00 ydat(5) = 3.0D+00 xdat(6) = 2.0D+00 ydat(6) = 1.0D+00 xdat(7) = 2.0D+00 ydat(7) = 3.0D+00 xdat(8) = 3.0D+00 ydat(8) = 1.0D+00 xdat(9) = 3.0D+00 ydat(9) = 2.0D+00 xdat(10) = 3.0D+00 ydat(10) = 3.0D+00 xdat(11) = 3.0D+00 ydat(11) = 4.0D+00 xdat(12) = 4.0D+00 ydat(12) = 3.0D+00 ! ! Compute matrix of basis polynomials (1,x,y,x**2,xy,y**2) evaluated ! at the data points. ! do i = 1,ndat ac(i,1) = 1.0D+00 ac(i,2) = xdat(i) ac(i,3) = ydat(i) ac(i,4) = xdat(i)**2 ac(i,5) = xdat(i)*ydat(i) ac(i,6) = ydat(i)**2 end do ! ! Compute the normal equations matrix. ! do i = 1,ncoef do j = 1,ncoef atac(i,j) = 0.0D+00 do k = 1,ndat atac(i,j) = atac(i,j)+ac(k,i)*ac(k,j) end do end do end do ! ! Factor the normal equations matrix. ! call sgetrf(ncoef,ncoef,atac,ncoef,ipivc,info) if ( info /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'GetDU2 - Fatal error!' write ( *, * ) ' SGETRF returned INFO = ',info write ( *, * ) ' while factoring matrix ATAC.' stop end if ! ! Each computation has its "center" at a pressure node, so we ! can organize our do loops by counting through each presure node. ! do icol = 1,nx do irow = 1,ny ic = 2*(icol-1)*(2*ny-1)+(2*irow-1) isww = ic-2*(2*ny-1)-1 iww = isww+1 issw = ic-(2*ny-1)-2 isw = issw+1 iw = isw+1 inw = iw+1 iss = ic-2 is = iss+1 in = ic+1 ise = ic+(2*ny-1)-1 ie = ise+1 ine = ie+1 inne = ine+1 iee = ic+2*(2*ny-1) inee = iee+1 ! ! CASE: CORNER NODES, SKIP EM ! if ( icol == 1.and.irow== 1 ) then else if (icol == 1.and.irow==ny ) then else if (icol ==nx.and.irow== 1 ) then else if (icol ==nx.and.irow==ny ) then ! ! CASE: THE NORTHERN BOUNDARY ! else if (irow ==ny ) then nodat(1) = isww nodat(2) = issw nodat(3) = isw nodat(4) = iw nodat(5) = is nodat(6) = ise nodat(7) = ie ndat = 7 ! ! Copy out the value of dUdY at the data points. ! do i = 1,ndat dat(i) = dudyn(nodat(i)) end do ! ! Compute right hand side of normal equations. ! do i = 1,ncoef b(i) = 0.0D+00 do j = 1,ndat b(i) = b(i)+an(j,i)*dat(j) end do end do ! ! Solve the normal equations. ! nrhs = 1 call sgetrs('N',ncoef,nrhs,atan,ncoef,ipivn,b,ncoef,info) ! ! Now evaluate the least squares polynomial at the nodes. ! if ( icol == 2 ) then x = 0.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(iww)) end if x = 2.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(ic)) if ( icol ==nx-1 ) then x = 4.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(iee)) end if ! ! CASE: THE EASTERN BOUNDARY ! else if (icol ==nx ) then nodat(1) = isww nodat(2) = issw nodat(3) = isw nodat(4) = iw nodat(5) = inw nodat(6) = is nodat(7) = in ndat = 7 ! ! Copy out the value of dUdY at the data points. ! do i = 1,ndat dat(i) = dudyn(nodat(i)) end do ! ! Compute right hand side of normal equations. ! do i = 1,ncoef b(i) = 0.0D+00 do j = 1,ndat b(i) = b(i)+ae(j,i)*dat(j) end do end do ! ! Solve the normal equations. ! nrhs = 1 call sgetrs('N',ncoef,nrhs,atae,ncoef,ipive,b,ncoef,info) ! ! Now evaluate the least squares polynomial at the nodes. ! x = 2.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(ic)) ! ! CASE: THE WESTERN BOUNDARY ! else if (icol == 1 ) then nodat(1) = is nodat(2) = in nodat(3) = ise nodat(4) = ie nodat(5) = ine nodat(6) = inne nodat(7) = inee ndat = 7 ! ! Copy out the value of dUdY at the data points. ! do i = 1,ndat dat(i) = dudyn(nodat(i)) end do ! ! Compute right hand side of normal equations. ! do i = 1,ncoef b(i) = 0.0D+00 do j = 1,ndat b(i) = b(i)+aw(j,i)*dat(j) end do end do ! ! Solve the normal equations. ! nrhs = 1 call sgetrs('N',ncoef,nrhs,ataw,ncoef,ipivw,b,ncoef,info) ! ! Now evaluate the least squares polynomial at the nodes. ! x = 2.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(ic)) ! ! CASE: THE SOUTHERN BOUNDARY ! else if (irow == 1 ) then nodat(1) = iw nodat(2) = inw nodat(3) = in nodat(4) = ie nodat(5) = ine nodat(6) = inne nodat(7) = inee ndat = 7 ! ! Copy out the value of dUdY at the data points. ! do i = 1,ndat dat(i) = dudyn(nodat(i)) end do ! ! Compute right hand side of normal equations. ! do i = 1,ncoef b(i) = 0.0D+00 do j = 1,ndat b(i) = b(i)+as(j,i)*dat(j) end do end do ! ! Solve the normal equations. ! nrhs = 1 call sgetrs('N',ncoef,nrhs,atas,ncoef,ipivs,b,ncoef,info) ! ! Now evaluate the least squares polynomial at the nodes. ! if ( icol == 2 ) then x = 0.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(iww)) end if x = 2.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(ic)) if ( icol ==nx-1 ) then x = 4.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(iee)) end if ! ! CASE: CENTRAL ! else nodat(1) = isww nodat(2) = issw nodat(3) = isw nodat(4) = iw nodat(5) = inw nodat(6) = is nodat(7) = in nodat(8) = ise nodat(9) = ie nodat(10) = ine nodat(11) = inne nodat(12) = inee ndat = 12 ! ! Copy out the value of dUdY at the data points. ! do i = 1,ndat dat(i) = dudyn(nodat(i)) end do ! ! Compute right hand side of normal equations. ! do i = 1,ncoef b(i) = 0.0D+00 do j = 1,ndat b(i) = b(i)+ac(j,i)*dat(j) end do end do ! ! Solve the normal equations. ! nrhs = 1 call sgetrs('N',ncoef,nrhs,atac,ncoef,ipivc,b,ncoef,info) ! ! Now evaluate the least squares polynomial at the nodes. ! x = 2.0D+00 y = 2.0D+00 call lspoly(b,ncoef,x,y,dudy2(ic)) end if end do end do return end subroutine getdu3(dudyn,dudy3,np,nx,ny) ! !******************************************************************************* ! !! GETDU3 uses interpolating serendipity quadrilaterals in an attempt ! to improve the accuracy of the computed values of dPdY, dUdY and ! dVdY at every node. ! ! Note a possible objection to this scheme: we are treating ! dUdY as primary data, instead of U. Perhaps, instead, we sh