Skip to content

Commit

Permalink
Merge pull request clawpack#174 from rjleveque/test_merge_fixProperNe…
Browse files Browse the repository at this point in the history 
…sting

Fix proper nesting merged with recent changes to master
  • Loading branch information
@rjleveque
rjleveque authored Sep 24, 2016
2 parents a592b7c + ab69bfe commit cab6449
Show file tree
Hide file tree
Showing 18 changed files with 951 additions and 820 deletions.
181 changes: 89 additions & 92 deletions src/2d/baseCheck.f
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,16 @@
c
c ----------------------------------------------------------------
c
logical function baseCheck(mnew,lbase,ilo,ihi,jlo,jhi)
logical function baseCheck(mnew,lbase,ilo,ihi,jlo,jhi,
. nvar,naux,thisBuff)

use amr_module
implicit double precision (a-h, o-z)

logical debug/.false./
logical debug/.true./
integer ist(3),iend(3),jst(3),jend(3),ishift(3),jshift(3)
logical borderx, bordery
integer thisBuff

c index into alloc from iclo:ichi and jclo:jchi, not 0..leni/j.
iadd(i,j) = locm + i - iclo + leni*(j-jclo)
Expand All @@ -28,19 +31,14 @@ logical function baseCheck(mnew,lbase,ilo,ihi,jlo,jhi)



if (debug) write(outunit,100) mnew,ilo,ihi,jlo,jhi,lbase
100 format("NESTCK2 testing grid ",i5," base level ",i5,/,
. " new grid from ilo:hi: ",2i12," to ",2i12)

levnew = node(nestlevel,mnew)
c
c initialize for potential periodicity
c each patch divided into 9 regions (some may be empty)
c e.g. i from (ilo,-1), (0,iregsz(level)-1),(iregsz(level),ihi)
c except using enlarged grid (ilo-1 to ihi+1)
c
call setIndices(ist,iend,jst,jend,ilo-1,ihi+1,jlo-1,jhi+1,
. ishift,jshift,levnew)
borderx = (ilo .eq. 0 .or. ihi .eq. iregsz(levnew)-1)
bordery = (jlo .eq. 0 .or. jhi .eq. jregsz(levnew)-1)


if (debug) write(outunit,100) mnew,lbase,ilo,ihi,jlo,jhi,levnew
100 format("NESTCK2 testing grid ",i5," base level ",i5,/,
. " new grid from ilo:hi: ",2i12," to ",2i12," at level ",i4)
c
c on to initializing for the given grid and its nest checking
levratx = 1
Expand All @@ -51,15 +49,21 @@ logical function baseCheck(mnew,lbase,ilo,ihi,jlo,jhi)
5 continue

c widen by 1 cell (proper nesting), then project to lbase
iclo = ilo-1
ichi = ihi+1
jclo = jlo-1
jchi = jhi+1
c this might stick out of domain, fix later
c figure out size for scratch storage on base grid for testing
iclo = ilo
ichi = ihi
jclo = jlo
jchi = jhi
do lev = levnew-1,lbase,-1
iclo = (dfloat(iclo)/intratx(lev))
ichi = (dfloat(ichi)/intratx(lev))
jclo = (dfloat(jclo)/intraty(lev))
jchi = (dfloat(jchi)/intraty(lev))
iclo = iclo/intratx(lev)
ichi = ichi/intratx(lev)
jclo = jclo/intraty(lev)
jchi = jchi/intraty(lev)
iclo = iclo - 1
ichi = ichi + 1
jclo = jclo - 1
jchi = jchi + 1
if (debug) then
write(outunit,111) lev, iclo,ichi,jclo,jchi
111 format(10x,"at level",i5," projected coords ilo:hi:",2i10,
Expand All @@ -73,29 +77,28 @@ logical function baseCheck(mnew,lbase,ilo,ihi,jlo,jhi)
if (debug) then
write(outunit,108) ilo-1,ihi+1,jlo-1,jhi+1
write(outunit,109) levratx,levraty
108 format(" enlarged fine grid from ilo:hi:",2i12,
108 format(" enlarged (by 1) fine grid from ilo:hi:",2i12,
. " to jlo:hi:", 2i12)
109 format(" refinement factors to base grid of ", 2i12)
write(outunit,101) iclo,ichi,jclo,jchi
101 format("coarsened to lbase, grid from ilo:hi: ",2i12,
. " to jlo:hi:",2i12)
101 format("coarsened to lbase, grid from iclo:hi: ",2i12,
. " to jclo:hi:",2i12)
endif

if (.not. (xperdom .and. borderx) .and.
. .not. (yperdom .and. bordery)) then
iclo = max(iclo,0) ! make sure in domain boundary when checking nesting
jclo = max(jclo,0)
ichi = min(ichi,iregsz(lbase)-1) ! subtract 1 since regsz is number of cells, so -1 is highest index
jchi = min(jchi,jregsz(lbase)-1)
endif


leni = ichi - iclo + 1
lenj = jchi - jclo + 1
lenrect = leni * lenj
locm = igetsp(lenrect)
c initialize on projected size of new grid. used to mark nesting
do k = 0, lenrect-1
alloc(locm + k) = 0.
end do
c
c corners dont count for nesting so mark as ok
c leave 0 for now so match older nestck
!-- alloc(iadd(iclo,jclo)) = 1.
!-- alloc(iadd(iclo,jchi)) = 1.
!-- alloc(iadd(ichi,jclo)) = 1.
!-- alloc(iadd(ichi,jchi)) = 1.
alloc(locm:locm+lenrect-1) = 0.
c
c if mnew on domain boundary fix flags so ok.
c fix extra border, and first/last real edge
Expand Down Expand Up @@ -126,73 +129,67 @@ logical function baseCheck(mnew,lbase,ilo,ihi,jlo,jhi)
endif

mptr = lstart(lbase)
20 iblo = node(ndilo, mptr)
ibhi = node(ndihi, mptr)
jblo = node(ndjlo, mptr)
jbhi = node(ndjhi, mptr)

if (debug) then
write(outunit,*)
write(outunit,102) mptr,lbase,iblo,ibhi,jblo,jbhi
102 format("TESTING AGAINST GRID ",i5," level ",i5,
. " from ilo:hi:",2i10," to jlo:hi:",2i10)
endif
20 iblo = node(ndilo, mptr) - thisBuff
ibhi = node(ndihi, mptr) + thisBuff
jblo = node(ndjlo, mptr) - thisbuff
jbhi = node(ndjhi, mptr) + thisBuff
c
! non periodic case, base level coordinates, just mark if nested.
if ((.not. (xperdom .and. borderx)) .and.
. .not. (yperdom .and. bordery)) then
ixlo = max(iclo,iblo)
ixhi = min(ichi,ibhi)
jxlo = max(jclo,jblo)
jxhi = min(jchi,jbhi)
if (.not.((ixlo.le.ixhi) .and. (jxlo.le.jxhi))) go to 30
do jx = jxlo, jxhi
do ix = ixlo, ixhi
alloc(iadd(ix,jx))=1.
end do
end do
go to 30
endif
c
c periodic case: initialize for potential periodicity
c each patch divided into 9 regions (some may be empty)
c e.g. i from (ilo,-1), (0,iregsz(level)-1),(iregsz(level),ihi)
c except using enlarged grid (ilo-1 to ihi+1)
c
call setIndices(ist,iend,jst,jend,iclo,ichi,jclo,jchi,
. ishift,jshift,lbase)

c compare all regions of patch with one lbase grid at a time
c compare all regions of coarsened patch with one lbase grid at a time
do 25 i = 1, 3
i1 = max(ilo-1,ist(i))
i2 = min(ihi+1, iend(i))
i1 = max(iclo,ist(i))
i2 = min(ichi, iend(i))
do 25 j = 1, 3
j1 = max(jlo-1, jst(j))
j2 = min(jhi+1, jend(j))
j1 = max(jclo, jst(j))
j2 = min(jchi, jend(j))

if (.not. ((i1 .le. i2) .and. (j1 .le. j2))) go to 25

if (debug) write(outunit,103)i,j,i1,i2,j1,j2
103 format("region ",2i10," from ilo:hi:",2i10,
. " to jlo:ji:",2i10)
c
c patch not empty. project coords to level lbase
iplo = i1+ishift(i)
iphi = i2+ishift(i)
jplo = j1+jshift(j)
jphi = j2+jshift(j)
do lev = levnew-1,lbase,-1
iplo = floor(dfloat(iplo)/intratx(lev))
iphi = ceiling(dfloat(iphi+1)/intratx(lev) - 1)
jplo = floor(dfloat(jplo)/intraty(lev))
jphi = ceiling(dfloat(jphi+1)/intraty(lev) - 1)
end do
if (debug) then
write(outunit,104) i,j,iplo,iphi,jplo,jphi
104 format("projected coords of region ",2i15," is ",2i12,
. " by ",2i12)
endif

ixlo = max(iplo,iblo)
ixhi = min(iphi,ibhi)
jxlo = max(jplo,jblo)
jxhi = min(jphi,jbhi)
c patch (possibly periodically wrapped) not empty.
c see if intersects base grid. wrap coords for periodicity
i1_shifted = i1 + ishift(i)
i2_shifted = i2 + ishift(i)
j1_shifted = j1 + jshift(j)
j2_shifted = j2 + jshift(j)

ixlo = max(i1_shifted,iblo)
ixhi = min(i2_shifted,ibhi)
jxlo = max(j1_shifted,jblo)
jxhi = min(j2_shifted,jbhi)

if (.not.((ixlo.le.ixhi) .and. (jxlo.le.jxhi))) go to 25

if (debug) write(outunit,105) ixlo,ixhi,jxlo,jxhi
105 format("intersected reg ",2i12," by ",2i12)
c FOR NOW DO NOT HANDLE PERIODIC NEST CHECKING

if (debug) then
write(outunit,106) ixlo,ixhi,jxlo,jxhi
write(outunit,107) ishift(i),jshift(j)
106 format("SETTING FROM ",4i5)
107 format("shifts are ",2i5)
endif
c mark intersected regions with 1
do jx = jxlo, jxhi
do ix = ixlo, ixhi
c need to mark nesting of orig coords, not shifted
alloc(iadd(ix-ishift(i)/levratx,jx-jshift(j)/levraty))=1.
end do
end do
c need to mark nesting of orig coords, not coarsened shifted indices
ix_unshifted = (ix - ishift(i)) ! back to unshifted coords
jx_unshifted = (jx - jshift(j)) ! to mark base grid nesting ok
alloc(iadd(ix_unshifted,jx_unshifted)) = 1.
end do
end do

25 continue

Expand Down
21 changes: 13 additions & 8 deletions src/2d/colate2.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ subroutine colate2 (badpts, len, lcheck, nUniquePts, lbase)
dimension badpts(2,len)
dimension ist(3), iend(3), jst(3), jend(3), ishift(3), jshift(3)
logical db/.false./
integer*8 largestIntEquiv

c
c index for flag array now based on integer index space, not 1:mibuff,1:mjbuff
Expand Down Expand Up @@ -162,15 +163,19 @@ subroutine colate2 (badpts, len, lcheck, nUniquePts, lbase)
c
c colate flagged points into single integer array for quicksorting
c
c call drivesort(npts,badpts,lcheck,nUniquePts,mbuff)
c ### bug found in drivesort- integer overflow
c ### temp bypass for rnady to run finer grids
if (lcheck .ge. 6) then
nUniquePts = npts
c sorting uses one dimensional packing of 2D indices
c check if domain will fit in integer*4.
c if not, just leave the duplicate, but rememer that efficiency
c of grids won't be correct (since divide by number of flaged pts in grid)
c If necessary, do whole process in integer*8 - then will have enough
c room, but will have to convert quicksort routine and drivesort
c the variable largestIntEquiv already declared integer*8 above.
largestIntEquiv = iregsz(lcheck)+mbuff +
. (iregsz(lcheck)+2*mbuff)*(jregsz(lcheck)+mbuff)
largestSingle = 2**30
if (largestSingle .le. largestIntEquiv) then
nUniquePts = npts ! bad name - they are not unique
else
c # assume points are unique
c # Cannot check for more than 6 levels due to integer overflow
c # bug that needs to be fixed!
call drivesort(npts,badpts,lcheck,nUniquePts,mbuff)
endif

Expand Down
42 changes: 26 additions & 16 deletions src/2d/filval.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,14 @@
!
! ------------------------------------------------------------------
!
subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
mjc, xleft, xright, ybot, ytop, nvar, mptr, ilo, ihi, &
subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
mjc, xleft, xright, ybot, ytop, nvar, mptr, ilo, ihi, &
jlo, jhi, aux, naux)

use amr_module, only: xlower, ylower, intratx, intraty, nghost, xperdom
use amr_module, only: yperdom, spheredom, xupper, yupper, alloc
use amr_module, only: outunit, NEEDS_TO_BE_SET, mcapa
use amr_module, only: newstl, iregsz, jregsz

implicit none

Expand All @@ -34,6 +35,9 @@ subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
real(kind=8) :: fliparray((mitot+mjtot)*nghost*(nvar+naux))
real(kind=8) :: setflags(mitot,mjtot),maxauxdif,aux2(naux,mitot,mjtot)
integer :: mjb
integer :: jm, im, nm
logical :: sticksoutxfine, sticksoutyfine,sticksoutxcrse,sticksoutycrse
logical :: DIAGONAL_CORNER

!for setaux timing
integer :: clock_start, clock_finish, clock_rate
Expand All @@ -43,6 +47,13 @@ subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
! External function definitions
real(kind=8) :: get_max_speed


DIAGONAL_CORNER(im,jm,mic,mjc) = (im .eq. 1 .and. jm .eq. mjc) .or. &
(im .eq. mic .and. jm .eq. mjc) .or. &
(im .eq. 1 .and. jm .eq. 1) .or. &
(im .eq. mic .and. jm .eq. 1)


refinement_ratio_x = intratx(level-1)
refinement_ratio_y = intraty(level-1)
dx_coarse = dx * refinement_ratio_x
Expand All @@ -64,16 +75,21 @@ subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
jchi = (jhi + 1) / refinement_ratio_y - 1 + 1
ng = 0

sticksoutxfine = ( (ilo .lt. 0) .or. (ihi .ge. iregsz(level)))
sticksoutyfine = ( (jlo .lt. 0) .or. (jhi .ge. jregsz(level)))
sticksoutxcrse = ((iclo .lt. 0) .or. (ichi .ge. iregsz(level-1)))
sticksoutycrse = ((jclo .lt. 0) .or. (jchi .ge. jregsz(level-1)))

if (naux == 0) then
if (xperdom .or. yperdom .or. spheredom) then
if ((xperdom .and. sticksoutxcrse).or. (yperdom.and.sticksoutycrse) .or. spheredom) then
call preintcopy(valc,mic,mjc,nvar,iclo,ichi,jclo,jchi,level-1,fliparray)
else
call intcopy(valc,mic,mjc,nvar,iclo,ichi,jclo,jchi,level-1,1,1)
endif
else
! intersect grids and copy all (soln and aux)
auxc(1,:,:) = NEEDS_TO_BE_SET
if (xperdom .or. yperdom .or. spheredom) then
if ((xperdom .and.sticksoutxcrse) .or. (yperdom.and.sticksoutycrse) .or. spheredom) then
call preicall(valc,auxc,mic,mjc,nvar,naux,iclo,ichi,jclo,jchi, &
level-1,fliparray)
else
Expand All @@ -90,14 +106,12 @@ subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
!!$ end do
! no ghost cells on coarse enlarged patch. set any remaining
! vals. should only be bcs that stick out of domain.
call system_clock(clock_start, clock_rate)
call cpu_time(cpu_start)
call setaux(ng,mic,mjc,xl,yb,dx_coarse,dy_coarse,naux,auxc)
call system_clock(clock_finish, clock_rate)
call cpu_time(cpu_finish)
endif
call bc2amr(valc,auxc,mic,mjc,nvar,naux,dx_coarse,dy_coarse,level-1, &
time,xl,xr,yb,yt)

call bc2amr(valc,auxc,mic,mjc,nvar,naux,dx_coarse,dy_coarse,level-1,time,xl,xr,yb, &
yt,xlower,ylower,xupper,yupper,xperdom,yperdom,spheredom)



! NOTE change in order of code. Since the interp from coarse to fine needs the aux
Expand All @@ -117,7 +131,7 @@ subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
! ## overwritten with coarse grid interpolation
aux(1,:,:) = NEEDS_TO_BE_SET ! indicates fine cells not yet set.

if (xperdom .or. yperdom .or. spheredom) then
if ((xperdom.and.sticksoutxfine) .or. (yperdom.and.sticksoutyfine) .or. spheredom) then
call preicall(val,aux,mitot,mjtot,nvar,naux,ilo-nghost,ihi+nghost, &
jlo-nghost,jhi+nghost,level,fliparray)
else
Expand All @@ -127,18 +141,14 @@ subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
setflags = aux(1,:,:) ! save since will overwrite in setaux when setting all aux vals
! need this so we know where to use coarse grid to set fine solution w/o overwriting
! set remaining aux vals not set by copying from prev existing grids
call system_clock(clock_start, clock_rate)
call cpu_time(cpu_start)
call setaux(nghost,nx,ny,xleft,ybot,dx,dy,naux,aux)
call system_clock(clock_finish, clock_rate)
call cpu_time(cpu_finish)
else ! either no aux exists, or cant reuse yet
! so only call intcopy (which copies soln) and not icall.
! in this case flag q(1,:) to NEEDS_TO_BE_SET flag so wont be overwritten
! by coarse grid interp. this is needed due to reversing order of
! work - first copy from fine grids, then interpolate from coarse grids
val(1,:,:) = NEEDS_TO_BE_SET
if (xperdom .or. yperdom .or. spheredom) then
if ((xperdom.and.sticksoutxfine) .or. (yperdom.and.sticksoutyfine) .or. spheredom) then
call preintcopy(val,mitot,mjtot,nvar,ilo-nghost,ihi+nghost, &
jlo-nghost,jhi+nghost,level,fliparray)
else
Expand Down
Loading

0 comments on commit cab6449

Please sign in to comment.