advanc.f.html | |
Source file: advanc.f | |
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/amrclaw/2d/lib | |
Converted: Mon Jan 21 2013 at 20:15:09 using clawcode2html | |
This documentation file will not reflect any later changes in the source file. |
c c -------------------------------------------------------------- c subroutine advanc (level,nvar,dtlevnew,vtime,naux) c implicit double precision (a-h,o-z) include "call.i" logical vtime integer omp_get_thread_num, omp_get_max_threads integer mythread/0/, maxthreads/1/ integer listgrids(numgrids(level)) c maxgr is maximum number of grids many things are c dimensioned at, so this is overall. only 1d array c though so should suffice. problem is c not being able to dimension at maxthreads integer locfp_save(maxgr) integer locgp_save(maxgr) c c ::::::::::::::; ADVANC ::::::::::::::::::::::::::::::::::::::::::: c integrate all grids at the input 'level' by one step of its delta(t) c this includes: setting the ghost cells c advancing the solution on the grid c adjusting fluxes for flux conservation step later c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c hx = hxposs(level) hy = hyposs(level) delt = possk(level) c this is linear alg. call prepgrids(listgrids,numgrids(level),level) c maxthreads initialized to 1 above in case no openmp !$ maxthreads = omp_get_max_threads() !$OMP PARALLEL DO PRIVATE(j,locnew, locaux, mptr,nx,ny,mitot !$OMP& ,mjtot,time), !$OMP& SHARED(level,nvar,naux,alloc,intrat,delt, !$OMP& nghost,node,rnode,numgrids,listgrids), !$OMP& SCHEDULE (dynamic,1) ccccccc!$OMP& DEFAULT(none), do j = 1, numgrids(level) c mget is n^2 alg. c mptr = mget(j,level) mptr = listgrids(j) nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost locnew = node(store1,mptr) locaux = node(storeaux,mptr) time = rnode(timemult,mptr) c call bound(time,nvar,nghost,alloc(locnew),mitot,mjtot,mptr, 1 alloc(locaux),naux) end do !$OMP END PARALLEL DO c c save coarse level values if there is a finer level for wave fixup if (level+1 .le. mxnest) then if (lstart(level+1) .ne. null) then call saveqc(level+1,nvar,naux) endif endif c dtlevnew = rinfinity cfl_level = 0.d0 !# to keep track of max cfl seen on each level c c each grid needs to get extra storage for integration serially (for now at least) in case c dynamic memory moves the alloc array during expansion. so get it once and for all for c biggest chunk it may need. That uses max1d on a side, which for parallel apps. c is generally much smaller than serial, though this could be a waste. c remember that max1d includes ghost cells c c if necessary can have each thread use max of grids it owns, but then c cant use dynamic grid assignments. c c next loop will get enough storage, it wont be dimensioned as called for c grids smaller than max1d on a side. do j = 1, maxthreads locfp_save(j) = igetsp(2*max1d*max1d*nvar) locgp_save(j) = igetsp(2*max1d*max1d*nvar) end do !$OMP PARALLEL DO PRIVATE(j,locold, locnew, mptr,nx,ny,mitot,mjtot) !$OMP& PRIVATE(locfp, locfm, locgp,locgm,ntot,xlow,ylow) !$OMP& PRIVATE(locaux,lenbc,locsvf,locsvq,locx1d,i) !$OMP& PRIVATE(dtnew,time,mythread) !$OMP& SHARED(rvol,rvoll,level,nvar,mxnest,alloc,intrat,delt) !$OMP& SHARED(nghost,intratx,intraty,hx,hy,naux,listsp) !$OMP& SHARED(node,rnode,dtlevnew,numgrids,listgrids) !$OMP& SHARED(locfp_save,locgp_save) !$OMP& SCHEDULE (DYNAMIC,1) !$OMP& DEFAULT(none) do j = 1, numgrids(level) c mptr = mget(j,level) mptr = listgrids(j) locold = node(store2, mptr) locnew = node(store1, mptr) nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 time = rnode(timemult,mptr) c mitot = nx + 2*nghost mjtot = ny + 2*nghost c ::: get scratch storage for fluxes and slopes !-- locfp = igetsp(mitot*mjtot*nvar) !-- locfm = igetsp(mitot*mjtot*nvar) !-- locgp = igetsp(mitot*mjtot*nvar) !-- locgm = igetsp(mitot*mjtot*nvar) c c next way so dont call igetsp so much, less parallel bottleneck in critical section !-- locfp = igetsp(2*mitot*mjtot*nvar) !-- locfm = locfp + mitot*mjtot*nvar !-- locgp = igetsp(2*mitot*mjtot*nvar) !-- locgm = locgp + mitot*mjtot*nvar c next way for dynamic memory enlargement safety !$ mythread = omp_get_thread_num() locfp = locfp_save(mythread+1) locfm = locfp + mitot*mjtot*nvar locgp = locgp_save(mythread+1) locgm = locgp + mitot*mjtot*nvar c c copy old soln. values into next time step's soln. values c since integrator will overwrite it. only for grids not at c the finest level. finest level grids do not maintain copies c of old and new time solution values. c if (level .lt. mxnest) then ntot = mitot * mjtot * nvar cdir$ ivdep do 10 i = 1, ntot 10 alloc(locold + i - 1) = alloc(locnew + i - 1) endif c xlow = rnode(cornxlo,mptr) - nghost*hx ylow = rnode(cornylo,mptr) - nghost*hy !$OMP CRITICAL(rv) rvol = rvol + nx * ny rvoll(level) = rvoll(level) + nx * ny !$OMP END CRITICAL(rv) locaux = node(storeaux,mptr) c if (node(ffluxptr,mptr) .ne. 0) then lenbc = 2*(nx/intratx(level-1)+ny/intraty(level-1)) locsvf = node(ffluxptr,mptr) locsvq = locsvf + nvar*lenbc locx1d = locsvq + nvar*lenbc call qad(alloc(locnew),mitot,mjtot,nvar, 1 alloc(locsvf),alloc(locsvq),lenbc, 2 intratx(level-1),intraty(level-1),hx,hy, 3 naux,alloc(locaux),alloc(locx1d),delt,mptr) endif c # see if the grid about to advanced has gauge data to output c # this corresponds to previous time step, but output done c # now to make linear interpolation easier, since grid c # now has boundary conditions filled in. c # no testing here for mgauges>0 so that do not c # need to use gauges.i. the only time advanc is c # called that isn't "real" is in the initial setting c # up of grids (setgrd), but source grids are 0 there so c # nothing will be output. c should change the way dumpguage does io - right now is critical section call dumpgauge(alloc(locnew),alloc(locaux),xlow,ylow, . nvar,mitot,mjtot,mptr) c call stepgrid(alloc(locnew),alloc(locfm),alloc(locfp), 1 alloc(locgm),alloc(locgp), 2 mitot,mjtot,nghost, 3 delt,dtnew,hx,hy,nvar, 4 xlow,ylow,time,mptr,naux,alloc(locaux)) if (node(cfluxptr,mptr) .ne. 0) 1 call fluxsv(mptr,alloc(locfm),alloc(locfp), 2 alloc(locgm),alloc(locgp), 3 alloc(node(cfluxptr,mptr)),mitot,mjtot, 4 nvar,listsp(level),delt,hx,hy) if (node(ffluxptr,mptr) .ne. 0) then lenbc = 2*(nx/intratx(level-1)+ny/intraty(level-1)) locsvf = node(ffluxptr,mptr) call fluxad(alloc(locfm),alloc(locfp), 1 alloc(locgm),alloc(locgp), 2 alloc(locsvf),mptr,mitot,mjtot,nvar, 4 lenbc,intratx(level-1),intraty(level-1), 5 nghost,delt,hx,hy) endif c !-- call reclam(locfp, mitot*mjtot*nvar) !-- call reclam(locfm, mitot*mjtot*nvar) !-- call reclam(locgp, mitot*mjtot*nvar) !-- call reclam(locgm, mitot*mjtot*nvar) c next way to reclaim was to minimize calls to c reclam, due to critical section and openmp !-- call reclam(locfp, 2*mitot*mjtot*nvar) !-- call reclam(locgp, 2*mitot*mjtot*nvar) !$OMP CRITICAL (newdt) dtlevnew = dmin1(dtlevnew,dtnew) !$OMP END CRITICAL (newdt) c rnode(timemult,mptr) = rnode(timemult,mptr)+delt end do !$OMP END PARALLEL DO c c debug statement: c write(*,*)" from advanc: level ",level," dtlevnew ",dtlevnew c new way to reclaim for safety with dynamic memory and openmp do j = 1, maxthreads call reclam(locfp_save(j),2*max1d*max1d*nvar) call reclam(locgp_save(j),2*max1d*max1d*nvar) end do c return end c c ------------------------------------------------------------- c integer function mget(j,level) implicit double precision (a-h,o-z) include "call.i" integer this_j mptr = lstart(level) this_j = 1 do while (this_j < j) mptr = node(levelptr, mptr) this_j = this_j + 1 if (this_j > numgrids(level)) then write(*,*)" *** ERROR: exceeded number of grids in mget *** " stop endif end do mget = mptr return end c c ------------------------------------------------------------- c subroutine prepgrids(listgrids,num, level) implicit double precision (a-h,o-z) include "call.i" integer listgrids(num) mptr = lstart(level) do j = 1, num listgrids(j) = mptr mptr = node(levelptr, mptr) end do if (mptr .ne. 0) then write(*,*)" Error in routine setting up grid array " stop endif return end