program main !*****************************************************************************80 ! !! MAIN is the main program for ARBY4. ! ! Discussion: ! ! ARBY4 solves a fluid flow problem which has several parameters. ! ! ARBY4 can solve the problem using the finite element method. The ! resulting solution will be called a "full" solution. ! ! Once a full solution has been calculated, ARBY4 can compute the ! first several derivatives of the solution with respect to the ! Reynolds number, derive a reduced basis, and find reduced ! solutions to the problem at a variety of Reynolds numbers. ! ! Diary: ! ! 12 December 2000 ! ! Double precision constants REQUIRE D+00 or something similar, otherwise ! they are computed as though they were REAL. So I tried to fix that, ! quickly. ! ! 14 September 1996 ! ! Mr Lee wants a code that will solve the driven cavity problem. ! I thought I would get FLOW4, but somehow ARBY4 may be the ! quicker solution. ! ! 12 September 1996 ! ! Worked on GFL2RB and GRB2FL. ! ! 11 September 1996 ! ! Max urged me to write a reduced basis paper, as Kazi Ito has already ! done so. ! ! I updated GFL2RB. ! ! 17 August 1996 ! ! OK, making transition BACK to representation of full solution as ! GFL(GRB) = GFLRB + Q * GRB. ! ! 16 August 1996 ! ! GFLRB is already stored. So now I have to reformulate the representation. ! ! 15 August 1996 ! ! I decided that I had to treat boundary conditions as follows: ! ! The reduced solution is represented as Y = Y0 + YBC*cbc + YFE*cfe ! ! Here, Y0 is the point at which the basis was generated. ! Each column of YBC corresponds to the flow solution of the problem ! at Y0 with boundary conditions differentiated with respect to ! parameter I. The YFE stuff is as usual. ! implicit none ! ! Set parameters that are independent. ! integer, parameter :: liv = 60 integer, parameter :: maxbcrb = 3 integer, parameter :: maxferb = 6 integer, parameter :: maxnx = 21 integer, parameter :: maxny = 21 integer, parameter :: maxparb = 5 integer, parameter :: maxparf = 5 ! ! Set parameters that are dependent on parameters. ! ! The assignment of LDAFL should really read (ldafl = 29*min(nx,ny)). ! integer, parameter :: ldafl = 29 * maxny integer, parameter :: maxcofrb = maxbcrb + maxferb integer, parameter :: maxelm = 2 * ( maxnx - 1 ) * ( maxny - 1 ) integer, parameter :: maxnfl = & 2 * ( 2 * maxnx - 1 ) * ( 2 * maxny - 1 ) + maxnx * maxny integer, parameter :: maxnp = ( 2 * maxnx - 1 ) * ( 2 * maxny - 1 ) integer, parameter :: maxpar = maxparb + maxparf + 1 ! ! Set parameters that are dependent on parameters that are dependent ! on parameters. ! integer, parameter :: lv = 78+maxpar*(maxpar+21)/2 double precision afl(ldafl,maxnfl) double precision arb(maxcofrb,maxcofrb) double precision area(3,maxelm) character ( len = 9 ) chtime character ( len = 80 ) command double precision cost double precision cost0 double precision costb double precision costp double precision costu double precision costv character ( len = 8 ) date double precision detlog double precision detman double precision difcof(maxcofrb) character ( len = 30 ) disfil double precision dopt(maxpar) double precision dpar double precision drey logical dvneq logical echo double precision epsdif character ( len = 2 ) eqn(maxnfl) real estart real estop double precision etaq(3) double precision factj double precision gfl(maxnfl) double precision gflafl(maxnfl) double precision gflnrm double precision gflopt(maxnfl) double precision gflrb(maxnfl) double precision gflsav(maxnfl) double precision gflsen(maxnfl) double precision gfltar(maxnfl) double precision gfltay(maxnfl) double precision gfltmp(maxnfl) double precision grb(maxcofrb) double precision grbarb(maxcofrb) double precision grbopt(maxcofrb) double precision grbsav(maxcofrb) double precision grbsen(maxcofrb) double precision grbtay(maxcofrb) character ( len = 20 ) gridx character ( len = 20 ) gridy double precision gsen(maxcofrb) double precision hx double precision hy integer i integer ibs integer ibump integer icolrb(maxcofrb) integer ierror integer ifs integer ihi integer ijac integer ilo integer indx(3,maxnp) integer info integer iopt(maxpar) integer ios integer ipar integer ipivfl(maxnfl) integer ipivrb(maxcofrb) integer iprint integer isotri(maxelm) integer itemp integer ival integer ival1 integer ival2 integer ivopt(liv) integer iwrite integer j integer jhi integer jlo integer jtay integer klo integer lchar integer lenc logical s_eqi integer maxnew integer maxopt integer maxsim integer mhi integer mlo integer nbcrb integer ncofrb integer nelem integer neqnfl integer nferb integer nhi integer nlband integer nlo integer node(6,maxelm) integer nodelm(maxnp) integer np integer npar integer nparb integer nparf integer npe integer nprof(2*maxny-1) integer nsenfl integer ntay integer numdif integer numnew integer numopt integer numsim integer nx integer ny double precision p(maxnp) double precision par(maxpar) double precision parafl(maxpar) double precision pararb(maxpar) double precision pardif ( maxpar) double precision paropt(maxpar) double precision parrb(maxpar) double precision parsav(maxpar) double precision parsen(maxpar) double precision partar(maxpar) double precision phifl(3,6,10,maxelm) double precision phirb(3,maxcofrb,15,maxelm) double precision rb(maxnfl,maxcofrb) character ( len = 20 ) region double precision resfl(maxnfl) double precision resflsav(maxnfl) double precision resfltmp(maxnfl) double precision resrb(maxcofrb) double precision reynld double precision reytay double precision rbase(maxcofrb,maxcofrb) double precision rmax double precision senfl(maxnfl,maxcofrb) double precision senrb(maxcofrb,maxcofrb) double precision splbmp(maxparb+2) double precision splflo(maxparf) real tarray(2) double precision taubmp(maxparb+2) double precision tauflo(maxparf) character ( len = 30 ) tecfil double precision temp character ( len = 10 ) time double precision tolnew double precision tolopt double precision tolsim character ( len = 10 ) tstart character ( len = 10 ) tstop double precision u(maxnp) double precision v(maxnp) double precision value double precision vopt(lv) double precision wateb double precision watep double precision wateu double precision watev double precision wquad(3) double precision xbl double precision xbr double precision xc(maxnp) double precision xmax double precision xmin double precision xopt(maxpar) double precision xprof double precision xquad(3,maxelm) double precision xrange double precision xsiq(3) double precision ybl double precision ybr double precision yc(maxnp) double precision ymax double precision ymin double precision yquad(3,maxelm) double precision yrange call timestamp ( ) iprint = 1 ! ! Get initial CPU clock reading. ! call etime ( tarray, estart ) echo = .false. call hello ( maxnx, maxny ) ! ! Get the starting time. ! call date_and_time ( date, time ) tstart = time ! ! Open the file in which we record the user input. ! open ( unit = 17, file = 'arby.in', status = 'replace' ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ARBY4 - Init: Initialize all data.' call init(afl,arb,area,command,cost,costb,costp,costu, & costv,difcof,disfil,drey,epsdif,eqn,etaq,gfl,gflafl, & gflrb,gflsav,gflsen,gfltar,gfltay,grb,grbarb,grbsav, & grbsen,grbtay,gridx,gridy,hx,hy,ibs,ibump,icolrb,ierror, & ifs,ihi,ijac,ilo,indx,iopt,ipar, & ipivfl,ipivrb,isotri,iwrite,jhi,jlo, & ldafl,maxcofrb,maxelm,maxnew,maxnfl,maxnp, & maxny,maxopt,maxpar,maxparb,maxparf,maxsim, & nbcrb,ncofrb,nelem,neqnfl,nferb, & nlband,node,nodelm,np,npar,nparb,nparf,npe,nprof,nsenfl,ntay, & numnew,numopt,numsim,nx, & ny,par,parafl,pararb,pardif,parrb,parsav,parsen,partar, & phifl,phirb,rbase,rb,region,resfl,resflsav,resrb,reynld, & reytay,senfl, & senrb,splbmp,splflo,taubmp,tauflo,tecfil,tolnew,tolopt, & tolsim,value,wateb,watep,wateu,watev,wquad,xbl,xbr, & xc,xmax,xmin,xprof,xquad,xrange,xsiq, & ybl,ybr,yc,ymax,ymin,yquad,yrange) ! ! Read the next command from the user ! command = '?' do if ( command(1:1) /= '#' .and. len_trim ( command ) /= 0 ) then if ( 0 < iprint ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter command:' end if end if read ( *, '(a)', iostat = ios ) command if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Fatal error!' write ( *, * ) ' I/O error for the input file.' stop end if if ( echo ) then write ( *, '(a)' ) trim ( command ) end if write ( 17, '(a)' ) trim ( command ) ! ! Check for output-suppressing semicolon at end. ! if ( 0 < lenc ) then if ( command(lenc:lenc) == ';' ) then iprint = 0 command(lenc:lenc) = ' ' lenc = lenc-1 else iprint = 1 end if end if 15 continue if ( command(1:1) == '#') then cycle else if ( command == ' ') then cycle end if if ( 0 < iprint ) then write ( *, '(a)' ) ' ' end if ! ! COMPARE ! if ( s_eqi ( command,'compare')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - Compare:' write ( *, '(a)' ) ' Compare full solutions GFL and GFLSAV.' end if call fxfl(area,eqn,gfl,ifs,indx,nelem,neqnfl,node,np,npar,nparf, & par,phifl,region,resfl,splflo,tauflo,xrange,yc,yrange) call fxfl(area,eqn,gflsav,ifs,indx,nelem,neqnfl,node,np,npar, & nparf,par,phifl,region,resflsav,splflo,tauflo,xrange,yc,yrange) gfltmp(1:neqnfl) = gfl(1:neqnfl) - gflsav(1:neqnfl) resfltmp(1:neqnfl) = resfl(1:neqnfl) - resflsav(1:neqnfl) call nrmflo ( gfltmp, indx, neqnfl, np, resfltmp ) ! ! COSTFL ! else if ( s_eqi ( command,'costfl')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - Cost FL:' write ( *, '(a)' ) ' Evaluate the cost function for GFL.' write ( *, '(a)' ) ' ' end if call getcst(cost,costb,costp,costu,costv,gfl,gfltar,indx,neqnfl,np, & nparb,nprof,ny,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) write ( *, '(a,g14.6)' ) 'Cost of GFL:',cost ! ! COSTRB ! else if ( s_eqi ( command,'costrb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - Cost RB:' write ( *, '(a)' ) ' Evaluate the cost function for GRB.' write ( *, '(a)' ) ' ' end if call grb2fl(gfl,gflrb,grb,maxnfl,ncofrb,neqnfl,rb) call getcst(cost,costb,costp,costu,costv,gfl,gfltar,indx,neqnfl,np, & nparb,nprof,ny,splbmp,taubmp,wateb,watep,wateu,watev,xbl,xbr,ybl,ybr,yc) write ( *, '(a,g14.6)' ) 'Cost of GRB:',cost ! ! DETFPFL ! else if ( s_eqi ( command,'detfpfl')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - DET FP FL' write ( *, '(a)' ) ' Compute determinant of full jacobian.' end if call fpfl(afl,area,eqn,gfl,indx,ldafl,maxelm,nelem,neqnfl, & nlband,node,np,npar,par,phifl) call dfacfl(afl,ldafl,neqnfl,nlband,nlband,ipivfl,info) call ddetfl(afl,detlog,detman,ipivfl,ldafl,neqnfl,nlband,nlband) write ( *, '(a,g14.6,a,g14.6)' ) 'Determinant = ', detman, & ' * 10 ** ', detlog ! ! DETFPRB ! else if ( s_eqi ( command,'detfprb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - DET FP RB' write ( *, '(a)' ) ' Compute the reduced jacobian,' write ( *, '(a)' ) ' and its determinant.' end if call fprb(arb,area,grb,indx,maxcofrb,maxelm,maxnfl,nbcrb, & ncofrb,nelem,nferb,node,np,nx,ny,phirb,rb,reynld,xc,xrange,yc,yrange) call dfacrb(arb,maxcofrb,ncofrb,ipivrb,info) call ddetrb(arb,detlog,detman,ipivrb,maxcofrb,ncofrb) write ( *, '(a,g14.6,a,g14.6)' ) 'Determinant = ', detman, & ' * 10 ** ', detlog ! ! DIFFPRB ! else if ( s_eqi ( command,'diffprb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - DIF FP RB' write ( *, '(a)' ) ' Estimate reduced jacobian by finite diff.' end if pararb(1:npar) = par(1:npar) grbarb(1:ncofrb) = grb(1:ncofrb) arb(1:maxcofrb,1:maxcofrb) = 0.0D+00 call diffprb(arb,area,epsdif,grb,indx,maxcofrb,maxelm, & maxnfl,nbcrb,ncofrb,nelem,nferb,node,np,npar,nparf,nx, & ny,par,phirb,rb,resrb,reynld,tauflo,xc,xrange,yc,yrange) ! ! DIFSENFL ! else if ( s_eqi ( command,'difsenfl')) then if ( ipar <= 0 .or. npar < ipar ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ARBY4 - Warning!' write ( *, '(a)' ) ' Cancelling the DIFSENFL command.' write ( *, '(a,i6)' ) ' IPAR = ',ipar write ( *, '(a)' ) ' but IPAR must be at least 0' write ( *, '(a,i6)' ) ' and no more than NPAR = ',npar cycle end if if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - DifSenFL:' write ( *, '(a)' ) ' Estimate full solution sensitivity ' write ( *, '(a)' ) ' with respect to parameter ',ipar write ( *, '(a)' ) ' via finite differences.' end if ipar = npar dpar = drey parsen(1:npar) = par(1:npar) gflsen(1:neqnfl) = gfl(1:neqnfl) call difsenfl(afl,area,difcof,dpar,eqn,gfl,gflafl,ifs,ijac,indx,ipar, & ipivfl,iwrite,ldafl,maxcofrb,maxelm,maxnew,maxnfl,ncofrb,nelem,neqnfl, & nlband,node,np,npar,nparf,par,parafl,phifl,region,resfl,senfl,splflo, & tauflo,tolnew,xrange,yc,yrange) ! ! DIFSENRB ! else if ( s_eqi ( command,'difsenrb')) then if ( ipar <= 0 .or. npar < ipar ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ARBY4 - Warning!' write ( *, '(a)' ) ' Cancelling the DIFSENRB command.' write ( *, '(a,i6)' ) ' IPAR = ', ipar write ( *, '(a)' ) ' but IPAR must be at least 0' write ( *, '(a,i6)' ) ' and no more than NPAR = ', npar cycle end if if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - Dif Sen RB:' write ( *, '(a)' ) ' Estimate reduced solution sensitivity ' write ( *, '(a,i6)' ) ' with respect to parameter ', ipar write ( *, '(a)' ) ' via finite differences.' end if ipar = npar dpar = drey parsen(1:npar) = par(1:npar) grbsen(1:ncofrb) = grb(1:ncofrb) call difsenrb(arb,area,difcof,dpar,grb,grbarb,indx,ipar,ipivrb,iwrite, & maxcofrb,maxelm,maxnew,maxnfl,nbcrb,ncofrb,nelem,nferb,node,np,npar, & nparf,nx,ny,par,pararb,phirb,rb,resrb,senrb,tauflo,tolnew,xc,xrange, & yc,yrange) ! ! DisPlot ! else if ( s_eqi ( command,'displot')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - DisPlot: Write data to DISPLAY plot file.' end if call intprs(gfl,indx,nelem,neqnfl,node,np,p) u(1:np) = gfl(indx(1,1:np)) v(1:np) = gfl(indx(2,1:np)) call wrdis(disfil,eqn,indx,isotri,maxcofrb,maxnfl,ncofrb,nelem,neqnfl, & node,np,npar,npe,nprof,nx,ny,p,par,rb,u,v,xc,xprof,yc) ! ! DREY = ! else if ( s_eqi ( command(1:4),'drey')) then if ( s_eqi ( command(1:5),'drey=')) then read(command(6:),*)drey else write ( *, '(a)' ) 'Enter value for DREY:' read(*,*)drey write(17,*)drey end if if ( 0 < iprint ) then write ( *, '(a,g14.6)' ) 'DREY set to ',drey end if ! ! EPSDIF = ! else if ( s_eqi ( command(1:6),'epsdif')) then if ( s_eqi ( command(1:7),'epsdif=')) then read(command(8:),*)epsdif else write ( *, '(a)' ) 'Enter value for EPSDIF:' read ( *, * ) epsdif write ( 17, '(g14.6)' ) epsdif end if if ( 0 < iprint ) then write ( *, '(a,g14.6)' ) 'EPSDIF set to ', epsdif end if ! ! ECHO ! else if ( s_eqi ( command,'echo')) then echo = .not.echo if ( echo) then write(*,'(a)')command if ( 0 < iprint ) then write ( *, '(a)' ) 'User commands will be echoed.' end if else if ( 0 < iprint ) then write ( *, '(a)' ) 'User commands will not be echoed.' end if end if ! ! EXPAND GRB ! else if ( s_eqi ( command,'expand grb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - Expand GRB:' write ( *, '(a)' ) ' Compute GFL = RB * GRB' end if call grb2fl(gfl,gflrb,grb,maxnfl,ncofrb,neqnfl,rb) ! ! DISFIL = ! else if ( s_eqi ( command(1:6),'disfil')) then if ( s_eqi ( command(1:7),'disfil=')) then disfil = command(8:) else write ( *, '(a)' ) 'Enter the DISPLAY output file name:' read(*,'(a)')disfil write(17,'(a)')disfil end if if ( 0 < iprint ) then write ( *, '(a)' ) 'DISPLAY output file name set to ' // trim ( disfil ) end if ! ! FILTEC = ! else if ( s_eqi ( command(1:6),'tecfil')) then if ( s_eqi ( command(1:7),'tecfil=')) then tecfil = command(8:) else write ( *, '(a)' ) 'Enter the TECPLOT output file name:' read(*,'(a)')tecfil write(17,'(a)')tecfil end if if ( 0 < iprint ) then write ( *, '(a)' ) 'TECPLOT output file name set to ' // trim ( tecfil ) end if ! ! FPFL ! else if ( s_eqi ( command,'fpfl')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FPFL - Evaluate full jacobian.' end if call fpfl(afl,area,eqn,gfl,indx,ldafl,maxelm,nelem,neqnfl, & nlband,node,np,npar,par,phifl) ! ! FPIRB ! else if ( s_eqi ( command,'fpirb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FPIRB:' write ( *, '(a)' ) ' Evaluate FP indirectly at reduced solution GRB,' write ( *, '(a)' ) ' by expanding GRB to GFL, evaluating FP(GFL),' write ( *, '(a)' ) ' and then reducing.' end if call fpirb(afl,arb,area,eqn,gflrb,grb,indx,ldafl,maxcofrb, & maxelm,maxnfl,nbcrb,ncofrb,nelem,neqnfl,nferb,nlband,node, & np,npar,nx,ny,par,phifl,rb,xc,xrange,yc,yrange) ! ! FPRB ! else if ( s_eqi ( command,'fprb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FPRB - Evaluate reduced jacobian.' end if pararb(1:npar) = par(1:npar) grbarb(1:ncofrb) = grb(1:ncofrb) call fprb(arb,area,grb,indx,maxcofrb,maxelm,maxnfl,nbcrb,ncofrb,nelem, & nferb,node,np,nx,ny,phirb,rb,reynld,xc,xrange,yc,yrange) ! ! FPRB = 0 ! else if ( s_eqi ( command,'fprb=0')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FPRB = 0 - Zero out the reduced jacobian.' end if arb(1:maxcofrb,1:maxcofrb) = 0.0D+00 ! ! FXFL ! else if ( s_eqi ( command,'fxfl')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FX FL:' write ( *, '(a)' ) ' Evaluate FXFL at full solution GFL.' end if call fxfl(area,eqn,gfl,ifs,indx,nelem,neqnfl,node,np,npar,nparf,par, & phifl,region,resfl,splflo,tauflo,xrange,yc,yrange) ! ! FXIRB ! else if ( s_eqi ( command,'fxirb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FX IRB:' write ( *, '(a)' ) ' Evaluate FX indirectly at reduced solution GRB,' write ( *, '(a)' ) ' by expanding GRB to GFL, evaluating F(GFL),' write ( *, '(a)' ) ' and then reducing.' end if call fxirb(area,eqn,gflrb,grb,ifs,indx,maxcofrb,maxnfl,nbcrb,ncofrb, & nelem,neqnfl,nferb,node,np,npar,nparf,nx,ny,par,phifl,rb,region,resrb, & splflo,tauflo,xc,xrange,yc,yrange) ! ! FXRB ! else if ( s_eqi ( command,'fxrb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FX RB:' write ( *, '(a)' ) ' Evaluate FX directly at reduced solution GRB.' end if reynld = par(npar) call fxrb(area,grb,indx,maxcofrb,maxelm,maxnfl,nbcrb,ncofrb,nelem,nferb, & node,np,npar,nparf,nx,ny,par,phirb,rb,resrb,reynld,tauflo,xc,xrange, & yc,yrange) ! ! FXRB = 0 ! else if ( s_eqi ( command,'fxrb=0')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - FXRB = 0' end if resrb(1:ncofrb) = 0.0D+00 ! ! GETGSEN ! else if ( s_eqi ( command,'getgsen')) then call getgsen(grb,gsen,icolrb,maxcofrb,nbcrb,ncofrb,nsenfl,rbase) ! ! GETRB ! else if ( s_eqi ( command,'getrb')) then if ( 0 < iprint ) then write ( *, '(a)' ) 'ARBY4 - Get RB:' write ( *, '(a)' ) ' Compute reduced basis RB at current GFL.' end if if ( ncofrb < 0) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ARBY4 - Warning!' write ( *, '(a)' ) ' The GETRB command is being cancelled,' write ( *, '(a,i6)' ) ' since NCOFRB is ', ncofrb write ( *, '(a)' ) ' Please use the "NCOFRB = " command first!' end if if ( dvneq ( npar, par, parsen ) ) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Error!' write ( *, * ) ' Please compute the sensitivities first,' write ( *, * ) ' using the GETSENFL command!' cycle end if ! ! Save the affine displacement vector GFLRB. ! parrb(1:npar) = par(1:npar) gflrb(1:neqnfl) = gfl(1:neqnfl) ! ! Get boundary condition reduced basis vectors. ! call getbcrb(gflrb,maxcofrb,maxnfl,nbcrb,neqnfl,rb) ! ! Get finite element reduced basis vectors. ! call getferb(icolrb,maxcofrb,maxnfl,nbcrb,ncofrb,neqnfl, & nferb,nsenfl,rb,rbase,senfl,senrb) grb(1:ncofrb) = 0.0D+00 ! ! Compute PHIRB, the reduced basis functions at quadrature points. ! if ( 0 < iprint ) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Note:' write ( *, * ) ' Automatically evaluating PHIRB, the reduced' write ( *, * ) ' basis functions at quadrature points.' end if call setprb(eqn,indx,maxcofrb,maxelm,maxnfl,nelem,neqnfl, & ncofrb,node,np,phifl,phirb,rb) ! ! GETSENFL ! else if ( s_eqi ( command,'getsenfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Get Sen FL:' write ( *, * ) ' Compute full solution sensitivities with ' write ( *, * ) ' respect to the parameter REYNLD, up to order' write ( *, * ) ' NSENFL = ',nsenfl end if if ( nsenfl < 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' The GETSENFL command is being cancelled,' write ( *, * ) ' since NSENFL is ',nsenfl write ( *, * ) ' Please use the "NSENFL = " command first!' else parsen(1:npar) = par(1:npar) gflsen(1:neqnfl) = gfl(1:neqnfl) call getsenfl(afl,area,eqn,gfl,indx,ipivfl,ldafl,maxcofrb,maxnfl, & nelem,neqnfl,nlband,node,np,npar,nsenfl,par,phifl,resfl,senfl) end if ! ! GETSENRB ! else if ( s_eqi ( command,'getsenrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Get Sen RB' write ( *, * ) ' Compute reduced basis sensitivities.' end if call getsenrb(maxcofrb,maxnfl,ncofrb,neqnfl,nsenfl,rb,senfl,senrb) ! ! GFL = 0, GFLSAV, GFLTAY, TAYLOR ! else if ( s_eqi ( command(1:4),'gfl=')) then if ( command(5:5) == '0') then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFL = 0' write ( *, * ) ' Set full solution estimate GFL to zero.' end if gfl(1:neqnfl) = 0.0D+00 else if ( s_eqi ( command(5:10),'gflsav')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFL = GFLSAV' write ( *, * ) ' Set full solution estimate GFL to GFLSAV.' end if gfl(1:neqnfl) = gflsav(1:neqnfl) else if ( s_eqi ( command(5:10),'gfltay')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFL = GFLTAY' write ( *, * ) ' Set full solution estimate GFL to GFLTAY.' end if gfl(1:neqnfl) = gfltay(1:neqnfl) else if ( s_eqi ( command(5:10),'taylor')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFL = TAYLOR' write ( *, * ) ' Set full solution GFL to Taylor prediction,' write ( *, * ) ' based at REYTAY = ',reytay write ( *, * ) ' using NTAY = ',ntay,' terms.' end if gfl(1:neqnfl) = gfltay(1:neqnfl) do i = 1, neqnfl do j = 1, ntay call fact(j,factj) temp = ((reynld-reytay)**j)/factj gfl(i) = gfl(i)+temp*senfl(i,j) end do end do end if ! ! GFLSAV = GFL ! else if ( s_eqi ( command,'gflsav=gfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLSAV = GFL' write ( *, * ) ' Save current full solution.' end if parsav(1:npar) = par(1:npar) gflsav(1:neqnfl) = gfl(1:neqnfl) resflsav(1:neqnfl) = resfl(1:neqnfl) ! ! GFLTAY = 0, GFL, GFLSAV ! else if ( s_eqi ( command(1:7),'gfltay=')) then if ( command(8:8) == '0') then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTAY = 0' write ( *, * ) ' Set Taylor base full solution to zero.' end if gfltay(1:neqnfl) = 0.0D+00 else if ( s_eqi ( command(8:13),'gflsav')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTAY = GFLSAV' write ( *, * ) ' Set Taylor base full solution to GFLSAV.' end if gfltay(1:neqnfl) = gflsav(1:neqnfl) else if ( s_eqi ( command(8:10),'gfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTAY = GFL' write ( *, * ) ' Set Taylor base full solution to GFL.' end if gfltay(1:neqnfl) = gfl(1:neqnfl) end if ! ! GFLTMP = 0/GFL/GFL-GFLSAV/GFL-GFLTAR/GFLSAV/GFLSAV-GFLTAY ! else if ( s_eqi ( command(1:7),'gfltmp=')) then if ( command(8:8) == '0') then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTMP = 0' end if gfltmp(1:neqnfl) = 0.0D+00 else if ( s_eqi ( command(8:),'gfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTMP = GFL' end if gfltmp(1:neqnfl) = gfl(1:neqnfl) else if ( s_eqi ( command(8:),'gfl-gflsav')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTMP = GFL-GFLSAV' end if gfltmp(1:neqnfl) = gfl(1:neqnfl) - gflsav(1:neqnfl) else if ( s_eqi ( command(8:),'gfl-gfltar')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTMP = GFL-GFLTAR' end if gfltmp(1:neqnfl) = gfl(1:neqnfl) - gfltar(1:neqnfl) else if ( s_eqi ( command(8:),'gflsav')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTMP = GFLSAV' end if gfltmp(1:neqnfl) = gflsav(1:neqnfl) else if ( s_eqi ( command(8:),'gflsav-gfltay')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GFLTMP = GFLSAV-GFLTAY' end if gfltmp(1:neqnfl) = gflsav(1:neqnfl) - gfltay(1:neqnfl) else write ( *, * ) 'ARBY4 - Error' write ( *, * ) ' Did not understand your command!' cycle end if ! ! GRB(*) = * ! else if ( s_eqi ( command(1:4),'grb(')) then call chrcti(command(5:),ival,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' ChrCTI returned nonzero error flag!' cycle end if if ( ival < 1 .or. ncofrb < ival ) then write ( *, * ) ' ' write ( *, * ) 'INPUT - Warning!' write ( *, * ) ' Index IVAL of GRB is out of bounds!' write ( *, * ) ' IVAL = ',ival cycle end if call chrctd(command(5+lchar+2:),value,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' ChrCTD returned nonzero error flag!' cycle end if grb(ival) = value if ( 0 < iprint ) then write ( *, * ) 'GRB(',ival,') set to ',grb(ival) end if ! ! GRB = 0, GRBSAV, GRBTAY, TAYLOR ! else if ( s_eqi ( command(1:4),'grb=')) then if ( s_eqi ( command(1:5),'grb=0')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRB = 0' write ( *, * ) ' Set reduced solution to 0.' end if grb(1:ncofrb) = 0.0D+00 else if ( s_eqi ( command,'grb=grbsav')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRB = GRBSAV' write ( *, * ) ' Set reduced solution to saved value.' end if grb(1:ncofrb) = grbsav(1:ncofrb) else if ( s_eqi ( command,'grb=grbtay')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRB = GRBTAY' write ( *, * ) ' Set reduced solution to Taylor base.' end if grb(1:ncofrb) = grbtay(1:ncofrb) else if ( s_eqi ( command,'grb=taylor')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRB = TAYLOR' write ( *, * ) ' Set reduced solution to Taylor prediction,' write ( *, * ) ' based at REYTAY, GRBTAY plus sensitivities.' end if do i = 1, ncofrb grb(i) = 0.0D+00 do j = 1, ntay jtay = j-1 call fact(jtay,factj) temp = ((reynld-reytay)**jtay)/factj grb(i) = grb(i)+temp*senrb(i,j) end do end do else if ( s_eqi ( command(1:5),'grb=(')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRB = (v0,v1,...,vNeqnRB)' end if klo = 6 do i = 1, ncofrb call chrctd(command(klo:),temp,ierror,lchar) if ( ierror /= 0) then write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' There was an error reading your data.' grb(i) = 0.0D+00 else grb(i) = temp end if klo = klo+lchar end do end if ! ! GRBSAV = 0, GRB ! else if ( s_eqi ( command(1:7),'grbsav=')) then if ( command(8:8) == '0') then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRBSAV = 0' end if grbsav(1:ncofrb) = 0.0D+00 else if ( s_eqi ( command,'grbsav=grb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRBSAV = GRB' end if grbsav(1:ncofrb) = grb(1:ncofrb) end if ! ! GRBTAY = 0, GRB, GRBSAV ! else if ( s_eqi ( command(1:7),'grbtay=')) then if ( command(8:8) == '0') then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRBTAY = 0' write ( *, * ) ' Set Taylor base reduced solution to zero.' end if grbtay(1:ncofrb) = 0.0D+00 else if ( s_eqi ( command(8:13),'grbsav')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRBTAY = GRBSAV' write ( *, * ) ' Set Taylor base reduced solution to GRBSAV.' end if grbtay(1:ncofrb) = grbsav(1:ncofrb) else if ( s_eqi ( command(8:10),'grb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - GRBTAY = GRB' write ( *, * ) ' Set Taylor base reduced solution to GRB.' end if grbtay(1:ncofrb) = grb(1:ncofrb) end if ! ! GRIDX = ! else if ( s_eqi ( command(1:5),'gridx')) then if ( s_eqi ( command(1:6),'gridx=')) then gridx = command(7:) else write ( *, * ) 'Enter GRIDX option: UNIFORM, COS, SINSQ:' read(*,'(a)')gridx write(17,'(a)')gridx end if if ( 0 < iprint ) then write ( *, * ) 'The GRIDX option set to '//gridx write ( *, * ) 'Remember to use the SETGEO command' write ( *, * ) 'before trying to solve your system!' end if ! ! GRIDY = ! else if ( s_eqi ( command(1:5),'gridy')) then if ( s_eqi ( command(1:6),'gridy=')) then gridy = command(7:) else write ( *, * ) 'Enter GRIDY option: UNIFORM, COS, SINSQ:' read(*,'(a)')gridy write(17,'(a)')gridy end if if ( 0 < iprint ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'The GRIDY option set to ' // trim ( gridy ) write ( *, '(a)' ) 'Remember to use the SETGEO command' write ( *, '(a)' ) 'before trying to solve your system!' end if ! ! HELLO ! else if ( s_eqi ( command,'hello')) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ARBY4 - Hello' write ( *, '(a)' ) ' ' call hello ( maxnx, maxny ) ! ! HELP ! else if ( s_eqi ( command,'help')) then call help ! ! IBS = ! else if ( s_eqi ( command(1:3),'ibs')) then if ( s_eqi ( command(1:4),'ibs=')) then read(command(5:),*)ibs else write ( *, * ) 'Enter value for IBS:' read(*,*)ibs write(17,*)ibs end if if ( 0 < iprint ) then write ( *, * ) 'IBS set to ',ibs end if ! ! IBUMP = ! else if ( s_eqi ( command(1:5),'ibump')) then if ( s_eqi ( command(1:6),'ibump=')) then read(command(7:),*)ibump else write ( *, * ) 'Enter value for IBUMP:' read(*,*)ibump write(17,*)ibump end if if ( 0 < iprint ) then write ( *, * ) 'IBUMP set to ',ibump end if ! ! IFS = ! else if ( s_eqi ( command(1:3),'ifs')) then if ( s_eqi ( command(1:4),'ifs=')) then read(command(5:),*)ifs else write ( *, * ) 'Enter value for IFS:' read(*,*)ifs write(17,*)ifs end if if ( 0 < iprint ) then write ( *, * ) 'IFS set to ',ifs end if ! ! IHI = ! else if ( s_eqi ( command(1:3),'ihi')) then if ( s_eqi ( command(1:4),'ihi=')) then if ( s_eqi ( command,'ihi=ncofrb')) then ihi = ncofrb else if ( s_eqi ( command,'ihi=neqnfl')) then ihi = neqnfl else if ( s_eqi ( command,'ihi=np')) then ihi = np else read(command(5:),*)ihi end if else write ( *, * ) 'Enter value for IHI:' read(*,*)ihi write(17,*)ihi end if if ( 0 < iprint ) then write ( *, * ) 'IHI set to ',ihi end if ! ! IJAC = ! else if ( s_eqi ( command(1:4),'ijac')) then if ( s_eqi ( command(1:5),'ijac=')) then read(command(6:),*)ijac else write ( *, * ) 'Enter value for IJAC:' read(*,*)ijac write(17,*)ijac end if if ( 0 < iprint ) then write ( *, * ) 'IJAC set to ',ijac end if ! ! ILO = ! else if ( s_eqi ( command(1:3),'ilo')) then if ( s_eqi ( command(1:4),'ilo=')) then read(command(5:),*)ilo else write ( *, * ) 'Enter value for ILO:' read(*,*)ilo write(17,*)ilo end if if ( 0 < iprint ) then write ( *, * ) 'ILO set to ',ilo end if ! ! INIT ! else if ( s_eqi ( command,'init')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Init' write ( *, * ) ' Initialize all data to zero.' write ( *, * ) ' ' end if call init(afl,arb,area,command,cost,costb,costp,costu,costv,difcof, & disfil,drey,epsdif,eqn,etaq,gfl,gflafl,gflrb,gflsav,gflsen,gfltar, & gfltay,grb,grbarb,grbsav,grbsen,grbtay,gridx,gridy,hx,hy,ibs,ibump, & icolrb,ierror,ifs,ihi,ijac,ilo,indx,iopt,ipar,ipivfl,ipivrb,isotri, & iwrite,jhi,jlo,ldafl,maxcofrb,maxelm,maxnew,maxnfl,maxnp,maxny, & maxopt,maxpar,maxparb,maxparf,maxsim,nbcrb,ncofrb,nelem,neqnfl, & nferb,nlband,node,nodelm,np,npar,nparb,nparf,npe,nprof,nsenfl,ntay, & numnew,numopt,numsim,nx,ny,par,parafl,pararb,pardif,parrb,parsav, & parsen,partar,phifl,phirb,rbase,rb,region,resfl,resflsav,resrb, & reynld,reytay,senfl,senrb,splbmp,splflo,taubmp,tauflo,tecfil,tolnew, & tolopt,tolsim,value,wateb,watep,wateu,watev,wquad,xbl,xbr,xc, & xmax,xmin,xprof,xquad,xrange,xsiq,ybl,ybr,yc,ymax,ymin,yquad,yrange) ! ! IOPT(*) = ! else if ( s_eqi ( command(1:5),'iopt(')) then call chrcti(command(6:),ival1,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'INPUT - Fatal error!' write ( *, * ) ' ChrCTI returned nonzero error flag!' stop end if if ( ival1 < 1 .or. maxpar < ival1 ) then write ( *, * ) ' ' write ( *, * ) 'INPUT - Fatal error!' write ( *, * ) ' Index IVAL1 of IOPT is out of bounds!' write ( *, * ) ' IVAL1 = ', ival1 stop end if call chrcti(command(6+lchar+2:),ival2,ierror,lchar) iopt(ival1) = ival2 if ( 0 < iprint ) then write ( *, * ) 'IOPT(',ival1,' ) set to ',ival2 end if ! ! IWRITE = ! else if ( s_eqi ( command(1:6),'iwrite')) then if ( s_eqi ( command(1:7),'iwrite=')) then read(command(8:),*)iwrite else write ( *, * ) 'Enter value for IWRITE:' read(*,*)iwrite write(17,*)iwrite end if if ( 0 < iprint ) then write ( *, * ) 'IWRITE set to ',iwrite end if ! ! JHI = ! else if ( s_eqi ( command(1:3),'jhi')) then if ( s_eqi ( command(1:4),'jhi=')) then if ( s_eqi ( command,'jhi=ncofrb')) then jhi = ncofrb else if ( s_eqi ( command,'jhi=neqnfl')) then jhi = neqnfl else if ( s_eqi ( command,'jhi=nsenfl')) then jhi = nsenfl else read(command(5:),*)jhi end if else write ( *, * ) 'Enter value for JHI:' read(*,*)jhi write(17,*)jhi end if if ( 0 < iprint ) then write ( *, * ) 'JHI set to ',jhi end if ! ! JLO = ! else if ( s_eqi ( command(1:3),'jlo')) then if ( s_eqi ( command(1:4),'jlo=')) then read(command(5:),*)jlo else write ( *, * ) 'Enter value for JLO:' read(*,*)jlo write(17,*)jlo end if if ( 0 < iprint ) then write ( *, * ) 'JLO set to ',jlo end if ! ! L2NORM GFL/GFLSAV/GFLTAR/GFLTAY/GFLTMP ! else if ( s_eqi ( command(1:6),'l2norm')) then if ( s_eqi ( command(8:),'gfl')) then call l2norm(gfl,gflnrm,indx,nelem,neqnfl,node,np,xc,yc) write ( *, * ) 'ARBY4 - L2Norm of GFL = ',gflnrm else if ( s_eqi ( command(8:),'gflsav')) then call l2norm(gflsav,gflnrm,indx,nelem,neqnfl,node,np,xc,yc) write ( *, * ) 'ARBY4 - L2Norm of GFLSAV = ',gflnrm else if ( s_eqi ( command(8:),'gfltar')) then call l2norm(gfltar,gflnrm,indx,nelem,neqnfl,node,np,xc,yc) write ( *, * ) 'ARBY4 - L2Norm of GFLTAR = ',gflnrm else if ( s_eqi ( command(8:),'gfltay')) then call l2norm(gfltay,gflnrm,indx,nelem,neqnfl,node,np,xc,yc) write ( *, * ) 'ARBY4 - L2Norm of GFLTAY = ',gflnrm else if ( s_eqi ( command(8:),'gfltmp')) then call l2norm(gfltmp,gflnrm,indx,nelem,neqnfl,node,np,xc,yc) write ( *, * ) 'ARBY4 - L2Norm of GFLTMP = ',gflnrm else write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Error!' write ( *, * ) ' Legal choices were GFL/GFLSAV/GFLTAY/GFLTMP.' write ( *, * ) ' Your choice was '//command(8:) cycle end if ! ! MAXNEW = ! else if ( s_eqi ( command(1:6),'maxnew')) then if ( s_eqi ( command(1:7),'maxnew=')) then read(command(8:),*)maxnew else write ( *, * ) 'Enter value for MAXNEW:' read(*,*)maxnew write(17,*)maxnew end if if ( 0 < iprint ) then write ( *, * ) 'MAXNEW set to ',maxnew end if ! ! MAXOPT = ! else if ( s_eqi ( command(1:6),'maxopt')) then if ( s_eqi ( command(1:7),'maxopt=')) then read(command(8:),*)maxopt else write ( *, * ) 'Enter value for MAXOPT:' read(*,*)maxopt write(17,*)maxopt end if if ( 0 < iprint ) then write ( *, * ) 'MAXOPT set to ',maxopt end if ! ! MAXSIM = ! else if ( s_eqi ( command(1:6),'maxsim')) then if ( s_eqi ( command(1:7),'maxsim=')) then read(command(8:),*)maxsim else write ( *, * ) 'Enter value for MAXSIM:' read(*,*)maxsim write(17,*)maxsim end if if ( 0 < iprint ) then write ( *, * ) 'MAXSIM set to ',maxsim end if ! ! NBCRB = ! else if ( s_eqi ( command(1:5),'nbcrb')) then if ( s_eqi ( command(1:6),'nbcrb=')) then read(command(7:),*)nbcrb else write ( *, * ) 'Enter value for NBCRB:' read(*,*)nbcrb write(17,*)nbcrb end if if ( 0 < iprint ) then write ( *, * ) 'NBCRB set to ',nbcrb end if ! ! NEWTFL ! Apply Newton's method to full solution estimate. ! else if ( s_eqi ( command,'newtfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - NewtFL' write ( *, * ) ' Apply Newton to full solution estimate GFL.' end if call newtfl(afl,area,eqn,gfl,gflafl,ierror,ifs,ijac,indx,ipivfl,iwrite, & ldafl,maxelm,maxnew,nelem,neqnfl,nlband,node,np,npar,nparf,numnew,par, & parafl,phifl,region,resfl,rmax,splflo,tauflo,tolnew,xrange,yc,yrange) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Fatal error!' write ( *, * ) ' NEWTFL failed!' write ( *, * ) ' The parameters at which failure occurred:' write ( *, * ) ' ' call prpar(iopt,npar,nparb,nparf,par) ierror = 1 cycle else if ( iwrite <= 1) then write ( *, * ) ' Newton step ',numsim,' residual norm = ',rmax end if end if ! ! NEWTRB ! Apply Newton's method to reduced solution estimate. ! else if ( s_eqi ( command,'newtrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - NewtRB' write ( *, * ) ' Apply Newton to reduced solution estimate GRB.' end if if ( ncofrb <= 0 ) then cycle end if call newtrb(arb,area,grb,grbarb,ierror,indx,ipivrb, & iwrite,maxcofrb,maxelm,maxnew,maxnfl,nbcrb,ncofrb,nelem, & nferb,node,np,npar,nparf,nx,ny,par,pararb,phirb, & rb,resrb,rmax,tauflo,tolnew,xc,xrange,yc,yrange) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Fatal error!' write ( *, * ) ' NEWTRB failed!' write ( *, * ) ' The parameters at which failure occurred:' call prpar(iopt,npar,nparb,nparf,par) ierror = 1 cycle else if ( iwrite <= 1) then write ( *, * ) ' Final Newton residual was MxNorm(FXRB) = ',rmax end if end if ! ! NPARB = ! else if ( s_eqi ( command(1:5),'nparb')) then if ( s_eqi ( command(1:6),'nparb=')) then read(command(7:),*)nparb else write ( *, * ) 'Enter value for NPARB:' read(*,*)nparb write(17,*)nparb end if if ( 0 < iprint ) then write ( *, * ) 'NPARB set to ',nparb end if ! ! NPARF = ! else if ( s_eqi ( command(1:5),'nparf')) then if ( s_eqi ( command(1:6),'nparf=')) then read(command(7:),*)nparf else write ( *, * ) 'Enter value for NPARF:' read(*,*)nparf write(17,*)nparf end if if ( 0 < iprint ) then write ( *, * ) 'NPARF set to ',nparf end if ! ! NSENFL = # ! else if ( s_eqi ( command(1:6),'nsenfl')) then if ( s_eqi ( command(1:7),'nsenfl=')) then read(command(8:),*)itemp else write ( *, * ) 'Enter value for NSENFL:' read(*,*)itemp write(17,*)itemp end if if ( itemp < 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' NSENFL must be at least 0!' write ( *, * ) ' but your value was ',itemp else nsenfl = itemp if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - NSENFL set to ',nsenfl end if end if ! ! NTAY = # ! NTAY = NCOFRB ! else if ( s_eqi ( command(1:4),'ntay')) then if ( s_eqi ( command(1:5),'ntay=')) then if ( s_eqi ( command(6:11),'ncofrb')) then itemp = ncofrb else read(command(6:),*)itemp end if else write ( *, * ) 'Enter value for NTAY:' read(*,*)itemp write(17,*)itemp end if if ( itemp < 0 .or. ncofrb < itemp ) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' NTAY must be between 0 and NCOFRB = ',ncofrb write ( *, * ) ' but your value was ',itemp else ntay = itemp if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - NTAY set to ',ntay end if end if ! ! NX = ! else if ( s_eqi ( command(1:2),'nx')) then if ( s_eqi ( command(1:3),'nx=')) then read(command(4:),*)itemp else write ( *, * ) 'Enter value for NX:' read(*,*)itemp write(17,*)itemp end if if ( itemp < 2 ) then write ( *, * ) 'ARBY4 - Unacceptable input.' write ( *, * ) ' NX must be at least 2.' write ( *, * ) ' Your value was ',itemp else if ( maxnx < itemp ) then write ( *, * ) 'ARBY4 - Unacceptable input.' write ( *, * ) ' NX must be no more than MAXNX = ',maxnx write ( *, * ) ' Your value was ',itemp else nx = itemp if ( 0 < iprint ) then write ( *, * ) 'NX set to ',nx write ( *, * ) 'Remember to use the SETLOG and SETGEO commands' write ( *, * ) 'before trying to solve your systems!' end if end if ! ! NY = ! else if ( s_eqi ( command(1:2),'ny')) then if ( s_eqi ( command(1:3),'ny=')) then read(command(4:),*)itemp else write ( *, * ) 'Enter value for NY:' read(*,*)itemp write(17,*)itemp end if if ( itemp < 2) then write ( *, * ) 'ARBY4 - Unacceptable input.' write ( *, * ) ' NY must be at least 2.' write ( *, * ) ' Your value was ',itemp else if ( maxny < itemp ) then write ( *, * ) 'ARBY4 - Unacceptable input.' write ( *, * ) ' NY must be no more than MAXNY = ',maxny write ( *, * ) ' Your value was ',itemp else ny = itemp if ( 0 < iprint ) then write ( *, * ) 'NY set to ',ny write ( *, * ) 'Remember to use the SETLOG and SETGEO commands' write ( *, * ) 'before trying to solve your system!' end if end if ! ! OPTFL ! else if ( s_eqi ( command,'optfl')) then write ( *, * ) 'ARBY4 - OPT FL:' write ( *, * ) ' Optimize the full system.' write ( *, * ) ' NO WAY, JOSE! NOT READY YET!' ! ! OPTDIFFL ! else if ( s_eqi ( command,'optdiffl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - OptDifFl:' write ( *, * ) ' Optimize the cost of the full system;' write ( *, * ) ' The optimization code will approximate cost' write ( *, * ) ' gradients by finite differences.' write ( *, * ) ' Initial estimate is (PAR,GFL,COST).' write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Note!' write ( *, * ) ' You must already have issued the TARGET command!' end if call optdiffl(afl,area,cost,dopt,eqn,etaq,gfl,gflafl, & gflopt,gfltar,gridx,gridy,ibs,ierror,ifs,ijac,indx,iopt, & ipivfl,isotri,ivopt,iwrite,ldafl,liv,lv,maxelm,maxnew,maxnfl, & maxnp,maxny,maxopt,maxpar,maxparb,maxparf,maxsim,nelem, & neqnfl,nlband,node,nodelm,np,npar,nparb,nparf,nprof,numdif, & numopt,nx,ny,par,parafl,paropt,phifl,region,resfl,splbmp, & splflo,taubmp,tauflo,tolnew,tolopt,tolsim,vopt,wateb,watep, & wateu,watev,wquad,xbl,xbr,xc,xopt,xquad,xrange,xsiq,ybl, & ybr,yc,yquad,yrange) if ( ierror == 0) then par(1:npar) = paropt(1:npar) gfl(1:neqnfl) = gflopt(1:neqnfl) write ( *, * ) ' ' write ( *, * ) 'Optimizing parameters:' call prpar(iopt,npar,nparb,nparf,par) write ( *, * ) ' ' write ( *, * ) 'Optimal cost = ',cost write ( *, * ) ' ' write ( *, * ) 'Number of standard full solutions:',numopt write ( *, * ) 'Number of auxilliary solutions: ',numdif else write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' The optimization was unsuccessful.' end if ! ! OPTDIFRB ! else if ( s_eqi ( command,'optdifrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - OptDifRB:' write ( *, * ) ' Optimize the cost of the reduced system;' write ( *, * ) ' The optimization code will approximate cost' write ( *, * ) ' gradients by finite differences.' write ( *, * ) ' Initial estimate is (PAR,GRB,COST).' write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Note!' write ( *, * ) ' You must already have issued the commands:' write ( *, * ) ' SETLOG, SETGEO, NEWTFL, GETSEN, GETRB, TARGET!' end if call optdifrb(arb,area,cost,dopt,gflrb,gfltar,gfltmp, & grb,grbarb,grbopt,ierror,indx,iopt,ipivrb, & ivopt,iwrite,liv,lv,maxcofrb,maxelm,maxnew,maxnfl,maxnp, & maxny,maxopt,maxpar,maxparb,maxsim,nbcrb,ncofrb,nelem,neqnfl, & nferb,node,np,npar,nparb,nparf,nprof,numdif,numopt,nx,ny,par, & pararb,paropt,phirb,rb,resrb,splbmp,tauflo,taubmp,tolnew, & tolopt,tolsim,vopt,wateb,watep,wateu,watev,xbl,xbr,xc,xopt, & xrange,ybl,ybr,yc,yrange) if ( ierror == 0) then par(1:npar) = paropt(1:npar) grb(1:ncofrb) = grbopt(1:ncofrb) write ( *, * ) ' ' write ( *, * ) 'Optimizing parameters:' call prpar(iopt,npar,nparb,nparf,par) write ( *, * ) ' ' write ( *, * ) 'Optimal cost = ',cost write ( *, * ) ' ' write ( *, * ) 'Number of standard full solutions:',numopt write ( *, * ) 'Number of auxilliary solutions: ',numdif else write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' The optimization was unsuccessful.' end if ! ! PAR(*) = * ! else if ( s_eqi ( command(1:4),'par(')) then call chrcti(command(5:),ival,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' ChrCTI returned nonzero error flag!' cycle end if if ( ival < 1 .or. maxpar < ival ) then write ( *, * ) ' ' write ( *, * ) 'INPUT - Warning!' write ( *, * ) ' Index IVAL of PAR is out of bounds!' write ( *, * ) ' IVAL = ',ival cycle end if call chrctd(command(5+lchar+2:),value,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' ChrCTD returned nonzero error flag!' cycle end if par(ival) = value if ( 0 < iprint ) then write ( *, * ) 'PAR(',ival,') set to ',par(ival) end if ! ! PARTAR(*) = * ! else if ( s_eqi ( command(1:7),'partar(')) then call chrcti(command(8:),ival,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' ChrCTI returned nonzero error flag!' cycle end if if ( ival < 1 .or. maxpar < ival ) then write ( *, * ) ' ' write ( *, * ) 'INPUT - Warning!' write ( *, * ) ' Index IVAL of PARTAR is out of bounds!' write ( *, * ) ' IVAL = ', ival cycle end if call chrctd(command(8+lchar+2:),value,ierror,lchar) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' ChrCTD returned nonzero error flag!' cycle end if partar(ival) = value par(ival) = value if ( 0 < iprint ) then write ( *, * ) 'PARTAR(',ival,') set to ',partar(ival) end if ! ! PICFL ! else if ( s_eqi ( command,'picfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - PicFL:' write ( *, * ) ' Apply Picard to full solution estimate GFL.' end if call picfl(afl,area,eqn,gfl,ierror,ifs,indx,ipivfl,iwrite,ldafl,maxsim, & nelem,neqnfl,nlband,node,np,npar,nparf,numsim,par,phifl,region,resfl, & rmax,splflo,tauflo,tolsim,xc,xrange,yc,yrange) if ( iwrite <= 1) then write ( *, * ) ' Picard step ',numsim,' residual norm = ',rmax end if ! ! PICRB ! else if ( s_eqi ( command,'picrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - PicRB:' write ( *, * ) ' Apply Picard to reduced solution estimate GRB.' end if if ( ncofrb <= 0) then cycle end if call picrb(arb,area,grb,ierror,indx,ipivrb,iwrite,maxcofrb,maxelm,maxnfl, & maxsim,nbcrb,ncofrb,nelem,nferb,node,np,npar,nparf,nx,ny,par,phirb,rb, & resrb,rmax,tauflo,tolsim,xc,xrange,yc,yrange) if ( iwrite <= 1) then write ( *, * ) ' Final Picard residual was MxNorm(FXRB) = ',rmax end if ! ! PRFPFL ! else if ( s_eqi ( command,'prfpfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr FP FL:' write ( *, * ) ' Print full jacobian.' write ( *, * ) ' Rows ILO = ',ilo,' to IHI=',ihi write ( *, * ) ' Cols JLO = ',jlo,' to JHI=',jhi write ( *, * ) ' ' write ( *, * ) ' Parameters for matrix, PARAFL:' call prpar(iopt,npar,nparb,nparf,parafl) end if call prbmat(afl,ihi,ilo,jhi,jlo,ldafl,neqnfl,nlband) ! ! PRFPRB ! else if ( s_eqi ( command,'prfprb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr FP RB:' write ( *, * ) ' Print reduced jacobian.' write ( *, * ) ' Rows ILO = ',ilo,' to IHI=',ihi write ( *, * ) ' Cols JLO = ',jlo,' to JHI=',jhi write ( *, * ) ' ' write ( *, * ) ' Parameters for matrix, PARARB:' call prpar(iopt,npar,nparb,nparf,pararb) end if mlo = 1 mhi = maxcofrb nlo = 1 nhi = maxcofrb call prdmat(arb,ihi,ilo,jhi,jlo,mhi,mlo,nhi,nlo) ! ! PRDAT ! else if ( s_eqi ( command,'prdat')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr Dat' write ( *, * ) ' Print current problem data.' end if call prdat(disfil,drey,epsdif,gridx,gridy,hx,hy,ibs,ibump,ifs,ijac,iopt, & maxnew,maxopt,maxsim,nbcrb,ncofrb,nelem,nferb,neqnfl,np,npar,nparb, & nparf,ntay,nx,ny,region,reytay,tecfil,tolnew,tolopt,tolsim,wateb, & watep,wateu,watev,xbl,xbr,xprof,xrange,ybl,ybr,yrange) ! ! PRELEM ! else if ( s_eqi ( command,'prelem')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr Elem' write ( *, * ) ' Print element data.' end if call prelem(ihi,ilo,nelem,node,np,xc,yc) ! ! PRFXFL ! else if ( s_eqi ( command,'prfxfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr FX FL' write ( *, * ) ' Print full residual FXFL.' end if call prvecfl(eqn,ihi,ilo,indx,neqnfl,np,resfl) ! ! PRFXFLNRM ! else if ( s_eqi ( command,'prfxflnrm')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr FX FL Nrm' write ( *, * ) ' Print norm of full residual FXFL.' end if call prfxfln(neqnfl,resfl) ! ! PRFXRB ! else if ( s_eqi ( command,'prfxrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr FX RB' write ( *, * ) ' Print reduced residual FXRB.' end if nlo = 1 nhi = ncofrb call prvecrb(ihi,ilo,nhi,nlo,resrb) ! ! PRGFL ! else if ( s_eqi ( command,'prgfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr G FL:' write ( *, * ) ' Print full solution GFL.' write ( *, * ) ' ' write ( *, * ) ' Flow parameters, PAR:' call prpar(iopt,npar,nparb,nparf,par) end if call prvecfl(eqn,ihi,ilo,indx,neqnfl,np,gfl) ! ! PRGFLNRM ! else if ( s_eqi ( command,'prgflnrm')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr GFL Nrm:' write ( *, * ) ' Print norms of full solution GFL.' end if call fxfl(area,eqn,gfl,ifs,indx,nelem,neqnfl,node,np,npar,nparf,par, & phifl,region,resfl,splflo,tauflo,xrange,yc,yrange) call nrmflo(gfl,indx,neqnfl,np,resfl) ! ! PRGRB ! else if ( s_eqi ( command,'prgrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr G RB:' write ( *, * ) ' Print reduced solution GRB.' end if call prgrb(grb,ncofrb) ! ! PRGSEN ! else if ( s_eqi ( command,'prgsen')) then call prgrb(gsen,nbcrb+nsenfl) ! ! PRINDX ! else if ( s_eqi ( command,'prindx')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr INDX' write ( *, * ) ' Print node/equation table,' write ( *, * ) ' for nodes ILO = ',ilo,' to IHI=',ihi end if call prindx(ihi,ilo,indx,np,xc,yc) ! ! PRPAR ! else if ( s_eqi ( command,'prpar')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr PAR: Print current parameters PAR.' end if call prpar(iopt,npar,nparb,nparf,par) ! ! PRPARTAR ! else if ( s_eqi ( command,'prpartar')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr PAR TAR:' write ( *, * ) ' Print target parameters PARTAR.' end if call prpar(iopt,npar,nparb,nparf,partar) ! ! PRRBASE ! else if ( s_eqi ( command,'prrbase')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr RBase' write ( *, * ) ' Print the "R" factor of the reduced basis.' end if ilo = 1 ihi = ncofrb jlo = 1 jhi = ncofrb mlo = 1 mhi = maxcofrb nlo = 1 nhi = maxcofrb call prdmat(rbase,ihi,ilo,jhi,jlo,mhi,mlo,nhi,nlo) ! ! PRRB ! else if ( s_eqi ( command,'prrb')) then if ( ncofrb <= 0) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' PRRB command cancelled, NCOFRB = ',ncofrb write ( *, * ) ' Use the GETRB command first!' cycle end if if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr RB' write ( *, * ) ' Print the reduced basis,' write ( *, * ) ' nodes ILO = ',ilo,' to IHI=',ihi write ( *, * ) ' columns JLO = ',jlo,' to JHI=',jhi write ( *, * ) ' ' write ( *, * ) ' Parameters at reduced basis, PARRB:' call prpar(iopt,npar,nparb,nparf,parrb) end if call prmatfl(rb,eqn,ihi,ilo,indx,jhi,jlo,maxnfl,ncofrb,neqnfl,np) ! ! PRSENFL: Print full sensitivities. ! else if ( s_eqi ( command,'prsenfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr Sen FL' write ( *, * ) ' Print full sensitivities.' write ( *, * ) ' ' write ( *, * ) ' Parameters at sensitivity, PARSEN:' call prpar(iopt,npar,nparb,nparf,parsen) end if call prmatfl(senfl,eqn,ihi,ilo,indx,jhi,jlo,maxnfl,nsenfl,neqnfl,np) ! ! PRSENNRM: Print full sensitivity norms. ! else if ( s_eqi ( command,'prsennrm')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr Sen Nrm' write ( *, * ) ' Print sensitivity norms.' end if call prsenn(maxcofrb,maxnfl,ncofrb,neqnfl,senfl) ! ! PRSENRB: Print reduced sensitivities. ! else if ( s_eqi ( command,'prsenrb')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr Sen RB' write ( *, * ) ' Print matrix of reduced sensitivities,' write ( *, * ) ' rows ILO = ',ilo,' to IHI=',ihi write ( *, * ) ' columns JLO = ',jlo,' to JHI=',jhi write ( *, * ) ' ' write ( *, * ) ' Parameters at sensitivity, PARSEN:' call prpar(iopt,npar,nparb,nparf,parsen) end if mlo = 1 mhi = maxcofrb nlo = 1 nhi = ncofrb call prdmat(senrb,ihi,ilo,jhi,jlo,mhi,mlo,nhi,nlo) ! ! PRUVPGFL ! else if ( s_eqi ( command,'pruvpgfl')) then call pruvpfl(gfl,indx,neqnfl,np,xc,xmax,xmin,yc,ymax,ymin) ! ! PRUVPRB ! else if ( s_eqi ( command,'pruvprb')) then do j = 1, ncofrb write ( *, * ) ' ' write ( *, * ) 'Reduced basis vector ',j write ( *, * ) ' ' call pruvpfl(rb(1,j),indx,neqnfl,np,xc,xmax,xmin,yc,ymax,ymin) end do ! ! PRUVPSENFL ! else if ( s_eqi ( command,'pruvpsenfl')) then do j = 1, nsenfl write ( *, * ) ' ' write ( *, * ) 'Sensitivity vector ',j write ( *, * ) ' ' call pruvpfl(senfl(1,j),indx,neqnfl,np,xc,xmax,xmin,yc,ymax,ymin) end do ! ! PRUVPGRB ! else if ( s_eqi ( command,'pruvpgrb')) then call pruvprb(grb,indx,maxnfl,ncofrb,nelem,node,nodelm,np, & rb,xc,xmax,xmin,yc,ymax,ymin) ! ! PRXY ! else if ( s_eqi ( command,'prxy')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Pr XY' write ( *, * ) ' Print out X and Y nodal coordinates.' end if call prxy(ihi,ilo,np,ny,xc,yc) ! ! QUIT ! else if ( s_eqi ( command, 'quit' ) ) then exit ! ! REDUCE GFL ! else if ( s_eqi ( command,'reduce gfl')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Reduce GFL:' write ( *, * ) ' Given a reduced basis RB computed at the' write ( *, * ) ' full solution GFLRB, and an arbitrary full' write ( *, * ) ' solution GFL, compute the reduced basis ' write ( *, * ) ' coefficients of GFL:' write ( *, * ) ' GRB = RB^T * GFL.' end if call gfl2rb(gfl,gflrb,grb,maxnfl,ncofrb,neqnfl,rb) ! ! REGION = ! else if ( s_eqi ( command(1:6),'region')) then if ( s_eqi ( command(1:7),'region=')) then region = command(8:) else write ( *, * ) 'Enter the region, CAVITY, CHANNEL or STEP:' read(*,'(a)')region write(17,'(a)')region end if if ( s_eqi ( region,'cavity')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Cavity:' write ( *, * ) ' Set user input values to cavity defaults.' end if call cavity(ibs,ibump,ifs,iopt,maxopt,maxpar,nbcrb,npar,nparb,nparf, & npe,nx,ny,par,region,reynld,tolnew,tolopt,tolsim,wateb,watep,wateu, & watev,xbl,xbr,xprof,xrange,ybl,ybr,yrange) else if ( s_eqi ( region,'cavity2')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Cavity2:' write ( *, * ) ' Set H C Lee cavity defaults.' end if call cavity2(ibs,ibump,ifs,iopt,maxopt,maxpar,nbcrb,npar,nparb,nparf, & npe,nx,ny,par,region,reynld,tolnew,tolopt,tolsim,wateb,watep,wateu, & watev,xbl,xbr,xprof,xrange,ybl,ybr,yrange) else if ( s_eqi ( region,'channel')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Channel:' write ( *, * ) ' Set user input values to channel defaults.' end if call channl(ibs,ibump,ifs,iopt,maxopt,maxpar,nbcrb,npar,nparb,nparf, & npe,nx,ny,par,region,reynld,tolnew,tolopt,tolsim,wateb,watep,wateu, & watev,xbl,xbr,xprof,xrange,ybl,ybr,yrange) else if ( s_eqi ( region,'step')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Step`:' write ( *, * ) ' Set user input values to step defaults.' end if call step(ibs,ibump,ifs,iopt,maxopt,maxpar,nbcrb,npar,nparb,nparf,npe, & nx,ny,par,region,reynld,tolnew,tolopt,tolsim,wateb,watep,wateu,watev, & xbl,xbr,xprof,xrange,ybl,ybr,yrange) end if ! ! REYNLD = ! else if ( s_eqi ( command(1:6),'reynld')) then if ( s_eqi ( command(1:7),'reynld=')) then read(command(8:),*)reynld else write ( *, * ) 'Enter value for REYNLD:' read(*,*)reynld write(17,*)reynld end if par(nparf+nparb+1) = reynld if ( 0 < iprint ) then write ( *, * ) 'REYNLD parameter set to ',reynld end if ! ! REYTAY = ! else if ( s_eqi ( command(1:6),'reytay')) then if ( s_eqi ( command(1:7),'reytay=')) then if ( s_eqi ( command,'reytay=reynld')) then reytay = reynld else read(command(8:),*)reytay end if else write ( *, * ) 'Enter value for REYTAY:' read(*,*)reytay write(17,*)reytay end if if ( 0 < iprint ) then write ( *, * ) 'REYTAY parameter set to ',reytay end if ! ! SETGEO ! else if ( s_eqi ( command,'setgeo')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - SetGeo: Set problem geometry.' end if call setgeo(area,etaq,gridx,gridy,ibs,isotri,nelem,node,nodelm,np,npar, & nparb,nparf,nx,ny,par,phifl,region,splbmp,taubmp,wquad,xbl,xbr,xc, & xquad,xrange,xsiq,ybl,ybr,yc,yquad,yrange) ! ! SETLOG ! else if ( s_eqi ( command,'setlog')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - SetLog: Set problem logical data.' end if call setlog(eqn,hx,hy,ibump,indx,isotri,ldafl,maxelm,maxnfl,maxnp,nelem, & neqnfl,nlband,node,np,nprof,nx,ny,region,xbl,xbr,xprof,xrange,ybr,yrange) ! ! STOP ! else if ( s_eqi ( command,'stop') ) then exit ! ! TARGET ! else if ( s_eqi ( command,'target')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Target:' write ( *, * ) ' Save current GFL as GTAR.' end if call target(cost0,gfl,gfltar,indx,maxnfl,maxny,maxparb,neqnfl,np,npar, & nparb,nprof,ny,par,partar,splbmp,taubmp,wateb,watep,wateu,watev,xbl, & xbr,ybl,ybr,yc) ! ! TEST2 ! else if ( s_eqi ( command,'test2')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Test2:' write ( *, * ) ' Compare full and reduced state variables' write ( *, * ) ' in elements ILO through IHI.' end if call test2(gfl,grb,ihi,ilo,indx,maxcofrb,maxelm,ncofrb, & nelem,neqnfl,node,np,phifl,phirb) ! ! TEST3 ! else if ( s_eqi ( command,'test3')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Test3:' write ( *, * ) ' Compare RB*Rfact and SenFL.' end if call test3(maxcofrb,maxnfl,ncofrb,neqnfl,rb,senfl,senrb) ! ! TEST4 ! else if ( s_eqi ( command,'test4')) then if ( ipar <= 0 .or. npar < ipar ) then write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' Cancelling the TEST4 command.' write ( *, * ) ' IPAR = ',ipar write ( *, * ) ' but IPAR must be at least 0' write ( *, * ) ' and no more than NPAR = ',npar cycle end if if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Test4' write ( *, * ) ' Compare full sensitivities computed directly' write ( *, * ) ' and via finite differences.' end if call test4(afl,area,difcof,dpar,drey,eqn,gfl,gflafl, & ifs,ijac,indx,ipar,ipivfl,iwrite,ldafl,maxcofrb,maxelm, & maxnew,maxnfl,ncofrb,nelem,neqnfl,nlband,node,np,npar, & nparf,nsenfl,par,parafl,phifl,region,resfl, & senfl,splflo,tauflo,tolnew,xrange,yc,yrange) ! ! TEST5 ! else if ( s_eqi ( command,'test5')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Test5:' write ( *, * ) ' Compare RB*Rfact and basis vectors.' end if call test5(maxcofrb,maxnfl,ncofrb,neqnfl,rb,rbase) ! ! TIME ! else if ( s_eqi ( command,'time')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - Time:' write ( *, * ) ' Report current time,' write ( *, * ) ' time elapsed since last TIME call,' write ( *, * ) ' time elapsed since the program began.' end if write ( *, * ) ' The (real) start time was '// trim ( tstart ) call date_and_time ( date, time ) write ( *, * ) ' The current (real) time is ' // trim ( time ) ! ! TOLNEW = ! else if ( s_eqi ( command(1:6),'tolnew')) then if ( s_eqi ( command(1:7),'tolnew=')) then read(command(8:),*)tolnew else write ( *, * ) 'Enter value for TOLNEW:' read(*,*)tolnew write(17,*)tolnew end if if ( 0 < iprint ) then write ( *, * ) 'TOLNEW set to ',tolnew end if ! ! TOLOPT = ! else if ( s_eqi ( command(1:6),'tolopt')) then if ( s_eqi ( command(1:7),'tolopt=')) then read(command(8:),*)tolopt else write ( *, * ) 'Enter value for TOLOPT:' read(*,*)tolopt write(17,*)tolopt end if if ( 0 < iprint ) then write ( *, * ) 'TOLOPT set to ',tolopt end if ! ! TOLSIM = ! else if ( s_eqi ( command(1:6),'tolsim')) then if ( s_eqi ( command(1:7),'tolsim=')) then read(command(8:),*)tolsim else write ( *, * ) 'Enter value for TOLSIM:' read(*,*)tolsim write(17,*)tolsim end if if ( 0 < iprint ) then write ( *, * ) 'TOLSIM set to ',tolsim end if ! ! TecPlot ! else if ( s_eqi ( command,'tecplot')) then if ( 0 < iprint ) then write ( *, * ) 'ARBY4 - TecPlot:' write ( *, * ) ' Write data to TECPLOT plot file.' end if call intprs(gfl,indx,nelem,neqnfl,node,np,p) u(1:np) = gfl(indx(1,1:np)) v(1:np) = gfl(indx(2,1:np)) call wrtec(nelem,node,np,p,tecfil,u,v,xc,yc) ! ! WATEB = ! else if ( s_eqi ( command(1:5),'wateb')) then if ( s_eqi ( command(1:6),'wateb=')) then read(command(7:),*)wateb else write ( *, * ) 'Enter value for WATEB:' read(*,*)wateb write(17,*)wateb end if if ( 0 < iprint ) then write ( *, * ) 'WATEB set to ',wateb end if ! ! WATEP = ! else if ( s_eqi ( command(1:5),'watep')) then if ( s_eqi ( command(1:6),'watep=')) then read(command(7:),*)watep else write ( *, * ) 'Enter value for WATEP:' read(*,*)watep write(17,*)watep end if if ( 0 < iprint ) then write ( *, * ) 'WATEP set to ',watep end if ! ! WATEU = ! else if ( s_eqi ( command(1:5),'wateu')) then if ( s_eqi ( command(1:6),'wateu=')) then read(command(7:),*)wateu else write ( *, * ) 'Enter value for WATEU:' read(*,*)wateu write(17,*)wateu end if if ( 0 < iprint ) then write ( *, * ) 'WATEU set to ',wateu end if ! ! WATEV = ! else if ( s_eqi ( command(1:5),'watev')) then if ( s_eqi ( command(1:6),'watev=')) then read(command(7:),*)watev else write ( *, * ) 'Enter value for WATEV:' read(*,*)watev write(17,*)watev end if if ( 0 < iprint ) then write ( *, * ) 'WATEV set to ',watev end if ! ! XBL = ! else if ( s_eqi ( command(1:3),'xbl')) then if ( s_eqi ( command(1:4),'xbl=')) then read(command(5:),*)xbl else write ( *, * ) 'Enter value for XBL:' read(*,*)xbl write(17,*)xbl end if if ( 0 < iprint ) then write ( *, * ) 'XBL set to ',xbl end if ! ! XBR = ! else if ( s_eqi ( command(1:3),'xbr')) then if ( s_eqi ( command(1:4),'xbr=')) then read(command(5:),*)xbr else write ( *, * ) 'Enter value for XBR:' read(*,*)xbr write(17,*)xbr end if if ( 0 < iprint ) then write ( *, * ) 'XBR set to ',xbr end if ! ! XMAX = ! else if ( s_eqi ( command(1:4),'xmax')) then if ( s_eqi ( command(1:5),'xmax=')) then read(command(6:),*)xmax else write ( *, * ) 'Enter value for XMAX:' read(*,*)xmax write(17,*)xmax end if if ( 0 < iprint ) then write ( *, * ) 'XMAX set to ',xmax end if ! ! XMIN = ! else if ( s_eqi ( command(1:4),'xmin')) then if ( s_eqi ( command(1:5),'xmin=')) then read(command(6:),*)xmin else write ( *, * ) 'Enter value for XMIN:' read(*,*)xmin write(17,*)xmin end if if ( 0 < iprint ) then write ( *, * ) 'XMIN set to ',xmin end if ! ! XPROF = ! else if ( s_eqi ( command(1:5),'xprof')) then if ( s_eqi ( command(1:6),'xprof=')) then read(command(7:),*)xprof else write ( *, * ) 'Enter value for XPROF:' read(*,*)xprof write(17,*)xprof end if if ( 0 < iprint ) then write ( *, * ) 'XPROF set to ',xprof end if ! ! XRANGE = ! else if ( s_eqi ( command(1:6),'xrange')) then if ( s_eqi ( command(1:7),'xrange=')) then read(command(8:),*)xrange else write ( *, * ) 'Enter value for XRANGE:' read(*,*)xrange write(17,*)xrange end if if ( 0 < iprint ) then write ( *, * ) 'XRANGE set to ',xrange end if ! ! YBL = ! else if ( s_eqi ( command(1:3),'ybl')) then if ( s_eqi ( command(1:4),'ybl=')) then read(command(5:),*)ybl else write ( *, * ) 'Enter value for YBL:' read(*,*)ybl write(17,*)ybl end if if ( 0 < iprint ) then write ( *, * ) 'YBL set to ',ybl end if ! ! YBR = ! else if ( s_eqi ( command(1:3),'ybr')) then if ( s_eqi ( command(1:4),'ybr=')) then read(command(5:),*)ybr else write ( *, * ) 'Enter value for YBR:' read(*,*)ybr write(17,*)ybr end if if ( 0 < iprint ) then write ( *, * ) 'YBR set to ',ybr end if ! ! YMAX = ! else if ( s_eqi ( command(1:4),'ymax')) then if ( s_eqi ( command(1:5),'ymax=')) then read(command(6:),*)ymax else write ( *, * ) 'Enter value for YMAX:' read(*,*)ymax write(17,*)ymax end if if ( 0 < iprint ) then write ( *, * ) 'YMAX set to ',ymax end if ! ! YMIN = ! else if ( s_eqi ( command(1:4),'ymin')) then if ( s_eqi ( command(1:5),'ymin=')) then read(command(6:),*)ymin else write ( *, * ) 'Enter value for YMIN:' read(*,*)ymin write(17,*)ymin end if if ( 0 < iprint ) then write ( *, * ) 'YMIN set to ',ymin end if ! ! YRANGE = ! else if ( s_eqi ( command(1:6),'yrange')) then if ( s_eqi ( command(1:7),'yrange=')) then read(command(8:),*)yrange else write ( *, * ) 'Enter value for YRANGE:' read(*,*)yrange write(17,*)yrange end if if ( 0 < iprint ) then write ( *, * ) 'YRANGE set to ',yrange end if ! ! Unrecognized command ! else write ( *, * ) ' ' write ( *, * ) 'ARBY4 - Warning!' write ( *, * ) ' Unrecognized command: ' // trim ( command ) end if end do write ( *, * ) ' ' write ( *, * ) 'ARBY4 - STOP command:' write ( *, * ) ' Halt the program!' write ( *, * ) ' ' close(unit = 17) write ( *, * ) ' Closing the user input file ARBY.IN.' write ( *, * ) ' ' write ( *, * ) ' The (real) start time was '// trim ( tstart ) call date_and_time ( date, time ) tstop = time write ( *, * ) ' The (real) stopping time was '// trim ( tstop ) call delhms ( tstart, tstop, itemp ) write ( *, * ) ' The (real) elapsed time in seconds is ',itemp write ( *, * ) ' The real elapsed time in minutes is ', & real(itemp) / 60.0D+00 call etime ( tarray, estop ) write ( *, * ) ' ' write ( *, * ) ' CPU in seconds = ',estop-estart write ( *, * ) ' CPU in minutes = ',(estop-estart)/60.0D+00 write ( *, * ) ' ' write ( *, * ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine diffprb ( arb, area, epsdif, grb, indx, maxcofrb, maxelm, & maxnfl, nbcrb, ncofrb, nelem, nferb, node, np, npar, nparf, nx, & ny, par, phirb, rb, resrb, reynld, tauflo, xc, xrange, yc, yrange ) !*****************************************************************************80 ! !! DIFFPRB estimates the jacobian of the reduced function. ! ! Discussion: ! ! Finite differences are used. ! ! Modified: ! ! 31 July 1996 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! double precision ARB(MAXCOFRB,MAXCOFRB). ! ARB contains the Jacobian or Picard matrix for the reduced ! Navier Stokes system, stored as an NCOFRB by NCOFRB array. ! ! double precision AREA(3,MAXELM). ! 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. ! ! double precision EPSDIF. ! EPSDIF is a small quantity, which is used to compute the ! perturbations for the finite difference approximations. ! ! double precision GRB(NCOFRB). ! GRB contains the reduced basis coefficients of the current ! estimate of the state solution. ! ! integer INDX(3,NP). ! INDX(I,J) contains, for each node J, the global index of U, ! V and P at that node, or 0 or a negative value. The global ! index of U, V, or P is the index of the coefficient vector ! that contains the value of the finite element coefficient ! associated with the corresponding basis function at the ! given node. ! ! If K = INDX(I,J) is positive, then the value of the degree ! of freedom is stored in the solution vector entry GFL(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. ! ! integer MAXCOFRB. ! MAXCOFRB is the maximum legal value for NCOFRB, the number ! of coefficients used to specify a particular reduced basis ! solution. ! ! integer MAXELM. ! MAXELM is the maximum number of elements. ! ! integer MAXNFL. ! MAXNFL is the maximum number of equations or coefficients allowed ! for the full system. MAXNFL must be used instead of NEQNFL as ! the leading dimension of certain multi-dimensional arrays. ! ! integer NBCRB. ! NBCRB is the number of independent boundary condition ! vectors used for the reduced basis. NBCRB is normally ! at least 1, and must be no more than MAXBCRB. ! ! integer NCOFRB. ! NCOFRB is the number of coefficients needed to determine ! a particular reduced basis function. ! NCOFRB is the sum of NBCRB and NFERB. ! ! integer NELEM. ! NELEM is the number of elements. ! NELEM can be determined as 2*(NX-1)*(NY-1). ! ! integer NFERB. ! NFERB is the number of reduced basis coefficients that will ! be determined via the finite element method. ! ! integer NODE(6,MAXELM) or NODE(6,NELEM). ! NODE(I,J) contains, for an element J, the global index of ! the node whose local number in J is I. ! ! The local ordering of the nodes is suggested by this diagram: ! ! Global nodes Elements NODE ! 1 2 3 4 5 6 ! 74 84 94 3-6-1 2 Left element = (94,72,74,83,73,84) ! | / /| ! 73 83 93 5 4 4 5 Right element = (72,94,92,83,93,82) ! |/ / | ! 72 82 92 2 1-6-3 ! ! integer NP. ! NP is the number of nodes used to define the finite element mesh. ! Typically, the mesh is generated as a rectangular array, with ! an odd number of nodes in the horizontal and vertical directions. ! The formula for NP is NP = (2*NX-1)*(2*NY-1). ! ! integer NPAR. ! NPAR is the number of parameters. ! NPAR = NPARF + NPARB + 1. ! The parameters control the shape and strength of the inflow, ! the shape of the bump, and the value of the Reynolds number. ! ! integer NPARF. ! NPARF is the number of parameters associated with the ! inflow. NPARF must be at least 1. ! ! integer NX. ! NX controls the spacing of nodes and elements in ! the X direction. There are 2*NX-1 nodes along various ! lines in the X direction. ! ! The number of elements along a line in the X direction is ! NX-1 (or 2*(NX-1) to make a full rectangular strip). ! ! integer NY. ! NY controls the spacing of nodes and elements in ! the Y direction. There are 2*NY-1 nodes along various ! lines in the Y direction. ! ! The number of elements along a line in the Y direction is ! NY-1 (or 2*(NY-1) to make a full vertical strip). ! ! double precision PAR(NPAR). ! PAR contains the values of the problem parameters. ! ! PAR(1:NPARF) = inflow controls. ! PAR(NPARF+1:NPARF+NPARB) = bump controls. ! PAR(NPARF+NPARB+1) = the REYNLD parameter. ! ! double precision PHIRB(3,MAXCOFRB,15,MAXELM). ! PHIRB contains the values of a finite element basis function ! or its X or Y derivative, in a given element, at a given ! quadrature point, for a particular reduced basis function. ! ! For PHIRB(I,J,K,L), index J refers to the reduced basis ! basis functions, for J = 0 to NCOFRB. ! ! The meaning of the K index of PHIRB(I,J,K,L) is as follows: ! ! For the quadrature point I, and reduced basis function J, ! in element L, PHIRB(I,J,K,L) represents the value of: ! ! K = 1, WUrb, the finite element U velocity basis function; ! K = 2, dWUrbdX, the X derivative of WUrb; ! K = 3, dWUrbdY, the Y derivative of WUrb; ! K = 4, WVrb, the finite element V velocity basis function; ! K = 5, dWVrbdX, the X derivative of WVrb; ! K = 6, dWVrbdY, the Y derivative of WVrb; ! K = 7, Q, the finite element pressure basis function. ! K = 8, dQrbdX, the X derivative of Qrb; ! K = 9, dQrbdY, the Y derivative of Qrb. ! K = 10, WU0rb, same as WUrb, with zero BC. ! K = 11, dWU0rbdX, same as dWUrbdX, with zero BC. ! K = 12, dWU0rbdY, same as dWUrbdY, with zero BC. ! K = 13, WV0rb, same as WVrb, with zero BC. ! K = 14, dWV0rbdX, same as dWVrbdX, with zero BC. ! K = 15, dWV0rbdY, same as dWVrbdY, with zero BC. ! ! double precision RB(MAXNFL,MAXCOFRB). ! RB is the NEQNFL by NCOFRB array of reduced basis vectors. ! ! RB is based on a particular solution of the full system, ! which is saved as GFLRB. ! ! We compute the first NCOFRB derivatives of GFLRB with ! respect to a parameter. The first derivative ! is stored in column 1 of RB, and so on. ! ! double precision RESRB(NCOFRB). ! RESRB contains the residual in the reduced basis equations, ! for the parameter values PAR and reduced basis coefficients GRB. ! ! double precision REYNLD. ! REYNLD is the current value of the Reynolds number. ! Normally, REYNLD is stored as PARA(NPARF+NPARB+1). ! ! double precision TAUFLO(NPARF). ! TAUFLO contains the location of the spline abscissas for ! the inflow. ! ! A recent code change was made. For a channel flow, where ! NPARF = 1 meant a fit through 3 points, now NPARF=3 means ! a fit through 3 points. The endpoints must be explicitly ! counted. ! ! double precision XC(NP). ! XC contains the X coordinates of the nodes. ! ! double precision XRANGE. ! XRANGE is the total width of the region. ! ! double precision YC(NP). ! YC contains the Y coordinates of the nodes. ! ! double precision YRANGE. ! YRANGE is the total height of the region. ! implicit none ! integer maxcofrb integer maxelm integer maxnfl integer nelem integer ncofrb integer np integer npar integer nparf ! double precision arb(maxcofrb,ncofrb) double precision area(3,nelem) double precision delta double precision epsdif double precision grb(ncofrb) double precision grbtmp(ncofrb) integer i integer indx(3,np) integer j integer nbcrb integer nferb integer node(6,maxelm) integer nx integer ny double precision par(npar) double precision phirb(3,maxcofrb,15,maxelm) double precision rb(maxnfl,maxcofrb) double precision resrb(ncofrb) double precision reynld double precision tauflo(nparf) double precision xc(np) double precision xrange double precision yc(np) double precision yrange ! ! Get the function value at the base value. ! call fxrb(area,grb,indx,maxcofrb,maxelm,maxnfl,nbcrb,ncofrb, & nelem,nferb,node,np,npar,nparf,nx,ny,par,phirb,rb, & resrb,reynld,tauflo,xc,xrange,yc,yrange) ! ! Start each column of the jacobian equal to F(GRB). ! do j = 1, ncofrb arb(1:ncofrb,j) = resrb(1:ncofrb) end do do j = 1, ncofrb ! ! Perturb G(J). ! grbtmp(1:ncofrb) = grb(1:ncofrb) delta = epsdif * ( 1.0D+00 + abs ( grb(j) ) ) grbtmp(j) = grb(j) + delta ! ! Evaluate F(I) at the perturbed value G(J). ! call fxrb(area,grbtmp,indx,maxcofrb,maxelm,maxnfl,nbcrb,ncofrb, & nelem,nferb,node,np,npar,nparf,nx,ny,par,phirb,rb, & resrb,reynld,tauflo,xc,xrange,yc,yrange) ! ! Estimate the dependence ARB(I,J) = dF(I)/dG(J) ! arb(1:ncofrb,j) = ( resrb(1:ncofrb) - arb(1:ncofrb,j) ) / delta end do return end subroutine difsenfl ( afl, area, difcof, dpar, eqn, gfl, gflafl, ifs, & ijac, indx, ipar, ipivfl, iwrite, ldafl, maxcofrb, maxelm, maxnew, & maxnfl, ncofrb, nelem, neqnfl, nlband, node, np, npar, nparf, par, & parafl, phifl, region, resfl, senfl, splflo, tauflo, tolnew, & xrange, yc, yrange ) !*****************************************************************************80 ! !! DIFSENFL estimates full solution derivatives with respect to parameters. ! ! Discussion: ! ! The routine computes a central difference estimate for the first ! NCOFRB derivatives of the full solution GFL with respect to the IPAR-th ! parameter. ! ! DIFSENFL is rather inefficient. ALTHOUGH SOME SOLUTIONS ! ARE USED SEVERAL TIMES, DIFSENFL RECOMPUTES THEM EACH TIME. ! A CORRECTION OF THIS PROBLEM WOULD BE TO COMPUTE THE ENTIRE ! TRIANGLE OF COEFFICIENTS FIRST, AND THEN COMPUTE JUST THE ! SOLUTIONS NEEDED ONCE. ! ! Modified: ! ! 01 July 1996. ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! double precision AFL(LDAFL,MAXNFL). ! If Newton iteration is being carried out, AFL contains the ! Jacobian matrix for the full system. ! If Picard iteration is being carried out, AFL contains the ! Picard matrix for the full system. ! AFL is stored in LINPACK general band storage mode, with ! logical dimensions (3*NLBAND+1, NEQNFL). ! Where is the (I,J) entry of AFL actually stored? ! AFL has actual storage for such an entry only if ! -NLBAND <= I-J <= NLBAND. ! In such a case, the (I,J) entry is actually stored in ! AFL(I-J+2*NLBAND+1,J) ! ! double precision AREA(3,MAXELM). ! 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. ! ! double precision DIFCOF(NDIF). ! DIFCOF contains the coefficients needed to approximate ! the 0-th through (NDIF-1)-th derivatives of a function F. ! ! double precision DREY. ! DREY is the suggested increment in the REYNLD value, ! to be used during the finite difference estimations. ! ! double precision DOPT(NPAR). ! DOPT contains a set of scale factors for the parameters, used ! by the optimization code. The suggestion is that DOPT(I) be ! chosen so that DOPT(I)*PAR(I) is roughly the same order of ! magnitude for I from 1 to NPAR. ! ! double precision DPAR. ! DPAR is the suggested increment in the parameter value, ! to be used during the finite difference estimations. ! ! character ( len = 2 ) EQN(MAXNFL). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. ! ! 'U' A 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. ! 'U0' A dummy value of U = 0 should be set. ! ! 'V' A 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. ! 'V0' A dummy value of V = 0 should be set. ! ! 'P' A continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! 'P0' A dummy value of P = 0 should be set. ! ! double precision GFL(NEQNFL). ! GFL contains the current solution estimate for the full problem, ! containing the pressure and velocity coefficients. ! The vector INDX must be used to index this data. ! ! double precision GFLAFL(NEQNFL). ! GFLAFL stores the value of GFL at which the Jacobian ! was generated. ! ! integer IFS. ! IFS is the inflow shape option. ! 0, piecewise constant function. ! 1, piecewise linear function. ! 2, piecewise quadratic function. ! ! integer IJAC. ! IJAC determines the frequency for evaluating and factoring ! the Jacobian matrix during any particular Newton process. ! 1, evaluate the Jacobian on every step of the Newton ! iteration. ! n, evaluate the Jacobian only at steps 0, n, 2*n, and so on. ! ! integer INDX(3,NP). ! INDX(I,J) contains, for each node J, the global index of U, ! V and P at that node, or 0 or a negative value. The global ! index of U, V, or P is the index of the coefficient vector ! that contains the value of the finite element coefficient ! associated with the corresponding basis function at the ! given node. ! ! If K = INDX(I,J) is positive, then the value of the degree ! of freedom is stored in the solution vector entry GFL(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. ! ! integer IPAR. ! IPAR is the index of the parameter to be varied. ! ! integer IPIVFL(NEQNFL). ! IPIVFL is a pivot vector for the solution of the full ! linear system. ! ! integer IWRITE. ! IWRITE controls the amount of output printed. ! 0, print out the least amount. ! 1, print out some. ! 2, print out a lot. ! ! integer LDAFL. ! LDAFL is the first dimension of the matrix AFL as declared in ! the main program. LDAFL must be at least 3*NLBAND+1. ! ! integer MAXCOFRB. ! MAXCOFRB is the maximum legal value for NCOFRB, the number ! of coefficients used to specify a particular reduced basis ! solution. ! ! integer MAXELM. ! MAXELM is the maximum number of elements. ! ! integer MAXNEW. ! MAXNEW is the maximum number of steps to take in one Newton ! iteration. A typical value is 20. ! ! integer MAXNFL. ! MAXNFL is the maximum number of equations or coefficients allowed ! for the full system. MAXNFL must be used instead of NEQNFL as ! the leading dimension of certain multi-dimensional arrays. ! ! integer NCOFRB. ! NCOFRB is the number of coefficients needed to determine ! a particular reduced basis function. ! NCOFRB is the sum of NBCRB and NFERB. ! ! integer NELEM. ! NELEM is the number of elements. ! NELEM can be determined as 2*(NX-1)*(NY-1). ! ! integer NEQNFL. ! NEQNFL is the number of equations (and coefficients) in the full ! finite element system. ! ! integer NLBAND. ! NLBAND is the lower bandwidth of the matrix AFL. ! The zero structure of AFL is assumed to be symmetric, and so ! NLBAND is also the upper bandwidth of AFL. ! ! integer NODE(6,MAXELM) or NODE(6,NELEM). ! NODE(I,J) contains, for an element J, the global index of ! the node whose local number in J is I. ! ! The local ordering of the nodes is suggested by this diagram: ! ! Global nodes Elements NODE ! 1 2 3 4 5 6 ! 74 84 94 3-6-1 2 Left element = (94,72,74,83,73,84) ! | / /| ! 73 83 93 5 4 4 5 Right element = (72,94,92,83,93,82) ! |/ / | ! 72 82 92 2 1-6-3 ! ! integer NP. ! NP is the number of nodes used to define the finite element mesh. ! Typically, the mesh is generated as a rectangular array, with ! an odd number of nodes in the horizontal and vertical directions. ! The formula for NP is NP = (2*NX-1)*(2*NY-1). ! ! integer NPAR. ! NPAR is the number of parameters. ! NPAR = NPARF + NPARB + 1. ! The parameters control the shape and strength of the inflow, ! the shape of the bump, and the value of the Reynolds number. ! ! integer NPARF. ! NPARF is the number of parameters associated with the ! inflow. NPARF must be at least 1. ! ! double precision PAR(NPAR). ! PAR contains the values of the problem parameters. ! ! PAR(1:NPARF) = inflow controls. ! PAR(NPARF+1:NPARF+NPARB) = bump controls. ! PAR(NPARF+NPARB+1) = the REYNLD parameter. ! ! double precision PARAFL(NPAR). ! PARAFL contains the parameters where the Picard matrix or ! Jacobian of the full system was generated. ! ! double precision PHIFL(3,6,10,NELEM). ! PHIFL contains the value of a finite element basis function, its ! derivative, or other information, evaluated at the quadrature ! points (which are the element midside nodes). ! ! The meaning of the entry PHIFL(I,J,K,L) is as follows. ! For the quadrature point I, and basis function J, in element L, ! PHIFL(I,J,K,L) represents the value of: ! ! K = 1, W, the finite element basis function for velocities; ! K = 2, dWdX, the X derivative of W; ! K = 3, dWdY, the Y derivative of W; ! K = 4, Q, the finite element basis function for pressures; ! K = 5, dQdX, the X derivative of Q; ! K = 6, dQdY, the Y derivative of Q; ! K = 7, dXsidX, the X derivative of the mapping (X,Y)->XSI; ! K = 8, dXsidY, the Y derivative of the mapping (X,Y)->XSI; ! K = 9, dEtadX, the X derivative of the mapping (X,Y)->ETA; ! K = 10, dEtadY, the Y derivative of the mapping (X,Y)->ETA; ! ! In particular, PHIFL(I,J,K,L) is the value of the quadratic ! basis function W associated with local node J in element L, ! evaluated at quadrature point I. ! ! Note that PHIFL(I,J,K,L) = 0 whenever J=4, 5, or 6 and K=4, 5, or 6, ! since there are only three linear basis functions. ! ! character ( len = 20 ) REGION. ! REGION specifies the flow region. ! ! 'cavity', a driven cavity, 1 unit on each side, open on ! the top with a tangential velocity specification there. ! ! 'cavity2', a driven cavity, 1 unit on each side, open on ! the top and bottome, with tangential velocity specifications ! there. ! ! 'channel', a channel, 10 units long by 3 high, inflow on ! the left, outflow on the right, with a bump on the bottom. ! ! 'step', a channel, 12 units long by 3 high, inflow on the ! left, outflow on the right, with a step on the bottom. ! ! double precision RESFL(NEQNFL). ! RESFL contains the residual in the full basis equations. ! ! double precision SENFL(MAXNFL,MAXCOFRB). ! Columns 1 through NSENFL of SENFL contain the sensitivities ! of the full solution with respect to the REYNLD parameter, for ! orders 0 through NSENFL-1. ! ! SENFL(I,J) contains the (J-1)-th sensitivity of the I-th full unknown ! with respect to REYNLD. ! ! double precision SPLFLO(NPARF). ! SPLFLO contains the spline coefficients for the inflow. ! ! double precision TAUFLO(NPARF). ! TAUFLO contains the location of the spline abscissas for ! the inflow. ! ! A recent code change was made. For a channel flow, where ! NPARF = 1 meant a fit through 3 points, now NPARF=3 means ! a fit through 3 points. The endpoints must be explicitly ! counted. ! ! double precision TOLNEW. ! TOLNEW is the convergence tolerance for the Newton iteration. ! ! double precision XRANGE. ! XRANGE is the total width of the region. ! ! double precision YRANGE. ! YRANGE is the total height of the region. ! implicit none ! integer ldafl integer maxcofrb integer maxelm integer maxnfl integer ncofrb integer nelem integer neqnfl integer np integer npar integer nparf ! double precision afl(ldafl,neqnfl) double precision area(3,nelem) double precision difcof(ncofrb) double precision dpar character ( len = 2 ) eqn(neqnfl) double precision gfl(neqnfl) double precision gflafl(neqnfl) double precision gfltmp(neqnfl) integer i integer ierror integer ifs integer ijac integer indx(3,np) integer ipar integer ipivfl(neqnfl) integer iwrite integer j integer maxnew integer ndif integer nlband integer node(6,nelem) integer numnew double precision par(npar) double precision parafl(npar) double precision partmp(npar) double precision phifl(3,6,10,nelem) character ( len = 20 ) region double precision resfl(neqnfl) double precision rmax double precision senfl(maxnfl,maxcofrb) double precision splflo(nparf) double precision tauflo(nparf) double precision tolnew double precision xrange double precision yc(np) double precision yrange ! ! Zero out the SENFL array. ! senfl(1:neqnfl,1:ncofrb) = 0.0D+00 write ( *, * ) ' ' write ( *, * ) ' DIFSENFL: DPAR = ',dpar ! ! Compute difference NDIF, for NDIF = 0 to NCOFRB. ! do ndif = 1, ncofrb if ( ndif == 0 ) then senfl(1:neqnfl,ndif) = gfl(1:neqnfl) else if ( 2 <= iwrite ) then write ( *, '(a)' ) ' ' write ( *, * ) 'DIFSENFL - Computing difference NDIF = ',ndif end if ! ! Get the NDIF-1 order difference coefficients. ! call difset(difcof,dpar,iwrite,ndif) ! ! Evaluate the solution at several values of the parameter. ! do i = 0, ndif ! ! Copy the parameters, but reset the IPAR-th parameter value. ! partmp(1:npar) = par(1:npar) partmp(ipar) = par(ipar)+(2*i-ndif)*dpar write ( *, * ) 'J = ',j,' PAR(IPAR)=',partmp(ipar) ! ! Estimate the solution GTMP at parameters PARTMP. ! gfltmp(1:neqnfl) = gfl(1:neqnfl) ! ! Call NEWTFL to get the solution more closely. ! call newtfl(afl,area,eqn,gfltmp,gflafl,ierror,ifs,ijac,indx,ipivfl, & iwrite,ldafl,maxelm,maxnew,nelem,neqnfl,nlband,node,np,npar,nparf, & numnew,partmp,parafl,phifl,region,resfl,rmax,splflo,tauflo,tolnew, & xrange,yc,yrange) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'DIFSENFL - Fatal error!' write ( *, * ) ' NEWTFL failed, with IERROR = ',ierror stop end if ! ! Add the term associated with this solution to the estimate ! of the NDIF-th derivative. ! senfl(1:neqnfl,ndif) = senfl(1:neqnfl,ndif) & + difcof(i+1) * gfltmp(1:neqnfl) end do end if end do return end subroutine difsenrb(arb,area,difcof,dpar,grb,grbarb,indx,ipar,ipivrb,iwrite, & maxcofrb,maxelm,maxnew,maxnfl,nbcrb,ncofrb,nelem,nferb,node,np,npar,nparf, & nx,ny,par,pararb,phirb,rb,resrb,senrb,tauflo,tolnew,xc,xrange,yc,yrange) !*****************************************************************************80 ! !! DIFSENRB estimates the reduced sensitivities using finite differences. ! ! Discussion: ! ! The denominators used in the difference calculations ! get very small if the original increment is smaller than 1. ! It is strongly suggested that the parameter increment not be ! made too small! ! ! Modified: ! ! 01 August 1996 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Workspace, double precision ARB(MAXNRB,MAXNRB). ! ARB contains the Jacobian or Picard matrix for the reduced ! Navier Stokes system, stored as a dense NCOFRB by NCOFRB ! array. ! ! Input, double precision AREA(3,MAXELM). ! 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. ! ! Workspace, double precision DCOF(0:NDIF). ! DCOF contains the coefficients needed to approximate ! the NDIF-th derivative of a function F. ! ! Input, double precision DPAR. ! DPAR is the suggested increment in the parameter value, ! to be used during the finite difference estimations. ! ! Input, double precision GRB(NCOFRB). ! GRB contains the reduced basis coefficients of the current ! estimate of the state solution. ! ! Output, double precision GRBARB(NCOFRB). ! GRBARB contains the reduced basis coefficients at which ! the matrix ARB was last evaluated. ! ! Workspace, double precision GRBTMP(NCOFRB). ! ! Input, integer IPAR. ! The index of the parameter to be varied. ! ! Workspace, integer IPIVRB(NCOFRB). ! IPIVRB is a pivot vector for the solution of the reduced ! linear system. ! ! Input, integer IWRITE. ! IWRITE controls the amount of output printed. ! 0, print out the least amount. ! 1, print out some. ! 2, print out a lot. ! ! Input, integer MAXNEW. ! MAXNEW is the maximum number of steps to take in one Newton ! iteration. A typical value is 20. ! ! Input, integer MAXNRB. ! MAXNRB is the maximum number of equations allowed for the ! reduced basis system. ! ! Input, integer NELEM. ! NELEM is the number of elements. ! NELEM can be determined as 2*(NX-1)*(NY-1). ! ! Input, integer NCOFRB. ! NCOFRB is the number of basis functions used for the ! reduced basis method. (The first basis vector is labeled ! "0"). In this program, that amounts to the number of columns ! in the matrix RB. NCOFRB is also the number of reduced basis ! state equations, and reduced basis coefficients GRB. ! ! Input, integer NPAR. ! NPAR is the number of parameters. ! NPAR = NPARF + NPARB + 1. ! The parameters control the shape of the inflow, ! the shape of the bump obstacle, and the strength of the ! flow. ! ! Input, double precision PAR(NPAR). ! PAR is the current estimate for the parameters. ! PAR(1:NPARF) = inflow controls. ! PAR(NPARF+1:NPARF+NPARB) = bump controls. ! PAR(NPARF+NPARB+1) = the REYNLD parameter. ! ! Output, double precision PARARB(NPAR). ! PARARB contains the parameters where the Picard matrix or ! Jacobian of the reduced system was generated. ! ! Workspace, double precision PARTMP(NPAR). ! ! double precision PHIRB(3,NCOFRB,15,NELEM). ! PHIRB contains the values of a finite element basis function ! or its X or Y derivative, in a given element, at a given ! quadrature point, for a particular reduced basis function. ! ! For PHIRB(I,J,K,L), index J refers to the reduced basis ! basis functions, for J = 0 to NCOFRB. ! ! The meaning of the K index of PHIRB(I,J,K,L) is as follows: ! ! For the quadrature point I, and reduced basis function J, ! in element L, PHIRB(I,J,K,L) represents the value of: ! ! K = 1, WUrb, the finite element U velocity basis function; ! K = 2, dWUrbdX, the X derivative of WUrb; ! K = 3, dWUrbdY, the Y derivative of WUrb; ! K = 4, WVrb, the finite element V velocity basis function; ! K = 5, dWVrbdX, the X derivative of WVrb; ! K = 6, dWVrbdY, the Y derivative of WVrb; ! K = 7, Q, the finite element pressure basis function. ! K = 8, dQrbdX, the X derivative of Qrb; ! K = 9, dQrbdY, the Y derivative of Qrb. ! K = 10, WU0rb, same as WUrb, with zero BC. ! K = 11, dWU0rbdX, same as dWUrbdX, with zero BC. ! K = 12, dWU0rbdY, same as dWUrbdY, with zero BC. ! K = 13, WV0rb, same as WVrb, with zero BC. ! K = 14, dWV0rbdX, same as dWVrbdX, with zero BC. ! K = 15, dWV0rbdY, same as dWVrbdY, with zero BC. ! ! Workspace, double precision RESRB(NCOFRB). ! RESRB contains the residual in the reduced basis equations, ! for the parameter values PAR and reduced basis coefficients GRB. ! ! Output, double precision SENRB(MAXNRB,NCOFRB). ! SENRB contains the first several order sensitivities of the ! reduced solution with respect to the REYNLD parameter. ! SENRB(I,J) contains the J-th sensitivity of the I-th reduced unknown ! with respect to REYNLD. ! ! Input, double precision TOLNEW. ! TOLNEW is the convergence tolerance for the Newton iteration. ! implicit none ! integer maxcofrb integer maxelm integer maxnfl integer ncofrb integer nelem integer np integer npar integer nparf ! double precision arb(maxcofrb,maxcofrb) double precision area(3,nelem) double precision difcof(maxcofrb) double precision dpar double precision grb(maxcofrb) double precision grbarb(maxcofrb) double precision grbtmp(maxcofrb) integer i integer icof integer ierror integer indx(3,np) integer ipar integer ipivrb(maxcofrb) integer iwrite integer j integer jcof integer jdif integer maxnew integer nbcrb integer nferb integer node(6,nelem) integer nx integer ny double precision par(npar) double precision pararb(npar) double precision partmp(npar) double precision phirb(3,maxcofrb,15,maxelm) double precision rb(maxnfl,maxcofrb) double precision rmax double precision resrb(maxcofrb) double precision senrb(maxcofrb,maxcofrb) double precision tauflo(nparf) double precision tolnew double precision xc(np) double precision xrange double precision yc(np) double precision yrange ! ! Zero out the SENRB array. ! senrb(1:maxcofrb,1:maxcofrb) = 0.0D+00 ! ! JCOF counts the number of coefficients we will compute on ! each pass. We're done on the last pass. ! do jcof = 1, ncofrb jdif = jcof-1 if ( jdif == 0) then senrb(1:ncofrb,jcof) = grb(1:ncofrb) else write ( *, * ) ' ' write ( *, * ) 'Computing difference order JDIF = ',jdif ! ! Get the JCOF difference coefficients DIFCOF. ! call difset(difcof,dpar,iwrite,jcof) ! ! Evaluate the solution at JCOF values of the parameter. ! do icof = 1, jcof ! ! Copy the parameters, but reset the IPAR-th parameter value. ! partmp(1:npar) = par(1:npar) partmp(ipar) = par(ipar)+(2*icof-jcof-1)*dpar write ( *, * ) 'ICOF = ',ICOF,' PAR(IPAR)=',partmp(ipar) ! ! Estimate the solution GRBTMP at parameters PARTMP. ! grbtmp(1:ncofrb) = grb(1:ncofrb) ! ! Call NEWTRB to get the solution more closely. ! write ( *, * ) 'About to call NEWTRB' call newtrb(arb,area,grbtmp,grbarb,ierror,indx,ipivrb,iwrite, & maxcofrb,maxelm,maxnew,maxnfl,nbcrb,ncofrb,nelem,nferb,node,np, & npar,nparf,nx,ny,partmp,pararb,phirb,rb,resrb,rmax,tauflo,tolnew, & xc,xrange,yc,yrange) if ( ierror /= 0) then write ( *, * ) ' ' write ( *, * ) 'DIFSENRB - Fatal error!' write ( *, * ) ' NEWTRB failed, with IERROR = ',ierror stop end if ! ! Add the term associated with this solution to the estimate ! of the JDIF-th derivative. ! write ( *, * ) 'ABOUT TO ADD BLEEDING TERM' do j = 1, ncofrb senrb(j,jcof) = senrb(j,jcof)+difcof(icof)*grbtmp(j) end do end do end if end do return end subroutine flowbc(ifs,npar,nparf,par,region,splflo,tauflo,ubc,vbc, & xrange,xval,yrange,yval) !*****************************************************************************80 ! !! FLOWBC computes the specified boundary values at a given position. ! ! Modified: ! ! 07 October 1996 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IFS. ! IFS is the inflow shape option. ! 0, piecewise constant function. ! 1, piecewise linear function. ! 2, piecewise quadratic function. ! ! Input, integer NPAR. ! NPAR is the number of parameters. ! NPAR = NPARF + NPARB + 1. ! The parameters control the shape and strength of the inflow, the ! shape of the bump, and the value of the Reynolds number. ! ! Input, integer NPARF. ! NPARF is the number of parameters associated with the ! inflow. NPARF must be at least 1. ! ! Input, double precision PAR(NPAR). ! PAR contains the values of the problem parameters. ! PAR(1:NPARF) = inflow controls. ! PAR(NPARF+1:NPARF+NPARB) = bump controls. ! PAR(NPARF+NPARB+1) = the REYNLD parameter. ! ! Input, character ( len = 20 ) REGION. ! REGION specifies the flow region. ! 'cavity', a driven cavity, 1 unit on each side, open on ! the top with a tangential velocity specification there. ! 'cavity2', a driven cavity, 1 unit on each side, open on ! the top and bottom, with tangential velocity specifications ! there. ! 'channel', a channel, 10 units long by 3 high, inflow on ! the left, outflow on the right, with a bump on the bottom. ! 'step', a channel, 12 units long by 3 high, inflow on the ! left, outflow on the right, with a step on the bottom. ! ! Workspace, double precision SPLFLO(NPARF). ! SPLFLO contains the spline coefficients for the inflow. ! ! Workspace, double precision TAUFLO(NPARF). ! TAUFLO contains the location of the spline abscissas for ! the inflow. There are NPARF of them, because the end ! values of the spline are constrained to have particular ! values. ! ! Output, double precision UBC. ! UBC is the value of the horizontal velocity specified ! at (XVAL,YVAL). ! ! Output, double precision VBC. ! VBC is the value of the vertical velocity specified at (XVAL,YVAL). ! ! Input, double precision XRANGE. ! XRANGE is the total width of the region. ! ! Input, double precision XVAL. ! XVAL is the X coordinate of the point on the inflow boundary at ! which the specified velocity is desired. ! ! Input, double precision YRANGE. ! YRANGE is the total width of the region. ! ! Input, double precision YVAL. ! YVAL is the Y coordinate of the point on the inflow boundary at ! which the specified velocity is desired. ! implicit none ! integer npar integer nparf ! integer i integer ifs logical s_eqi double precision par(npar) character ( len = 20 ) region double precision splflo(nparf) double precision tauflo(nparf) double precision ubc double precision vbc double precision xrange double precision xval double precision yrange double precision yval ! ! Inflow points for the cavity have the form (X,YRANGE). ! NPARF must be at least 1. ! if ( s_eqi ( region,'cavity')) then if ( nparf == 1) then tauflo(1) = xrange/2.0D+00 else do i = 1, nparf tauflo(i) = xrange*dble((i-1))/dble(nparf-1) end do end if splflo(1:nparf) = par(1:nparf) if ( ifs == 0) then call pcval(nparf,xval,tauflo,ubc,splflo) else if ( ifs == 1) then call plval(nparf,xval,tauflo,ubc,splflo) else if ( ifs == 2) then call pqval(nparf,xval,tauflo,ubc,splflo) else write ( *, * ) ' ' write ( *, * ) 'FlowBC - Fatal error!' write ( *, * ) ' Illegal value of IFS = ',ifs stop end if vbc = 0.0D+00 ! ! Inflow points for cavity2 have the form (X,0) or (X,YRANGE). ! NPARF must be at least 2. ! else if ( s_eqi ( region,'cavity2')) then if ( nparf == 2) then tauflo(1) = xrange/2.0D+00 else do i = 1, nparf/2 tauflo(i) = xrange*dble((i-1))/dble(nparf/2-1) end do end if if ( yval == 0.0D+00 ) then do i = 1, nparf/2 splflo(i) = par(i) end do else if ( yval == yrange ) then do i = 1, nparf/2 splflo(i) = par(i+nparf/2) end do end if if ( ifs == 0) then call pcval(nparf/2,xval,tauflo,ubc,splflo) else if ( ifs == 1) then call plval(nparf/2,xval,tauflo,ubc,splflo) else if ( ifs == 2) then call pqval(nparf/2,xval,tauflo,ubc,splflo) else write ( *, * ) ' ' write ( *, * ) 'FlowBC - Fatal error!' write ( *, * ) ' Illegal value of IFS = ',ifs stop end if vbc = 0.0D+00 ! ! Inflow points for the channel have the form (0,Y). ! ! NPARF must be at least 2, which specifies values at ! (0,0) and (0,YRANGE). ! else if ( s_eqi ( region,'channel')) then do i = 1, nparf tauflo(i) = yrange*dble((i-1))/dble(nparf-1) end do splflo(1:nparf) = par(1:nparf) if ( ifs == 0) then call pcval(nparf,yval,tauflo,ubc,splflo) else if ( ifs == 1) then call plval(nparf,yval,tauflo,ubc,splflo) else if ( ifs == 2) then call pqval(nparf,yval,tauflo,ubc,splflo) else write ( *, * ) ' ' write ( *, * ) 'FlowBC - Fatal error!' write ( *, * ) ' Illegal value of IFS = ',ifs stop end if vbc = 0.0D+00 ! ! Inflow points for the step have the coordinates (0,Y). ! ! NPARF must be at least 2, which specifies values at ! (0,0) and (0,YRANGE). ! else if ( s_eqi ( region,'step')) then do i = 1, nparf tauflo(i) = yrange*dble((i-1))/dble(nparf-1) end do splflo(1:nparf) = par(1:nparf) if ( ifs == 0) then call pcval(nparf,yval,tauflo,ubc,splflo) else if ( ifs == 1) then call plval(nparf,yval,tauflo,ubc,splflo) else if ( ifs == 2) then call pqval(nparf,yval,tauflo,ubc,splflo) else write ( *, * ) ' ' write ( *, * ) 'FlowBC - Fatal error!' write ( *, * ) ' Illegal value of IFS = ',ifs stop end if vbc = 0.0D+00 end if return end subroutine fpbcrb ( arb, indx, maxcofrb, maxnfl, nbcrb, ncofrb, & nelem, node, np, nx, ny, rb, xc, xrange, yc, yrange ) !*****************************************************************************80 ! !! FPBCRB evaluates the jacobian of the reduced boundary conditions. ! ! Modified: ! ! 01 August 1996 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! double precision ARB(MAXCOFRB,MAXCOFRB). ! ARB contains the Jacobian or Picard matrix for the reduced ! Navier Stokes system, stored as an NCOFRB by NCOFRB array. ! ! integer INDX(3,NP). ! INDX(I,J) contains, for each node J, the global index of U, ! V and P at that node, or 0 or a negative value. The global ! index of U, V, or P is the index of the coefficient vector ! that contains the value of the finite element coefficient ! associated with the corresponding basis function at the ! given node. ! If K = INDX(I,J) is positive, then the value of the degree ! of freedom is stored in the solution vector entry GFL(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. ! ! Integer MAXCOFRB. ! MAXCOFRB is the maximum legal value for NCOFRB, the number ! of coefficients used to specify a particular reduced basis ! solution. ! ! integer MAXNFL. ! MAXNFL is the maximum number of equations or coefficients allowed ! for the full system. MAXNFL must be used instead of NEQNFL as ! the leading dimension of certain multi-dimensional arrays. ! ! integer NBCRB. ! NBCRB is the number of independent boundary condition ! vectors used for the reduced basis. NBCRB is normally ! at least 1, and must be no more than MAXBCRB. ! ! integer NCOFRB. ! NCOFRB is the number of coefficients needed to determine ! a particular reduced basis function. ! NCOFRB is the sum of NBCRB and NFERB. ! ! integer NELEM. ! NELEM is the number of elements. ! NELEM can be determined as 2*(NX-1)*(NY-1). ! ! integer NODE(6,MAXELM) or NODE(6,NELEM). ! NODE(I,J) contains, for an element J, the global index of ! the node whose local number in J is I. ! The local ordering of the nodes is suggested by this diagram: ! ! Global nodes Elements NODE ! 1 2 3 4 5 6 ! 74 84 94 3-6-1 2 Left element = (94,72,74,83,73,84) ! | / /| ! 73 83 93 5 4 4 5 Right element = (72,94,92,83,93,82) ! |/ / | ! 72 82 92 2 1-6-3 ! ! integer NP. ! NP is the number of nodes used to define the finite element mesh. ! Typically, the mesh is generated as a rectangular array, with ! an odd number of nodes in the horizontal and vertical directions. ! The formula for NP is NP = (2*NX-1)*(2*NY-1). ! ! integer NX. ! NX controls the spacing of nodes and elements in ! the X direction. There are 2*NX-1 nodes along various ! lines in the X direction. ! The number of elements along a line in the X direction is ! NX-1 (or 2*(NX-1) to make a full rectangular strip). ! ! integer NY. ! NY controls the spacing of nodes and elements in ! the Y direction. There are 2*NY-1 nodes along various ! lines in the Y direction. ! The number of elements along a line in the Y direction is ! NY-1 (or 2*(NY-1) to make a full vertical strip). ! ! double precision RB(MAXNFL,MAXCOFRB). ! RB is the NEQNFL by NCOFRB array of reduced basis vectors. ! RB is generated by computing a finite element solution GFL, ! which is saved for later reference as "GFLRB". ! GFLRB is copied into the first column of RB. ! Then, we compute the first NCOFRB derivatives of GFLRB with ! respect to a parameter. The first derivative ! is stored in column 1 of RB, and so on. ! Now we compute the QR factorization of this matrix. ! We intend that NEQNFL >> NCOFRB, and RB is a matrix with orthogonal ! columns, so that: ! Transpose(RB) * RB = Identity(1+NCOFRB) ! If GFL is any set of finite element coefficients, the corresponding ! set of reduced basis coefficients can be computed as: ! GRB = Transpose(RB) * GFL ! If GRB is a set of reduced basis coefficients, a corresponding ! set of finite element coefficients can be computed as: ! GFL = RB * GRB. ! While it is the case that you can expand and then reduce, ! and always get the same result, it is not the case that ! when you reduce and then expand you get the same result! ! It is true, for ANY GRB, that ! GRB = Transpose(RB) * RB * GRB ! which follows from Transpose(RB) * RB = Identity(1+NCOFRB). ! However, for a general GFL, it is the case that ! GFL = /= RB * Transpose(RB) * GFL. ! Only if GFL was generated from a reduced basis coefficient ! vector will equality apply. In other words, if GFL was generated ! from a reduced basis coefficient: ! GFL = RB * GRB ! then ! RB * Transpose(RB) * GFL = RB * Transpose(RB) * (RB * GRB) ! = RB * GRB = GFL ! so in this strictly limited case, ! RB * Transpose(RB) = Identity(NEQNFL). ! ! double precision XRANGE. ! XRANGE is the total width of the region. ! ! double precision YRANGE. ! YRANGE is the total height of the region. ! implicit none ! integer maxcofrb integer maxnfl integer ncofrb integer nelem integer np ! double precision arb(maxcofrb,maxcofrb) double precision dwdx double precision dwdy integer ibcrb integer icoffl integer icofrb integer icol integer ielem integer indx(3,np) integer inode integer iq integer nbcrb integer node(6,nelem) integer nx integer ny double precision rb(maxnfl,maxcofrb) double precision w double precision wurb double precision xbc double precision xc(np) double precision xrange double precision ybc double precision yc(np) double precision yrange ! ! Zero out the BC rows of the matrix. ! arb(1:nbcrb,1:ncofrb) = 0.0D+00 do ibcrb = 1, nbcrb ! ! For the driven cavity, the boundary collocation points are evenly ! spaced between the ends of the upper boundary. ! xbc = xrange * dble(ibcrb)/dble(nbcrb+1) ybc = yrange icol = 1 + int ( xbc * dble ( nx - 1 ) / xrange ) if ( nx - 1 < icol ) then icol = nx-1 end if ielem = icol*(2*ny-2)-1 ! ! Evaluate the reduced solution UBCRB at (XBC,YBC). ! do icofrb = 1, ncofrb wurb = 0.0D+00 do iq = 1, 6 call qbf(ielem,iq,w,dwdx,dwdy,nelem,node,np,xc,xbc,yc,ybc) inode = node(iq,ielem) icoffl = indx(1,inode) wurb = wurb+rb(icoffl,icofrb)*w end do arb(ibcrb,icofrb) = wurb end do end do return end subroutine fpferb ( arb, area, grb, maxcofrb, maxelm, nbcrb, ncofrb, & nelem, nferb, phirb, reynld ) !*****************************************************************************80 ! !! FPFERB evaluates the reduced basis jacobian directly. ! ! Discussion: ! ! FPFERB computes the reduced basis jacobian without computing the ! full basis jacobian first. ! ! FPFERB is given ! ! GRB, the reduced basis coefficients of an approximate solution, ! PHIRB, the reduced basis functions, evaluated at the quadrature ! points, ! REYNLD, the current Reynolds number, ! ! and computes ! ! ARB, the reduced basis jacobian of the Navier Stokes equations. ! ! Modified: ! ! 01 August 1996 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! double precision ARB(MAXCOFRB,MAXCOFRB). ! ARB contains the Jacobian or Picard matrix for the reduced ! Navier Stokes system, stored as an NCOFRB by NCOFRB array. ! ! double precision AREA(3,MAXELM). ! 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. ! ! double precision GRB(NCOFRB). ! GRB contains the reduced basis coefficients of the current ! estimate of the state solution. ! ! integer MAXCOFRB. ! MAXCOFRB is the maximum legal value for NCOFRB, the number ! of coefficients used to specify a particular reduced basis ! solution. ! ! integer MAXELM. ! MAXELM is the maximum number of elements. ! ! integer NBCRB. ! NBCRB is the number of independent boundary condition ! vectors used for the reduced basis. NBCRB is normally ! at least 1, and must be no more than MAXBCRB. ! ! integer NCOFRB. ! NCOFRB is the number of coefficients needed to determine ! a particular reduced basis function. ! NCOFRB is the sum of NBCRB and NFERB. ! ! integer NELEM. ! NELEM is the number of elements. ! NELEM can be determined as 2*(NX-1)*(NY-1). ! ! integer NFERB. ! NFERB is the number of reduced basis coefficients that will ! be determined via the finite element method. ! ! double precision PHIRB(3,MAXCOFRB,15,MAXELM). ! PHIRB contains the values of a finite element basis function ! or its X or Y derivative, in a given element, at a given ! quadrature point, for a particular reduced basis function. ! ! For PHIRB(I,J,K,L), index J refers to the reduced basis ! basis functions, for J = 0 to NCOFRB. ! ! The meaning of the K index of PHIRB(I,J,K,L) is as follows: ! ! For the quadrature point I, and reduced basis function J, ! in element L, PHIRB(I,J,K,L) represents the value of: ! ! K = 1, WUrb, the finite element U velocity basis function; ! K = 2, dWUrbdX, the X derivative of WUrb; ! K = 3, dWUrbdY, the Y derivative of WUrb; ! K = 4, WVrb, the finite element V velocity basis function; ! K = 5, dWVrbdX, the X derivative of WVrb; ! K = 6, dWVrbdY, the Y derivative of WVrb; ! K = 7, Q, the finite element pressure basis function. ! K = 8, dQrbdX, the X derivative of Qrb; ! K = 9, dQrbdY, the Y derivative of Qrb. ! K = 10, WU0rb, same as WUrb, with zero BC. ! K = 11, dWU0rbdX, same as dWUrbdX, with zero BC. ! K = 12, dWU0rbdY, same as dWUrbdY, with zero BC. ! K = 13, WV0rb, same as WVrb, with zero BC. ! K = 14, dWV0rbdX, same as dWVrbdX, with zero BC. ! K = 15, dWV0rbdY, same as dWVrbdY, with zero BC. ! ! double precision REYNLD. ! REYNLD is the current value of the Reynolds number. ! Normally, REYNLD is stored as PARA(NPARF+NPARB+1). ! implicit none ! integer maxcofrb integer maxelm integer ncofrb integer nelem ! double precision ar double precision arb(maxcofrb,maxcofrb) double precision area(3,maxelm) double precision dqjdx double precision dqjdy double precision dprbdx double precision dprbdy double precision durbdx double precision durbdy double precision dvrbdx double precision dvrbdy double precision dwu0dx double precision dwujdx double precision dwu0dy double precision dwujdy double precision dwv0dx double precision dwvjdx double precision dwv0dy double precision dwvjdy double precision grb(ncofrb) integer icofrb integer ielem integer iquad integer jcofrb logical s_eqi integer nbcrb integer nferb double precision prb double precision phirb(3,maxcofrb,15,maxelm) double precision reynld double precision urb double precision vrb double precision wu0 double precision wuj double precision wv0 double precision wvj ! ! Zero out the FE rows of the matrix. ! do icofrb = nbcrb+1, nbcrb+nferb arb(icofrb,1:ncofrb) = 0.0D+00 end do ! ! Consider an element IELEM... ! do ielem = 1, nelem ! ! ...and a quadrature point IQUAD... ! do iquad = 1, 3 ar = area(iquad,ielem) ! ! For the given reduced coefficients GRB, and basis functions ! PHIRB, evaluate U, V, and P, and their spatial derivatives. ! call uvpqrb(dprbdx,dprbdy,durbdx,durbdy,dvrbdx,dvrbdy,grb, & ielem,iquad,maxcofrb,maxelm,ncofrb,phirb,prb,urb,vrb) ! ! Consider FE reduced basis function ICOFRB. ! do icofrb = nbcrb+1, nbcrb+nferb wu0 = phirb(iquad,icofrb,10,ielem) dwu0dx = phirb(iquad,icofrb,11,ielem) dwu0dy = phirb(iquad,icofrb,12,ielem) wv0 = phirb(iquad,icofrb,13,ielem) dwv0dx = phirb(iquad,icofrb,14,ielem) dwv0dy = phirb(iquad,icofrb,15,ielem) ! ! Take the derivative with respect to basis function JCOFRB. ! do jcofrb = 1, ncofrb wuj = phirb(iquad,jcofrb,1,ielem) dwujdx = phirb(iquad,jcofrb,2,ielem) dwujdy = phirb(iquad,jcofrb,3,ielem) wvj = phirb(iquad,jcofrb,4,ielem) dwvjdx = phirb(iquad,jcofrb,5,ielem) dwvjdy = phirb(iquad,jcofrb,6,ielem) dqjdx = phirb(iquad,jcofrb,8,ielem) dqjdy = phirb(iquad,jcofrb,9,ielem) ! ! The horizontal momentum equations. ! arb(icofrb,jcofrb) = arb(icofrb,jcofrb)+ar* & (dwujdx*dwu0dx + dwujdy*dwu0dy+reynld & *(wuj*durbdx+urb*dwujdx+wvj*durbdy+vrb*dwujdy+dqjdx)*wu0) ! ! The vertical momentum equations. ! arb(icofrb,jcofrb) = arb(icofrb,jcofrb)+ar* & (dwvjdx*dwv0dx + dwvjdy*dwv0dy +reynld & *(wuj*dvrbdx+urb*dwvjdx+wvj*dvrbdy+vrb*dwvjdy+dqjdy)*wv0) end do end do end do end do return end subroutine fpfl ( afl, area, eqn, gfl, indx, ldafl, maxelm, nelem, neqnfl, & nlband, node, np, npar, par, phifl ) !*****************************************************************************80 ! !! FPFL computes the jacobian of the residual function of the full solution. ! ! Discussion: ! ! The differentiated Navier Stokes functions have the form: ! ! ! d U-Eqn/d U-Coef: ! ! Integral ! ! dWj/dx * dWi/dx + dWj/dy * dWi/dy ! + reynld * (Wj*dUold/dx + Uold*dWj/dx+ Vold*dWj/dy) * Wi dx dy ! ! d U-Eqn/d V-Coef: ! ! Integral ! ! reynld * Wj*dUold/dy * Wi dx dy ! ! d U-Eqn/d P-Coef: ! ! Integral ! ! reynld * dQj/dx * Wi dx dy ! ! d V-Eqn/d U-Coef: ! ! Integral ! ! reynld * Wj*dVold/dx * Wi dx dy ! ! d V-Eqn/d V-Coef: ! ! Integral ! ! dWj/dx * dWi/dx + dWj/dy * dWi/dy ! + reynld * (Uold*dWj/dx + Wj*dVold/dy + Vold*dWj/dy) * Wi dx dy ! ! d V-Eqn/d P-Coef: ! ! Integral ! ! reynld * 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: ! ! 21 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! double precision AFL(LDAFL,MAXNFL). ! If Newton iteration is being carried out, AFL contains the ! Jacobian matrix for the full system. ! If Picard iteration is being carried out, AFL contains the ! Picard matrix for the full system. ! ! AFL is stored in LINPACK general band storage mode, with ! logical dimensions (3*NLBAND+1, NEQNFL). ! ! Where is the (I,J) entry of AFL actually stored? ! AFL has actual storage for such an entry only if ! -NLBAND <= I-J <= NLBAND. ! In such a case, the (I,J) entry is actually stored in ! AFL(I-J+2*NLBAND+1,J) ! ! double precision AREA(3,MAXELM). ! 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. ! ! character ( len = 2 ) EQN(MAXNFL). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. ! ! 'U' A 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. ! 'U0' A dummy value of U = 0 should be set. ! ! 'V' A 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. ! 'V0' A dummy value of V = 0 should be set. ! ! 'P' A continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! 'P0' A dummy value of P = 0 should be set. ! ! double precision GFL(NEQNFL). ! GFL contains the current solution estimate for the full problem, ! containing the pressure and velocity coefficients. ! The vector INDX must be used to index this data. ! ! integer INDX(3,NP). ! INDX(I,J) contains, for each node J, the global index of U, ! V and P at that node, or 0 or a negative value. The global ! index of U, V, or P is the index of the coefficient vector ! that contains the value of the finite element coefficient ! associated with the corresponding basis function at the ! given node. ! ! integer LDAFL. ! LDAFL is the first dimension of the matrix AFL as declared in ! the main program. LDAFL must be at least 3*NLBAND+1. ! ! integer MAXELM. ! MAXELM is the maximum number of elements allowed. ! ! integer NELEM. ! NELEM is the number of elements. ! NELEM can be determined as 2*(NX-1)*(NY-1). ! ! integer NEQNFL. ! NEQNFL is the number of equations (and coefficients) in the full ! finite element system. ! ! integer NLBAND. ! NLBAND is the lower bandwidth of the matrix AFL. ! The zero structure of AFL is assumed to be symmetric, and so ! NLBAND is also the upper bandwidth of AFL. ! ! integer NODE(6,MAXELM) or NODE(6,NELEM). ! NODE(I,J) contains, for an element J, the global index of ! the node whose local number in J is I. ! ! integer NP. ! NP is the number of nodes used to define the finite element mesh. ! Typically, the mesh is generated as a rectangular array, with ! an odd number of nodes in the horizontal and vertical directions. !