program main !*****************************************************************************80 ! !! MAIN is the main program for BUMP_ORIGINAL. ! ! Discussion: ! ! BUMP_ORIGINAL is the main program for the original version of the ! flow-with-bump code. ! ! Diary: ! ! 8 September 1992 ! ! The PVAL routine does not handle isoparametric elements ! correctly. ! ! G is now always the current solution. F is just an auxilliary ! vector. I also added SENS for the sensitivities. ! ! I have made GDUMP, UVDUMP and XYDUMP active on every iteration ! of the loop now, although I don't think PLOT3D can handle ! multiple snapshots. ! ! Options: I could start working on continuation now. ! ! 4 September 1992 ! ! Added LONG to GDUMP output. ! ! 3 September 1992 ! ! Replaced NLBUMP and NRBUMP by XBLEFT and XBRITE, ! and compute the nodes based on the X's. ! ! Added IBUMP, which allows me to use no, some, or all isoparametric ! elements. ! ! 2 September 1992 ! ! I moved the profile line from X = 9 to X=4. So now is ! the time to make that a variable, defined by XPROF. ! ! Now, set NX = 11, NY=4. Results are not good. ! ! One explanation for this could be that my elements are ! not all isoparametric above the bump. ! ! Fixed the problem. It occurred because I was treating the ! area of isometric elements incorrectly! ! ! I'm happy now. I use 0 for the initial guess in the secant ! iteration. And I copy f into g properly...I'm back to ! 4 secant, 12 newton. ! ! 1 September 1992 ! ! Unrolled loops in LINSYS that ran over unknowns, 1 to 3. ! The good news: convergence in 40 seconds, with 3 secant ! steps, and 10 Newton steps. ! The bad news: I can tell, from NSTOKE's difference report, ! that the solutions differ. ! ! I did a profile analysis today. Remember, this is the ! "buggy" code. Here are the results: ! ! % time procedure ! ! 54.4 lapack ! 24.6 linsys ! 3.7 lapack ! 2.8 lapack ! 2.5 qbf ! 1.5 uval ! 0.7 setbas ! 0.4 trans ! 0.3 bsp ! 0.2 setban ! ! 31 August 1992 ! ! ! Replacing calls to QBF/REFQBF and BSP/REFBSP by lookups in ! PHI and PSI. Code now runs in .976 minutes = 58 seconds. ! And there's some more optimization I can do. ! ! 27 August 1992 ! ! I FOUND THE ERROR! The quadrature points change from step to step, ! in the NON isoparametric elements! ! ! 24 August 1992 ! ! With the tighter variable ordering, BUMP carries out 4 secant ! iterations, 12 Newton iterations, in 1.38 minutes. ! ! 21 August 1992 ! ! Code written 10 July takes 555 seconds (9.25 minutes) to run ! 21 by 7 problem. ! ! Converting LINPACK calls to LAPACK calls brings the time ! down to 262 seconds (4.38 minutes). ! ! ! List of routines: ! ! BSP is given a point (X,Y), and evaluates, for a particular ! triangle, the linear basis functions associated with ! pressure. ! ! BUMP is the main program. It initializes data, and controls ! an iteration which is intended to produce a value of the ! parameter A which produces a desired internal flow profile. ! ! DELETE delete old copies of files before a new copy is opened. ! ! GDUMP writes information about a particular time step to a file for ! use by a graphics display program. ! ! GETG returns a vector containing the values of a quantity at each ! node where the internal profile is taken, by extracting those ! values from a global list. ! ! GRAM computes the Gram matrix, GR(I,J) = INTEGRAL PHI(I)*PHI(J), ! and the line integral of the internal velocity profile, ! INTEGRAL UI*PHI(J). ! ! IGETL returns the local unknown number, given the global unknown ! number for any node on the line X = XZERO, where the internal ! velocity profile is taken. ! ! LINSYS sets up and solves the Navier Stokes linear system. ! ! NSTOKE solves the Navier Stokes equations for a given set of ! boundary conditions, including the inflow determined by ! the parameter A. The solution is done via Newton's method. ! ! PVAL computes a table of pressures for UVDUMP. ! ! QBF is given a point (X,Y), and evaluates ! a particular basis function for a particular quadratic ! triangular element at that point. ! ! REFBSP evaluates the linear basis functions in a reference ! triangle. ! ! REFQBF evaluates the quadratic basis functions in a reference ! triangle. ! ! RESID computes the residual. ! ! SETBAN computes the bandwidth of the matrix. ! ! SETGRD computes the numbering of nodes. ! Constructs and numbers the elements, and assigns them nodes. ! Calculates the number of unknowns. ! ! SETLIN records the horizontal velocity indices that lie along ! the profile line. ! ! SETQUD computes the location of the quadrature points. ! ! SETXY computes the location of the grid points. ! ! TRANS computes the transformation from the reference element ! to any isoparametric element. ! ! UBDRY produces the value of the inflow at a given point as ! determined by the value of the parameter A. ! ! UBUMP calculates the sensitivity on the bump. ! ! UVAL evaluates the horizontal and vertical velocities, and their ! derivatives with respect to X and Y, at a given point in an ! element. ! ! UVDUMP dump velocites and pressures for PLOT3D. ! ! XYDUMP dump X and Y locations for PLOT3D. ! ! ! The Navier Stokes equations: ! ! The primitive variable formulation involves horizontal velocity U, ! vertical velocity V, and pressure P. The equations are: ! ! U dUdx + V dUdy + dPdx - mu*(ddU/dxdx + ddU/dydy) = F1 ! ! U dVdx + V dVdy + dPdy - mu*(ddV/dxdx + ddV/dydy) = F2 ! ! dUdx + dVdy = 0 ! ! ! a double precision a(maxrow,neqn), contains the matrix ! in LINPACK general band storage mode. ! ! aprof double precision aprof, the value of the parameter for ! which the target profile was generated. ! ! area double precision area(nelemn), contains the area of the ! element. ! ! dcda double precision dcda(my), the sensitivities. ! ! f double precision f(neqn), is used as a right hand side ! vector in LINSYS, and is overwritten there by the ! new solution (which is then copied into g). ! ! fileg character fileg*30, the name of the file into ! which the DISPLAY graphics information will be stored. ! ! fileu character fileu*30, the name of the file into ! which the PLOT3D graphics velocity information will be ! stored. ! ! filex character filex*30, the name of the file into ! which the PLOT3D graphics coordinate information will be ! stored. ! ! g double precision g(neqn), the current solution vector, ! in which are stored pressures and velocities. ! ! ibump integer ibump. IBUMP determines where isoparametric elements ! will be used. ! ! 0, no isoparametric elements will be used. ! ! 1, isoparametric elements will be used only for the ! elements which directly impinge on the bump. ! ! 2, isoparametric elements will be used for all elements which ! are above the bump. ! ! 3, isoparametric elements will be used for all elements. ! ! indx integer indx(maxnp,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! insc integer insc(np), contains, for each node I, the ! index of the pressure P at that node, or 0. ! ! iounit integer iounit, the FORTRAN output unit to which graphics ! output is written. ! ! ipivot integer ipivot(neqn), pivoting space. ! ! isotri integer isotri(nelemn), contains, for each element, ! a 1 if the element is isoparametric, or 0 otherwise. ! ! iwrite integer iwrite, controls the amount of output printed. ! ! ivunit integer ivunit, the unit to which UV graphics is written. ! ! ixunit integer ixunit, the unit to which XY graphics is written. ! ! long logical long. ! ! .TRUE. if region is "long and thin", and ! .FALSE. if region is "tall and skinny". ! ! maxelm integer maxelm, the maximum number of elements. ! This is actually simply the first dimension of the array NODE. ! ! maxeqn integer maxeqn, the maximum number of equations, and functions. ! ! maxnew integer maxnew, the maximum number of Newton steps to take ! in one iteration. ! ! maxnp integer maxnp, the maximum number of nodes. ! This is actually simply the first dimension of the array INDX. ! ! maxrow integer maxrow, the first dimension of the matrix A. ! ! maxsec integer maxsec, the maximum number of secant steps to take. ! ! mx integer mx, the number of nodes in the X direction. ! ! my integer my, the number of nodes in the Y direction. ! ! nband integer nband, the bandwidth of the linear system. ! ! nbleft integer nbleft, the column at which the left corner of the ! bump lies. ! ! nbrite integer nbrite, the column at which the right corner of the ! bump lies. ! ! nelemn integer nelemn, the number of elements. ! ! neqn integer neqn, the number of equations, and functions. ! ! nlband integer nlband, the lower bandwidth of the matrix. ! ! nnodes integer nnodes, the number of nodes per element. ! ! node integer node(maxelm,nnodes), contains, for each element, ! the global node indices of the element nodes. ! ! nodex0 integer nodex0, the node whose Y coordinate is zero, ! and whose X coordinate is XPROF. This is the first ! node along the line where the velocity profile is ! measured. ! ! np integer np, the number of nodes. ! ! nquad integer nquad, the number of quadrature points per ! element. ! ! nrow integer nrow, the number of rows need to store the ! matrix A. ! ! numnew integer numnew, total number of Newton steps taken. ! ! nx integer nx, controls the spacing of nodes and elements in ! the X direction. There are 2*NX+1 nodes in the X ! direction. ! ! ny integer ny, controls the spacing of nodes and elements in ! the Y direction. There are 2*NY+1 nodes in the Y ! direction. ! ! phi double precision phi(nelemn,nquad,nnodes,3). Each entry of PHI ! contains the value of a quadratic basis function or its derivative, ! evaluated at a quadrature point. ! ! In particular, PHI(I,J,K,1) is the value of the quadratic basis ! function associated with local node K in element I, evaluatated ! at quadrature point J. ! ! PHI(I,J,K,2) is the X derivative of that same basis function, ! PHI(I,J,K,3) is the Y derivative of that same basis function. ! ! psi double precision psi(nelemn,nquad,nnodes). Each entry of PSI contains ! the value of a linear basis function evaluated at a ! quadrature point. ! ! PSI(I,J,K) is the value of the linear basis function associated ! with local node K in element I, evaluated at quadrature point J. ! ! res double precision res(neqn), holds the residual. ! ! reynld double precision reynld, the value of the Reynolds ! number. ! ! sens double precision sens(neqn), the sensitivities. ! ! tolnew double precision tolnew, convergence tolerance for Newton ! iteration. ! ! tolsec double precision tolsec, convergence tolerance for secant ! iteration. ! ! xbleft double precision xbleft, the left X coordinate of the bump. ! ! xbrite double precision xbrite, the right X coordinate of the bump. ! ! xc double precision xc(np), the X coordinates of the nodes. ! ! xlngth double precision xlngth, the length of the region. ! ! xm double precision xm(nelemn,nquad), contains the X coordinates ! of the quadrature points for each element. ! ! xprof double precision xprof, the X coordinate at which the ! profile is measured. This value should be a grid value! ! ! yc double precision yc(np), the Y coordinates of the nodes. ! ! ylngth double precision ylngth, the height of the region. ! ! ym double precision ym(nelemn,nquad), contains the Y coordinates ! of the quadrature points for each element. ! implicit double precision(a-h,o-z) external ddot integer, parameter :: maxnew = 4 integer, parameter :: maxsec = 10 integer, parameter :: nx = 21 integer, parameter :: ny = 7 ! ! This assignment should really read (maxrow = 27*min(nx,ny)) ! integer, parameter :: maxrow = 27*ny integer, parameter :: nelemn = 2*(nx-1)*(ny-1) integer, parameter :: mx = 2*nx-1 integer, parameter :: my = 2*ny-1 integer, parameter :: maxeqn = 2*mx*my+nx*ny integer, parameter :: np = mx*my integer, parameter :: nnodes = 6 integer, parameter :: nquad = 3 double precision a(maxrow,maxeqn) double precision anew double precision anext double precision aold double precision aprof double precision area(nelemn) double precision dcda(my) double precision ddot double precision f(maxeqn) character fileg*30 character fileu*30 character filex*30 double precision g(maxeqn) double precision gr(my,my) integer i integer ibump integer iline(my) integer indx(np,2) integer insc(np) integer iounit integer ipivot(maxeqn) integer isotri(nelemn) integer ivunit integer iwrite integer ixunit integer j logical long double precision minutes integer nband integer neqn integer nlband integer node(nelemn,nnodes) integer nrow double precision phi(nelemn,nquad,nnodes,3) double precision psi(nelemn,nquad,nnodes) double precision r(my) double precision res(maxeqn) double precision reynld double precision rtemp(my) double precision seconds double precision sens(maxeqn) real tarray(2) double precision temp double precision tolnew double precision tolsec double precision uprof(my) double precision xbleft double precision xbrite double precision xc(np) double precision xlngth double precision xm(nelemn,nquad) double precision xprof double precision yc(np) double precision ylngth double precision ym(nelemn,nquad) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BUMP_ORIGINAL' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' 8 September 1992' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Control problem for channel flow over a bump.' ! ! Set input data ! anew = 0.0D+00 anext = 0.3 aold = 0.0D+00 aprof = 0.25 fileg = 'display.dat' fileu = 'uv.dat' filex = 'xy.dat' ibump = 2 iounit = 2 ivunit = 4 iwrite = 1 ixunit = 3 numnew = 0 numsec = 0 reynld = 1.0D+00 rjpnew = 0.0D+00 rjpold = 0.0D+00 tolnew = 0.0001 tolsec = 0.0001 xbleft = 1.0D+00 xbrite = 3.0D+00 xlngth = 10.0D+00 xprof = 4.0D+00 ylngth = 3.0D+00 write(*,*)' ' write(*,*)'Bump generated with height of ',aprof write(*,*)' ' write(*,*)'NX = ',nx write(*,*)'NY = ',ny write(*,*)'Number of elements = ',nelemn write(*,*)'Reynolds number = ',reynld write(*,*)'Secant tolerance = ',tolsec write(*,*)'Newton tolerance = ',tolnew ! ! SETGRD constructs grid, numbers unknowns, calculates areas, ! and points for midpoint quadrature rule. ! call setgrd(ibump,indx,insc,isotri,iwrite,long,maxeqn,mx,my, & nelemn,neqn,nnodes,node,np,nx,ny,xbleft,xbrite,xlngth) ! ! Set X Y coordinates of grid points ! ypert = aprof call setxy(iwrite,long,mx,my,np,nx,ny,xc,xlngth,yc,ylngth,ypert) ! ! Set quadrature points ! call setqud(area,isotri,iwrite,nelemn,nnodes,node,np,nquad, & xc,xm,yc,ym) ! ! Set basis functions values at grid points ! call setbas(isotri,nelemn,nnodes,node,np,nquad,phi,psi,xc,xm,yc,ym) ! ! Find points on velocity profile sampling line ! call setlin(iline,indx,iwrite,long,mx,my,np,nx,ny,xlngth,xprof) ! ! Compute bandwidth ! call setban(indx,insc,maxrow,nband,nelemn,nlband,nnodes,node,np,nrow) ! ! Generate solution of Navier Stokes equation, using ! an initial estimate of g = 0. ! g(1:neqn) = 0.0D+00 call nstoke(a,area,f,g,indx,insc,ipivot,isotri,maxnew, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad, & nrow,numnew,phi,psi,reynld,tolnew,xc,xm,yc,ym) ! ! Compute the residual for the solution g. ! call resid(area,g,indx,insc,isotri,iwrite,nelemn,neqn,nnodes, & node,np,nquad,phi,psi,res,reynld,xc,xm,yc,ym) ! ! Output flow along line in array UPROF ! call getg(g,iline,my,neqn,uprof) if ( iwrite >= 1 ) then write(*,*)' ' write(*,*)'Velocity profile:' write(*,*)' ' write(*,'(5g14.6)') (uprof(i),i = 1,my) end if ! ! Calculate Gram matrix GR and vector R = line integral of UPROF*phi ! call gram(gr,iline,indx,iwrite,my,nelemn,nnodes,node, & np,r,uprof,xc,xprof,yc) ! ! Write data to file for display by DISPLAY ! if ( iwrite >= 1 ) then npara = 1 para = aprof call delete(fileg) open(unit = iounit,file=fileg,form='formatted',status='new') call gdump(g,indx,insc,iounit,isotri,long,nelemn,neqn, & nnodes,node,np,npara,nx,ny,para,reynld,rjpnew,xc,yc) end if ! ! XYDUMP creates a PLOT3D "XY" file ! if ( iwrite >= 1 ) then call delete(filex) open(unit = ixunit,file=filex,form='formatted',status='new') call xydump(ixunit,long,np,nx,ny,xc,yc) end if ! ! UVDUMP creates a PLOT3D "Q" file ! if ( iwrite >= 1 ) then call delete(fileu) open(unit = ivunit,file=fileu,form='formatted',status='new') call uvdump(g,indx,insc,ivunit,long,mx,my, & nelemn,neqn,nnodes,node,np,a,reynld,yc) end if ! ! Destroy information about true solution before beginning ! secant iteration. ! g(1:neqn) = 0.0D+00 ! ! Secant iteration loop: ! do iter = 1,maxsec write(*,*)' ' write(*,*)'Secant iteration ',iter numsec = numsec+1 ! ! Update the grid ! ypert = anew call setxy(iwrite,long,mx,my,np,nx,ny,xc,xlngth,yc,ylngth,ypert) ! ! Set quadrature points ! call setqud(area,isotri,iwrite,nelemn,nnodes,node,np,nquad,xc,xm,yc,ym) ! ! Set basis functions values at grid points ! call setbas(isotri,nelemn,nnodes,node,np,nquad,phi,psi,xc,xm,yc,ym) ! ! Solve for flow at new value of parameter ! call nstoke(a,area,f,g,indx,insc,ipivot,isotri,maxnew, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad, & nrow,numnew,phi,psi,reynld,tolnew,xc,xm,yc,ym) ! ! Get velocity profile along sampling line. ! call getg(g,iline,my,neqn,uprof) if ( iwrite >= 1 ) then write(*,*)' ' write(*,*)'Velocity profile:' write(*,*)' ' write(*,'(5g14.6)') (uprof(i),i = 1,my) end if ! ! Solve linear system for sensitivities. ! itype = -2 call linsys(a,area,sens,g,indx,insc,ipivot,isotri,itype, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad,nrow, & phi,psi,reynld,xc,xm,yc,ym) ! ! Get sensitivities DCDA along sampling line. ! call getg(sens,iline,my,neqn,dcda) if ( iwrite >= 2 ) then write(*,*)' ' write(*,*)'Sensitivities:' write(*,*)' ' write(*,'(5g14.6)')(dcda(i),i = 1,my) end if ! ! Evaluate J prime at current value of parameter where J is ! the functional to be minimized. ! ! JPRIME = 2.0D+00 * DCDA(I) * ( GR(I,J) * UPROF(J) - R(I) ) ! rjpnew = 0.0D+00 do i = 1, my temp = - r(i) do j = 1, my temp = temp + gr(i,j) * uprof(j) end do rjpnew = rjpnew + 2.0D+00 * dcda(i) * temp end do ! ! Write data to graphics files. ! if ( iwrite >= 1 ) then npara = 1 para = anew call gdump(g,indx,insc,iounit,isotri,long,nelemn,neqn, & nnodes,node,np,npara,nx,ny,para,reynld,rjpnew,xc,yc) call xydump(ixunit,long,np,nx,ny,xc,yc) call uvdump(g,indx,insc,ivunit,long,mx,my, & nelemn,neqn,nnodes,node,np,a,reynld,yc) end if write(*,*)' ' write(*,*)'Parameter = ',anew,' J prime=',rjpnew ! ! Update the estimate of the parameter using the secant step. ! if ( iter > 1 ) then anext = aold-rjpold*(anew-aold)/(rjpnew-rjpold) end if aold = anew anew = anext rjpold = rjpnew if ( anew /= 0.0D+00 ) then test = abs(anew-aold)/anew else test = 0.0D+00 end if write(*,*)'New value of parameter = ',anew write(*,*)'Convergence test = ',test if ( abs((anew-aold)) <= abs(anew)*tolsec.and.iter > 1 ) then write(*,*)'Secant iteration converged.' go to 40 end if end do write(*,*)'Secant iteration failed to converge.' 40 continue ! ! Produce total CPU time used. ! ! tarray(2) = 0.0D+00 ! tarray(1) = second() ! seconds = etime ( tarray ) minutes = seconds / 60.0D+00 write(*,*) ' ' write(*,*) 'Total execution time = ', seconds ,' seconds.' write(*,*) ' = ', minutes ,' minutes.' write(*,*)'Number of secant steps = ',numsec write(*,*)'Number of Newton steps = ',numnew ! ! Close graphics data files. ! if ( iwrite >= 1 ) then close(unit = iounit) close(unit = ixunit) close(unit = ivunit) end if write ( *, * ) ' ' write ( *, * ) 'BUMP_ORIGINAL' write ( *, * ) ' Normal end of execution.' stop end function bsp(it,iq,id,nelemn,nnodes,node,np,xc,xq,yc,yq) !*****************************************************************************80 ! !! BSP evaluates the linear basis function associated with pressure. ! implicit double precision(a-h,o-z) integer nelemn integer nnodes integer np integer node(nelemn,nnodes) double precision xc(np) double precision yc(np) iq1 = iq iq2 = mod(iq ,3)+1 iq3 = mod(iq+1,3)+1 i1 = node(it,iq1) i2 = node(it,iq2) i3 = node(it,iq3) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) if ( id == 1 ) then bsp = 1+((yc(i2)-yc(i3))*(xq-xc(i1))+(xc(i3)-xc(i2))*(yq-yc(i1)))/d else if ( id == 2 ) then bsp = (yc(i2)-yc(i3))/d else if ( id == 3 ) then bsp = (xc(i3)-xc(i2))/d else write(*,*)'BSP - Illegal value of ID = ',id stop end if return end subroutine dcopy ( n, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DCOPY copies a vector, x, to a vector, y. ! ! uses unrolled loops for increments equal to one. ! jack dongarra, linpack, 3/11/78. ! modified 12/3/93, array(1) declarations changed to array(*) ! double precision dx(*),dy(*) integer i,incx,incy,ix,iy,m,n if ( n <= 0)return if ( incx == 1.and.incy == 1)go to 20 ! ! code for unequal increments or equal increments ! not equal to 1 ! ix = 1 iy = 1 if ( incx < 0)ix = (-n+1)*incx + 1 if ( incy < 0)iy = (-n+1)*incy + 1 do i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy end do return ! ! code for both increments equal to 1 ! ! ! clean-up loop ! 20 m = mod(n,7) if ( m == 0 ) go to 40 dy(1:m) = dx(1:m) if ( n < 7 ) return 40 continue do i = m+1, n ,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) end do return end function ddot ( n, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DDOT forms the dot product of two vectors. ! ! uses unrolled loops for increments equal to one. ! jack dongarra, linpack, 3/11/78. ! modified 12/3/93, array(1) declarations changed to array(*) ! double precision ddot double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,n ddot = 0.0d+00 dtemp = 0.0d+00 if ( n <= 0)return if ( incx == 1.and.incy == 1)go to 20 ! ! code for unequal increments or equal increments ! not equal to 1 ! ix = 1 iy = 1 if ( incx < 0)ix = (-n+1)*incx + 1 if ( incy < 0)iy = (-n+1)*incy + 1 do i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy end do ddot = dtemp return ! ! code for both increments equal to 1 ! ! ! clean-up loop ! 20 m = mod(n,5) if ( m == 0 ) go to 40 do i = 1,m dtemp = dtemp + dx(i)*dy(i) end do if ( n < 5 ) go to 60 40 continue do i = m+1, n, 5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + & dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) end do 60 ddot = dtemp return end subroutine delete(filnam) !*****************************************************************************80 ! !! DELETE is used to dispose of old copies of files before opening a new copy ! character filnam*(*) open(unit = 99,file=filnam,status='old',err=10) close(unit = 99,status='delete',err=10) 10 continue ! return end subroutine gdump(f,indx,insc,iounit,isotri,long,nelemn,neqn, & nnodes,node,np,npara,nx,ny,para,reynld,rjpnew,xc,yc) !*****************************************************************************80 ! !! GDUMP writes graphics information to a file. ! ! In order to keep things simple, exactly one ! value, real or integer, is written per record. ! ! NPARA Input, INTEGER NPARA, the number of parameters. Fixed at 1 ! for now. ! ! PARA Input, DOUBLE PRECISION PARA(MAXPAR), the parameters. ! implicit double precision (a-h,o-z) integer nelemn integer neqn integer nnodes integer np double precision f(neqn) double precision fval integer i integer indx(np,2) integer insc(np) integer iounit integer iset integer isotri(nelemn) integer j logical long integer node(nelemn,nnodes) integer npara integer nx integer ny double precision para double precision reynld double precision rjpnew double precision ubdry double precision xc(np) double precision yc(np) save iset data iset /0/ iset = iset+1 write(iounit,*)long write (iounit,*) nelemn write (iounit,*) np write (iounit,*) npara write (iounit,*) nx write (iounit,*) ny ! ! Pressures ! do 10 i = 1, np j = insc(i) if (j <= 0) then fval = 0.0D+00 else fval = f(j) end if write (iounit,*) fval 10 continue ! ! Horizontal velocities, U ! do 20 i = 1, np j = indx(i,1) if (j == 0) then fval = 0.0D+00 else if (j < 0) then fval = ubdry(1,yc(i)) else fval = f(j) end if write (iounit,*) fval 20 continue ! ! Vertical velocities, V ! do i = 1, np j = indx(i,2) if (j <= 0) then fval = 0.0D+00 else fval = f(j) end if write (iounit,*) fval end do do 40 i = 1, np write (iounit,*) indx(i,1) write (iounit,*) indx(i,2) 40 continue do 50 i = 1, np write (iounit,*) insc(i) 50 continue do 60 i = 1, nelemn write (iounit,*) isotri(i) 60 continue do 80 i = 1, nelemn do 70 j = 1, 6 write (iounit,*) node(i,j) 70 continue 80 continue write (iounit,*) para write (iounit,*) reynld write (iounit,*) rjpnew do 100 i = 1, np write (iounit,*) xc(i) 100 continue do 110 i = 1, np write (iounit,*) yc(i) 110 continue write (*,*) 'GDUMP wrote data set ',iset,' to file.' return end subroutine getg(f,iline,my,neqn,u) !*****************************************************************************80 ! !! GETG outputs the array u along the line x = xprof given the ! entire field in vector f ! implicit double precision(a-h,o-z) integer my integer neqn double precision f(neqn) integer i integer iline(my) integer j double precision u(my) do i = 1,my u(i) = 0.0D+00 j = iline(i) if ( j <= 0 ) then u(i) = 0.0D+00 else u(i) = f(j) end if end do return end subroutine gram(gr,iline,indx,iwrite,my,nelemn,nnodes,node, & np,r,uprof,xc,xprof,yc) !*****************************************************************************80 ! !! GRAM computes and stores the gram matrix and the vector ! whose components are the line integral of UI*phi(j) ! 3 point Gauss quadrature rule is used for line integral ! implicit double precision(a-h,o-z) ! integer my integer nelemn integer nnodes integer np ! double precision gr(my,my) integer iline(my) integer indx(np,2) integer node(nelemn,nnodes) double precision r(my) double precision uprof(my) double precision wt(3) double precision xc(np) double precision xprof double precision yc(np) double precision yq(3) ! ! Input values for 3 point Gauss quadrature ! wt(1) = 5.0/9.0D+00 wt(2) = 8.0/9.0D+00 wt(3) = wt(1) yq(1) = -0.7745966692 yq(2) = 0.0D+00 yq(3) = -yq(1) ! ! zero arrays ! do i = 1,my r(i) = 0.0D+00 do 10 j = 1,my gr(i,j) = 0.0D+00 10 continue end do ! ! Compute line integral by looping over intervals along line ! using three point Gauss quadrature ! do 70 it = 1,nelemn ! ! Check to see if we are in a triangle with a side along line ! x = xprof ! k = node(it,1) kk = node(it,2) if ( abs(xc(k)-xprof) > 1.0e-4.or. & abs(xc(kk)-xprof) > 1.0e-4) go to 70 do 60 iquad = 1,3 bma2 = (yc(kk)-yc(k))/2.0D+00 ar = bma2*wt(iquad) x = xprof y = yc(k)+bma2*(yq(iquad)+1.0) ! ! Compute u internal at quadrature points ! uiqdpt = 0 do iq = 1,nnodes if ( iq == 1.or.iq.eq.2.or.iq.eq.4 ) then call qbf(x,y,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) ip = node(it,iq) iun = indx(ip,1) if ( iun > 0) then ii = igetl(iun,iline,my) uiqdpt = uiqdpt+bb*uprof(ii) else if ( iun == -1) then ubc = ubdry(1,yc(ip)) uiqdpt = uiqdpt+bb*ubc end if end if end do ! ! Only loop over nodes lying on line x = xprof ! do 50 iq = 1,nnodes if ( iq == 1.or.iq.eq.2.or.iq.eq.4 ) then ip = node(it,iq) call qbf(x,y,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) i = indx(ip,1) if ( i > 0 ) then ii = igetl(i,iline,my) r(ii) = r(ii)+bb*uiqdpt*ar do 40 iqq = 1,nnodes if ( iqq == 1.or.iqq.eq.2.or.iqq.eq.4 ) then ipp = node(it,iqq) call qbf(x,y,it,iqq,bbb,bbx,bby,nelemn,nnodes,node,np,xc,yc) j = indx(ipp,1) if ( j /= 0 ) then jj = igetl(j,iline,my) gr(ii,jj) = gr(ii,jj)+bb*bbb*ar end if end if 40 continue end if end if 50 continue 60 continue 70 continue ! if ( iwrite >= 3 ) then write(*,*)' ' write(*,*)'Gram matrix:' write(*,*)' ' do 90 i = 1,my do 80 j = 1,my write(*,*)i,j,gr(i,j) 80 continue 90 continue write(*,*)' ' write(*,*)'R vector:' write(*,*)' ' do 100 i = 1,my write(*,*)r(i) 100 continue end if ! return end function idamax ( n, dx, incx ) !*****************************************************************************80 ! !! IDAMAX finds the index of element having max. absolute value. ! ! jack dongarra, linpack, 3/11/78. ! modified 3/93 to return if incx <= 0. ! modified 12/3/93, array(1) declarations changed to array(*) ! integer idamax double precision dx(*),dmax integer i,incx,ix,n ! idamax = 0 if ( n < 1 .or. incx <= 0 ) return idamax = 1 if ( n == 1)return if ( incx == 1)go to 20 ! ! code for increment not equal to 1 ! ix = 1 dmax = dabs(dx(1)) ix = ix + incx do i = 2,n if ( dmax < dabs(dx(ix)) ) then idamax = i dmax = dabs(dx(ix)) end if ix = ix + incx end do return ! ! code for increment equal to 1 ! 20 continue dmax = abs ( dx(1) ) do i = 2,n if ( abs ( dx(i) ) > dmax ) then idamax = i dmax = abs ( dx(i) ) end if end do return end function igetl(i,iline,my) !*****************************************************************************80 ! !! IGETL gets the local unknown number along line x = xprof given the ! global unknown number ! implicit double precision(a-h,o-z) ! integer my ! integer i integer igetl integer iline(my) integer j do j = 1,my if ( iline(j) == i ) then igetl = j return end if end do write(*,*)'IGETL - Fatal error' write(*,*)'Unable to get local unknown number for ',i stop end subroutine linsys(a,area,f,g,indx,insc,ipivot,isotri,itype, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad,nrow, & phi,psi,reynld,xc,xm,yc,ym) !*****************************************************************************80 ! !! LINSYS solves the linearized Navier Stokes equation. ! ! ITYPE = -1 for solutions of the Navier Stokes equation. ! ITYPE = -2 for solutions of the sensitivity equation. ! implicit double precision(a-h,o-z) integer maxrow integer nelemn integer neqn integer nnodes integer np integer nquad double precision a(maxrow,neqn) double precision area(nelemn) double precision det double precision etax double precision etay double precision f(neqn) double precision g(neqn) integer i integer indx(np,2) integer insc(np) integer ioff integer ipivot(neqn) integer iq integer isotri(nelemn) integer itype integer j integer nband integer nlband integer node(nelemn,nnodes) integer nrow double precision phi(nelemn,nquad,nnodes,3) double precision psi(nelemn,nquad,nnodes) double precision reynld double precision un(2) double precision unx(2) double precision uny(2) double precision visc double precision xc(np) double precision xix double precision xiy double precision xm(nelemn,nquad) double precision xq double precision yc(np) double precision ym(nelemn,nquad) double precision yq ioff = nlband+nlband+1 visc = 1.0D+00 / reynld f(1:neqn) = 0.0D+00 a(1:nrow,1:neqn) = 0.0D+00 do 90 it = 1,nelemn ar = area(it)/3.0D+00 do 80 iquad = 1,nquad yq = ym(it,iquad) xq = xm(it,iquad) if ( isotri(it) == 1 ) then call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) ar = det*area(it)/3.0D+00 end if call uval(etax,etay,g,indx,isotri,it,nelemn,neqn, & nnodes,node,np,un,uny,unx,xc,xix,xiy,xq,yc,yq) ! ! For each basis function: ! do 70 iq = 1,nnodes ip = node(it,iq) bb = phi(it,iquad,iq,1) bx = phi(it,iquad,iq,2) by = phi(it,iquad,iq,3) bbl = psi(it,iquad,iq) ihor = indx(ip,1) iver = indx(ip,2) iprs = insc(ip) if ( ihor > 0 ) then f(ihor) = f(ihor)+ar*bb*(un(1)*unx(1)+un(2)*uny(1)) end if if ( iver > 0 ) then f(iver) = f(iver)+ar*bb*(un(1)*unx(2)+un(2)*uny(2)) end if ! ! For another basis function, ! do 50 iqq = 1,nnodes ipp = node(it,iqq) bbb = phi(it,iquad,iqq,1) bbx = phi(it,iquad,iqq,2) bby = phi(it,iquad,iqq,3) bbbl = psi(it,iquad,iqq) ju = indx(ipp,1) jv = indx(ipp,2) jp = insc(ipp) ! ! Horizontal velocity variable ! if ( ju > 0 ) then if ( ihor > 0 ) then iuse = ihor-ju+ioff a(iuse,ju) = a(iuse,ju)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))) end if if ( iver > 0 ) then iuse = iver-ju+ioff a(iuse,ju) = a(iuse,ju)+ar*bb*bbb*unx(2) end if if ( iprs > 0 ) then iuse = iprs-ju+ioff a(iuse,ju) = a(iuse,ju)+ar*bbx*bbl end if else if ( ju == itype ) then if ( ju == -1 ) then ubc = ubdry(1,yc(ipp)) else if ( ju == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,1,nelemn, & neqn,nnodes,node,np,xc,yc) end if if ( ihor > 0 ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))) f(ihor) = f(ihor)-ubc*aij end if if ( iver > 0 ) then aij = ar*bb*bbb*unx(2) f(iver) = f(iver)-ubc*aij end if if ( iprs > 0 ) then aij = ar*bbx*bbl f(iprs) = f(iprs)-ubc*aij end if end if ! ! Vertical velocity variable ! if ( jv > 0 ) then if ( ihor > 0 ) then iuse = ihor-jv+ioff a(iuse,jv) = a(iuse,jv)+ar*bb*bbb*uny(1) end if if ( iver > 0 ) then iuse = iver-jv+ioff a(iuse,jv) = a(iuse,jv)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))) end if if ( iprs > 0 ) then iuse = iprs-jv+ioff a(iuse,jv) = a(iuse,jv)+ar*bby*bbl end if else if ( jv == itype ) then if ( jv == -1 ) then ubc = ubdry(2,yc(ipp)) else if ( jv == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,2,nelemn, & neqn,nnodes,node,np,xc,yc) end if if ( ihor > 0 ) then aij = ar*bb*bbb*uny(1) f(ihor) = f(ihor)-ubc*aij end if if ( iver > 0 ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))) f(iver) = f(iver)-ubc*aij end if if ( iprs > 0 ) then aij = ar*bby*bbl f(iprs) = f(iprs)-ubc*aij end if end if ! ! Pressure variable ! if ( jp > 0 ) then if ( ihor > 0 ) then iuse = ihor-jp+ioff a(iuse,jp) = a(iuse,jp)-ar*bx*bbbl end if if ( iver > 0 ) then iuse = iver-jp+ioff a(iuse,jp) = a(iuse,jp)-ar*by*bbbl end if end if 50 continue 70 continue 80 continue 90 continue ! ! The last equation is "reset" to require that the last pressure ! be zero. ! f(neqn) = 0.0D+00 do j = neqn-nlband,neqn-1 i = neqn-j+ioff a(i,j) = 0.0D+00 end do a(nband,neqn) = 1.0D+00 ! ! Factor the matrix ! call dgbfa ( a, maxrow, neqn, nlband, nlband, ipivot, info ) if (info /= 0) then write (*,*) ' ' write (*,*) 'LINSYS - fatal error!' write (*,*) 'DGBFA returns INFO = ',info stop end if ! ! Solve the linear system ! call dgbsl ( a, maxrow, neqn, nlband, nlband, ipivot, f, 0 ) if (info /= 0) then write (*,*) ' ' write (*,*) 'LINSYS - fatal error!' write (*,*) 'DGBTRS returns INFO = ',info stop end if return end subroutine nstoke(a,area,f,g,indx,insc,ipivot,isotri,maxnew, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad, & nrow,numnew,phi,psi,reynld,tolnew,xc,xm,yc,ym) !*****************************************************************************80 ! !! NSTOKE solves the Navier Stokes equation using Taylor-Hood elements. ! implicit double precision(a-h,o-z) integer maxrow integer nelemn integer neqn integer nnodes integer np integer nquad double precision a(maxrow,neqn) double precision area(nelemn) double precision f(neqn) double precision g(neqn) integer indx(np,2) integer insc(np) integer ipivot(neqn) integer isotri(nelemn) integer maxnew integer nband integer nlband integer node(nelemn,nnodes) integer nrow integer numnew double precision phi(nelemn,nquad,nnodes,3) double precision psi(nelemn,nquad,nnodes) double precision reynld double precision tolnew double precision xc(np) double precision xm(nelemn,nquad) double precision yc(np) double precision ym(nelemn,nquad) ! ! G contains an initial estimate of the solution. ! do 20 iter = 1,maxnew numnew = numnew+1 itype = -1 call linsys(a,area,f,g,indx,insc,ipivot,isotri,itype, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad,nrow, & phi,psi,reynld,xc,xm,yc,ym) ! ! Check for convergence ! do i = 1,neqn g(i) = g(i)-f(i) end do diff = abs(g(idamax(neqn,g,1))) write(*,*)'NSTOKE: Iteration ',iter,' MaxNorm(diff) = ',diff call dcopy(neqn,f,1,g,1) if ( diff <= tolnew ) then write(*,*)'NSTOKE converged.' return end if if ( iter == maxnew ) then write(*,*)'NSTOKE failed!' stop end if 20 continue return end subroutine pval (g,insc,long,mx,my,nelemn,neqn,nnodes,node,np,press) !*****************************************************************************80 ! !! PVAL computes a table of pressures at all the nodes. ! ! ! This is needed for use with PLOT3D. ! ! This routine does not handle general isoparametric elements ! accurately. ! implicit double precision (a-h,o-z) integer mx integer my integer nelemn integer neqn integer nnodes integer np double precision g(neqn) integer i integer insc(np) integer ip integer iq integer it integer ivar integer j logical long integer node(nelemn,nnodes) double precision press(mx,my) do 20 i = 1,mx do 10 j = 1,my press(i,j) = 0.0D+00 10 continue 20 continue ! ! Read the pressures where they are computed. ! These are "(odd, odd)" points. ! do 40 it = 1,nelemn do 30 iq = 1,3 ip = node(it,iq) ivar = insc(ip) if ( long ) then i = ((ip-1)/my)+1 j = mod(ip-1,my)+1 else i = mod(ip-1,mx)+1 j = ((ip-1)/mx)+1 end if if ( ivar > 0 ) then press(i,j) = g(ivar) end if 30 continue 40 continue ! ! Interpolate the pressures at points (even, odd) and (odd, even). ! do i = 2,mx-1,2 do j = 1,my,2 press(i,j) = 0.5*(press(i-1,j)+press(i+1,j)) end do end do do j = 2,my-1,2 do i = 1,mx,2 press(i,j) = 0.5*(press(i,j-1)+press(i,j+1)) end do end do ! ! Interpolate the pressures at points (even,even). ! do j = 2,my-1,2 do i = 2,mx-1,2 press(i,j) = 0.5*(press(i-1,j-1)+press(i+1,j+1)) end do end do return end subroutine qbf(xq,yq,it,in,bb,bx,by,nelemn,nnodes,node,np,xc,yc) !*****************************************************************************80 ! !! QBF evaluates the quadratic basis functions. ! implicit double precision(a-h,o-z) integer nelemn integer nnodes integer np double precision bb double precision bx double precision by integer node(nelemn,nnodes) double precision xc(np) double precision xq double precision yc(np) double precision yq if ( in <= 3 ) then in1 = in in2 = mod(in,3)+1 in3 = mod(in+1,3)+1 i1 = node(it,in1) i2 = node(it,in2) i3 = node(it,in3) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) t = 1+((yc(i2)-yc(i3))*(xq-xc(i1))+(xc(i3)-xc(i2))*(yq-yc(i1)))/d bb = t*(2*t-1) bx = (yc(i2)-yc(i3))*(4*t-1)/d by = (xc(i3)-xc(i2))*(4*t-1)/d else inn = in-3 in1 = inn in2 = mod(inn,3)+1 in3 = mod(inn+1,3)+1 i1 = node(it,in1) i2 = node(it,in2) i3 = node(it,in3) j1 = i2 j2 = i3 j3 = i1 d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) c = (xc(j2)-xc(j1))*(yc(j3)-yc(j1))-(xc(j3)-xc(j1))*(yc(j2)-yc(j1)) t = 1+((yc(i2)-yc(i3))*(xq-xc(i1))+(xc(i3)-xc(i2))*(yq-yc(i1)))/d s = 1+((yc(j2)-yc(j3))*(xq-xc(j1))+(xc(j3)-xc(j2))*(yq-yc(j1)))/c bb = 4*s*t bx = 4*(t*(yc(j2)-yc(j3))/c+s*(yc(i2)-yc(i3))/d) by = 4*(t*(xc(j3)-xc(j2))/c+s*(xc(i3)-xc(i2))/d) end if return end function refbsp(xq,yq,iq) !*****************************************************************************80 ! !! REFBSP evaluates the linear basis functions in a reference triangle. ! implicit double precision(a-h,o-z) integer iq double precision refbsp double precision xq double precision yq if ( iq == 1 ) then refbsp = 1.0-xq else if ( iq == 2 ) then refbsp = yq else if ( iq == 3 ) then refbsp = xq-yq end if return end subroutine refqbf(x,y,in,bb,bx,by,etax,etay,xix,xiy) ! !*****************************************************************************80 ! !! REFQBF evaluates the quadratic basis functions on reference triangle ! implicit double precision(a-h,o-z) double precision bb double precision bx double precision by double precision etax double precision etay integer in double precision tbx double precision tby double precision x double precision xix double precision xiy double precision y ! if ( in == 1) then bb = 1-3*x+2*x*x tbx = -3+4*x tby = 0 else if ( in == 2) then bb = -y+2*y*y tbx = 0 tby = -1+4*y else if (in == 3) then bb = -x+2*x*x+y-4*x*y+2*y*y tbx = -1+4*x-4*y tby = 1-4*x+4*y else if ( in == 4) then bb = 4*y-4*x*y tbx = -4*y tby = 4-4*x else if ( in == 5) then bb = 4*x*y-4*y*y tbx = 4*y tby = 4*x-8*y else if ( in == 6) then bb = 4*x-4*x*x-4*y+4*x*y tbx = 4-8*x+4*y tby = -4+4*x else write(*,*)'REFQBF - Illegal value of IN = ',in stop end if bx = tbx*xix+tby*etax by = tbx*xiy+tby*etay return end subroutine resid(area,g,indx,insc,isotri,iwrite,nelemn,neqn, & nnodes,node,np,nquad,phi,psi,res,reynld,xc,xm,yc,ym) !*****************************************************************************80 ! !! RESID computes the residual. ! implicit double precision(a-h,o-z) integer nelemn integer neqn integer nnodes integer np integer nquad double precision area(nelemn) double precision det double precision etax double precision etay double precision g(neqn) integer i integer indx(np,2) integer insc(np) integer isotri(nelemn) integer it integer itype integer iwrite integer node(nelemn,nnodes) double precision phi(nelemn,nquad,nnodes,3) double precision psi(nelemn,nquad,nnodes) double precision res(neqn) double precision reynld double precision un(2) double precision unx(2) double precision uny(2) double precision visc double precision xc(np) double precision xix double precision xiy double precision xm(nelemn,nquad) double precision xq double precision yc(np) double precision ym(nelemn,nquad) double precision yq itype = -1 visc = 1.0/reynld res(1:neqn) = 0.0D+00 do 90 it = 1,nelemn ar = area(it)/3.0D+00 do 80 iquad = 1,nquad yq = ym(it,iquad) xq = xm(it,iquad) if ( isotri(it) == 1 ) then call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) ar = det*area(it)/3.0D+00 end if call uval(etax,etay,g,indx,isotri,it,nelemn,neqn, & nnodes,node,np,un,uny,unx,xc,xix,xiy,xq,yc,yq) ! ! For each basis function: ! do 70 iq = 1,nnodes ip = node(it,iq) bb = phi(it,iquad,iq,1) bx = phi(it,iquad,iq,2) by = phi(it,iquad,iq,3) bbl = psi(it,iquad,iq) iprs = insc(ip) ihor = indx(ip,1) iver = indx(ip,2) ! if ( ihor > 0 ) then res(ihor) = res(ihor)-ar*bb*(un(1)*unx(1)+un(2)*uny(1)) end if if ( iver > 0 ) then res(iver) = res(iver)-ar*bb*(un(1)*unx(2)+un(2)*uny(2)) end if ! ! For another basis function, ! do 50 iqq = 1,nnodes ipp = node(it,iqq) bbb = phi(it,iquad,iqq,1) bbx = phi(it,iquad,iqq,2) bby = phi(it,iquad,iqq,3) bbbl = psi(it,iquad,iqq) ju = indx(ipp,1) jv = indx(ipp,2) jp = insc(ipp) ! ! Horizontal velocity variable ! if ( ju > 0 ) then if ( ihor > 0 ) then res(ihor) = res(ihor)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2)))*g(ju) end if if ( iver > 0 ) then res(iver) = res(iver)+ar*bb*bbb*unx(2)*g(ju) end if if ( iprs > 0 ) then res(iprs) = res(iprs)+ar*bbx*bbl*g(ju) end if else if ( ju == itype ) then if ( ju == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,1,nelemn, & neqn,nnodes,node,np,xc,yc) else if ( ju == -1 ) then ubc = ubdry(1,yc(ipp)) end if if ( ihor > 0 ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))) res(ihor) = res(ihor)+ubc*aij end if if ( iver > 0 ) then aij = ar*bb*bbb*unx(2) res(iver) = res(iver)+ubc*aij end if if ( iprs > 0 ) then aij = ar*bbx*bbl res(iprs) = res(iprs)+ubc*aij end if end if ! ! Vertical velocity variable ! if ( jv > 0 ) then if ( ihor > 0 ) then res(ihor) = res(ihor)+ar*bb*bbb*uny(1)*g(jv) end if if ( iver > 0 ) then res(iver) = res(iver)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1)))*g(jv) end if if ( iprs > 0 ) then res(iprs) = res(iprs)+ar*bby*bbl*g(jv) end if else if ( jv == itype ) then if ( jv == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,2,nelemn, & neqn,nnodes,node,np,xc,yc) else if ( jv == -1 ) then ubc = ubdry(2,yc(ipp)) end if if ( ihor > 0 ) then aij = ar*bb*bbb*uny(1) res(ihor) = res(ihor)+ubc*aij end if if ( iver > 0 ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))) res(iver) = res(iver)+ubc*aij end if if ( iprs > 0 ) then aij = ar*bby*bbl res(iprs) = res(iprs)+ubc*aij end if end if ! ! Pressure variable ! if ( jp > 0 ) then if ( ihor > 0 ) then res(ihor) = res(ihor)-ar*bx*bbbl*g(jp) end if if ( iver > 0 ) then res(iver) = res(iver)-ar*by*bbbl*g(jp) end if end if 50 continue 70 continue 80 continue 90 continue ! ! The last equation is "reset" to require that the last pressure ! be zero. ! res(neqn) = g(neqn) rmax = 0.0D+00 imax = 0 ibad = 0 do i = 1,neqn test = abs(res(i)) if ( test > rmax ) then rmax = test imax = i end if if ( test > 1.0e-3 ) then ibad = ibad+1 end if end do if ( iwrite >= 1 ) then write(*,*)' ' write(*,*)'RESIDUAL INFORMATION:' write(*,*)' ' write(*,*)'Worst residual is number ',IMAX write(*,*)'of magnitude ',RMAX write(*,*)' ' write(*,*)'Number of "bad" residuals is ',IBAD,' out of ',NEQN write(*,*)' ' end if if ( iwrite >= 2 ) then write(*,*)'Raw residuals:' write(*,*)' ' i = 0 do 120 j = 1,np if ( indx(j,1) > 0 ) then i = i+1 if ( abs(res(i)) <= 1.0e-3 ) then write(*,'(1x,a1,2i5,g14.6)')'U',i,j,res(i) else write(*,'(a1,a1,2i5,g14.6)')'*','U',i,j,res(i) end if end if if ( indx(j,2) > 0 ) then i = i+1 if ( abs(res(i)) <= 1.0e-3 ) then write(*,'(1x,a1,2i5,g14.6)')'V',i,j,res(i) else write(*,'(a1,a1,2i5,g14.6)')'*','V',i,j,res(i) end if end if if ( insc(j) > 0 ) then i = i+1 if ( abs(res(i)) <= 1.0e-3 ) then write(*,'(1x,a1,2i5,g14.6)')'P',i,j,res(i) else write(*,'(a1,a1,2i5,g14.6)')'*','P',i,j,res(i) end if end if 120 continue end if return end subroutine setban(indx,insc,maxrow,nband,nelemn,nlband,nnodes, & node,np,nrow) !*****************************************************************************80 ! !! SETBAN computes the half band width. ! implicit double precision (a-h,o-z) integer nelemn integer nnodes integer np integer i integer indx(np,2) integer insc(np) integer ip integer ipp integer iq integer iqq integer it integer iuk integer iukk integer j integer maxrow integer nband integer nlband integer node(nelemn,nnodes) integer nrow nlband = 0 do 50 it = 1, nelemn do 40 iq = 1, nnodes ip = node(it,iq) do 30 iuk = 1, 3 if (iuk == 3) then i = insc(ip) else i = indx(ip,iuk) end if if (i > 0 ) then do 20 iqq = 1, nnodes ipp = node(it,iqq) do iukk = 1, 3 if (iukk == 3) then j = insc(ipp) else j = indx(ipp,iukk) end if if ( j > 0 ) then nlband = max(nlband,j-i) end if end do 20 continue end if 30 continue 40 continue 50 continue nband = nlband+nlband+1 nrow = nlband+nlband+nlband+1 ! write(*,*)' ' write(*,*)'SETBAN:' write(*,*)' ' write(*,*)' Lower bandwidth = ',nlband write(*,*)' Total bandwidth = ',nband write(*,*)' Required matrix rows = ',nrow if ( nrow > maxrow ) then write(*,*)'SETBAN - NROW is too large!' write(*,*)'The maximum allowed is ',maxrow stop end if return end subroutine setbas(isotri,nelemn,nnodes,node,np,nquad,phi,psi, & xc,xm,yc,ym) !*****************************************************************************80 ! !! SETBAS computes the value of the basis functions at each ! integration point. ! implicit double precision (a-h,o-z) integer nelemn integer nnodes integer np integer nquad double precision bb double precision bx double precision by double precision det double precision etax double precision etay integer iq integer isotri(nelemn) integer it integer j integer node(nelemn,nnodes) double precision phi(nelemn,nquad,nnodes,3) double precision psi(nelemn,nquad,nnodes) double precision xc(np) double precision xix double precision xiy double precision xm(nelemn,nquad) double precision xq double precision yc(np) double precision ym(nelemn,nquad) double precision yq do it = 1,nelemn do j = 1,nquad xq = xm(it,j) yq = ym(it,j) call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) do iq = 1,nnodes if ( isotri(it) == 0 ) then psi(it,j,iq) = bsp(it,iq,1,nelemn,nnodes,node,np,xc,xq,yc,yq) call qbf(xq,yq,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) else call refqbf(xq,yq,iq,bb,bx,by,etax,etay,xix,xiy) psi(it,j,iq) = refbsp(xq,yq,iq) end if phi(it,j,iq,1) = bb phi(it,j,iq,2) = bx phi(it,j,iq,3) = by end do end do end do return end subroutine setgrd(ibump,indx,insc,isotri,iwrite,long,maxeqn,mx,my, & nelemn,neqn,nnodes,node,np,nx,ny,xbleft,xbrite,xlngth) !*****************************************************************************80 ! !! SETGRD sets up the grid for the problem assuming we are using quadratics ! for the velocity and linears for the pressure ! implicit double precision(a-h,o-z) integer my integer nelemn integer nnodes integer np integer ibump integer indx(np,2) integer insc(np) integer isotri(nelemn) integer iwrite logical long integer maxeqn integer mx integer nbleft integer nbrite integer neqn integer node(nelemn,nnodes) integer nx integer ny double precision xbleft double precision xbrite double precision xlngth write(*,*)' ' write(*,*)'SETGRD:' write(*,*)' ' ! ! Determine whether region is long or skinny. This will determine ! how we number the nodes and elements. ! if ( nx > ny ) then long = .true. write(*,*)'Using vertical ordering.' else long = .false. write(*,*)'Using horizontal ordering.' end if ! ! Report on isoparametric strategy: ! if ( ibump == 0 ) then write(*,*)'No isoparametric elements will be used.' else if ( ibump == 1 ) then write(*,*)'Isoparametric elements directly on bump.' else if ( ibump == 2 ) then write(*,*)'All elements above bump are isoparametric.' else if ( ibump == 3 ) then write(*,*)'All elements are isoparametric.' else write(*,*)'Unexpected value of IBUMP = ',ibump stop end if ! ! Compute node locations of bump corners based on X coordinates ! nbleft = nint(xbleft*(mx-1)/xlngth)+1 nbrite = nint(xbrite*(mx-1)/xlngth)+1 write(*,*)'Bump extends from ',xbleft,' at node ',nbleft write(*,*)' to ',xbrite,' at node ',nbrite ! ! Assign nodes to elements. ! neqn = 0 ielemn = 0 do 10 ip = 1,np if ( long ) then ic = ((ip-1)/my)+1 jc = mod((ip-1),my)+1 else ic = mod((ip-1),mx)+1 jc = ((ip-1)/mx)+1 end if icnt = mod(ic,2) jcnt = mod(jc,2) ! ! If both the row count and the column count are odd, ! and we're not in the last row or top column, ! then we can define two new triangular elements based at the node. ! ! For horizontal ordering, ! given the following arrangement of nodes, for instance: ! ! 21 22 23 24 25 ! 16 17 18 19 20 ! 11 12 13 14 15 ! 06 07 08 09 10 ! 01 02 03 04 05 ! ! when we arrive at node 13, we will define ! ! element 7: (13, 23, 25, 18, 24, 19) ! element 8: (13, 25, 15, 19, 20, 14) ! ! ! For vertical ordering, ! given the following arrangement of nodes, for instance: ! ! 05 10 15 20 25 ! 04 09 14 19 24 ! 03 08 13 18 23 ! 02 07 12 17 22 ! 01 06 11 16 21 ! ! when we arrive at node 13, we will define ! ! element 7: (13, 25, 23, 19, 24, 18) ! element 8: (13, 15, 25, 14, 20, 19) ! if ( (icnt == 1.and.jcnt == 1).and.(ic /= mx).and.(jc /= my) ) then if ( long ) then ip1 = ip+my ip2 = ip+my+my ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip+2 node(ielemn,3) = ip2+2 node(ielemn,4) = ip+1 node(ielemn,5) = ip1+2 node(ielemn,6) = ip1+1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then isotri(ielemn) = 0 else if ( ibump == 2 ) then if ( ic >= nbleft.and.ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip2+2 node(ielemn,3) = ip2 node(ielemn,4) = ip1+1 node(ielemn,5) = ip2+1 node(ielemn,6) = ip1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then if ( jc == 1.and.ic >= nbleft.and.ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else if ( ibump == 2 ) then if ( ic >= nbleft.and.ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if else ip1 = ip+mx ip2 = ip+mx+mx ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip2 node(ielemn,3) = ip2+2 node(ielemn,4) = ip1 node(ielemn,5) = ip2+1 node(ielemn,6) = ip1+1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then isotri(ielemn) = 0 else if ( ibump == 2 ) then if ( ic >= nbleft.and.ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip2+2 node(ielemn,3) = ip+2 node(ielemn,4) = ip1+1 node(ielemn,5) = ip1+2 node(ielemn,6) = ip+1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then if ( jc == 1.and.ic >= nbleft.and.ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else if ( ibump == 2 ) then if ( ic >= nbleft.and.ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if end if end if ! ! Left hand boundary, horizontal and vertical velocities specified. ! if ( ic == 1.and.1 < jc.and.jc.lt.my ) then indx(ip,1) = -1 indx(ip,2) = -1 ! ! Right hand boundary, horizontal velocities unknown, vertical ! velocities specified ! else if ( ic == mx.and.1 < jc.and.jc.lt.my) then neqn = neqn+1 indx(ip,1) = neqn indx(ip,2) = 0 ! ! Lower boundary, with isoperimetric triangle ! else if ( jc == 1.and.isotri(ielemn).eq.1 ) then indx(ip,1) = -2 indx(ip,2) = -2 ! ! Otherwise, just a wall ! else if ( ic == 1.or.ic.eq.mx.or.jc.eq.1.or.jc.eq.my ) then indx(ip,1) = 0 indx(ip,2) = 0 ! ! Otherwise, a normal interior node where both velocities are unknown ! else neqn = neqn+2 indx(ip,1) = neqn-1 indx(ip,2) = neqn end if if ( jcnt == 1.and.icnt.eq.1 ) then neqn = neqn+1 insc(ip) = neqn else insc(ip) = 0 end if 10 continue if (iwrite >= 1 ) then write(*,*)' ' write(*,*)' I INDX 1, INDX 2, INSC' write(*,*)' ' do i = 1,np write(*,'(4i5)')i,indx(i,1),indx(i,2),insc(i) end do write(*,*)' ' write(*,*)'Isoparametric triangles:' write(*,*)' ' do i = 1,nelemn if ( isotri(i) == 1)write(*,*)i end do write(*,*)' ' write(*,*)' IT NODE(IT,*)' write(*,*)' ' do it = 1,nelemn write(*,'(7i6)') it,(node(it,i),i = 1,6) end do end if write(*,*)' ' write(*,*)'SETGRD: Number of unknowns = ',neqn if ( neqn > maxeqn ) then write(*,*)'SETGRD - Too many unknowns!' write(*,*)'The maximum allowed is MAXEQN = ',maxeqn write(*,*)'This problem requires NEQN = ',neqn stop end if return end subroutine setlin(iline,indx,iwrite,long,mx,my,np,nx,ny,xlngth,xprof) !*****************************************************************************80 ! !! SETLIN determines unknown numbers along line x = xprof. ! implicit double precision (a-h,o-z) integer my integer np integer i integer iline(my) integer indx(np,2) integer iwrite logical long integer mx integer nodex0 integer nx integer ny double precision xlngth double precision xprof ! ! Determine the number of a node on the profile line ! itemp = nint((2.0*(nx-1)*xprof)/xlngth) if ( long ) then nodex0 = itemp*(2*ny-1)+1 else nodex0 = itemp+1 end if write(*,*)' ' write(*,*)'SETLIN:' write(*,*)' ' write(*,*)' Profile generated at X = ',xprof write(*,*)' which is above node = ',nodex0 do i = 1,my if ( long ) then ip = nodex0+(i-1) else ip = nodex0+mx*(i-1) end if iline(i) = indx(ip,1) end do if ( iwrite >= 1 ) then write(*,*)' ' write(*,*)' Unknown numbers along line:' write(*,*)' ' write(*,'(5i5)') (iline(i),i = 1,my) end if return end subroutine setqud(area,isotri,iwrite,nelemn,nnodes,node,np,nquad, & xc,xm,yc,ym) !*****************************************************************************80 ! !! SETQUD sets midpoint quadrature rule information. ! implicit double precision (a-h,o-z) integer nelemn integer nnodes integer np integer nquad double precision area(nelemn) integer isotri(nelemn) integer iwrite integer node(nelemn,nnodes) double precision xc(np) double precision xm(nelemn,nquad) double precision yc(np) double precision ym(nelemn,nquad) do it = 1,nelemn ip1 = node(it,1) ip2 = node(it,2) ip3 = node(it,3) x1 = xc(ip1) x2 = xc(ip2) x3 = xc(ip3) y1 = yc(ip1) y2 = yc(ip2) y3 = yc(ip3) if ( isotri(it) == 0 ) then xm(it,1) = 0.5*(x1+x2) xm(it,2) = 0.5*(x2+x3) xm(it,3) = 0.5*(x3+x1) ym(it,1) = 0.5*(y1+y2) ym(it,2) = 0.5*(y2+y3) ym(it,3) = 0.5*(y3+y1) area(it) = 0.5*abs((y1+y2)*(x2-x1)+(y2+y3)*(x3-x2)+(y3+y1)*(x1-x3)) else xm(it,1) = 0.5 ym(it,1) = 0.5 xm(it,2) = 1.0D+00 ym(it,2) = 0.5 xm(it,3) = 0.5 ym(it,3) = 0.0D+00 area(it) = 0.5 end if end do if ( iwrite >= 3 ) then write(*,*)' ' write(*,*)'SETQUD: Element Areas and Quadrature points:' write(*,*)' ' do i = 1,nelemn write(*,*)i,area(i) do j = 1,nquad write(*,*)i,j,xm(i,j),ym(i,j) end do end do end if return end subroutine setxy(iwrite,long,mx,my,np,nx,ny,xc,xlngth,yc,ylngth,ypert) !*****************************************************************************80 ! !! SETXY sets the grid coordinates based on the value of the parameter. ! implicit double precision(a-h,o-z) integer my integer np integer iwrite logical long integer mx integer nx integer ny double precision xc(np) double precision xlngth double precision yc(np) double precision ylngth double precision ypert ! ! set function and endpoints for bottom bump boundary ! ybot(x) = -ypert*(x-3.0)*(x-1.0) do ip = 1,np if ( long ) then ic = ((ip-1)/my)+1 jc = mod((ip-1),my)+1 else ic = mod((ip-1),mx)+1 jc = ((ip-1)/mx)+1 end if xc(ip) = (ic-1)*xlngth/(2*nx-2) if ( ybot(xc(ip)) > 0.0D+00 ) then ylo = ybot(xc(ip)) else ylo = 0.0D+00 end if yc(ip) = ( (my-jc)*ylo + (jc-1)*ylngth ) / (2*ny-2) end do if (iwrite >= 2 ) then write(*,*)' ' write(*,*)'SETGRD:' write(*,*)' ' write(*,*)' I XC YC' write(*,*)' ' do i = 1,np write(*,'(i5,2f12.5)')i,xc(i),yc(i) end do end if return end subroutine trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) !*****************************************************************************80* ! !! TRANS calculates the transformation which maps reference element ! into isoparametric element. ! Outputs the determinant and partial derivatives for transformation ! implicit double precision(a-h,o-z) integer nelemn integer nnodes integer np double precision det double precision etax double precision etay integer it integer node(nelemn,nnodes) double precision xc(np) double precision xix double precision xiy double precision xq double precision yc(np) double precision yq i1 = node(it,1) i2 = node(it,2) i3 = node(it,3) i4 = node(it,4) i5 = node(it,5) i6 = node(it,6) x1 = xc(i1) y1 = yc(i1) x2 = xc(i2) y2 = yc(i2) x3 = xc(i3) y3 = yc(i3) x4 = xc(i4) y4 = yc(i4) x5 = xc(i5) y5 = yc(i5) x6 = xc(i6) y6 = yc(i6) ! ! Set coefficients in transformation ! a1 = 2*x3-4*x6+2*x1 b1 = -4*x3-4*x4+4*x5+4*x6 c1 = 2*x2+2*x3-4*x5 d1 = -3*x1-x3+4*x6 e1 = -x2+x3+4*x4-4*x6 a2 = 2*y3-4*y6+2*y1 b2 = -4*y3-4*y4+4*y5+4*y6 c2 = 2*y2+2*y3-4*y5 d2 = -3*y1-y3+4*y6 e2 = -y2+y3+4*y4-4*y6 ! ! Compute partial derivatives at point (xq,yq) ! etax = 2*a2*xq+b2*yq+d2 etay = 2*a1*xq+b1*yq+d1 xix = b2*xq+2*c2*yq+e2 xiy = b1*xq+2*c1*yq+e1 ! ! Compute determinant of transformation evaluated at point (xq,yq) ! det = (2*a1*b2-2*a2*b1)*xq*xq+(4*a1*c2-4*a2*c1)*xq*yq & +(2*b1*c2-2*b2*c1)*yq*yq+(2*a1*e2+b2*d1-b1*d2-2*a2*e1)*xq & +(2*c2*d1+b1*e2-b2*e1-2*c1*d2)*yq+d1*e2-d2*e1 ! ! Compute xi_x, xi_y, eta_x, eta_y ! etax = -etax/det etay = etay/det xix = xix/det xiy = -xiy/det return end function ubdry(iuk,yy) !*****************************************************************************80 ! !! UBDRY sets the parabolic inflow in terms of the value of the parameter ! implicit double precision(a-h,o-z) integer iuk double precision ubdry double precision yy if ( iuk == 1 ) then ubdry = (-2.0*yy+6.0)*yy/9.0D+00 else ubdry = 0.0D+00 end if return end function ubump(g,indx,ip,iqq,isotri,it,iukk,nelemn,neqn,nnodes,node,np,xc,yc) !*****************************************************************************80 ! !! UBUMP calculates the value of the sensitivity du/da on the bump ! Sets du/da = -uy*phi and dv/da=-vy*phi where phi is shape of bump ! implicit double precision(a-h,o-z) integer nelemn integer neqn integer nnodes integer np double precision g(neqn) integer indx(np,2) integer isotri(nelemn) integer node(nelemn,nnodes) double precision ubump double precision un(2) double precision unx(2) double precision uny(2) double precision xc(np) double precision xq double precision yc(np) double precision yq if ( isotri(it) == 0 ) then xq = xc(ip) yq = yc(ip) else if ( iqq == 1 ) then xq = 0.0D+00 yq = 0.0D+00 else if ( iqq == 2 ) then xq = 1.0D+00 yq = 1.0D+00 else if ( iqq == 3 ) then xq = 1.0D+00 yq = 0.0D+00 else if ( iqq == 4 ) then xq = 0.5 yq = 0.5 else if ( iqq == 5 ) then xq = 1.0D+00 yq = 0.5 else if ( iqq == 6 ) then xq = 0.5 yq = 0.0D+00 end if call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) end if ! ! Calculate value of uy and vy (old solutions) at node ! call uval(etax,etay,g,indx,isotri,it,nelemn,neqn, & nnodes,node,np,un,uny,unx,xc,xix,xiy,xq,yc,yq) if ( iukk == 1 ) then ubump = uny(1)*(xc(ip)-1.0)*(xc(ip)-3.0) else if ( iukk == 2 ) then ubump = uny(2)*(xc(ip)-1.0)*(xc(ip)-3.0) else write(*,*)'UBUMP called for iukk = ',iukk stop end if return end subroutine uval(etax,etay,g,indx,isotri,it,nelemn,neqn, & nnodes,node,np,un,uny,unx,xc,xix,xiy,xq,yc,yq) !*****************************************************************************80 ! !! UVAL evaluates the velocities and their X, Y derivatives at ! a given quadrature point. ! implicit double precision (a-h,o-z) external ubdry integer nelemn integer neqn integer nnodes integer np double precision bb double precision bx double precision by double precision etax double precision etay double precision g(neqn) integer i integer indx(np,2) integer ip integer iq integer isotri(nelemn) integer it integer iuk integer iun integer node(nelemn,nnodes) double precision ubc double precision ubdry double precision un(2) double precision unx(2) double precision uny(2) double precision xc(np) double precision xix double precision xiy double precision xq double precision yc(np) double precision yq un(1:2) = 0.0D+00 unx(1:2) = 0.0D+00 uny(1:2) = 0.0D+00 do iq = 1,nnodes if ( isotri(it) == 1 ) then call refqbf(xq,yq,iq,bb,bx,by,etax,etay,xix,xiy) else call qbf(xq,yq,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) end if ip = node(it,iq) do iuk = 1,2 iun = indx(ip,iuk) if ( iun > 0 ) then un(iuk) = un(iuk)+bb*g(iun) unx(iuk) = unx(iuk)+bx*g(iun) uny(iuk) = uny(iuk)+by*g(iun) else if ( iun == -1 ) then ubc = ubdry(iuk,yc(ip)) un(iuk) = un(iuk)+bb*ubc unx(iuk) = unx(iuk)+bx*ubc uny(iuk) = uny(iuk)+by*ubc end if end do end do return end subroutine uvdump (f,indx,insc,ivunit,long,mx,my, & nelemn,neqn,nnodes,node,np,press,reynld,yc) !*****************************************************************************80 ! !! UVDUMP creates a velocity file for use by PLOT3D. ! ! Given the following set of nodes: ! ! A B C ! D E F ! G H I ! ! the file will have the form: ! ! D, U(G), V(G), P ! D, U(H), V(H), P ! D, U(I), V(I), P ! D, U(D), V(D), P ! D, U(E), V(E), P ! D, U(F), V(F), P ! D, U(A), V(A), P ! D, U(B), V(B), P ! D, U(C), V(C), P ! ! Here both D and P are set to 1 for now, representing dummy values ! of density and pressure. ! ! PRESS Workspace, DOUBLE PRECISION PRESS(MX,MY), used to hold the ! computed pressures. ! implicit double precision (a-h,o-z) integer nelemn integer neqn integer nnodes integer np double precision alpha double precision dval double precision f(neqn) double precision fsmach integer i integer indx(np,2) integer insc(np) integer ip integer iset integer ivunit integer j integer k logical long integer mx integer my integer node(nelemn,nnodes) double precision press(mx,my) double precision reynld double precision time double precision ubdry double precision uval double precision vval double precision yc(np) save iset data iset /0/ iset = iset+1 dval = 1.0D+00 fsmach = 1.0D+00 alpha = 1.0D+00 time = 1.0D+00 call pval (f,insc,long,mx,my,nelemn,neqn,nnodes,node,np,press) ! ! If NX > NY, then nodes with a constant Y value are numbered consecutively. ! if ( long ) then write(ivunit,'(2I5)')mx,my write(ivunit,'(4G15.5)')fsmach,alpha,reynld,time do 30 ii = 1,4 do 20 j = 1,my do 10 i = 1,mx ip = (i-1)*my+j if ( ii == 1 ) then write(ivunit,'(G15.5)')dval else if ( ii == 2 ) then k = indx(ip,1) if (k == 0) then uval = 0.0D+00 else if (k < 0) then uval = ubdry(1,yc(ip)) else uval = f(k) end if write(ivunit,'(G15.5)')uval else if ( ii == 3 ) then k = indx(ip,2) if (k == 0) then vval = 0.0D+00 else vval = f(k) end if write(ivunit,'(G15.5)')vval else write(ivunit,'(G15.5)')press(i,j) end if 10 continue 20 continue 30 continue ! ! If NX < NY, then nodes with a constant X value are numbered consecutively. ! else write(ivunit,'(2I5)')mx,my write(ivunit,'(4G15.5)')fsmach,alpha,reynld,time do ii = 1,4 do i = 1,mx do j = 1,my if ( ii == 1 ) then write(ivunit,'(G15.5)')dval else if ( ii == 2 ) then ip = (i-1)*my+j k = indx(ip,1) if (k == 0) then uval = 0.0D+00 else if (k < 0) then uval = ubdry(1,yc(i)) else uval = f(k) end if write(ivunit,'(G15.5)')uval else if ( ii == 3 ) then k = indx(ip,2) if (k == 0) then vval = 0.0D+00 else vval = f(k) end if write(ivunit,'(G15.5)')vval else write(ivunit,'(G15.5)')press(i,j) end if end do end do end do end if write (*,*) 'UVDUMP wrote data set ',iset,' to file.' return end subroutine xydump (ixunit,long,np,nx,ny,xc,yc) !*****************************************************************************80 ! !! XYDUMP creates a grid file for use by PLOT3D. ! ! Given the following set of nodes: ! ! A B C ! D E F ! G H I ! ! the file will have the form: ! ! X(G), X(H), X(I), X(D), X(E), X(F), X(A), X(B), X(C), ! Y(G), Y(H), Y(I), Y(D), Y(E), Y(F), Y(A), Y(B), Y(C). ! implicit double precision (a-h,o-z) integer np integer i integer ip integer iset integer ixunit integer j logical long integer mx integer my integer nx integer ny double precision xc(np) double precision yc(np) save iset data iset /0/ iset = iset+1 mx = 2*nx-1 my = 2*ny-1 ! ! If NX > NY, then nodes with a constant Y value are numbered consecutively. ! if ( long ) then write(ixunit,'(2I15)')mx,my do i = 1,my do j = 1,mx ip = (j-1)*my+i write(ixunit,'(G15.5)')xc(ip) end do end do do i = 1,my do j = 1,mx ip = (j-1)*my+i write(ixunit,'(G15.5)')yc(ip) end do end do ! ! If NX < NY, then nodes with a constant X value are numbered consecutively. ! else write(ixunit,'(2I15)')my,mx do j = 1,mx do i = 1,my ip = (j-1)*my+i write(ixunit,'(G15.5)')xc(ip) end do end do do j = 1,mx do i = 1,my ip = (j-1)*my+i write(ixunit,'(G15.5)')yc(ip) end do end do end if write (*,*) 'XYDUMP wrote data set ',iset,' to file.' return end subroutine dgbfa ( abd, lda, n, ml, mu, ipvt, info ) !*****************************************************************************80 ! !! DGBFA factors a double precision band matrix by elimination. ! ! Discussion: ! ! dgbfa is usually called by dgbco, but it can be called ! directly with a saving in time if rcond is not needed. ! ! Parameters: ! ! on entry ! ! abd double precision(lda, n) ! contains the matrix in band storage. the columns ! of the matrix are stored in the columns of abd and ! the diagonals of the matrix are stored in rows ! ml+1 through 2*ml+mu+1 of abd . ! see the comments below for details. ! ! lda integer ! the leading dimension of the array abd . ! lda must be >= 2*ml + mu + 1 . ! ! n integer ! the order of the original matrix. ! ! ml integer ! number of diagonals below the main diagonal. ! 0 <= ml < n . ! ! mu integer ! number of diagonals above the main diagonal. ! 0 <= mu < n . ! more efficient if ml <= mu . ! on return ! ! abd an upper triangular matrix in band storage and ! the multipliers which were used to obtain it. ! the factorization can be written a = l*u where ! l is a product of permutation and unit lower ! triangular matrices and u is upper triangular. ! ! ipvt integer(n) ! an integer vector of pivot indices. ! ! info integer ! = 0 normal value. ! = k if u(k,k) == 0.0D+00 . this is not an error ! condition for this subroutine, but it does ! indicate that dgbsl will divide by zero if ! called. use rcond in dgbco for a reliable ! indication of singularity. ! ! band storage ! ! if a is a band matrix, the following program segment ! will set up the input. ! ! ml = (band width below the diagonal) ! mu = (band width above the diagonal) ! m = ml + mu + 1 ! do j = 1, n ! i1 = max ( 1, j-mu ) ! i2 = min ( n, j+ml ) ! do i = i1, i2 ! k = i - j + m ! abd(k,j) = a(i,j) ! end do ! end do ! ! this uses rows ml+1 through 2*ml+mu+1 of abd . ! in addition, the first ml rows in abd are used for ! elements generated during the triangularization. ! the total number of rows needed in abd is 2*ml+mu+1 . ! the ml+mu by ml+mu upper left triangle and the ! ml by ml lower right triangle are not referenced. ! ! linpack. this version dated 08/14/78 . ! cleve moler, university of new mexico, argonne national lab. ! integer lda integer n double precision abd(lda,n) integer i integer i0 integer info integer ipvt(n) integer idamax integer j integer j0 integer j1 integer ju integer jz integer k integer l integer lm integer m integer ml integer mm integer mu double precision t m = ml + mu + 1 info = 0 ! ! Zero initial fill-in columns. ! j0 = mu + 2 j1 = min ( n, m ) - 1 do jz = j0, j1 i0 = m + 1 - jz do i = i0, ml abd(i,jz) = 0.0D+00 end do end do jz = j1 ju = 0 ! ! Gaussian elimination with partial pivoting. ! do k = 1, n-1 ! ! Zero next fill-in column. ! jz = jz + 1 if ( jz <= n ) then abd(1:ml,jz) = 0.0D+00 end if ! ! Find L = pivot index. ! lm = min ( ml, n-k ) l = idamax ( lm+1, abd(m,k), 1 ) + m - 1 ipvt(k) = l + k - m ! ! Zero pivot implies this column already triangularized. ! if ( abd(l,k) == 0.0D+00 ) then info = k ! ! Interchange if necessary. ! else if ( l /= m ) then t = abd(l,k) abd(l,k) = abd(m,k) abd(m,k) = t end if ! ! Compute multipliers. ! t = -1.0D+00 / abd(m,k) call dscal ( lm, t, abd(m+1,k), 1 ) ! ! Row elimination with column indexing. ! ju = min ( max ( ju, mu+ipvt(k) ), n ) mm = m do j = k+1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if ( l /= mm ) then abd(l,j) = abd(mm,j) abd(mm,j) = t end if call daxpy ( lm, t, abd(m+1,k), 1, abd(mm+1,j), 1 ) end do end if end do ipvt(n) = n if ( abd(m,n) == 0.0D+00 ) then info = n end if return end subroutine dgbsl ( abd, lda, n, ml, mu, ipvt, b, job ) !*****************************************************************************80 ! !! DGBSL solves a double precision banded system factored by SGBCO or SGBFA. ! ! Discussion: ! ! SGBSL can solve either a * x = b or trans(a) * x = b. ! ! Parameters: ! ! on entry ! ! abd double precision(lda, n) ! the output from dgbco or dgbfa. ! ! lda integer ! the leading dimension of the array abd . ! ! n integer ! the order of the original matrix. ! ! ml integer ! number of diagonals below the main diagonal. ! ! mu integer ! number of diagonals above the main diagonal. ! ! ipvt integer(n) ! the pivot vector from dgbco or dgbfa. ! ! b double precision(n) ! the right hand side vector. ! ! job integer ! = 0 to solve a*x = b , ! = nonzero to solve trans(a)*x = b , where ! trans(a) is the transpose. ! ! on return ! ! b the solution vector x . ! ! error condition ! ! a division by zero will occur if the input factor contains a ! zero on the diagonal. technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of lda . it will not occur if the subroutines are ! called correctly and if dgbco has set rcond > 0.0D+00 ! or dgbfa has set info == 0 . ! ! to compute inverse(a) * c where c is a matrix ! with p columns ! call dgbco ( abd, lda, n, ml, mu, ipvt, rcond, z ) ! if (rcond is too small) go to ... ! do j = 1, p ! call dgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 ) ! end do ! ! linpack. this version dated 08/14/78 . ! cleve moler, university of new mexico, argonne national lab. ! integer lda integer n double precision abd(lda,n) double precision b(n) integer ipvt(n) integer job integer k integer l integer la integer lb integer lm integer m integer ml integer mu double precision ddot double precision t m = mu + ml + 1 ! ! JOB = 0, Solve a * x = b. ! ! First solve l*y = b. ! if ( job == 0 ) then if ( ml > 0 ) then do k = 1, n-1 lm = min ( ml, n-k ) l = ipvt(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if call daxpy ( lm, t, abd(m+1,k), 1, b(k+1), 1 ) end do end if ! ! Now solve u*x = y. ! do k = n, 1, -1 b(k) = b(k) / abd(m,k) lm = min ( k, m ) - 1 la = m - lm lb = k - lm t = -b(k) call daxpy ( lm, t, abd(la,k), 1, b(lb), 1 ) end do ! ! JOB nonzero, solve trans(a) * x = b. ! ! First solve trans(u)*y = b. ! else do k = 1, n lm = min ( k, m ) - 1 la = m - lm lb = k - lm t = ddot ( lm, abd(la,k), 1, b(lb), 1 ) b(k) = ( b(k) - t ) / abd(m,k) end do ! ! Now solve trans(l)*x = y ! if ( ml > 0 ) then do k = n-1, 1, -1 lm = min ( ml, n-k ) b(k) = b(k) + ddot ( lm, abd(m+1,k), 1, b(k+1), 1 ) l = ipvt(k) if ( l /= k ) then t = b(l) b(l) = b(k) b(k) = t end if end do end if end if return end subroutine daxpy ( n, da, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DAXPY computes a constant times a vector plus a vector. ! ! uses unrolled loops for increments equal to one. ! jack dongarra, linpack, 3/11/78. ! modified 12/3/93, array(1) declarations changed to array(*) ! double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,n ! if ( n <= 0)return if (da == 0.0d+00 ) return if ( incx == 1.and.incy == 1)go to 20 ! ! code for unequal increments or equal increments ! not equal to 1 ! ix = 1 iy = 1 if ( incx < 0)ix = (-n+1)*incx + 1 if ( incy < 0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return ! ! code for both increments equal to 1 ! ! ! clean-up loop ! 20 m = mod(n,4) do i = 1,m dy(i) = dy(i) + da*dx(i) end do if ( n < 4 ) return 40 continue do i = m+1, n, 4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) end do return end subroutine dscal ( n, da, dx, incx ) ! !*****************************************************************************80 ! !! DSCAL scales a vector by a constant. ! ! uses unrolled loops for increment equal to one. ! jack dongarra, linpack, 3/11/78. ! modified 3/93 to return if incx <= 0. ! modified 12/3/93, array(1) declarations changed to array(*) ! double precision da,dx(*) integer i,incx,m,n,nincx ! if ( n <= 0 .or. incx <= 0 )return if ( incx == 1)go to 20 ! ! code for increment not equal to 1 ! nincx = n*incx do i = 1,nincx,incx dx(i) = da*dx(i) end do return ! ! code for increment equal to 1 ! ! ! clean-up loop ! 20 continue m = mod(n,5) if ( m == 0 ) go to 40 dx(1:m) = da*dx(1:m) if ( n < 5 ) return 40 continue do i = m+1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) end do return end