c c -------------------------------------------------------------- c program main implicit double precision (a-h,o-z) parameter(n1d=2100,memsize=333333) common /coord/ xgrid(0:n1d) common /domain/ corn1,corn2 common /fluid/ gamma(5),zf(5),qchem(5),gmole(5),mphase common /frontin/ xloop(10),varloop(2,6,10), & ibdloop(10),mwloop(10),istloop(10), & newloop(10),irrt0(10),irrtn,nloops common /hypmeth/ mthslp(3),mthlim(3) common /init0/ rhol,ul,pl,rhor,ur,pr common /phypar/ den0,sqrm common /space/ lfree(50,2),lhead(6),lenf,idimf,lentot,lenmax common /sterms/ nsource,isource(10),igeom common /trackin/ qjtol(3),itrack,idetec(3),isurgy dimension alloc(memsize),work(6) dimension ibc(2),ql(n1d,6),qr(n1d,6) character*8 pltfile,infile,outfile logical vtime c c # main program for solving Euler's equations of gas dynamics c # numerical methods: front-tracking methods and relates c infile = 'tape5' outfile = 'tape6' pltfile = 'plot.3' open(3,file=pltfile,status='unknown',form='formatted') open(5,file=infile,status='old',form='formatted') open(6,file=outfile,status='unknown',form='formatted') open(97,file='front.d',status='unknown',form='formatted') c c # input data read(5,*) iprob ! type of problem read(5,*) meqn,mwave ! meqn, wave families read(5,*) corn1,corn2 ! computational domain read(5,*) ibc(1),ibc(2) ! boundary conditions read(5,*) nsource,(isource(i),i=1,nsource) ! source terms read(5,*) igeom ! type of geometry read(5,*) hx,dt0 ! mesh spacing, fixed timestep read(5,*) vtime,cfl ! variable timestep, Courant # read(5,*) nstop,tstop ! stopping step and time read(5,*) io2d ! 2D output read(5,*) iorder ! order of hyperbolic integrator read(5,*) (mthslp(mw),mw=1,mwave) ! method of slopes read(5,*) (mthlim(mw),mw=1,mwave) ! limiter c read(5,*) itrack ! tracking option read(5,*) nloops ! number of tracked points do i=1,nloops read(5,*) istloop(i),newloop(i),ibdloop(i),xloop(i) enddo c read(5,*) (idetec(mw),mw=1,mwave) ! shock detection criterion read(5,*) (qjtol(mw),mw=1,mwave) ! shock detection tolerance read(5,*) isurgy ! wave interaction c read(5,*) iriem1 ! regular cell Riemann solver read(5,*) iriem2 ! irregular cell Riemann solver read(5,*) mphase ! multi-phase model read(5,*) nlayer,(gamma(i),i=1,nlayer) read(5,*) (zf(i),i=1,nlayer) read(5,*) (gmole(i),i=1,nlayer) read(5,*) den0,sqrm ! density ratio, Mach # c write(6,1001) corn1,corn2,hx,dt0,vtime,cfl,nstop,tstop, & iorder,itrack,iprob,ibc(1),ibc(2), & iriem1,iriem2 c time = 0.d0 dt = dt0 nstart = 0 c c # numerical boundary zones lwidth = iorder c c # total number of cells mitot = idnint((corn2-corn1)/hx)+2*lwidth c c # regular mesh points xlow = corn1-dfloat(lwidth)*hx do i=0,mitot+2 xgrid(i) = xlow+dfloat(i-1)*hx enddo c c # computational domain ix1 = lwidth+1 ixn = mitot-lwidth imin = ix1-lwidth+1 imax = ixn+lwidth c c # set data structure call setdoc(alloc,memsize) c c # allocate memory loc = maloc(mitot*meqn) locold = maloc(mitot*meqn) locirr = maloc(mitot) lociro = maloc(mitot) loca = maloc(mitot*mwave) locb = maloc(mitot*meqn*mwave) locc = maloc(mitot*mwave) locw = maloc(mitot) c c # set initial condition call initc(alloc(loc),mitot,meqn,hx,xlow,iprob) c c # set grid structure call setgrd(alloc(loc),alloc(locirr),mitot,meqn,hx, & lstgrd,lstfrt) c if (io2d .eq. 1) then write(91,1003) time call outvar(alloc(loc),alloc(locirr),mitot,meqn, & hx,ix1,ixn,work,3,lstgrd,1) endif c c # main integration loop do 900 nstep=nstart+1,nstop c # set boundary states call physbd(alloc(loc),alloc(locirr),mitot,meqn,lwidth,ibc, & work,lstgrd,iprob) c c # compute variable time step if (vtime) call vartme(alloc(loc),alloc(locirr),mitot, & meqn,hx,time,nstep,dt,cfl, & work,lstgrd) c if(time+dt .gt. tstop) dt = tstop-time c c # handle source term if (iorder .eq. 2) call source(alloc(loc),alloc(locirr), & mitot,meqn,hx,dt,iorder,work,lstgrd) c c # prepare Riemann data at regular cell interfaces call rpdata(ql,qr,alloc(loc),alloc(locirr),alloc(locw), & mitot,meqn,imin,imax,lstgrd) c c # solve Riemann problems at the regular cell interfaces call vrmjmp(ql,qr,alloc(loca),alloc(locb),alloc(locc), & mitot,meqn,mwave,imin,imax,iriem1,napprx1) c c # solve Riemann problems at the tracked interfaces call rfront(alloc(locirr),mitot,meqn,mwave,work,lstfrt, & iriem2,napprx2) c c # flag tracked points call ffront(alloc(loca),alloc(locb),alloc(locc),mitot,meqn, & mwave,imin,imax,time,work,lstfrt) c c # check time step to avoid wave interactions occurred c # in the current time step call fsurgy(dt,alloc(locw),mitot) c c # initialize grid structure at the next time step call mkgrid(alloc(loc),alloc(locold),alloc(locirr), & alloc(lociro),mitot,meqn,work,lstgrd,lstold) c c # advance tracked fronts to the next time step c # and insert them to the grid call ifront(alloc(locirr),alloc(lociro),mitot,hx,time,dt, & lstgrd,lstfrt) c c # set values for the irregular cells call setqir(alloc(locold),alloc(lociro),mitot,meqn,hx, & work,lstgrd,lstold,irrtn) c c # update cell averages: wave propagation method c # for regular cells call godunv(alloc(loc),alloc(locirr),alloc(lociro), & alloc(loca),alloc(locb),mitot,meqn,mwave, & hx,dt,ix1,ixn,iorder,work,lstgrd,lstold) c c # for irregular cells call pfront(alloc(loc),alloc(locirr),alloc(lociro), & alloc(loca),alloc(locb),alloc(locw), & mitot,meqn,mwave,hx,dt,iorder,work, & lstgrd,lstold) c c # delete old fronts call dfront(lstold,lstfrt) c c # handle source terms call source(alloc(loc),alloc(locirr),mitot,meqn,hx,dt, & iorder,work,lstgrd) c c # update time time = time+dt c if (io2d .eq. 1) then write(91,1003) time call outvar(alloc(loc),alloc(locirr),mitot,meqn, & hx,ix1,ixn,work,3,lstgrd,1) endif c write(6,1002) nstep,time,dt c c # check for stopping the program if ((time .ge. tstop) .or. (nstep .eq. nstop)) & go to 999 900 continue c 999 continue c # output state variables iunit = 21 call outvar(alloc(loc),alloc(locirr),mitot,meqn, & hx,ix1,ixn,work,iunit,lstgrd,0) c write(6,1004) c 1001 format(' input parameters: ',/, & ' corn1 (left bdry) ',e12.5,/, & ' corn2 (right bdry) ',e12.5,/, & ' hx ',e12.5,/, & ' fixed dt ',e12.5,/, & ' variable time step ',l9,/, & ' cfl # (if var. delt) ',e12.5,/, & ' nstop ',i9,/, & ' tstop ',e12.5,/, & ' iorder ',i9,/, & ' itrack ',i9,/, & ' iprob ',i9,/, & ' ibc(1) ',i9,/, & ' ibc(2) ',i9,/, & ' iriem1 ',i9,/, & ' iriem2 ',i9,/) 1002 format('nstep=',i4, ', time=',f10.6,', dt=',f10.6) 1003 format(f16.8) 1004 format(/,' ------ end of integration -------- ') c stop end c c ------------------------------------------------------------------- c subroutine initc(q,mitot,meqn,hx,xlow,iprob) implicit double precision (a-h,o-z) double precision mach common /domain/ corn1,corn2 common /fluid/ gamma(5),zf(5),qchem(5),gmole(5),mphase common /frontin/ xloop(10),varloop(2,6,10), & ibdloop(10),mwloop(10),istloop(10), & newloop(10),irrt0(10),irrtn,nloops common /fshockc/ xp common /gravity/ gravc(2) common /phypar/ den0,sqrm dimension q(mitot,meqn) dimension ql(6),qm(6),qr(6) c c # set initial condition c sloc1 = .5d0 sloc2 = sloc1 irrt0(1) = 1 irrtn = 2 c go to (1,2,3,4,5,6,7) iprob c 1 continue c # Sod's Riemann data rhol = 1.d0 ul = 0.d0 pl = 1.d0 c rhor = .125d0 ur = 0.d0 pr = 0.1d0 go to 900 c 2 continue c # Sod's Riemann data c # level set: distance function rhol = 1.d0 ul = 0.d0 pl = 1.d0 c rhor = .125d0 ur = 0.d0 pr = 0.1d0 c do i=1,mitot x = xlow+(dfloat(i-1)+0.5d0)*hx z0 = x-sloc1 if (z0 .lt. 0.d0) then call prmtoc(ql,rhol,ul,pl,gamma(1)) do m=1,meqn q(i,m) = ql(m) enddo else call prmtoc(qr,rhor,ur,pr,gamma(2)) do m=1,meqn q(i,m) = qr(m) enddo endif enddo return c 3 continue c # gravity-driven flow c # (1D) Rayleigh-Taylor instability c # heavy fluid rhol = 1.d0 cl = 1.d0 ul = 0.0d0 p00 = cl**2*rhol/gamma(1) pl = p00 templ = pl*gmole(1)/rhol c c # light fluid rhor = rhol/den0 ur = 0.0d0 pr = p00 cr = dsqrt(gamma(2)*pr/rhor) tempr = pr*gmole(2)/rhor c atwood = (rhol-rhor)/(rhol+rhor) c write(6,1992) rhol,ul,pl,templ write(6,1994) rhor,ur,pr,tempr go to 900 c 4 continue c # shock-contact interaction c # (1D) Richtmyer-Meshkov instability c # middle: cm = 1.d0 um = 0.d0 pm = 1.d0 rhom = gamma(2)*pm/cm**2 tempm = pm*gmole(2)/rhom c c # bottom: light or heavy gas rhol = rhom/den0 ul = 0.d0 pl = 1.d0 templ = pl*gmole(1)/rhol c write(6,1992) rhol,ul,pl,templ write(6,1994) rhom,um,pm,tempm c c # Atwood number atwood = (rhol-rhom)/(rhol+rhom) c call prmtoc(ql,rhol,ul,pl,gamma(1),zf(1)) call prmtoc(qm,rhom,um,pm,gamma(2),zf(2)) c irrt0(1) = 1 xp = xloop(1) do i=1,mitot call celave(i,hx,wl) do m=1,meqn q(i,m) = wl*ql(m)+(1.d0-wl)*qm(m) enddo enddo c do m=1,meqn varloop(1,m,1) = ql(m) varloop(2,m,1) = qm(m) enddo c c # given pre-shock state and Mach number c # compute post-shock state mach = (um-sqrm)/cm call shock(rhor,ur,pr,rhom,um,pm,mach,gamma(2)) ur = dsign(1.d0,mach)*ur call prmtoc(qr,rhor,ur,pr,gamma(2),zf(2)) c irrt0(2) = 1 xp = xloop(2) do i=1,mitot call celave(i,hx,wl) if (wl .ne. 1.d0) then do m=1,meqn q(i,m) = wl*qm(m)+(1.d0-wl)*qr(m) enddo endif enddo c do m=1,meqn varloop(1,m,2) = qm(m) varloop(2,m,2) = qr(m) enddo return c 5 continue c # radially symmetric problem: expanding flow irrt0(1) = 3 xp = xloop(1) rs = xp c ur = 0.d0 pr = 1.d0 cr = 1.d0 rhor = gamma(1)*pr/cr**2 tempr = pr*gmole(1)/rhor c c # given pre-shock state and Mach number c # compute post-shock state cr = dsqrt(gamma(1)*pr/rhor) mach = (ur-sqrm)/cr call shock(rhol,ul,pl,rhor,ur,pr,mach,gamma(1)) call prmtoc(ql,rhol,ul,pl,gamma(1),zf(1)) call prmtoc(qr,rhor,ur,pr,gamma(2),zf(2)) c do i=1,mitot r = xlow+(dfloat(i-1)+0.5d0)*hx u = (r/rs)**2*ul call prmtoc(qm,rhol,u,pl,gamma(1),zf(1)) c call celave(i,hx,wl) do m=1,meqn q(i,m) = wl*qm(m)+(1.d0-wl)*qr(m) enddo enddo c do m=1,meqn varloop(1,m,1) = ql(m) varloop(2,m,1) = qr(m) enddo return c 6 continue c # radially symmetric problem: converging flow irrt0(1) = 3 xp = xloop(1) rs = xp c ur = 0.d0 pr = 1.d0 cr = 1.d0 rhor = gamma(1)*pr/cr**2 tempr = pr*gmole(1)/rhor c c # given pre-shock state and Mach number c # compute post-shock state cr = dsqrt(gamma(1)*pr/rhor) mach = (ur-sqrm)/cr call shock(rhol,ul,pl,rhor,ur,pr,mach,gamma(1)) call prmtoc(ql,rhol,ul,pl,gamma(1),zf(1)) call prmtoc(qr,rhor,ur,pr,gamma(2),zf(2)) c do i=1,mitot r = xlow+(dfloat(i-1)+0.5d0)*hx u = (r/rs)**2*ul call prmtoc(qm,rhol,u,pl,gamma(1),zf(1)) c call celave(i,hx,wl) do m=1,meqn q(i,m) = wl*qm(m)+(1.d0-wl)*qr(m) enddo enddo c do m=1,meqn varloop(1,m,1) = ql(m) varloop(2,m,1) = qr(m) enddo return c 7 continue c # steady state nozzle flow xp = xloop(1) rhol = .502d0 ul = 1.299d0 pl = .3809d0 zl = 0.d0 rhor = .776d0 c htot = corn2-corn1 b = (corn2*rhol-corn1*rhor)/htot slope = (rhor-rhol)/htot c do i=1,mitot x = xlow+(dfloat(i-1)+0.5d0)*hx xnext = x+hx rho = 0.5d0*(x+xnext)*slope+b call prmtoc(qm,rho,ul,pl,gamma(1),zf(1)) c do m=1,meqn q(i,m) = qm(m) enddo enddo return c 900 continue call prmtoc(ql,rhol,ul,pl,gamma(1),zf(1)) call prmtoc(qr,rhor,ur,pr,gamma(2),zf(2)) c do i=1,mitot x = xlow + (dfloat(i-1)+0.5d0)*hx if (x .lt. sloc1) then do m=1,meqn q(i,m) = ql(m) enddo else do m=1,meqn q(i,m) = qr(m) enddo endif enddo c 1992 format('contact (lf): rho=',f9.5,', u=',f9.5, & ', p=',f9.5, ', T=',f9.5) 1994 format('contact (rg): rho=',f9.5,', u=',f9.5, & ', p=',f9.5, ', T=',f9.5) 1996 format('sqmach=',f10.6,', gravc(1)=',f10.6,', D =',f10.6, & ', atwood=',f10.6) c return end c c ---------------------------------------------------------------------- c subroutine physbd(q,irr,mitot,meqn,lwidth,ibc,work,lstgrd,iprob) implicit double precision (a-h,o-z) parameter(maxbd=126) common /init0/ rhol,ul,pl,rhor,ur,pr common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),work(meqn),ibc(2) c c # set up boundary conditions c ix1 = lwidth+1 ixn = mitot-lwidth c c # left hand side go to (1,2,3,4) ibc(1) c 1 continue c # outflow bc do i=1,lwidth do m=1,meqn q(i,m) = q(ix1,m) enddo enddo go to 101 c 2 continue c # solid wall ico = 0 do i=lwidth,1,-1 do m=1,meqn q(i,m) = q(ix1+ico,m) enddo q(i,2) = -q(i,2) ico = ico+1 enddo go to 101 c 3 continue c # periodic bc ico = 0 do i=lwidth,1,-1 do m=1,meqn q(i,m) = q(ixn-ico,m) enddo ico = ico+1 enddo go to 101 c 4 continue c # characteristic bc rho0 = rhol if (iprob .eq. 15) rho0 = .56d0 c if (irr(ix1) .ne. lstgrd) then iptr = irr(ix1) do m=1,meqn work(m) = qirr(1,m,iptr) enddo else do m=1,meqn work(m) = q(ix1,m) enddo endif call ctoprm(work,rho,u,p,gamma0,zf0) gasc1 = gamma0-1.d0 c c = dsqrt(gamma0*p/rho) s0 = p/rho**gamma0 p0 = rho0**gamma0*s0 c0 = dsqrt(gamma0*p0/rho0) rim = u-2.d0/gasc1*c rip = rim+4.d0/gasc1*c0 u0 = .5d0*(rim+rip) c call prmtoc(work,rho0,u0,p0,gamma0,zf0) c do i=1,lwidth do m=1,meqn q(i,m) = work(m) enddo enddo go to 101 c 101 continue c # right hand side go to (110,120,130,140) ibc(2) c 110 continue c # outflow bc do i=ixn+1,mitot do m=1,meqn q(i,m) = q(ixn,m) enddo enddo return c 120 continue c # solid wall ico = 0 do i=ixn+1,mitot do m=1,meqn q(i,m) = q(ixn-ico,m) enddo q(i,2) = -q(i,2) ico = ico+1 enddo return c 130 continue c # periodic bc ico = 0 do i=ixn+1,mitot do m=1,meqn q(i,m) = q(ix1+ico,m) enddo ico = ico+1 enddo return c 140 continue c # characteristic bc rho0 = rhol if (iprob .eq. 15) rho0 = .56d0 c if (irr(ix1) .ne. lstgrd) then iptr = irr(ix1) do m=1,meqn work(m) = qirr(1,m,iptr) enddo else do m=1,meqn work(m) = q(ix1,m) enddo endif call ctoprm(work,rho,u,p,gamma0,zf0) gasc1 = gamma0-1.d0 c c = dsqrt(gamma0*p/rho) s0 = p/rho**gamma0 p0 = rho0**gamma0*s0 c0 = dsqrt(gamma0*p0/rho0) rim = u-2.d0/gasc1*c rip = rim+4.d0/gasc1*c0 u0 = .5d0*(rim+rip) c call prmtoc(work,rho0,u0,p0,gamma0,zf0) c do i=1,lwidth do m=1,meqn q(i,m) = work(m) enddo enddo return c end c c ------------------------------------------------------------------------ c subroutine vartme(q,irr,mitot,meqn,hx,time,nstep,dt,cfl,work, & lstgrd) implicit double precision (a-h,o-z) parameter(maxbd=126) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),work(meqn) c c # compute variable time step c cmax = 0.d0 do 10 i=1,mitot if (iabs(irr(i)) .ne. lstgrd) then iptr = irr(i) do icell=1,ncells(iptr) do m=1,meqn work(m) = qirr(icell,m,iptr) enddo call ctoprm(work,rho,u,p,gamma,zf) c if ((p .le. 0.d0) .or. (rho .le. 0.d0)) then write(6,1002) iptr,ixirr(iptr),icell write(6,1001) nstep,time,rho,p stop endif c c = dsqrt(gamma*p/rho) cmax = dmax1(cmax,dabs(u)+c) enddo else if (irr(i) .eq. -lstgrd) go to 10 do m=1,meqn work(m) = q(i,m) enddo call ctoprm(work,rho,u,p,gamma,zf) c if ((p .lt. 0.d0) .or. (rho .le. 0.d0)) then write(6,*) 'in vartime: rho or p < 0: i=',i write(6,1001) nstep,time,rho,p stop endif endif c c = dsqrt(gamma*p/rho) cmax = dmax1(cmax,dabs(u)+c) 10 continue c dt = cfl*hx/cmax c 1001 format('at step=',i4,', time=',f10.6,', rho=',f12.6,', p=',f12.6) 1002 format('in irreg cell: at iptr=',i3,',i=',i3, & ',is=',i3,', rho or p .le. 0') c return end c c ------------------------------------------------------------------------ c subroutine godunv(q,irr,irrold,s,qjump,mitot,meqn,mwave, & hx,dt,ix1,ixn,iorder,work,lstgrd,lstold) implicit double precision (a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /hypmeth/ mthslp(3),mthlim(3) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),irrold(mitot) dimension s(mitot,mwave),qjump(mitot,meqn,mwave) dimension work(meqn),qja(6),qjb(6),pwave(3) data idebug/0/ c c # Godunov-type method c # update cell averages by wave propagation c do 19 i=ix1,ixn+1 c # ignore irregular cells here if ((irrold(i) .ne. lstold) .or. (irrold(i-1) .ne. lstold)) & go to 19 c x0 = xgrid(i) c c # handle each wave in turn do 10 mw=1,mwave if (idebug .eq. 1) write(6,1002) mw,s(i,mw), & (qjump(i,m,mw),m=1,meqn) c if (s(i,mw) .eq. 0.d0) go to 10 if (s(i,mw) .gt. 0.d0) then ix0 = i do m=1,meqn work(m) = -qjump(i,m,mw) enddo else ix0 = i-1 do m=1,meqn work(m) = qjump(i,m,mw) enddo endif c x1 = x0+s(i,mw)*dt c if (iabs(irr(ix0)) .eq. lstgrd) then c # regular cell case frac = dabs(s(i,mw)*dt)/hx if (frac .gt. 1.d0) then write(6,*) 'CFL violation: cfl=',frac stop endif c do m=1,meqn q(ix0,m) = q(ix0,m)+frac*work(m) enddo else c # irregular cell case pwave(1) = dmin1(x0,x1) pwave(2) = dmax1(x0,x1) call irrwv1(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) endif c if (iorder .eq. 1) go to 10 c c # second order correction: ha = hx do m=1,meqn qja(m) = qjump(i,m,mw) enddo c if (s(i,mw) .gt. 0.d0) then is = i-1 js = is-1 hlen = -hx else is = i+1 js = is hlen = hx endif c hb = hx if (iabs(irrold(js)) .ne. lstold) then iptr = irrold(js) icell = irrcls(xgrid(is),iptr) hb = .5d0*(hx+ar(icell,iptr)) endif c do m=1,meqn qjb(m) = qjump(is,m,mw) enddo call limter(work,qja,qjb,ha,hb,meqn, & mthslp(mw),mthlim(mw)) c c # piecewise linear wave polygon pwave(1) = dmin1(x1,x1+hlen) pwave(3) = dmax1(x1,x1+hlen) pwave(2) = .5d0*(pwave(1)+pwave(3)) c ss = dabs(s(i,mw)) frac = .5d0*ss*dt*(hx-ss*dt)/hx if (frac .gt. 1.d0) then write(6,*) 'CFL violation: cfl=',frac stop endif c do ix0 =i,i-1,-1 if (iabs(irr(ix0)) .eq. lstgrd) then c # regular cell case if (ix0 .eq. i-1) frac = -frac do m=1,meqn q(ix0,m) = q(ix0,m)+frac*work(m) enddo else c # irregular cell case call irrwv2(ix0,q,irr,work,pwave,mitot, & meqn,lstgrd) endif enddo 10 continue 19 continue c 1001 format('q=',5f10.5) 1002 format('mw=',i3,', s=',f10.5,', qjump=',5f10.5) c return end c c ------------------------------------------------------------------ c subroutine vrmjmp(ql,qr,s,qjump,ltrack,mitot,meqn,mwave, & imin,imax,iriemn,napprx) implicit double precision (a-h,o-z) parameter(n1d=2100) dimension s(mitot,mwave),qjump(mitot,meqn,mwave) logical ltrack(mitot,mwave) dimension ql(n1d,6),qr(n1d,6) dimension rhol(n1d),ul(n1d),pl(n1d) dimension rhor(n1d),ur(n1d),pr(n1d) dimension clsql(n1d),clsqr(n1d),pstar(n1d) dimension gam1(n1d),gam2(n1d),zf1(n1d),zf2(n1d) dimension delta(6),qmidl(6),qmidr(6) data itno/10/, small/1.d-8/ c c # vector Riemann solver for Euler equations of gas dynamics c do i=imin,imax do m=1,meqn delta(m) = ql(i,m) enddo call ctoprm(delta,rhol(i),ul(i),pl(i),gam1(i),zf1(i)) c do m=1,meqn delta(m) = qr(i,m) enddo call ctoprm(delta,rhor(i),ur(i),pr(i),gam2(i),zf2(i)) enddo c do i=imin,imax if((rhol(i) .le. 0.d0) .or. (rhor(i) .le. 0.d0) .or. & (pl(i) .le. 0.d0) .or. (pr(i) .le. 0.d0)) then write(6,*) '> negative density or pressure in vrmjmp:', & ' i=',i write(6,*) 'rhol=',rhol(i),', pl=',pl(i),gam1(i) write(6,*) 'rhor=',rhor(i),', pr=',pr(i),gam2(i) stop endif enddo c go to (10,20,30) iriemn c 10 continue c # "exact" Riemann solver- Netwon's iteration version c # see van Leer's 1979 jcp paper, or Colella's 1982 sissc paper c # for more detail c c # iterate to find pstar, ustar: do i=imin,imax clsql(i) = gam1(i)*rhol(i)*pl(i) clsqr(i) = gam2(i)*rhor(i)*pr(i) wl = dsqrt(clsql(i)) wr = dsqrt(clsqr(i)) pstar(i) = (wl*pr(i)+wr*pl(i)-wr*wl* & (ur(i)-ul(i)))/(wl+wr) pstar(i) = dmax1(pstar(i),small) enddo c do 12 i=imin,imax gasc1 = 0.5d0*(gam1(i)+1.d0)/gam1(i) gasc2 = 0.5d0*(gam2(i)+1.d0)/gam2(i) do iter=1,itno wsql = clsql(i)*(1.d0+gasc1*(pstar(i)/pl(i)-1.d0)) wl = dsqrt(wsql) c wsqr = clsqr(i)*(1.d0+gasc2*(pstar(i)/pr(i)-1.d0)) wr = dsqrt(wsqr) c zl = 2.d0*wl*wsql/(wsql+clsql(i)) zr = 2.d0*wr*wsqr/(wsqr+clsqr(i)) c ustarl = ul(i)+(pl(i)-pstar(i))/wl ustarr = ur(i)-(pr(i)-pstar(i))/wr epsp = zl*zr*(ustarr-ustarl)/(zl+zr) pstar(i) = pstar(i)-epsp pstar(i) = dmax1(pstar(i),small) if(dabs(epsp) .le. small) go to 12 enddo c write(6,*) 'in vrmjmp:' write(6,*) 'iteration in Riemann problem does not converge' write(6,*) 'i=', i, ', epsp=',epsp write(6,*) rhol(i),ul(i),pl(i),gam1(i) write(6,*) rhor(i),ur(i),pr(i),gam2(i) stop 12 continue c do 15 i=imin,imax gasc1 = 0.5d0*(gam1(i)+1.d0)/gam1(i) wsql = clsql(i)*(1.d0+gasc1*(pstar(i)/pl(i)-1.d0)) wl = dsqrt(wsql) rho1 = rhol(i)/(1.d0-rhol(i)*(pstar(i)-pl(i))/wsql) c gasc2 = 0.5d0*(gam2(i)+1.d0)/gam2(i) wsqr = clsqr(i)*(1.d0+gasc2*(pstar(i)/pr(i)-1.d0)) wr = dsqrt(wsqr) rho2 = rhor(i)/(1.d0-rhor(i)*(pstar(i)-pr(i))/wsqr) c ustar = (pl(i)-pr(i)+ul(i)*wl+ur(i)*wr)/(wl+wr) c c # left wave: c ============ ltrack(i,1) = pstar(i) .gt. pl(i) s(i,1) = ustar-wl/rho1 call prmtoc(qmidl,rho1,ustar,pstar(i),gam1(i),zf1(i)) do m=1,meqn qjump(i,m,1) = qmidl(m)-ql(i,m) enddo c c # right wave: c ============= ltrack(i,3) = pstar(i) .gt. pr(i) s(i,3) = ustar+wr/rho2 call prmtoc(qmidr,rho2,ustar,pstar(i),gam2(i),zf2(i)) do m=1,meqn qjump(i,m,3) = qr(i,m)-qmidr(m) enddo c c # slip line: c ============ ltrack(i,2) = .true. s(i,2) = ustar do m=1,meqn qjump(i,m,2) = qmidr(m)-qmidl(m) enddo 15 continue return c end c c ------------------------------------------------------------- c subroutine outvar(q,irr,mitot,meqn,hx,ix1,ixn,work, & iunit,lstgrd,iave) implicit double precision (a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),work(meqn) c c # output primitive variables c do i=ix1,ixn x = xgrid(i)+0.5d0*hx if (iabs(irr(i)) .ne. lstgrd) then c # irregular cell if (iave .eq. 1) then write(iunit,*) x, ' NaN NaN NaN NaN NaN' else iptr = irr(i) c do icell=1,ncells(iptr) do m=1,meqn work(m) = qirr(icell,m,iptr) enddo call ctoprm(work,rho,u,p,gamma,zf) c x = xcen(icell,iptr) write(iunit,1001) x,rho,u,p,gamma,zf enddo endif else c # regular cell do m=1,meqn work(m) = q(i,m) enddo call ctoprm(work,rho,u,p,gamma,zf) c write(iunit,1001) x,rho,u,p,gamma,zf endif enddo c 1001 format(f9.5,5f15.8) return end c c ------------------------------------------------------------------- c double precision function philim(a,b,meth) implicit double precision (a-h,o-z) c c # Compute a limiter based on wave strengths a and b. c # meth determines what limiter is used. c # a is assumed to be nonzero. c r = b/a c go to (1,2,3,4,5) meth c 1 continue c # Lax-Wendroff limiter philim = 1.d0 return c 2 continue c # minmod limiter philim = dmax1(0.d0, dmin1(1.d0, r)) return c 3 continue c # superbee limiter philim = dmax1(0.d0, dmin1(1.d0, 2.d0*r), dmin1(2.d0, r)) return c 4 continue c # monotonized centered philim = dmax1(0.d0, dmin1(2.d0, 2.d0*r, & .5d0*(1.d0+r))) return c 5 continue c # van Leer philim = (r+ dabs(r))/(1.d0+dabs(r)) return c end c c -------------------------------------------------------------- c subroutine prmtoc(q,rho,u,p,gamma,z) implicit double precision (a-h,o-z) common /fluid/ geos(5),zf(5),qchem(5),gmole(5),mphase dimension q(6) c c # change primitive variables to conservative variables c q(1) = rho q(2) = rho*u q(3) = p/(gamma-1.d0)+.5d0*q(2)*q(2)/q(1) return c end c c -------------------------------------------------------------- c subroutine ctoprm(q,rho,u,p,gamma) implicit double precision (a-h,o-z) common /fluid/ geos(5),zf(5),qchem(5),gmole(5),mphase dimension q(6) c c # change conservative variables to primitive variables c rho = q(1) u = q(2)/q(1) gamma = geos(1) p = (gamma-1.d0)*(q(3)-.5d0*q(2)*q(2)/q(1)) return c end c c ------------------------------------------------------------------ c integer function maloc(nwords) implicit double precision (a-h,o-z) common /space/ lfree(50,2),lhead(6),lenf,idimf,lentot,lenmax logical sprint data sprint/.false./ c c # allocate contiguous space of length nword in main storage array c # alloc. user code (or pointer to the owner of this storage) c # is mptr. lenf = current length of lfree list. c c # find first fit from free space list itake = 0 do 20 i=1,lenf if (lfree(i,2) .lt. nwords) go to 20 itake = i go to 25 20 continue go to 900 c c # put remaining words back on list. 25 left = lfree(itake,2)-nwords maloc = lfree(itake,1) c c # the following code which is ignored for now adds the new if (left .le. 0) go to 30 lfree(itake,2) = left lfree(itake,1) = lfree(itake,1)+nwords go to 99 c c # item is totally removed. move next items in list up one. 30 lenf = lenf-1 do i=itake,lenf lfree(i,1) = lfree(i+1,1) lfree(i,2) = lfree(i+1,2) enddo go to 99 c 900 write(6,901) nwords write(6,902) ((lfree(i,j),j=1,2),i=1,lenf) stop c 99 lentot = lentot+nwords if (lenmax .lt. lentot) lenmax = lentot if (sprint) write(6,100) nwords, maloc, lentot, lenmax c 100 format(' allocating ',i5,' words in location ',i5, & ' lentot, lenmax ', 2i10) 901 format(' require ',i10,' words - either none left or not big', & ' enough space') 902 format(' free list: ',//,2x,50(i10,4x,i10,/,2x)) c return end c c -------------------------------------------------------------- c subroutine shock(rhol,ul,pl,rhor,ur,pr,mach,gamma) implicit double precision (a-h,o-z) double precision mach c c # given pre-shock state and Mach number c # compute post-shock state c c # Care should be taken in using the following formulas c gasc1 = 2.d0*gamma/(gamma+1.d0) gasc2 = (gamma-1.d0)/(gamma+1.d0) sqma = mach**2 urel = mach*dsqrt(gamma*pr/rhor) sksp = ur-urel c pl = pr*(gasc1*sqma-gasc2) rhol = rhor*(0.5d0*(gamma+1.d0)*sqma)/ & (1.d0+0.5d0*(gamma-1.d0)*sqma) c ul = sksp+urel*(gasc1/gamma/sqma+gasc2) c write(6,*)'Mach number=',mach write(6,1001) rhor,ur,pr write(6,1002) rhol,ul,pl c 1001 format('pre-shock state: rho=',f12.6,', u=',f12.6, & ', p=',f12.6) 1002 format('post-shock state:rho=',f12.6,', u=',f12.6, & ', p=',f12.6) return c end c c --------------------------------------------------------------- c subroutine source(q,irr,mitot,meqn,hx,dt,iorder,work,lstgrd) implicit double precision (a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /gravity/ gravc(2) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) common /sterms/ nsource,isource(10),igeom dimension q(mitot,meqn),irr(mitot),work(meqn) dimension g(6),gprev(6) c c # handle source terms c dt2 = dt if (iorder .eq. 2) dt2 = 0.5d0*dt c do 900 isour=1,nsource go to (1,2,3,4) isource(isour) c 1 continue c # gravitational source terms do i=1,mitot do m=1,meqn work(m) = q(i,m) enddo c call ctoprm(work,rho,u,p,gamma,z) u = u+dt2*gravc(1) c call prmtoc(work,rho,u,p,gamma,z) c do m=1,meqn q(i,m) = work(m) enddo enddo go to 900 c 2 continue c # geometric source term do i=1,mitot if (irr(i) .ne. lstgrd) then iptr = irr(i) do icell=1,ncells(iptr) x = xcen(icell,iptr) gfrac = -fpgeom(x)/fgeom(x) do m=1,meqn work(m) = qirr(icell,m,iptr) enddo call ctoprm(work,rho,u,p,gamma,z) c if (iorder .eq. 1) then c # first order ODE solver g(1) = gfrac*work(2) g(2) = gfrac*work(2)**2/work(1) g(3) = gfrac*(work(3)+p)*u g(4) = gfrac*work(4)*u do m=1,meqn qirr(icell,m,iptr) = qirr(icell,m,iptr)+ & dt*g(m) enddo else c # second order ODE solver gprev(1) = gfrac*work(2) gprev(2) = gfrac*work(2)**2/work(1) gprev(3) = gfrac*(work(3)+p)*u gprev(4) = gfrac*work(4)*u do m=1,meqn work(m) = work(m)+dt2*gprev(m) enddo call ctoprm(work,rho,u,p,gamma,z) c g(1) = gfrac*work(2) g(2) = gfrac*work(2)**2/work(1) g(3) = gfrac*(work(3)+p)*u g(4) = gfrac*work(4)*u do m=1,meqn qirr(icell,m,iptr) = qirr(icell,m,iptr)+ & .5d0*dt2*(g(m)+gprev(m)) enddo endif enddo else x = xgrid(i)+0.5d0*hx gfrac = -fpgeom(x)/fgeom(x) do m=1,meqn work(m) = q(i,m) enddo call ctoprm(work,rho,u,p,gamma,z) c if (iorder .eq. 1) then c # first order ODE solver g(1) = gfrac*work(2) g(2) = gfrac*work(2)**2/work(1) g(3) = gfrac*(work(3)+p)*u g(4) = gfrac*work(4)*u do m=1,meqn q(i,m) = q(i,m)+dt*g(m) enddo else c # second order ODE solver gprev(1) = gfrac*work(2) gprev(2) = gfrac*work(2)**2/work(1) gprev(3) = gfrac*(work(3)+p)*u gprev(4) = gfrac*work(4)*u c do m=1,meqn work(m) = work(m)+dt2*gprev(m) enddo call ctoprm(work,rho,u,p,gamma,z) c g(1) = gfrac*work(2) g(2) = gfrac*work(2)**2/work(1) g(3) = gfrac*(work(3)+p)*u g(4) = gfrac*work(4)*u do m=1,meqn q(i,m) = q(i,m)+.5d0*dt2*(g(m)+gprev(m)) enddo endif endif enddo go to 900 c 3 continue c # reacting source terms c # the ODE is solved exactly go to 900 c 4 continue c # radiation source terms go to 900 900 continue c return end c c ----------------------------------------------------------------------- c subroutine limter(work,qja,qjb,ha,hb,meqn,islope,meth) implicit double precision (a-h,o-z) dimension work(meqn),qja(6),qjb(6) c c # compute limited slope c do m=1,meqn work(m) = 0.d0 enddo c go to (1,2) islope c 1 continue c # limit over each component of jumps separately do m=1,meqn if (qja(m) .ne. 0.d0) & work(m) = philim(qja(m)/ha,qjb(m)/hb,meth)* & qja(m)/ha enddo return c 2 continue c # limit over each characteristic mode of jumps separately wnorm = 0.d0 wvdot = 0.d0 do m=1,meqn wnorm = wnorm+qja(m)*qja(m) wvdot = wvdot+qja(m)*qjb(m) enddo wnorm = dsqrt(wnorm)/ha wvdot = wvdot/(ha*hb) if (wnorm .eq. 0.d0) return c wlimitr = philim(wnorm,wvdot/wnorm,meth) do m=1,meqn work(m) = wlimitr*qja(m)/ha enddo return c end c c ----------------------------------------------------------------- c subroutine rfront(irr,mitot,meqn,mwave,work,lstfrt,iriem) implicit double precision (a-h,o-z) parameter(maxfr=10,maxbd=126) common /frpsoln/ sfront(maxbd,3),qjfront(maxbd,6,3), & mwtrack(maxbd,3) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension irr(mitot),work(meqn) dimension ql(6),qr(6),s(3),qjump(6,3) logical ltrack(3),mwtrack data idebug/0/ c c # solve Riemann problems at tracked fronts c if(idebug .eq. 1) write(6,*) 'entering rfront routine' c kfront = lstfrt do 10 ifr=1,maxfr kfront = iabs(nxtfront(kfront)) if(kfront .eq. 0) go to 19 c if(idebug .eq. 1) & write(6,1008) kfront,ityfront(kfront) c c # handle each tracked point in turn kbond = lstbond(kfront) do ibd=1,maxbd kbond = iabs(nxtbond(kbond)) if(kbond .eq. 0) go to 10 ix = ixbond(kbond) iptr = irr(ix) c c # Riemann data at the front: irrl = irrcell(1,kbond) irrr = irrcell(2,kbond) c if (itybond(kbond) .eq. 2) then c # boundary interface do m=1,meqn ql(m) = qirr(irrr,m,iptr) qr(m) = qirr(irrr,m,iptr) enddo ql(2) = -qr(2) else c # ordinary tracked interface (eg, for shock) do m=1,meqn ql(m) = qirr(irrl,m,iptr) qr(m) = qirr(irrr,m,iptr) enddo endif c call rmnjmp(ql,qr,s,qjump,ltrack,meqn,iriem) c do mw=1,mwave sfront(kbond,mw) = s(mw) mwtrack(kbond,mw) = ltrack(mw) do m=1,meqn qjfront(kbond,m,mw) = qjump(m,mw) enddo enddo c if (idebug .eq. 1) then write(6,1003) iptr,ix,kbond, & irrcell(1,kbond),irrcell(2,kbond) write(6,1002) (ql(m),m=1,meqn) write(6,1002) (qr(m),m=1,meqn) do mw=1,mwave write(6,1001) mw,s(mw),(qjump(m,mw),m=1,meqn) enddo endif enddo 10 continue 19 continue c if (idebug .eq. 1) & write(6,*)'leaveing rfront:---------------------' c 1001 format('mw=',i3,', s=',f10.5,', qj=',5f10.5) 1002 format(5f12.6) 1003 format('iptr=',i3,', ix=',i3,', kbond=',i3, & ', icl=',i3,', icr=',i3) 1008 format('Riemann solutions at front',i3,', type=',i3) return c end c c ---------------------------------------------------------------- c subroutine mkgrid(q,qold,irr,irrold,mitot,meqn,work,lstgrd, & lstold) implicit double precision (a-h,o-z) parameter(maxbd=126) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),qold(mitot,meqn) dimension irr(mitot),irrold(mitot),work(meqn) c c # make new grid c c # save a copy of the old grid states on uniform grid call irrst0(q,mitot,meqn,work,lstgrd) c c # copy the state values do m=1,meqn do i=1,mitot qold(i,m) = q(i,m) enddo enddo c c # copy irregular cell pointers do i=1,mitot irrold(i) = irr(i) enddo lstold = lstgrd c c # set new grid pointer lstgrd = lstptr(1) nxtirr(lstgrd) = 0 do i=1,mitot if (irrold(i) .eq. -lstold) then irr(i) = -lstgrd else irr(i) = lstgrd endif enddo ar(1,lstgrd) = ar(1,lstold) c return end c c --------------------------------------------------------------- c subroutine irrst0(q,mitot,meqn,work,lstgrd) implicit double precision (a-h,o-z) parameter(maxbd=126) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),work(meqn) c c # compute irregular cells' cell averages over the uniform grid c iptr = lstgrd do 100 ik=1,maxbd iptr = iabs(nxtirr(iptr)) if (iptr .eq. 0) go to 199 ix = ixirr(iptr) c area = 0.d0 do m=1,meqn work(m) = 0.d0 enddo c ibeg = 1 if (ityirr(iptr) .eq. 2) ibeg = 2 c c # compute cell average do icell=ibeg,ncells(iptr) do m=1,meqn work(m) = work(m)+ar(icell,iptr)*qirr(icell,m,iptr) enddo area = area+ar(icell,iptr) enddo c do m=1,meqn q(ix,m) = work(m)/area enddo 100 continue 199 continue c return end c c ----------------------------------------------------------------------- c subroutine setgrd(q,irr,mitot,meqn,hx,lstgrd,lstfrt) implicit double precision (a-h,o-z) parameter(n1d=2100,maxfr=10,maxbd=126) common /coord/ xgrid(0:n1d) common /frontin/ xloop(10),varloop(2,6,10), & ibdloop(10),mwloop(10),istloop(10), & newloop(10),irrt0(10),irrtn,nloops common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) common /local/ iloops dimension q(mitot,meqn),irr(mitot) c c # set irregular cells c c # data structure c # starting grid pointer lstgrd = lstptr(1) nxtirr(lstgrd) = 0 c c # starting front pointer lstfrt = lstptr(2) nxtfront(lstfrt) = 0 c c # initialize irr and area arrays do i=1,mitot irr(i) = lstgrd enddo c do icell=1,6 do k=1,maxbd ar(icell,k)= hx enddo enddo c do 10 iloops=1,nloops if (istloop(iloops) .eq. 1) then c # tracking case if (newloop(iloops) .eq. 1) then c # initialize new front structure newfr = maptr(2) newfront(newfr) = 0 ityfront(newfr) = ibdloop(iloops) c newbd = lstptr(6) lstbond(newfr) = newbd nxtbond(newbd) = 0 endif endif c xp = xloop(iloops) ix = int((xp-xgrid(0))/hx) c if((ix .lt. 0) .or. (ix .gt. mitot+1)) go to 10 c c # make irregular cell call mkirrc(irr,ix,xp,mitot,lstgrd) c c # set values for the irregular cells call setqir(q,irr,mitot,meqn,hx,lstgrd,lstgrd,irrt0(iloops)) c c # set bond type kbond = lstbond(newfr) do ibd=1,maxbd kbond = nxtbond(kbond) if (kbond .eq. 0) go to 10 itybond(kbond) = ibdloop(iloops) enddo 10 continue c return end c c ------------------------------------------------------------------ c subroutine rmnjmp(ql,qr,s,qjump,ltrack,meqn,iriemn) implicit double precision (a-h,o-z) logical ltrack(3) dimension s(3),qjump(6,3) dimension ql(6),qr(6),delta(6),qmidl(6),qmidr(6) data itno/10/, small/1.d-8/ c c # single Riemann solver for Euler equations of gas dynamics c do m=1,meqn delta(m) = ql(m) enddo call ctoprm(delta,rhol,ul,pl,gam1,zf1) c do m=1,meqn delta(m) = qr(m) enddo call ctoprm(delta,rhor,ur,pr,gam2,zf2) c if((rhol .le. 0.d0) .or. (rhor .le. 0.d0) .or. & (pl .le. 0.d0) .or. (pr .le. 0.d0)) then write(6,*) '> negative density or pressure in rmnjmp:' write(6,*) 'rhol=',rhol,', pl=',pl,', gam1=',gam1 write(6,*) 'rhor=',rhor,', pr=',pr,', gam2=',gam2 stop endif c go to (10,20,30) iriemn c 10 continue c # "exact" Riemann solver - version of Netwon's iteration c # see van Leer's 1979 jcp paper, or Colella's 1982 sissc paper c # for more detail c c # iterate to find pstar, ustar: clsql = gam1*rhol*pl clsqr = gam2*rhor*pr wl = dsqrt(clsql) wr = dsqrt(clsqr) pstar = (wl*pr+wr*pl-wr*wl*(ur-ul))/(wl+wr) pstar = dmax1(pstar,small) c gasc1 = 0.5d0*(gam1+1.d0)/gam1 gasc2 = 0.5d0*(gam2+1.d0)/gam2 do iter=1,itno wsql = clsql*(1.d0+gasc1*(pstar/pl-1.d0)) wl = dsqrt(wsql) c wsqr = clsqr*(1.d0+gasc2*(pstar/pr-1.d0)) wr = dsqrt(wsqr) c zl = 2.d0*wl*wsql/(wsql+clsql) zr = 2.d0*wr*wsqr/(wsqr+clsqr) c ustarl = ul+(pl-pstar)/wl ustarr = ur-(pr-pstar)/wr epsp = zl*zr*(ustarr-ustarl)/(zl+zr) pstar = pstar-epsp pstar = dmax1(pstar,small) if(dabs(epsp) .le. small) go to 12 enddo c write(6,*) 'iteration in rmnjmp does not converge: epsp=', epsp write(6,*) rhol,ul,pl,gam1 write(6,*) rhor,ur,pr,gam2 stop 12 continue c wsql = clsql*(1.d0+gasc1*(pstar/pl-1.d0)) wl = dsqrt(wsql) rho1 = rhol/(1.d0-rhol*(pstar-pl)/wsql) c wsqr = clsqr*(1.d0+gasc2*(pstar/pr-1.d0)) wr = dsqrt(wsqr) rho2 = rhor/(1.d0-rhor*(pstar-pr)/wsqr) c ustar = (pl-pr+ul*wl+ur*wr)/(wl+wr) c c # left wave: c ============ ltrack(1) = pstar .gt. pl s(1) = ustar-wl/rho1 call prmtoc(qmidl,rho1,ustar,pstar,gam1,zf1) do m=1,meqn qjump(m,1) = qmidl(m)-ql(m) enddo c c # right wave: c ============= ltrack(3) = pstar .gt. pr s(3) = ustar+wr/rho2 call prmtoc(qmidr,rho2,ustar,pstar,gam2,zf2) do m=1,meqn qjump(m,3) = qr(m)-qmidr(m) enddo c c # slip line: c ============ ltrack(2) = .true. s(2) = ustar do m=1,meqn qjump(m,2) = qmidr(m)-qmidl(m) enddo return c c ------------------------------------------------------------- c integer function lstptr(itype) implicit double precision (a-h,o-z) parameter(maxfr=10,maxbd=126) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) common /space/ lfree(50,2),lhead(6),lenf,idimf,lentot,lenmax c c # first free spot in linked list of irregular info. pointed to c # by lhead. c if (lhead(itype) .ne. 0) go to 10 write(6,1001) itype stop c 10 lstptr = lhead(itype) if (itype .eq. 1) then lhead(itype) = iabs(nxtirr(lhead(itype))) elseif (itype .eq. 2) then lhead(itype) = iabs(nxtfront(lhead(itype))) nxtfront(lstptr) = 0 elseif (itype .eq. 6) then lhead(itype) = iabs(nxtbond(lhead(itype))) else write(6,*) 'no pointer available for itype=',itype stop endif c 1001 format('out of irregular list space: type=',i3) return end c c --------------------------------------------------------------- c subroutine rpdata(ql,qr,q,irr,rwork,mitot,meqn,imin,imax, & lstgrd) implicit double precision(a-h,o-z) parameter(maxbd=126,n1d=2100) common /coord/ xgrid(0:n1d) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),rwork(mitot) dimension ql(n1d,6),qr(n1d,6) data idebug/0/ c c # Prepare Riemann data at regular cell interfaces c do i=imin,imax do m=1,meqn ql(i,m) = q(i-1,m) qr(i,m) = q(i,m) enddo enddo c c # re-adjust Riemann data ql and qr c # if there are irregular cells next to it c do k=1,mitot rwork(k) = 0.d0 enddo c iptr = lstgrd do 100 ik=1,maxbd iptr = iabs(nxtirr(iptr)) if (iptr .eq. 0) go to 199 ix = ixirr(iptr) c c # left hand side if (rwork(ix) .eq. 1.d0) go to 70 c if (idebug .eq. 1) then write(6,*) 'in rpdata: ix= ',ix write(6,1001) (ql(ix,m),m=1,meqn) write(6,1001) (qr(ix,m),m=1,meqn) endif c rwork(ix) = 1.d0 iptr1 = irr(ix-1) if (iabs(iptr1) .eq. lstgrd) then icell = irrcls(xgrid(ix),iptr) do m=1,meqn qr(ix,m) = qirr(icell,m,iptr) enddo else icell = irrcls(xgrid(ix),iptr1) jcell = irrcls(xgrid(ix),iptr) do m=1,meqn ql(ix,m) = qirr(icell,m,iptr1) qr(ix,m) = qirr(jcell,m,iptr) enddo endif c if (idebug .eq. 1) then write(6,*) 'out rpdata:' write(6,1001) (ql(ix,m),m=1,meqn) write(6,1001) (qr(ix,m),m=1,meqn) endif c 70 continue c # right hand side if (rwork(ix+1) .eq. 1.d0) go to 100 c if (idebug .eq. 1) then write(6,*) 'in rpdata: ix+1=',ix+1 write(6,1001) (ql(ix+1,m),m=1,meqn) write(6,1001) (qr(ix+1,m),m=1,meqn) endif c rwork(ix+1) = 1.d0 iptr1 = irr(ix+1) if (iabs(iptr1) .eq. lstgrd) then icell = irrcls(xgrid(ix+1),iptr) do m=1,meqn ql(ix+1,m) = qirr(icell,m,iptr) enddo else icell = irrcls(xgrid(ix+1),iptr) jcell = irrcls(xgrid(ix+1),iptr1) do m=1,meqn ql(ix+1,m) = qirr(icell,m,iptr) qr(ix+1,m) = qirr(jcell,m,iptr1) enddo endif c if (idebug .eq. 1) then write(6,*) 'out rpdata:' write(6,1001) (ql(ix+1,m),m=1,meqn) write(6,1001) (qr(ix+1,m),m=1,meqn) endif c 100 continue 199 continue c 1001 format('q =',6f10.5) return end c c ----------------------------------------------------------------- c subroutine ffront(s,qjump,ltrack,mitot,meqn,mwave,imin,imax, & time,work,lstfrt) implicit double precision (a-h,o-z) parameter(maxfr=10,maxbd=126,n1d=2100) common /coord/ xgrid(0:n1d) common /frpsoln/ sfront(maxbd,3),qjfront(maxbd,6,3), & mwtrack(maxbd,3) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /trackin/ qjtol(3),itrack,idetec(3),isurgy common /stags/ oldpt(maxbd),sppt(maxbd),mwpt(maxbd),ntags dimension s(mitot,mwave),qjump(mitot,meqn,mwave),work(meqn) logical ltrack(mitot,mwave),mwtrack integer sdetec data idebug/1/ c c # flag strong waves c # store their locations, wave number, and speeds c c # criterion for choosing moving points: (for shock or contact) c # if either of the following quantity is greater than c # some prescribed tolerance, c # (1) the maximum norm of the 'qjump' c # (2) the first component of the 'qjump' c k = 0 c c # tracking points finding options go to (1,2,3) itrack+1 c 1 continue c # no tracking return c 2 continue c # employ shock detection algorithm to determining c # tracked points c c # check tracked cell interfaces kfront = lstfrt do 200 ifr=1,maxfr kfront = iabs(nxtfront(kfront)) if(kfront .eq. 0) go to 209 c kbond = lstbond(kfront) do ibd=1,maxbd kbond = iabs(nxtbond(kbond)) if(kbond .eq. 0) go to 200 c do mw=1,mwave if (mwtrack(kbond,mw) .eq. .true.) then do m=1,meqn work(m) = qjfront(kbond,m,mw) enddo c if (sdetec(work,meqn,idetec(mw),qjtol(mw)) & .eq. 1) then k = k+1 oldpt(k) = points(kbond) sppt(k) = sfront(kbond,mw) mwpt(k) = mw else mwtrack(kbond,mw) = .false. endif endif enddo enddo 200 continue 209 continue c c # check fixed cell boundaries do i=imin,imax do mw=1,mwave if (ltrack(i,mw) .eq. .true.) then do m=1,meqn work(m) = qjump(i,m,mw) enddo c if (sdetec(work,meqn,idetec(mw),qjtol(mw)) & .eq. 1) then k = k+1 oldpt(k) = xgrid(i) sppt(k) = s(i,mw) mwpt(k) = mw else ltrack(i,mw) = .false. endif endif enddo enddo go to 900 c 3 continue c # track specified moving points kfront = lstfrt do 300 ifr=1,maxfr kfront = iabs(nxtfront(kfront)) if (kfront .eq. 0) go to 309 c if (idebug .eq. 1) & write(6,1001) kfront,ityfront(kfront), & mwfront(kfront) c kbond = lstbond(kfront) do ibd=1,maxbd kbond = iabs(nxtbond(kbond)) if (kbond .eq. 0) go to 300 c do mw=1,mwave mwtrack(kbond,mw) = .false. enddo mwtrack(kbond,mwfront(kfront)) = .true. c k = k+1 oldpt(k) = points(kbond) sppt(k) = sfront(kbond,mwfront(kfront)) mwpt(k) = mwfront(kfront) enddo 300 continue 309 continue go to 900 c 900 continue ntags = k if (ntags .eq. 0) write(6,1002) time c if (idebug .eq. 1) then write(6,*) 'entering ffront:' do k=1,ntags write(6,1003) k,oldpt(k),sppt(k),mwpt(k) enddo endif c 1001 format('front',i3,', type=',i3,', wave-type=',i3) 1002 format('no moving point at time t=',f13.6) 1003 format('k=',i3,', xloc=',f12.6,', s=',f12.6, & ', mw=',i3) c return end c c ------------------------------------------------------------------- c integer function sdetec(work,meqn,idetec,tol) implicit double precision (a-h,o-z) dimension work(meqn) c c # shock detection criterion c sdetec = 0 c go to (1,2) idetec c 1 continue c # density strength qmax = dabs(work(1)) go to 99 c 2 continue c # maximum norm of the wave strength "work" qmax = 0.d0 do m=1,meqn qmax = dmax1(dabs(qmax),dabs(work(m))) enddo go to 99 c 99 continue if (qmax .ge. tol) sdetec = 1 return c end c c -------------------------------------------------------------------- c subroutine celave(ix,hx,wl) implicit double precision (a-h,o-z) parameter(n1d=2100) common /coord/ xgrid(0:n1d) common /fshockc/ xp dimension poly(2) logical alll,allr,fl(2) c c # compute wl, fraction of cell that lies in left state. c # for initial data with two states ql and qr separated by a c # discontinuity at the point "xp" c poly(1) = xgrid(ix) poly(2) = xgrid(ix+1) alll = .true. allr = .true. c do i=1,2 fl(i) = poly(i)-xp .lt. 0.d0 alll = alll .and. fl(i) allr = allr .and. (.not. fl(i)) enddo c if (alll .eq. .true.) then wl = 1.d0 elseif (allr .eq. .true.) then wl = 0.d0 else area = xp-poly(1) wl = area/hx endif c return end c c ----------------------------------------------------------------- c subroutine pfront(q,irr,irrold,s,qjump,rwork,mitot,meqn, & mwave,hx,dt,iorder,work,lstgrd,lstold) implicit double precision(a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /frpsoln/ sfront(maxbd,3),qjfront(maxbd,6,3), & mwtrack(maxbd,3) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),irrold(mitot) dimension s(mitot,mwave),qjump(mitot,meqn,mwave) dimension rwork(mitot),work(meqn),pwave(3) data idebug/0/ c c # Update cell values based on solutions of Riemann problems c # from irregular cell interfaces c if (idebug .ne. 0) then write(6,*) 'entering pfront:' iptr = lstgrd do ik=1,maxbd iptr = nxtirr(iptr) if(iptr .eq. 0) go to 919 c write(6,1001) iptr,ixirr(iptr),ityirr(iptr) c do icell=1,ncells(iptr) write(6,1003) ar(icell,iptr), & (qirr(icell,m,iptr),m=1,meqn) enddo enddo 919 continue write(6,*) '-----------------------------------------' endif c do k=1,mitot rwork(k) = 0.d0 enddo c iptr = lstold do 100 ik=1,maxbd iptr = iabs(nxtirr(iptr)) if(iptr .eq. 0) go to 199 c ix = ixirr(iptr) c if (idebug .eq. 1) write(6,*) 'for ix=',ix c c # handle tracked interface first do 40 ibd=1,nbonds(iptr) kbond = nbondc(ibd,iptr) x0 = points(kbond) c if (idebug .eq. 1) write(6,*) 'tracked interface' c do 30 mw=1,mwave if (dabs(sfront(kbond,mw)) .lt. 1.d-12) go to 30 c if (sfront(kbond,mw) .gt. 0.d0) then ix1 = ix do m=1,meqn work(m) = -qjfront(kbond,m,mw) enddo else ix1 = ix-1 do m=1,meqn work(m) = qjfront(kbond,m,mw) enddo endif ix2 = ix1+1 c x1 = x0+sfront(kbond,mw)*dt pwave(1) = dmin1(x0,x1) pwave(2) = dmax1(x0,x1) c c # update cell values do ix0=ix1,ix2 call irrwv1(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) enddo c if (iorder .eq. 1) go to 30 c c # second order correction: call irrslp(work,hlen,irrold,sfront(kbond,mw), & qjump,mitot,meqn,mwave,ix,ix,mw,kbond, & lstold,.false.) c pwave(1) = dmin1(x1,x1+hlen) pwave(3) = dmax1(x1,x1+hlen) pwave(2) = .5d0*(pwave(1)+pwave(3)) c do ix0=ix1,ix2 call irrwv2(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) enddo 30 continue 40 continue c 49 continue c # handle fixed cell bdry c # left hand side if (rwork(ix) .eq. 1.d0) go to 60 rwork(ix) = 1.d0 c x0 = xgrid(ix) ix1 = ix-1 ix2 = ix c if (idebug .eq. 1) write(6,*) 'left interface' c do 59 mw=1,mwave if (dabs(s(ix,mw)) .lt. 1.d-12) go to 59 c if (s(ix,mw) .gt. 0.d0) then ix0 = ix2 do m=1,meqn work(m) = -qjump(ix,m,mw) enddo else ix0 = ix1 do m=1,meqn work(m) = qjump(ix,m,mw) enddo endif c x1 = x0+s(ix,mw)*dt pwave(1) = dmin1(x0,x1) pwave(2) = dmax1(x0,x1) c c # update cell values: first order call irrwv1(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) c if(iorder .eq. 1) go to 59 c c # second order correction call irrslp(work,hlen,irrold,s(ix,mw),qjump, & mitot,meqn,mwave,ix1,ix2,mw,idummy, & lstold,.true.) c pwave(1) = dmin1(x1,x1+hlen) pwave(3) = dmax1(x1,x1+hlen) pwave(2) = .5d0*(pwave(1)+pwave(3)) c c # update cell values: second order do ix0=ix1,ix2 call irrwv2(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) enddo c c # boundary treatment 59 continue c 60 continue c # right hand side if (rwork(ix+1) .eq. 1.d0) go to 100 rwork(ix+1) = 1.d0 c x0 = xgrid(ix+1) ix1 = ix ix2 = ix+1 c if (idebug .eq. 1) write(6,*) 'right interface' c do 79 mw=1,mwave if (dabs(s(ix+1,mw)) .lt. 1.d-12) go to 79 c if (s(ix+1,mw) .gt. 0.d0) then ix0 = ix2 do m=1,meqn work(m) = -qjump(ix+1,m,mw) enddo else ix0 = ix1 do m=1,meqn work(m) = qjump(ix+1,m,mw) enddo endif c x1 = x0+s(ix+1,mw)*dt pwave(1) = dmin1(x0,x1) pwave(2) = dmax1(x0,x1) c c # update cell values: first order call irrwv1(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) c c # boundary treatment c if(iorder .eq. 1) go to 79 c c # second order correction call irrslp(work,hlen,irrold,s(ix+1,mw),qjump, & mitot,meqn,mwave,ix1,ix2,mw,idummy, & lstold,.true.) c pwave(1) = dmin1(x1,x1+hlen) pwave(3) = dmax1(x1,x1+hlen) pwave(2) = .5d0*(pwave(1)+pwave(3)) c c # update cell values: second order do ix0=ix1,ix2 call irrwv2(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) enddo c c # boundary treatment 79 continue 100 continue 199 continue c if (idebug .ne. 0) then write(6,*) 'leaving pfront:' iptr = lstgrd do ik=1,maxbd iptr = nxtirr(iptr) if(iptr .eq. 0) go to 929 c write(6,1001) iptr,ixirr(iptr),ityirr(iptr) c do icell=1,ncells(iptr) write(6,1003) ar(icell,iptr), & (qirr(icell,m,iptr),m=1,meqn) enddo enddo 929 continue write(6,*) '-----------------------------------------' endif c 1001 format('iptr=',i4,', ix=',i3,', type=',i3) 1003 format(e12.6,2x,6f10.5) c return end c c --------------------------------------------------------------- c subroutine ifront(irr,irrold,mitot,hx,time,dt,lstgrd,lstfrt) implicit double precision(a-h,o-z) parameter(n1d=2100,maxfr=10,maxbd=126) common /coord/ xgrid(0:n1d) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) common /stags/ oldpt(maxbd),sppt(maxbd),mwpt(maxbd),ntags dimension irr(mitot),irrold(mitot) c c # advance tracked fronts to the next time step c # and insert them to the new grid c c # handle stationary fronts first kfront = lstfrt do 20 ifr=1,maxfr kfront = nxtfront(kfront) if (kfront .eq. 0) go to 29 if ((newfront(kfront) .eq. 1) .or. & (ityfront(kfront) .eq. 1)) go to 20 c c # set new front structure newfr = maptr(2) newbd = lstptr(6) lstbond(newfr) = newbd nxtbond(newbd) = 0 c c # copy stationary front and create new grid kbond = lstbond(kfront) do ibd=1,maxbd kbond = nxtbond(kbond) if (kbond .eq. 0) go to 19 ix = ixbond(kbond) xp = points(kbond) c c # make irregular cell call mkirrc(irr,ix,xp,mitot,lstgrd) c ityirr(irr(ix)) = ityirr(irrold(ix)) enddo 19 continue c c # set front type ityfront(newfr) = ityfront(kfront) kbond = lstbond(newfr) do ibd=1,maxbd kbond = nxtbond(kbond) if (kbond .eq. 0) go to 20 itybond(kbond) = ityfront(kfront) c c # output front write(97,1001) points(kbond),time, & points(kbond),time+dt enddo 20 continue 29 continue c c # handle tracked moving fronts if (ntags .ne. 0) then c # new front structure newfr = maptr(2) newbd = lstptr(6) lstbond(newfr) = newbd nxtbond(newbd) = 0 endif c do k=1,ntags xp = oldpt(k)+sppt(k)*dt ix = int((xp-xgrid(0))/hx) c if ((ix .gt. 0) .or. (ix .lt. mitot+1)) then write(97,1001) oldpt(k),time,xp,time+dt c c # make irregular cell call mkirrc(irr,ix,xp,mitot,lstgrd) endif enddo c 1001 format(4f16.8) return end c c ----------------------------------------------------------------- c subroutine dfront(lstold,lstfrt) implicit double precision (a-h,o-z) parameter(maxfr=10) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) c c # delete old front and reset new front c c # set irregular cells link list c # remove old grid from the list call putptr(lstold,1) c c # remove old front from the front list kfront = lstfrt do ifr=1,maxfr kfront = nxtfront(kfront) if (kfront .eq. 0) go to 19 if (newfront(kfront) .eq. 0) then call putptr(lstbond(kfront),6) call putptr(kfront,2) endif enddo 19 continue c c # set new front flag to old kfront = lstfrt do ifr=1,maxfr kfront = nxtfront(kfront) if (kfront .eq. 0) go to 29 newfront(kfront) = 0 enddo 29 continue c return end c c ------------------------------------------------------------- c subroutine putptr(iptr,itype) implicit double precision (a-h,o-z) parameter(maxfr=10,maxbd=126) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) common /space/ lfree(50,2),lhead(6),lenf,idimf,lentot,lenmax c c # put pointer "iptr" back to the free list c if(itype .eq. 1) then c # for irregular cell kirr = iptr 10 if(nxtirr(kirr) .eq. 0) go to 12 kirr = nxtirr(kirr) go to 10 c 12 nxtirr(kirr) = lhead(1) iprirr(lhead(1)) = kirr lhead(1) = iptr elseif(itype .eq. 2) then c # for tracked front kpre = iprfront(iptr) knxt = iabs(nxtfront(iptr)) nxtfront(kpre) = knxt if(knxt .eq. 0) go to 20 iprfront(knxt) = kpre c c # insert iptr to the free list klst = iprfront(lhead(2)) iprfront(iptr) = klst 20 continue nxtfront(iptr) = lhead(2) iprfront(lhead(2)) = iptr lhead(2) = iptr iptr = kpre elseif(itype .eq. 6) then c # for tracked bond kbond = iptr 60 if(nxtbond(kbond) .eq. 0) go to 62 kbond = nxtbond(kbond) go to 60 c 62 nxtbond(kbond) = lhead(6) iprbond(lhead(6)) = kbond lhead(6) = iptr else write(6,*) 'no pointer available for itype=',itype stop endif c return 1001 format(3(i5)) c end c c -------------------------------------------------------------------- c subroutine irrwv1(ix,q,irr,work,pwave,mitot,meqn,lstgrd) implicit double precision(a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),work(meqn) dimension pwave(3),p1(3),p2(3) data idebug/0/ c c # first order pwave propagation on the irregular cells c # given polygon of pwave and the index of the cells which may c # be affected by this pwave, c # then find the intersection area and use it to update c # the cell-averaged values c if((ix .lt. 1) .or. (ix .gt. mitot)) return c if (idebug .eq. 1) write(6,1001) (work(m), m=1,meqn) c iptr = irr(ix) c if (iptr .eq. -lstgrd) return c c # find the intersection of pwave with this polygon (ix,jy0) if (iptr .eq. lstgrd) then c # regular cell case p1(1) = xgrid(ix) p1(2) = xgrid(ix+1) c p2(1) = dmax1(pwave(1),p1(1)) p2(1) = dmin1(p2(1),p1(2)) p2(2) = dmin1(pwave(2),p1(2)) p2(2) = dmax1(p2(2),p1(1)) c area = p2(2)-p2(1) c if (idebug .eq. 1) then write(6,1002) lstgrd,ix,area,ar(1,lstgrd) write(6,1003) (q(ix,m),m=1,meqn) endif c frac = area/ar(1,lstgrd) do m=1,meqn q(ix,m) = q(ix,m)+frac*work(m) enddo c if (idebug .eq. 1) write(6,1003) (q(ix,m),m=1,meqn) else c # irregular cell case c # find the intersection area with each subcell do icell=1,ncells(iptr) p1(1) = poly(icell,1,iptr) p1(2) = poly(icell,2,iptr) c p2(1) = dmax1(pwave(1),p1(1)) p2(1) = dmin1(p2(1),p1(2)) p2(2) = dmin1(pwave(2),p1(2)) p2(2) = dmax1(p2(2),p1(1)) c area = p2(2)-p2(1) c if (idebug .eq. 1) then write(6,1002) iptr,icell,area,ar(icell,iptr) write(6,1003) (qirr(icell,m,iptr),m=1,meqn) endif c frac = area/ar(icell,iptr) do m=1,meqn qirr(icell,m,iptr) = qirr(icell,m,iptr)+ & frac*work(m) enddo c if (idebug .eq. 1) write(6,1003) & (qirr(icell,m,iptr),m=1,meqn) enddo endif c if (idebug .eq. 1) & write(6,*)'-----------------------------------------' c 1001 format('entering irrwv1: jump=',5f10.5) 1002 format('iptr=',i3,', icell=',i3,', area=',f12.6,', artot=',f12.6) 1003 format('q=',5f12.6) c return end c c ----------------------------------------------------------------- c double precision function fgeom(x) implicit double precision (a-h,o-z) common /sterms/ nsource,isource(10),igeom c c # igeom: c # 1 : divergent duct c # 2 : convergent duct c # 3 : convergent-divergent duct c # 4 : cylindrical symmetry c # 5 : spherical symmetry c go to (1,2,3,4,5) igeom c 1 continue fgeom = 1.398d0+.347d0*dtanh(8.d0*x-4.d0) return c 2 continue fgeom = 1.398d0+.347d0*dtanh(-8.d0*x+4.d0) return c 3 continue if(x .gt. 0.d0) then fgeom = 1.398d0+.347d0*dtanh(8.d0*x-4.d0) else x = x+1.d0 fgeom = 1.398d0+.347d0*dtanh(-8.d0*x+4.d0) endif return c 4 continue fgeom = x return c 5 continue fgeom = x*x return c end c c -------------------------------------------------------------------- c subroutine irrwv2(ix0,q,irr,work,pwave,mitot,meqn,lstgrd) implicit double precision(a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension q(mitot,meqn),irr(mitot),work(meqn) dimension pwave(3),p1(3),p2(3) data idebug/0/ c c # first order pwave propagation on irregular cells c # given wave polygon and index of the cell which may c # be affected, c # find intersection area of the wave polygon with grid cell, c # and use it to update cell average c if ((ix0 .lt. 1) .or. (ix0 .gt. mitot)) return c if (idebug .eq. 1) write(6,1001) (work(m),m=1,meqn) c iptr = irr(ix0) c if (iptr .eq. -lstgrd) return c c # find the intersection of pwave with this polygon (ix0,jy0) if (iptr .eq. lstgrd) then c # regular cell p1(1) = xgrid(ix0) p1(2) = xgrid(ix0+1) c p2(1) = dmax1(pwave(1),p1(1)) p2(1) = dmin1(p2(1),p1(2)) p2(2) = dmin1(pwave(3),p1(2)) p2(2) = dmax1(p2(2),p1(1)) c area = .5d0*(p2(2)-p2(1))*(p2(2)+p2(1)-2.d0*pwave(2)) c if (idebug .eq. 1) then write(6,1002) lstgrd,ix0,area,ar(1,lstgrd) write(6,1003) (q(ix0,m),m=1,meqn) endif c frac = area/ar(1,lstgrd) do m=1,meqn q(ix0,m) = q(ix0,m)+frac*work(m) enddo c if (idebug .eq. 1) write(6,1003) (q(ix0,m),m=1,meqn) else c # irregular cell c # find the intersection area with each subcell do icell=1,ncells(iptr) p1(1) = poly(icell,1,iptr) p1(2) = poly(icell,2,iptr) c p2(1) = dmax1(pwave(1),p1(1)) p2(1) = dmin1(p2(1),p1(2)) p2(2) = dmin1(pwave(3),p1(2)) p2(2) = dmax1(p2(2),p1(1)) c area = .5d0*(p2(2)-p2(1))*(p2(2)+p2(1)-2.d0*pwave(2)) c if (idebug .eq. 1) then write(6,1002) iptr,icell,area,ar(icell,iptr) write(6,1003) (qirr(icell,m,iptr),m=1,meqn) endif c frac = area/ar(icell,iptr) do m=1,meqn qirr(icell,m,iptr) = qirr(icell,m,iptr)+ & frac*work(m) enddo c if (idebug .eq. 1) write(6,1003) & (qirr(icell,m,iptr),m=1,meqn) enddo endif c 1001 format('entering irrwv2: sigma=', 5f10.5) 1002 format('iptr=',i3,', icell=',i3,', area=',f12.6,', artot=',f12.6) 1003 format('q=',5f12.6) c return end c c ------------------------------------------------------------------ c subroutine mkirrc(irr,ix,xp,mitot,lstgrd) implicit double precision (a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension irr(mitot),p1(3) data idebug/0/ c c # make irregular subcells c iptr = irr(ix) if (iabs(iptr) .eq. lstgrd) then c # regular cell case c if (iptr .eq. -lstgrd) return c icell = 1 jcell = 1 p1(1) = xgrid(ix) p1(2) = xgrid(ix+1) else c # irregular cell case icell = ncells(iptr) c c # locate the splitted subcell index do jcell=1,ncells(iptr) if ((xp-poly(jcell,1,iptr))* & (xp-poly(jcell,2,iptr)) .le. 0.d0) go to 10 enddo c 10 continue c # split this subcell into two parts do ii=1,3 p1(ii) = poly(jcell,ii,iptr) enddo endif c area1 = xp-p1(1) area2 = p1(2)-xp if ((area1 .le. 1.d-12) .or. & (area2 .le. 1.d-12)) return c c # allocate irregular cell pointer if (icell .eq. 1) iptr = maptr(1) c c # create subcell polygon poly(jcell,1,iptr) = p1(1) poly(jcell,2,iptr) = xp ar(jcell,iptr) = area1 xcen(jcell,iptr) = .5d0*(p1(1)+xp) c poly(icell+1,1,iptr) = xp poly(icell+1,2,iptr) = p1(2) ar(icell+1,iptr) = area2 xcen(icell+1,iptr) = .5d0*(xp+p1(2)) c c # make bond newbd = maptr(6) ixbond(newbd) = ix points(newbd) = xp c c # set other associated data sturcture irr(ix) = iptr ixirr(iptr) = ix ncells(iptr) = ncells(iptr)+1 nbonds(iptr) = nbonds(iptr)+1 nbondc(nbonds(iptr),iptr) = newbd c c # set bond-cell correspondence call setbdc(iptr) c if (idebug .eq. 1) then write(6,*) 'entering mkirrc:' write(6,1004) iptr,ix,newbd,ityirr(iptr) do i0=1,2 is = irrcell(i0,newbd) write(6,1003) is,ar(is,iptr) write(6,1002) (poly(is,ii,iptr),ii=1,2) enddo write(6,*) '---------------------------------------------' endif c 1002 format('poly=',2f12.6) 1003 format(i3,x,f14.10) 1004 format('iptr=',i4,', ix=',i3,', kbond=',i4, & ', ityirr=',i3) return end c c ------------------------------------------------------------------------ c subroutine setbdc(iptr) implicit double precision (a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension sign(6),ineg(6),ipos(6) c c # set bond-cells correspondence c if (nbonds(iptr) .ne. 1) go to 10 c c # single bond case kbond = nbondc(1,iptr) irrcell(1,kbond) = 1 irrcell(2,kbond) = 2 return c 10 continue c # multiple bond case do 50 ibd=1,nbonds(iptr) kbond = nbondc(ibd,iptr) xp = points(kbond) ip0 = 0 ip1 = 0 do icell=1,ncells(iptr) sign(icell) = dsign(1.d0,xcen(icell,iptr)-xp) if (sign(icell) .lt. 0.d0) then c # left-hand side subcell ip0 = ip0+1 ineg(ip0) = icell else c # right-hand side subcell ip1 = ip1+1 ipos(ip1) = icell endif enddo c c # set left-hand side bond-cell info do i=1,ip0 p1 = poly(ineg(i),1,iptr) p2 = poly(ineg(i),2,iptr) if ((dabs(xp-p1) .le. 1.d-12) .or. & (dabs(xp-p2) .le. 1.d-12)) go to 39 enddo go to 90 39 irrcell(1,kbond) = ineg(i) c c # set right-hand side bond-cell info do i=1,ip1 p1 = poly(ipos(i),1,iptr) p2 = poly(ipos(i),2,iptr) if ((dabs(xp-p1) .le. 1.d-12) .or. & (dabs(xp-p2) .le. 1.d-12)) go to 49 enddo go to 90 49 irrcell(2,kbond) = ipos(i) 50 continue return c 90 continue write(6,1004) iptr,ixirr(iptr),ityirr(iptr) do ibd=1,nbonds(iptr) kbond = nbondc(ibd,iptr) write(6,1001) kbond,points(kbond),irrcell(1,kbond), & irrcell(2,kbond),itybond(kbond) enddo c do 96 icell=1,ncells(iptr) write(6,1003) icell,ar(icell,iptr),xcen(icell,iptr), & sign(icell) write(6,1002) (poly(icell,ii,iptr),ii=1,2) 96 continue c 1001 format('kbond=',i3,', point=',f12.6, & ', icl=',i3,', icr=',i3,', itype=',i3) 1002 format('poly=',2f10.6) 1003 format('icell=',i3,', area=',e12.6,', xcen=',e12.6, & ', sign=',e12.6) 1004 format('ERROR in setbdc: iptr=',i3,', ix=',i3, & ', itype=',i3) stop c end c c -------------------------------------------------------------------- c integer function maptr(itype) implicit double precision (a-h,o-z) parameter(maxfr=10,maxbd=126) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) c c # get pointer for a irregular cell or a tracked front c nextk = lstptr(itype) c if(itype .eq. 1) then c # irregular cell k0 = iprirr(nextk) nxtirr(k0) = nextk k = nextk nxtirr(k) = 0 iprirr(k) = k0 ityirr(k) = 1 ncells(k) = 1 nbonds(k) = 0 elseif(itype .eq. 2) then c # tracked front k0 = iprfront(nextk) nxtfront(k0) = nextk k = nextk nxtfront(k) = 0 iprfront(k) = k0 ityfront(k) = 1 newfront(k) = 1 elseif(itype .eq. 6) then c # tracked bond k0 = iprbond(nextk) nxtbond(k0) = nextk k = nextk nxtbond(k) = 0 iprbond(k) = k0 itybond(k) = 0 else write(6,*) 'no pointer available for itype=',itype stop endif c maptr = k c return end c c ------------------------------------------------------------------- c subroutine setdoc(alloc,memsize) implicit double precision (a-h,o-z) parameter(maxfr=10,maxbd=126) common /space/ lfree(50,2),lhead(6),lenf,idimf,lentot,lenmax common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkfront/ iprfront(maxfr),nxtfront(maxfr), & ityfront(maxfr),newfront(maxfr), & mwfront(maxfr),lstbond(maxfr) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension alloc(memsize) c c # initialize main storage do i=1,memsize alloc(i) = 0.0d0 enddo c c # set free space array lenf = 3 idimf = 50 do i=1,idimf lfree(i,1) = 0 lfree(i,2) = 0 enddo lfree(3,1) = memsize+2 lfree(2,1) = 1 lfree(2,2) = memsize c c # set doubly linked list for irregular cells do i=1,maxbd-1 iprirr(i) = i-1 nxtirr(i) = i+1 enddo nxtirr(maxbd) = 0 lhead(1) = 1 c c # set doubly linked list for tracked fronts do i=1,maxfr-1 iprfront(i) = i-1 nxtfront(i) = i+1 enddo nxtfront(maxfr) = 0 lhead(2) = 1 c c # set doubly linked list for tracked bonds do i=1,maxbd-1 iprbond(i) = i-1 nxtbond(i) = i+1 enddo nxtbond(maxbd) = 0 lhead(6) = 1 return c end c c ----------------------------------------------------------------- c double precision function fpgeom(x) implicit double precision (a-h,o-z) common /sterms/ nsource,isource(10),igeom c go to (1,2,3,4,5) igeom c 1 continue fpgeom = 8.d0*.347d0/dcosh(8.d0*x-4.d0)**2 return c 2 continue fpgeom = -8.d0*.347d0/dcosh(-8.d0*x+4.d0)**2 return c 3 continue if(x .gt. 0.d0) then fpgeom = 8.d0*.347d0/dcosh(8.d0*x-4.d0)**2 else x = x+1.d0 fpgeom = -8.d0*.347d0/dcosh(-8.d0*x+4.d0)**2 endif return c 4 continue fpgeom = 1.d0 return c 5 continue fpgeom = 2.d0*x return c end c c ----------------------------------------------------------------------- c subroutine setqir(qold,irrold,mitot,meqn,hx,work,lstgrd, & lstold,irric) implicit double precision (a-h,o-z) parameter(maxbd=126) common /fluid/ gamma(5),zf(5),qchem(5),gmole(5),mphase common /frontin/ xloop(10),varloop(2,6,10), & ibdloop(10),mwloop(10),istloop(10), & newloop(10),irrt0(10),irrtn,nloops common /init0/ rhol,ul,pl,rhor,ur,pr common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) common /local/ iloops dimension qold(mitot,meqn),irrold(mitot),work(meqn) dimension po1(3),po2(3),p2(3) data idebug/1/ c c # set values for irregular cells c go to (1,2,3) irric c 1 continue c # set values based on the input states iptr = lstgrd do ik=1,maxbd iptr = nxtirr(iptr) if (iptr .eq. 0) go to 109 if (ityirr(iptr) .ne. 2) then kbond = nbondc(1,iptr) do m=1,meqn qirr(irrcell(1,kbond),m,iptr) = varloop(1,m,iloops) qirr(irrcell(2,kbond),m,iptr) = varloop(2,m,iloops) enddo endif enddo 109 continue go to 990 c 2 continue c # set values based on the old grid states iptr = lstgrd do ik=1,maxbd iptr = nxtirr(iptr) if(iptr .eq. 0) go to 219 c ix = ixirr(iptr) iptr1 = irrold(ix) if (iptr1 .ne. lstold) then c # irregular cell case do icell=1,ncells(iptr) do m=1,meqn qirr(icell,m,iptr) = 0.d0 enddo c po1(1) = poly(icell,1,iptr) po1(2) = poly(icell,2,iptr) c c # set values based on area weighted cell averages do icold=1,ncells(iptr1) po2(1) = poly(icold,1,iptr1) po2(2) = poly(icold,2,iptr1) c p2(1) = dmax1(po1(1),po2(1)) p2(1) = dmin1(p2(1),po2(2)) p2(2) = dmin1(po1(2),po2(2)) p2(2) = dmax1(p2(2),po2(1)) c area = p2(2)-p2(1) frac = area/ar(icell,iptr) do m=1,meqn qirr(icell,m,iptr) = qirr(icell,m,iptr)+ & frac*qirr(icold,m,iptr1) enddo enddo enddo else c # regular cell case do icell=1,ncells(iptr) do m=1,meqn qirr(icell,m,iptr) = qold(ix,m) enddo enddo endif enddo 219 continue go to 990 c 3 continue c # set values for radially symmetric flow rs = xloop(iloops) c iptr = lstgrd do ik=1,maxbd iptr = nxtirr(iptr) if (iptr .eq. 0) go to 309 c do m=1,meqn qirr(1,m,iptr) = varloop(1,m,iloops) qirr(2,m,iptr) = varloop(2,m,iloops) enddo c r = xcen(1,iptr) u = (r/rs)**2*ul call prmtoc(work,rhol,u,pl,gamma(1),zf(1)) do m=1,meqn qirr(1,m,iptr) = work(m) enddo enddo 309 continue go to 990 c 990 continue if (idebug .eq. 1) then write(6,*) 'entering setqir' iptr = lstgrd do ik=1,maxbd iptr = nxtirr(iptr) if (iptr .eq. 0) go to 199 c write(6,1001) iptr,ixirr(iptr),ityirr(iptr) c do ibond=1,nbonds(iptr) kbond = nbondc(ibond,iptr) write(6,1002) kbond,points(kbond),itybond(kbond) enddo c do icell=1,ncells(iptr) write(6,1003) ar(icell,iptr), & (qirr(icell,m,iptr),m=1,meqn) enddo enddo 199 continue write(6,*) '------------------------------------------' endif c 1001 format('iptr=',i4,', ix=',i3,', type=',i3) 1002 format('kbond=',i4,', point=',f10.5,', type=',i3) 1003 format(e12.6,2x,6f10.5) c return end c c ----------------------------------------------------------------- c subroutine irrslp(work,hlen,irrold,s,qjump,mitot,meqn, & mwave,ix1,ix2,mw,kbond,lstold,static) implicit double precision(a-h,o-z) parameter(n1d=2100,maxbd=126) common /coord/ xgrid(0:n1d) common /frpsoln/ sfront(maxbd,3),qjfront(maxbd,6,3), & mwtrack(maxbd,3) common /hypmeth/ mthslp(3),mthlim(3) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension irrold(mitot),qjump(mitot,meqn,mwave) dimension work(meqn),qja(6),qjb(6) logical mwtrack,static data idebug /0/ c c # compute slopes for irregular cells c # do one-sided difference for the strong wave point c if (idebug .eq. 1) write(6,1001) mw,s c if (.not. static) go to 200 c c # fixed cell boundary case iptr1 = irrold(ix1) if (iabs(iptr1) .eq. lstold) then hs1 = ar(1,lstold) else icell = irrcls(xgrid(ix2),iptr1) hs1 = ar(icell,iptr1) endif c iptr2 = irrold(ix2) if (iabs(iptr2) .eq. lstold) then hs2 = ar(1,lstold) else icell = irrcls(xgrid(ix2),iptr2) hs2 = ar(icell,iptr2) endif c ha = .5d0*(hs1+hs2) do m=1,meqn qja(m) = qjump(ix2,m,mw) enddo c if (s .gt. 0.d0) then hlen = -hs1 if (iabs(iptr1) .eq. lstold) then c # regular cell case iptr11 = irrold(ix1-1) if (iabs(iptr11) .eq. lstold) then hs2 = ar(1,lstold) else icell = irrcls(xgrid(ix1),iptr11) hs2 = ar(icell,iptr11) endif c hb = .5d0*(hs1+hs2) do m=1,meqn qjb(m) = qjump(ix1,m,mw) enddo else c # irregular cell case kbond = nbondc(kbcls(xgrid(ix2),iptr1),iptr1) hb = .5d0*(ar(irrcell(1,kbond),iptr1)+ & ar(irrcell(2,kbond),iptr1)) do m=1,meqn qjb(m) = qjfront(kbond,m,mw) enddo endif else hlen = hs2 if (iabs(iptr2) .eq. lstold) then c # regular cell case iptr21 = irrold(ix2+1) if (iabs(iptr21) .eq. lstold) then hs1 = ar(1,lstold) else icell = irrcls(xgrid(ix1),iptr21) hs1 = ar(icell,iptr21) endif c hb = .5d0*(hs1+hs2) do m=1,meqn qjb(m) = qjump(ix2+1,m,mw) enddo else c # irregular cell case kbond = nbondc(kbcls(xgrid(ix2),iptr2),iptr2) hb = .5d0*(ar(irrcell(1,kbond),iptr2)+ & ar(irrcell(2,kbond),iptr2)) do m=1,meqn qjb(m) = qjfront(kbond,m,mw) enddo endif endif go to 300 c 200 continue c # tracked cell boundary case iptr = irrold(ix1) hs1 = ar(irrcell(1,kbond),iptr) hs2 = ar(irrcell(2,kbond),iptr) c ha = .5d0*(hs1+hs2) do m=1,meqn qja(m) = qjfront(kbond,m,mw) enddo c if (s .gt. 0.d0) then hlen = -hs1 iptr1 = irrold(ix1-1) if (iabs(iptr1) .eq. lstold) then c # regular cell case hs2 = ar(1,lstold) else icell = irrcls(xgrid(ix1),iptr1) hs2 = ar(icell,iptr1) endif c hb = .5d0*(hs1+hs2) do m=1,meqn qjb(m) = qjump(ix1,m,mw) enddo else hlen = hs2 iptr2 = irrold(ix1+1) if (iabs(iptr2) .eq. lstold) then c # regular cell case hs1 = ar(1,lstold) else icell = irrcls(xgrid(ix1+1),iptr2) hs1 = ar(icell,iptr2) endif c hb = .5d0*(hs1+hs2) do m=1,meqn qjb(m) = qjump(ix1+1,m,mw) enddo endif c 300 continue c # compute slopes call limter(work,qja,qjb,ha,hb,meqn,mthslp(mw), & mthlim(mw)) 1001 format('entering irrslp: mw=',i3,', s=',f12.6) return end c c ------------------------------------------------------------------- c integer function irrcls(x,iptr) implicit double precision (a-h,o-z) parameter(maxbd=126) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension dar(6) c c # determine the closest subcell to point x c do ic=1,ncells(iptr) dar(ic) = x-xcen(ic,iptr) enddo c irrcls = 1 do ic=2,ncells(iptr) if (dabs(dar(ic)) .lt. dabs(dar(irrcls))) irrcls = ic enddo c return end c c ------------------------------------------------------------------- c integer function kbcls(x,iptr) implicit double precision (a-h,o-z) parameter(maxbd=126) common /lkbond/ iprbond(maxbd),nxtbond(maxbd),itybond(maxbd), & irrcell(2,maxbd),ixbond(maxbd),points(maxbd) common /lkirr/ iprirr(maxbd),nxtirr(maxbd),ityirr(maxbd), & ixirr(maxbd),ncells(maxbd),nbonds(maxbd), & nbondc(6,maxbd),ar(6,maxbd),xcen(6,maxbd), & poly(6,10,maxbd),qirr(6,6,maxbd) dimension dar(6) c c # determine the closest "bond" to point x c do ibd=1,nbonds(iptr) kbond = nbondc(ibd,iptr) dar(ibd) = x-points(kbond) enddo c kbcls = 1 do ibd=2,nbonds(iptr) if (dabs(dar(ibd)) .lt. dabs(dar(kbcls))) kbcls = ibd enddo c return end c c ----------------------------------------------------------------- c subroutine fsurgy(dt,rwork,mitot) implicit double precision (a-h,o-z) double precision kjdt,newpt parameter(maxbd=126) common /trackin/ qjtol(3),itrack,idetec(3),isurgy common /stags/ oldpt(maxbd),sppt(maxbd),mwpt(maxbd),ntags dimension rwork(mitot) dimension kjdt(maxbd),newpt(maxbd) c c # determine the appropriate time step to avoid c # strong wave interactions at the current time step c c # kjtot : total number of interaction points c # kjdt(i): the interaction time step of the ith interaction c kj = 0 c go to (1,2,3) isurgy+1 c 1 continue go to 900 c 2 continue c # handle strong wave interaction do k=1,ntags rwork(k) = 0.d0 newpt(k) = oldpt(k)+sppt(k)*dt enddo c do k=1,ntags rwork(k) = 1.d0 c do 20 k0=1,ntags if (rwork(k0) .eq. 1.d0) go to 20 if ((newpt(k)-newpt(k0))*(oldpt(k)-oldpt(k0)) .lt. & 0.d0) then c # wavet interaction kj = kj+1 kjdt(kj) = (oldpt(k0)-oldpt(k))/(sppt(k)-sppt(k0)) c write(6,*) 'interaction happened' write(6,*) 'x0=',oldpt(k), ', x1=',newpt(k) write(6,*) 'x0=',oldpt(k0),', x1=',newpt(k0) write(6,*) 'kj=',kj,', kjdt=',kjdt(kj) endif 20 continue enddo c kjtot = kj c c # choose the minimum time step if (kjtot .ne. 0) then do kj=1,kjtot dt = dmin1(dt,kjdt(kj)) enddo write(6,*) 'new dt = ', dt endif go to 900 c 3 continue go to 900 c 900 continue return end