Blockstructured Adaptive Mesh Refinement in object-oriented C++
c c c ===================================================== subroutine flux2(ixy,maxm,meqn,maux,mwaves,mbc,mx, & q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, & faddm,faddp,gaddm,gaddp,cfl1d,wave,s, & amdq,apdq,cqxx,bmasdq,bpasdq,work,mwork, & rpn2,rpt2) c ===================================================== c c Author: Randall J. LeVeque c Modified for AMROC: Ralf Deiterding c c # Compute the modification to fluxes f and g that are generated by c # all interfaces along a 1D slice of the 2D grid. c # ixy = 1 if it is a slice in x c # 2 if it is a slice in y c # This value is passed into the Riemann solvers. The flux modifications c # go into the arrays fadd and gadd. The notation is written assuming c # we are solving along a 1D slice in the x-direction. c c # fadd(i,.) modifies F to the left of cell i c # gadd(i,.,1) modifies G below cell i c # gadd(i,.,2) modifies G above cell i c c # The method used is specified by method(2:3): c c method(2) = 1 if only first order increment waves are to be used. c = 2 if second order correction terms are to be added, with c a flux-limiter as specified by mthlim. c = 3 ==> Slope-limiting of the conserved variables, with a c slope-limiter as specified by mthlim(1). c > 3 ==> User defined slope-limiter method. c Slope-limiting is intended to be used with dimensional c splitting! c c method(3) = -1 0 Gives dimensional splitting using Godunov c splitting, i.e. formally first order accurate. c On a single computational nodes =-1 and =-2 c will give identical results, but =-1 will be faster. c = -2 Dimensional splitting using Godunov splitting with c boundary update after each directional step. c The necessary ghost cell synchronization is done by c the surrounding AMROC framework. c This selection ensures that the solution of the c splitting method is independent of the number of c computational nodes. c = 0 Gives the Donor cell method. No transverse c propagation of neither the increment wave c nor the correction wave. c = 1 if transverse propagation of increment waves c (but not correction waves, if any) is to be applied. c = 2 if transverse propagation of correction waves is also c to be included. c c Note that if mcapa>0 then the capa array comes into the second c order correction terms, and is already included in dtdx1d: c If ixy = 1 then c dtdx1d(i) = dt/dx if mcapa= 0 c = dt/(dx*aux(i,jcom,mcapa)) if mcapa = 1 c If ixy = 2 then c dtdx1d(j) = dt/dy if mcapa = 0 c = dt/(dy*aux(icom,j,mcapa)) if mcapa = 1 c c Notation: c The jump in q (q1d(i,:)-q1d(i-1,:)) is split by rpn2 into c amdq = the left-going flux difference A^- Delta q c apdq = the right-going flux difference A^+ Delta q c Each of these is split by rpt2 into c bmasdq = the down-going transverse flux difference B^- A^* Delta q c bpasdq = the up-going transverse flux difference B^+ A^* Delta q c where A^* represents either A^- or A^+. c c implicit double precision (a-h,o-z) include "call.i" c dimension q1d(1-mbc:maxm+mbc, meqn) dimension amdq(1-mbc:maxm+mbc, meqn) dimension apdq(1-mbc:maxm+mbc, meqn) dimension bmasdq(1-mbc:maxm+mbc, meqn) dimension bpasdq(1-mbc:maxm+mbc, meqn) dimension cqxx(1-mbc:maxm+mbc, meqn) dimension faddm(1-mbc:maxm+mbc, meqn) dimension faddp(1-mbc:maxm+mbc, meqn) dimension gaddm(1-mbc:maxm+mbc, meqn, 2) dimension gaddp(1-mbc:maxm+mbc, meqn, 2) dimension dtdx1d(1-mbc:maxm+mbc) dimension aux1(1-mbc:maxm+mbc, maux) dimension aux2(1-mbc:maxm+mbc, maux) dimension aux3(1-mbc:maxm+mbc, maux) dimension s(1-mbc:maxm+mbc, mwaves) dimension wave(1-mbc:maxm+mbc, meqn, mwaves) dimension method(7),mthlim(mwaves) dimension work(mwork) external rpn2, rpt2 logical limit c i0ql = 1 i0qr = i0ql + (maxm+2*mbc)*meqn i0fl = i0qr + (maxm+2*mbc)*meqn i0fr = i0fl + (maxm+2*mbc)*meqn iused = i0fr + (maxm+2*mbc)*meqn - 1 c if (iused.gt.mwork) then write(6,*) '*** not enough work space in flux2ex' write(6,*) '*** iused = ', iused, ' mwork =',mwork stop endif c limit = .false. do 5 mw=1,mwaves if (mthlim(mw) .gt. 0) limit = .true. 5 continue c c # initialize flux increments: c ----------------------------- c do 30 jside=1,2 do 20 m=1,meqn do 10 i = 1-mbc, mx+mbc faddm(i,m) = 0.d0 faddp(i,m) = 0.d0 gaddm(i,m,jside) = 0.d0 gaddp(i,m,jside) = 0.d0 10 continue 20 continue 30 continue c c c # solve Riemann problem at each interface and compute Godunov updates c --------------------------------------------------------------------- c if (method(2).le.2) then call rpn2(ixy,maxm,meqn,mwaves,mbc,mx,q1d,q1d,maux,aux2,aux2, & wave,s,amdq,apdq) else call slope2(ixy,maxm,meqn,maux,mwaves,mbc,mx,q1d,aux2,dx,dt, & method,mthlim,wave,s,amdq,apdq,dtdx1d,rpn2, & work(i0ql),work(i0qr),work(i0fl),work(i0fr)) endif c c # Set fadd for the donor-cell upwind method (Godunov) do 40 m=1,meqn do 40 i=1,mx+1 faddp(i,m) = faddp(i,m) - apdq(i,m) faddm(i,m) = faddm(i,m) + amdq(i,m) 40 continue c c # compute maximum wave speed for checking Courant number: cfl1d = 0.d0 do 50 mw=1,mwaves do 50 i=1,mx+1 cfl1d = dmax1(cfl1d, dtdx1d(i)*dabs(s(i,mw))) 50 continue c if (method(2).le.1.or.method(2).ge.3) go to 130 c c # modify F fluxes for second order q_{xx} correction terms: c ----------------------------------------------------------- c c # apply limiter to waves: if (limit) call limiter(maxm,meqn,mwaves,mbc,mx,wave,s,mthlim) c do 120 i = 1, mx+1 c c # For correction terms below, need average of dtdx in cell c # i-1 and i. Compute these and overwrite dtdx1d: c dtdx1d(i-1) = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i)) c do 120 m=1,meqn cqxx(i,m) = 0.d0 do 119 mw=1,mwaves c c # second order corrections: cqxx(i,m) = cqxx(i,m) + dabs(s(i,mw)) & * (1.d0 - dabs(s(i,mw))*dtdx1d(i-1)) * wave(i,m,mw) c 119 continue faddm(i,m) = faddm(i,m) + 0.5d0 * cqxx(i,m) faddp(i,m) = faddp(i,m) + 0.5d0 * cqxx(i,m) 120 continue c c 130 continue c if (method(3).le.0) go to 999 !# no transverse propagation c if (method(3).eq.2) then c # incorporate cqxx into amdq and apdq so that it is split also. do 150 m=1,meqn do 150 i = 1, mx+1 amdq(i,m) = amdq(i,m) + cqxx(i,m) apdq(i,m) = apdq(i,m) - cqxx(i,m) 150 continue endif c c c # modify G fluxes for transverse propagation c -------------------------------------------- c c c # split the left-going flux difference into down-going and up-going: call rpt2(ixy,maxm,meqn,mwaves,mbc,mx, & q1d,q1d,maux,aux1,aux2,aux3, & 1,amdq,bmasdq,bpasdq) c c # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q: do 160 m=1,meqn do 160 i = 1, mx+1 gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(i,m) gaddm(i-1,m,1) = gaddm(i-1,m,1) - gupdate gaddp(i-1,m,1) = gaddp(i-1,m,1) - gupdate c gupdate = 0.5d0*dtdx1d(i-1) * bpasdq(i,m) gaddm(i-1,m,2) = gaddm(i-1,m,2) - gupdate gaddp(i-1,m,2) = gaddp(i-1,m,2) - gupdate 160 continue c c # split the right-going flux difference into down-going and up-going: call rpt2(ixy,maxm,meqn,mwaves,mbc,mx, & q1d,q1d,maux,aux1,aux2,aux3, & 2,apdq,bmasdq,bpasdq) c c # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q: do 180 m=1,meqn do 180 i = 1, mx+1 gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(i,m) gaddm(i,m,1) = gaddm(i,m,1) - gupdate gaddp(i,m,1) = gaddp(i,m,1) - gupdate c gupdate = 0.5d0*dtdx1d(i-1) * bpasdq(i,m) gaddm(i,m,2) = gaddm(i,m,2) - gupdate gaddp(i,m,2) = gaddp(i,m,2) - gupdate 180 continue c 999 continue return end c c c =================================================================== subroutine slope2(ixy,maxm,meqn,maux,mwaves,mbc,mx, & q,aux,dx,dt,method,mthlim, & wave,s,amdq,apdq,dtdx,rpn2, & ql,qr,fl,fr) c =================================================================== c c Author: Ralf Deiterding c c # Implements the standard MUSCL-Hancock method, which gives 2nd c # order accuracy for conventional finite-volume schemes that return c # the numerical fluxes instead of fluctuations. c # The variable reconstruction is one-dimensional and is intended c # to be used with dimensional splitting. c implicit double precision (a-h,o-z) 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 rpn2 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,m) = q(-1,m) qr(-1,m) = q(-1,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 rec2(ixy,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 # 2nd order spatial reconstuction: om!=0.d0 c # 3rd order spatial reconstruction: om=1.d0/3.d0 c # Higher order reconstructions 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 flx2(ixy,maxm,meqn,mbc,mx,ql,maux,aux,fl) call flx2(ixy,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 rpn2(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,aux,aux, & wave,s,amdq,apdq) c return end c
Quickstart Users Guide Programmers Reference Installation Examples Download
AMROC Main Home Contactlast update: 06/01/04