Blockstructured Adaptive Mesh Refinement in object-oriented C++
c c c ========================================================== subroutine unsp2(maxm,maxmx,maxmy,mvar,meqn,maux,mwaves,mbc, & mx,my,qold,aux,dx,dy,dt,method,mthlim,cflgrid, & fm,fp,gm,gp,faddm,faddp,gaddm,gaddp, & q1d,dtdx1d,dtdy1d,aux1,aux2,aux3, & work,mwork,rpn2,rpt2) c ========================================================== c c # Take one time step, updating q with the wave propagation method c # of Randall J. LeVeque. c c # fm, fp are fluxes to left and right of single cell edge, c # gm, gp are fluxes at bottom and top. c # fadd and gadd are used to return flux increments from flux2. c # See the flux2 documentation for more information. c implicit double precision (a-h,o-z) include "call.i" c dimension qold(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension fm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension fp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension gm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension gp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension q1d(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 aux(maux, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension aux1(1-mbc:maxm+mbc, maux) dimension aux2(1-mbc:maxm+mbc, maux) dimension aux3(1-mbc:maxm+mbc, maux) dimension dtdx1d(1-mbc:maxm+mbc) dimension dtdy1d(1-mbc:maxm+mbc) dimension method(7),mthlim(mwaves) dimension work(mwork) external rpn2,rpt2 c c # partition work array into pieces needed for local storage in c # flux2 routine. Find starting index of each piece: c i0wave = 1 i0s = i0wave + (maxm+2*mbc)*meqn*mwaves i0amdq = i0s + (maxm+2*mbc)*mwaves i0apdq = i0amdq + (maxm+2*mbc)*meqn i0cqxx = i0apdq + (maxm+2*mbc)*meqn i0bmadq = i0cqxx + (maxm+2*mbc)*meqn i0bpadq = i0bmadq + (maxm+2*mbc)*meqn i0slope = i0bpadq + (maxm+2*mbc)*meqn iused = i0slope + (maxm+2*mbc)*meqn*4 - 1 mslope = iused-i0slope+1 c if (iused.gt.mwork) then write(6,*) '*** not enough work space in unsp2' write(6,*) '*** iused = ', iused, ' mwork =',mwork stop endif c mcapa = method(6) c cflgrid = 0.d0 dtdx = dt/dx dtdy = dt/dy c do 10 j=1-mbc,my+mbc do 10 i=1-mbc,mx+mbc do 10 m=1,mvar fm(m,i,j) = 0.d0 fp(m,i,j) = 0.d0 gm(m,i,j) = 0.d0 gp(m,i,j) = 0.d0 10 continue c if (mcapa.eq.0) then c # no capa array: do 5 i=1-mbc,maxm+mbc dtdx1d(i) = dtdx dtdy1d(i) = dtdy 5 continue endif c c c # perform x-sweeps c ================== c do 50 j = 0,my+1 c c # copy data along a slice into 1d arrays: do 20 i = 1-mbc, mx+mbc do 20 m=1,meqn q1d(i,m) = qold(m,i,j) 20 continue c if (mcapa.gt.0) then do 21 i = 1-mbc, mx+mbc dtdx1d(i) = dtdx / aux(mcapa,i,j) 21 continue endif c if (maux .gt. 0) then do 22 i = 1-mbc, mx+mbc do 22 ma=1,maux aux1(i,ma) = aux(ma,i,j-1) aux2(i,ma) = aux(ma,i,j) aux3(i,ma) = aux(ma,i,j+1) 22 continue endif c c c # Store the value of j along this slice in the common block c # comxyt in case it is needed in the Riemann solver (for c # variable coefficient problems) jcom = j c c # compute modifications fadd and gadd to fluxes along this slice: call flux2(1,maxm,meqn,maux,mwaves,mbc,mx, & q1d,dtdx1d,aux1,aux2,aux3,method,mthlim, & faddm,faddp,gaddm,gaddp,cfl1d, & work(i0wave),work(i0s),work(i0amdq),work(i0apdq), & work(i0cqxx),work(i0bmadq),work(i0bpadq), & work(i0slope),mslope,rpn2,rpt2) c cflgrid = dmax1(cflgrid,cfl1d) c c # update fluxes for use in AMR: c do 25 i=1,mx+1 do 25 m=1,meqn fm(m,i,j) = fm(m,i,j) + faddm(i,m) fp(m,i,j) = fp(m,i,j) + faddp(i,m) gm(m,i,j) = gm(m,i,j) + gaddm(i,m,1) gp(m,i,j) = gp(m,i,j) + gaddp(i,m,1) gm(m,i,j+1) = gm(m,i,j+1) + gaddm(i,m,2) gp(m,i,j+1) = gp(m,i,j+1) + gaddp(i,m,2) 25 continue 50 continue c c c c # perform y sweeps c ================== c c do 100 i = 0, mx+1 c c # copy data along a slice into 1d arrays: do 70 j = 1-mbc, my+mbc do 70 m=1,meqn q1d(j,m) = qold(m,i,j) 70 continue c if (mcapa.gt.0) then do 71 j = 1-mbc, my+mbc dtdy1d(j) = dtdy / aux(mcapa,i,j) 71 continue endif c if (maux .gt. 0) then do 72 j = 1-mbc, my+mbc do 72 ma=1,maux aux1(j,ma) = aux(ma,i-1,j) aux2(j,ma) = aux(ma,i, j) aux3(j,ma) = aux(ma,i+1,j) 72 continue endif c c c # Store the value of i along this slice in the common block c # comxyt in case it is needed in the Riemann solver (for c # variable coefficient problems) icom = i c c # compute modifications fadd and gadd to fluxes along this slice: call flux2(2,maxm,meqn,maux,mwaves,mbc,my, & q1d,dtdy1d,aux1,aux2,aux3,method,mthlim, & faddm,faddp,gaddm,gaddp,cfl1d, & work(i0wave),work(i0s),work(i0amdq),work(i0apdq), & work(i0cqxx),work(i0bmadq),work(i0bpadq), & work(i0slope),mslope,rpn2,rpt2) c cflgrid = dmax1(cflgrid,cfl1d) c c # c # update fluxes for use in AMR: c do 75 j=1,my+1 do 75 m=1,meqn gm(m,i,j) = gm(m,i,j) + faddm(j,m) gp(m,i,j) = gp(m,i,j) + faddp(j,m) fm(m,i,j) = fm(m,i,j) + gaddm(j,m,1) fp(m,i,j) = fp(m,i,j) + gaddp(j,m,1) fm(m,i+1,j) = fm(m,i+1,j) + gaddm(j,m,2) fp(m,i+1,j) = fp(m,i+1,j) + gaddp(j,m,2) 75 continue 100 continue c do 110 j=1,my do 110 i=1,mx do 110 m=1,meqn if (mcapa.gt.0) then c # with capa array. qold(m,i,j) = qold(m,i,j) & - (dtdx * (fm(m,i+1,j) - fp(m,i,j)) & + dtdy * (gm(m,i,j+1) - gp(m,i,j))) / & aux(mcapa,i,j) c # no capa array. Standard flux differencing: else qold(m,i,j) = qold(m,i,j) & - dtdx * (fm(m,i+1,j) - fp(m,i,j)) & - dtdy * (gm(m,i,j+1) - gp(m,i,j)) endif 110 continue return end
Quickstart Users Guide Programmers Reference Installation Examples Download
AMROC Main Home Contactlast update: 06/01/04