c c c =================================================================== subroutine step1(maxmx,mvar,meqn,maux,mwaves,mbc,mx, & qold,aux,dx,t,dt,method,mthlim,cflgrid, & fm,fp,wave,s, & amdq,apdq,dtdx,aux1,q1,work,mwork,rp1) c =================================================================== c c Author: Randall J. LeVeque c Modified for AMROC: Ralf Deiterding c c # Main entry from AMROC. c # Take one time step, updating qold. c c # fm, fp are fluxes to left and right of single cell edge c c method(2) = 1 ==> Godunov method c method(2) = 2 ==> Wave limiter method c mthlim(p) controls what limiter is used in the pth family c method(2) = 3 ==> Slope-limiting for the conserved quantities c method(2) > 3 ==> User defined slope-limiter method c c amdq, apdq, wave, s, and f are used locally: c c amdq(1-mbc:maxmx+mbc, meqn) = left-going flux-differences c apdq(1-mbc:maxmx+mbc, meqn) = right-going flux-differences c e.g. amdq(i,m) = m'th component of A^- \Delta q from i'th Riemann c problem (between cells i-1 and i). c c wave(1-mbc:maxmx+mbc, meqn, mwaves) = waves from solution of c Riemann problems, c wave(i,m,mw) = mth component of jump in q across c wave in family mw in Riemann problem between c states i-1 and i. c c s(1-mbc:maxmx+mbc, mwaves) = wave speeds, c s(i,mw) = speed of wave in family mw in Riemann problem between c states i-1 and i. c c f(1-mbc:maxmx+mbc, meqn) = correction fluxes for second order method c f(i,m) = mth component of flux at left edge of ith cell c -------------------------------------------------------------------- c implicit double precision (a-h,o-z) include "call.i" c dimension qold(mvar, 1-mbc:maxmx+mbc) dimension fm(mvar, 1-mbc:maxmx+mbc) dimension fp(mvar, 1-mbc:maxmx+mbc) dimension wave(1-mbc:maxmx+mbc, meqn, mwaves) dimension s(1-mbc:maxmx+mbc, mwaves) dimension amdq(1-mbc:maxmx+mbc, meqn) dimension apdq(1-mbc:maxmx+mbc, meqn) dimension aux(maux, 1-mbc:maxmx+mbc) dimension dtdx(1-mbc:maxmx+mbc) dimension aux1(1-mbc:maxmx+mbc, maux) dimension q1(1-mbc:maxmx+mbc, meqn) dimension method(7),mthlim(mwaves) dimension work(mwork) external rp1 common /rpnflx/ mrpnflx c i0ql = 1 i0qr = i0ql + (maxmx+2*mbc)*meqn i0fl = i0qr + (maxmx+2*mbc)*meqn i0fr = i0fl + (maxmx+2*mbc)*meqn if (method(1).ge.1) then i0qls = i0fr + (maxmx+2*mbc)*meqn i0qrs = i0qls + (maxmx+2*mbc)*mvar iused = i0qrs + (maxmx+2*mbc)*mvar - 1 else iused = i0fr + (maxmx+2*mbc)*meqn - 1 endif c if (iused.gt.mwork) then write(6,*) '*** not enough work space in step1' write(6,*) '*** iused = ', iused, ' mwork =',mwork stop endif c tcom = t dtcom = dt dxcom = dx c mcapa = method(6) c c # check if any limiters are used: limit = 0 do 5 mw=1,mwaves if (mthlim(mw) .gt. 0) limit = 1 5 continue c do 10 i=1-mbc,mx+mbc if (mcapa.gt.0) then if (aux(mcapa,i) .le. 0.d0) then write(6,*) 'Error -- capa must be positive' stop endif dtdx(i) = dt / (dx*aux(mcapa,i)) else dtdx(i) = dt/dx endif 10 continue do 20 i = 1-mbc, mx+mbc do 20 m=1,meqn q1(i,m) = qold(m,i) 20 continue if (maux .gt. 0) then do 22 i = 1-mbc, mx+mbc do 22 ma=1,maux aux1(i,ma) = aux(ma,i) 22 continue endif do 23 i = 1-mbc, mx+mbc do 23 m=1,mvar fm(m,i) = 0.d0 fp(m,i) = 0.d0 23 continue c c c c # solve Riemann problem at each interface c ----------------------------------------- c if (method(2).le.2) then call rp1(maxmx,meqn,mwaves,mbc,mx,q1,q1,maux,aux1,aux1, & wave,s,amdq,apdq) else call slope1(maxmx,meqn,maux,mwaves,mbc,mx,q1,aux1,dx,dt, & method,mthlim,wave,s,amdq,apdq,dtdx,rp1, & work(i0ql),work(i0qr),work(i0fl),work(i0fr)) if (method(1).ge.1) & call saverec1(maxmx,mvar,meqn,mbc,mx,qold,work(i0qls), & work(i0qrs),work(i0ql),work(i0qr)) endif c do 40 i=1,mx+1 do 40 m=1,meqn fp(m,i) = fp(m,i) - apdq(i,m) fm(m,i) = fm(m,i) + amdq(i,m) 40 continue c c # compute maximum wave speed: cflgrid = 0.d0 do 50 mw=1,mwaves do 50 i=1,mx+1 cflgrid = dmax1(cflgrid, dtdx(i)*dabs(s(i,mw))) 50 continue c if (method(2).eq.1.or.method(2).ge.3) go to 900 c if (mrpnflx.ne.0) then write(6,*) '*** Riemann solver returns fluxes.' write(6,*) '*** Wave limiting not possible.' write(6,*) '*** Set method(2)>=3 for slope limiting.' stop endif c c # compute correction fluxes for second order q_{xx} terms: c ---------------------------------------------------------- c c # apply limiter to waves: if (limit.eq.1) call limiter(maxmx,meqn,mwaves,mbc,mx, & wave,s,mthlim) c do 120 i=1,mx+1 dtdxave = 0.5d0 * (dtdx(i-1) + dtdx(i)) do 120 m=1,meqn cqxx = 0.d0 do 110 mw=1,mwaves c # second order corrections: cqxx = cqxx + dabs(s(i,mw)) & * (1.d0 - dabs(s(i,mw))*dtdxave) * wave(i,m,mw) 110 continue fm(m,i) = fm(m,i) + 0.5d0 * cqxx fp(m,i) = fp(m,i) + 0.5d0 * cqxx c 120 continue c 900 continue c if (method(1).ge.1) then if (method(2).le.2) then call fmod1(maxmx,mvar,meqn,maux,mbc,mx,qold,qold,qold, & aux,dx,dt,method,cflgrid,fm,fp,q1,aux1, & amdq,apdq) else call fmod1(maxmx,mvar,meqn,maux,mbc,mx,qold,work(i0qls), & work(i0qrs),aux,dx,dt,method,cflgrid,fm,fp, & q1,aux1,amdq,apdq) endif endif c c # update q c ============================================ c do 150 i=1,mx do 150 m=1,meqn qold(m,i) = qold(m,i) - dtdx(i) * (fm(m,i+1) - fp(m,i)) 150 continue c c ============================================ c return end c c c =================================================================== subroutine slope1(maxm,meqn,maux,mwaves,mbc,mx, & q,aux,dx,dt,method,mthlim, & wave,s,amdq,apdq,dtdx,rp1, & ql,qr,fl,fr) c =================================================================== c c # Implements the standard MUSCL-Hancock method, which gives 2nd c # order accuracy for conventional finite-volume schemes that c # return the numerical fluxes instead of fluctuations. c c # Author: Ralf Deiterding, ralf@cacr.caltech.edu c implicit double precision (a-h,o-z) common /rpnflx/ mrpnflx c dimension q(1-mbc:maxm+mbc, meqn) dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) dimension aux(1-mbc:maxm+mbc, maux) dimension fl(1-mbc:maxm+mbc, meqn) dimension fr(1-mbc:maxm+mbc, meqn) dimension wave(1-mbc:maxm+mbc, meqn, mwaves) dimension s(1-mbc:maxm+mbc, mwaves) dimension amdq(1-mbc:maxm+mbc, meqn) dimension apdq(1-mbc:maxm+mbc, meqn) dimension dtdx(1-mbc:maxm+mbc) dimension method(7),mthlim(mwaves) external rp1 c mlim = 0 do 90 mw=1,mwaves if (mthlim(mw) .gt. 0) then mlim = mthlim(mw) goto 95 endif 90 continue 95 continue c do 100 m=1,meqn ql(1-mbc,m) = q(1-mbc,m) qr(1-mbc,m) = q(1-mbc,m) ql(mx+mbc,m) = q(mx+mbc,m) qr(mx+mbc,m) = q(mx+mbc,m) 100 continue c if (method(2).gt.3) then call rec1(maxm,meqn,mwaves,mbc,mx,q,method,mthlim,ql,qr) c c # MUSCL reconstruction with slope-limiting for conserved c # variables. Linear reconstruction: om=0.d0 c # Quadratic spatial reconstuction: om!=0.d0 c # 2nd order spatial reconstruction: om=1.d0/3.d0 c # Reconstructions with om!=0.d0 are not conservative! c else om = 0.d0 do 110 i=2-mbc,mx+mbc-1 do 110 m=1,meqn call reclim(q(i,m),q(i-1,m),q(i+1,m), & mlim,om,ql(i,m),qr(i,m)) 110 continue endif c call flx1(maxm,meqn,mbc,mx,ql,maux,aux,fl) call flx1(maxm,meqn,mbc,mx,qr,maux,aux,fr) c do 200 i=2-mbc,mx+mbc-1 do 200 m=1,meqn ql(i,m) = ql(i,m) + 0.5d0*dtdx(i)*(fl(i,m)-fr(i,m)) qr(i,m) = qr(i,m) + 0.5d0*dtdx(i)*(fl(i,m)-fr(i,m)) 200 continue c call rp1(maxm,meqn,mwaves,mbc,mx,ql,qr,maux,aux,aux, & wave,s,amdq,apdq) c c # Add differences of fluxes between original and reconstructed values c # to fluctuations, if a wave propagation scheme is applied. c if (mrpnflx.eq.0) then call flx1(maxm,meqn,mbc,mx,ql,maux,aux,fl) call flx1(maxm,meqn,mbc,mx,qr,maux,aux,fr) do 300 i=2-mbc,mx+mbc-1 do 300 m=1,meqn amdq(i,m) = amdq(i,m) + fr(i-1,m) apdq(i,m) = apdq(i,m) - fl(i ,m) 300 continue call flx1(maxm,meqn,mbc,mx,q,maux,aux,fl) do 310 i=2-mbc,mx+mbc-1 do 310 m=1,meqn amdq(i,m) = amdq(i,m) - fl(i-1,m) apdq(i,m) = apdq(i,m) + fl(i ,m) 310 continue endif c return end c c c c =================================================================== subroutine saverec1(maxm,mvar,meqn,mbc,mx,q,qls,qrs,ql,qr) c =================================================================== c c # Store reconstructed values for later use. c implicit double precision (a-h,o-z) dimension q(mvar, 1-mbc:maxm+mbc) dimension qls(mvar, 1-mbc:maxm+mbc) dimension qrs(mvar, 1-mbc:maxm+mbc) dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) c do 10 i = 1-mbc, mx+mbc do m=1,meqn qls(m,i) = ql(i,m) qrs(m,i) = qr(i,m) enddo do m=meqn+1,mvar qls(m,i) = q(m,i) qrs(m,i) = q(m,i) enddo 10 continue c return end c