diff --git a/src/3d/Makefile.amr_3d b/src/3d/Makefile.amr_3d index 0731cd658..dc148f87f 100644 --- a/src/3d/Makefile.amr_3d +++ b/src/3d/Makefile.amr_3d @@ -101,5 +101,10 @@ COMMON_SOURCES += \ $(AMRLIB)/resize_alloc.f90 \ $(AMRLIB)/restrt_alloc.f90 \ $(AMRLIB)/inlinelimiter.f \ + $(AMRLIB)/init_nodes.f90 \ + $(AMRLIB)/restrt_nodes.f90 \ + $(AMRLIB)/resize_nodes.f90 \ + $(AMRLIB)/init_bndryList.f \ + $(AMRLIB)/resize_bndryList.f \ diff --git a/src/3d/amr_module.f90 b/src/3d/amr_module.f90 index 33aebd163..1e0841f46 100644 --- a/src/3d/amr_module.f90 +++ b/src/3d/amr_module.f90 @@ -59,7 +59,7 @@ module amr_module integer, parameter :: iplane = 1 integer, parameter :: jplane = 2 integer, parameter :: kplane = 3 - integer, parameter :: maxgr = 5000 + !! integer, parameter :: maxgr = 5000 ! No longer fixed integer, parameter :: maxlv = 10 integer, parameter :: maxcl = 5000 @@ -77,17 +77,30 @@ module amr_module ! order is coarsest level to finest. is redone after regridding ! and on restarting. note use of sentinel in listStart - integer :: listOfGrids(maxgr),listStart(0:maxlv+1) - integer, parameter :: bndListSize = 8*maxgr - integer :: bndList(bndListSize,2) ! guess size, average # nbors 6? manage as linnked list + !! CHANGED mjb 6/15/18 to make maxgr a variable + !!integer :: listOfGrids(maxgr),listStart(0:maxlv+1) + integer :: listStart(0:maxlv+1) + !integer, parameter :: bndListSize = 8*maxgr + !integer :: bndList(bndListSize,2) ! guess size, average # nbors 6? manage as linked list + integer :: bndListSize + integer,allocatable, dimension(:,:) :: bndList ! new way is allocatable - real(kind=8) hxposs(maxlv),hyposs(maxlv),hzposs(maxlv), & - possk(maxlv),rnode(rsize, maxgr) + !real(kind=8) hxposs(maxlv),hyposs(maxlv),hzposs(maxlv), & + ! possk(maxlv),rnode(rsize, maxgr) + real(kind=8) hxposs(maxlv),hyposs(maxlv),hzposs(maxlv),possk(maxlv) + + ! start of dynamic allocation for maxgr and associated arrays + integer maxgr + real(kind=8), allocatable, dimension(:,:) :: rnode + ! new way, use allocatable, not pointer + integer, allocatable, dimension(:,:) :: node + integer, allocatable, dimension(:) :: listOfGrids real(kind=8) tol, tolsp - integer ibuff, mstart,ndfree,ndfree_bnd,lfine,node(nsize, maxgr), & + !integer ibuff, mstart,ndfree,ndfree_bnd,lfine,node(nsize, maxgr), & + integer ibuff, mstart,ndfree,ndfree_bnd,lfine, & icheck(maxlv),lstart(maxlv),newstl(maxlv), & listsp(maxlv),intratx(maxlv),intraty(maxlv), & intratz(maxlv),kratio(maxlv), & diff --git a/src/3d/check.f b/src/3d/check.f index 5c4d9a4f3..5478d9fe6 100644 --- a/src/3d/check.f +++ b/src/3d/check.f @@ -60,7 +60,9 @@ subroutine check(nsteps,time,nvar,naux) c c ### dump the data c - write(chkunit) lenmax,lendim,memsize + ! new version has variable node size arrays, so need to output + ! current size to read back in on restart + write(chkunit) lenmax,lendim,memsize,maxgr write(chkunit) (alloc(i),i=1,lendim) write(chkunit) hxposs,hyposs,hzposs,possk,icheck write(chkunit) lfree,lenf diff --git a/src/3d/domain.f b/src/3d/domain.f index a22bc5586..a6ea9742d 100644 --- a/src/3d/domain.f +++ b/src/3d/domain.f @@ -89,7 +89,7 @@ subroutine domain (nvar,vtime,nx,ny,nz,naux,start_time) dz = hzposs(1) 60 mitot = node(ndihi,mptr)-node(ndilo,mptr) + 1 + 2*nghost mjtot = node(ndjhi,mptr)-node(ndjlo,mptr) + 1 + 2*nghost - mktot = node(ndkhi,mptr)-node(ndklo,mptr) + 1 + 2*nghost + mktot = node(ndkhi,mptr)-node(ndklo,mptr) + 1 + 2*nghost locaux = node(storeaux,mptr) c 4/1/02 : Added cfl to call to estdt, so that we dont need call.i in estdt call estdt(alloc(node(store1,mptr)),mitot,mjtot,mktot, diff --git a/src/3d/gauges_module.f90 b/src/3d/gauges_module.f90 index 762ab46c7..82fbed29f 100644 --- a/src/3d/gauges_module.f90 +++ b/src/3d/gauges_module.f90 @@ -75,13 +75,14 @@ module gauges_module type(gauge_type), allocatable :: gauges(:) ! Gauge source info - integer, allocatable, dimension(:) :: mbestsrc, mbestorder, & - igauge, mbestg1, mbestg2 + !!integer, allocatable, dimension(:) :: mbestsrc, mbestorder, & + !! igauge, mbestg1, mbestg2 + integer, allocatable, dimension(:) :: mbestsrc contains subroutine set_gauges(restart, num_eqn, num_aux, fname) - use amr_module, only: maxgr + !use amr_module, only: maxgr use utility_module, only: get_value_count implicit none @@ -111,8 +112,9 @@ subroutine set_gauges(restart, num_eqn, num_aux, fname) allocate(gauges(num_gauges)) ! Initialize gauge source data - allocate(mbestsrc(num_gauges), mbestorder(num_gauges)) - allocate(mbestg1(maxgr), mbestg2(maxgr)) + allocate(mbestsrc(num_gauges)) + !! changed mjb 6/15/18, didn't scale right to be maxgr dependent + !!allocate(mbestg1(maxgr), mbestg2(maxgr)) mbestsrc = 0 ! Original gauge information @@ -267,71 +269,39 @@ subroutine setbestsrc() ! # for debugging, initialize sources to 0 then check that all set mbestsrc = 0 - do lev = 1, lfine + do 40 i = 1, num_gauges + + do 30 lev = lfine, 1, -1 mptr = lstart(lev) - do - do i = 1, num_gauges - if ((gauges(i)%x >= rnode(cornxlo,mptr)) .and. & - (gauges(i)%x <= rnode(cornxhi,mptr)) .and. & - (gauges(i)%y >= rnode(cornylo,mptr)) .and. & - (gauges(i)%y <= rnode(cornyhi,mptr)) .and. & - (gauges(i)%z >= rnode(cornzlo,mptr)) .and. & - (gauges(i)%z <= rnode(cornzhi,mptr)) ) then - mbestsrc(i) = mptr - end if - end do - mptr = node(levelptr, mptr) - if (mptr == 0) exit - end do - end do + 20 if ((gauges(i)%x >= rnode(cornxlo,mptr)) .and. & + (gauges(i)%x <= rnode(cornxhi,mptr)) .and. & + (gauges(i)%y >= rnode(cornylo,mptr)) .and. & + (gauges(i)%y <= rnode(cornyhi,mptr)) .and. & + (gauges(i)%z >= rnode(cornzlo,mptr)) .and. & + (gauges(i)%z <= rnode(cornzhi,mptr)) ) then + mbestsrc(i) = mptr + !! best source found for this gauge, go on to next one + !! we knwo its best because we started at finst level + go to 40 ! on to next gauge + else + mptr = node(levelptr, mptr) + if (mptr .ne. 0) then + go to 20 ! try another grid + else + go to 30 ! try next coarser level grids + endif + endif + 30 continue + if (mbestsrc(i) .eq. 0) & + print *, "ERROR in setting grid src for gauge data", i + 40 continue - do i = 1, num_gauges - if (mbestsrc(i) .eq. 0) & - print *, "ERROR in setting grid src for gauge data", i - end do - -! Sort the source arrays for easy testing during integration - call qsorti(mbestorder,num_gauges,mbestsrc) - -! After sorting, -! mbestsrc(mbestorder(i)) = grid index to be used for gauge i -! and mbestsrc(mbestorder(i)) is non-decreasing as i=1,2,..., num_gauges - -! write(6,*) '+++ mbestorder: ',mbestorder -! write(6,*) '+++ mbestsrc: ',mbestsrc - -! Figure out the set of gauges that should be handled on each grid: -! after loop below, grid k should handle gauges numbered -! mbestorder(i) for i = mbestg1(k), mbestg1(k)+1, ..., mbestg2(k) -! This will be used for looping in print_gauges subroutine. - - ! initialize arrays to default indicating grids that contain no gauges: - mbestg1 = 0 - mbestg2 = 0 - - k1 = 0 - do i=1,num_gauges - ki = mbestsrc(mbestorder(i)) - if (ki > k1) then - ! new grid number seen for first time in list - if (k1 > 0) then - ! mark end of gauges seen by previous grid - mbestg2(k1) = i-1 -! write(6,*) '+++ k1, mbestg2(k1): ',k1,mbestg2(k1) - endif - mbestg1(ki) = i -! write(6,*) '+++ ki, mbestg1(ki): ',ki,mbestg1(ki) - endif - k1 = ki - enddo - if (num_gauges > 0) then - ! finalize - mbestg2(ki) = num_gauges -! write(6,*) '+++ ki, mbestg2(ki): ',ki,mbestg2(ki) - endif +!!! NO MORE qsort and mbestg arrays. +!!! Each grid now loops over mbestsrc array to see which gauges it owns. - end subroutine setbestsrc + + end subroutine setbestsrc ! ! ------------------------------------------------------------------------- @@ -373,22 +343,12 @@ subroutine update_gauges(q, aux, xlow, ylow, zlow, num_eqn, mitot, mjtot, & integer :: i, j, n, i1, i2, iindex, jindex, kindex, ii, index, level integer :: var_index - ! REMOVE - integer :: ivar - real(kind=8) :: var1, var2 - + ! No gauges to record, exit if (num_gauges == 0) then return endif - i1 = mbestg1(mptr) - i2 = mbestg2(mptr) - - if (i1 == 0) then - ! no gauges to be handled by this grid - return - endif ! Grid info tgrid = rnode(timemult, mptr) @@ -397,28 +357,25 @@ subroutine update_gauges(q, aux, xlow, ylow, zlow, num_eqn, mitot, mjtot, & hy = hyposs(level) hz = hzposs(level) + !! this loop refactor to remove mbestg, which was way larger + !! than gauge array. faster to loop thru gauges this way ! Main Gauge Loop ====================================================== - do i = i1, i2 - ii = mbestorder(i) - if (mptr /= mbestsrc(ii)) then - print *, '*** should not happen... i, ii, mbestsrc(ii), mptr:' - print *, i, ii, mbestsrc(ii), mptr - stop - endif - if (tgrid < gauges(ii)%t_start .or. tgrid > gauges(ii)%t_end) then + do i = 1, num_gauges + if (mptr /= mbestsrc(i)) cycle + if (tgrid < gauges(i)%t_start .or. tgrid > gauges(i)%t_end) then cycle end if ! Minimum increment ! TODO Maybe always allow last time output recording? - if (tgrid - gauges(ii)%last_time < gauges(ii)%min_time_increment) then + if (tgrid - gauges(i)%last_time < gauges(i)%min_time_increment) then cycle end if ! compute indexing and bilinear interpolant weights ! Note: changes 0.5 to 0.5d0 - iindex = int(0.5d0 + (gauges(ii)%x - xlow) / hx) - jindex = int(0.5d0 + (gauges(ii)%y - ylow) / hy) - kindex = int(0.5d0 + (gauges(ii)%z - zlow) / hz) + iindex = int(0.5d0 + (gauges(i)%x - xlow) / hx) + jindex = int(0.5d0 + (gauges(i)%y - ylow) / hy) + kindex = int(0.5d0 + (gauges(i)%z - zlow) / hz) if ((iindex < nghost .or. iindex > mitot-nghost) .or. & (jindex < nghost .or. jindex > mjtot-nghost) .or. & (kindex < nghost .or. kindex > mktot-nghost)) then @@ -427,9 +384,9 @@ subroutine update_gauges(q, aux, xlow, ylow, zlow, num_eqn, mitot, mjtot, & xcent = xlow + (iindex - 0.5d0) * hx ycent = ylow + (jindex - 0.5d0) * hy zcent = zlow + (kindex - 0.5d0) * hz - xoff = (gauges(ii)%x - xcent) / hx - yoff = (gauges(ii)%y - ycent) / hy - zoff = (gauges(ii)%z - zcent) / hz + xoff = (gauges(i)%x - xcent) / hx + yoff = (gauges(i)%y - ycent) / hy + zoff = (gauges(i)%z - zcent) / hz ! Gauge interpolation seems to work, so error test is commented out. ! For debugging, use the code below... @@ -440,14 +397,14 @@ subroutine update_gauges(q, aux, xlow, ylow, zlow, num_eqn, mitot, mjtot, & ! yoff .lt. -1.d-4 .or. yoff .gt. 1.0001d0 .or. & ! zoff .lt. 0.d0 .or. zoff .gt. 1.d0) then ! write(6,*) "*** print_gauges: Interpolation problem at gauge ",& - ! igauge(ii) + ! igauge(i) ! write(6,*) " xoff,yoff,zoff: ", xoff,yoff,zoff !endif ! Bilinear interpolation var_index = 0 - do n = 1, size(gauges(ii)%q_out_vars, 1) - if (gauges(ii)%q_out_vars(n)) then + do n = 1, size(gauges(i)%q_out_vars, 1) + if (gauges(i)%q_out_vars(n)) then var_index = var_index + 1 var(var_index) = (1.d0 - xoff) * (1.d0 - yoff) * q(n, iindex, jindex, kindex) & + xoff * (1.d0 - yoff) * q(n, iindex + 1, jindex, kindex) & @@ -462,9 +419,9 @@ subroutine update_gauges(q, aux, xlow, ylow, zlow, num_eqn, mitot, mjtot, & end if end do - if (allocated(gauges(ii)%aux_out_vars)) then - do n = 1, size(gauges(ii)%aux_out_vars, 1) - if (gauges(ii)%aux_out_vars(n)) then + if (allocated(gauges(i)%aux_out_vars)) then + do n = 1, size(gauges(i)%aux_out_vars, 1) + if (gauges(i)%aux_out_vars(n)) then var_index = var_index + 1 var(var_index) = (1.d0 - xoff) * (1.d0 - yoff) * q(n, iindex, jindex, kindex) & + xoff * (1.d0 - yoff) * q(n, iindex + 1, jindex, kindex) & @@ -481,28 +438,28 @@ subroutine update_gauges(q, aux, xlow, ylow, zlow, num_eqn, mitot, mjtot, & end if ! Check to make sure we grabbed all the values - if (gauges(ii)%num_out_vars /= var_index) then - print *, gauges(ii)%num_out_vars, var_index - print *, gauges(ii)%q_out_vars - print *, gauges(ii)%aux_out_vars + if (gauges(i)%num_out_vars /= var_index) then + print *, gauges(i)%num_out_vars, var_index + print *, gauges(i)%q_out_vars + print *, gauges(i)%aux_out_vars stop "Somehow we did not grab all the values we wanted..." end if ! save info for this time - index = gauges(ii)%buffer_index + index = gauges(i)%buffer_index - gauges(ii)%level(index) = level - gauges(ii)%data(1,index) = tgrid - do j = 1, gauges(ii)%num_out_vars - gauges(ii)%data(1 + j, index) = var(j) + gauges(i)%level(index) = level + gauges(i)%data(1,index) = tgrid + do j = 1, gauges(i)%num_out_vars + gauges(i)%data(1 + j, index) = var(j) end do - gauges(ii)%buffer_index = index + 1 - if (gauges(ii)%buffer_index > MAX_BUFFER) then - call print_gauges_and_reset_nextLoc(ii) + gauges(i)%buffer_index = index + 1 + if (gauges(i)%buffer_index > MAX_BUFFER) then + call print_gauges_and_reset_nextLoc(i) endif - gauges(ii)%last_time = tgrid + gauges(i)%last_time = tgrid end do ! End of gauge loop ============================================= diff --git a/src/3d/init_bndryList.f b/src/3d/init_bndryList.f new file mode 100644 index 000000000..24685fbc8 --- /dev/null +++ b/src/3d/init_bndryList.f @@ -0,0 +1,135 @@ +c +c ----------------------------------------------------------------- +c + subroutine initBndryList() + + use amr_module + implicit none + + integer i +c +c need to manage the boundary List too +c + bndListSize = 8*maxgr ! guess at initial size + + ! allocate the bndList + if (.not. allocated(bndList)) then + allocate(bndList(bndListSize,2)) + print *, "bndList allocated..." + else + print *, "bndList already allocated!" + endif + + ! thread the bndry list + do i = 1, bndListSize + bndList(i,nextfree) = i+1 + end do + + bndList(bndListSize,nextfree) = null + ndfree_bnd = 1 + + end +c +c ----------------------------------------------------------------- +c +!> Preprocess each grid on level **level** to have a linked list of +!! other grids at the same level that supply ghost cells. +!! +!! The linked list is implemented with an array. +!! Each node in this linked list is represented by a row (with 2 elements) +!! in the array, named **bndList**. +!! +!! bndList(pos,gridNbor) is grid number stored at node **pos** +!! bndList(pos,nextfree) is pointer to next node in the linked list +!! +!! node(bndListSt,mptr) point to the head node in this linked list. +!! Thus bndList(node(bndListSt,mptr),gridNbor) is grid number +!! stored at first node of the linked list. +!! +!! + subroutine makeBndryList(level) +c + use amr_module + implicit none + + integer level, n, levSt, k, nborCount + integer nodget_bnd, nextSpot, prevNbor, msrc, mptr + integer imin, imax, jmin, jmax + integer imlo, imhi, jmlo, jmhi + integer ixlo, ixhi, jxlo, jxhi + +c :::::::::::::::::::::::::::: makeBndryList ::::::::::::::::::::::::: +c preprocess each grid to have linked list of other grids at +c same level that supply ghost cells. +c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + +c traverse linked list into array. list already sorted by arrangegrids + levSt = listStart(level) + do n = 1, numgrids(level) + mptr = listOfGrids(levSt+n-1) + imin = node(ndilo,mptr) - nghost ! ghost cells included since + imax = node(ndihi,mptr) + nghost ! this is what you want to fill + jmin = node(ndjlo,mptr) - nghost ! may also use for filval stuff, + jmax = node(ndjhi,mptr) + nghost ! change nghost to mbuff, etc + nborCount = 0 + + do k = 1, numgrids(level) ! loop over all other grids once to find touching ones + if (k .eq. n) cycle ! dont count yourself as source grid + msrc = listOfgrids(levSt+k-1) + + ! Check if grid mptr and patch intersect + imlo = node(ndilo, msrc) + jmlo = node(ndjlo, msrc) + imhi = node(ndihi, msrc) + jmhi = node(ndjhi, msrc) + + ixlo = max(imlo,imin) + ixhi = min(imhi,imax) + jxlo = max(jmlo,jmin) + jxhi = min(jmhi,jmax) + + if (ixlo .le. ixhi .and. jxlo .le. jxhi) then ! put on bnd list for mptr + nborCount = nborCount + 1 + nextSpot = nodget_bnd() + bndList(nextSpot,gridNbor) = msrc + ! get spot in bnd list. insert next grid at front to avoid traversing + bndList(nextSpot,nextfree) = node(bndListSt,mptr) + node(bndListSt,mptr) = nextSpot + endif + + end do + +! save final count + node(bndListNum,mptr) = nborcount + end do + + return + end +c +c ----------------------------------------------------------------- +c +!> Free the linked list of intersecting "boundary" grids for grid 'mold' +!! that is no longer active. +!! The linked list starts at node(bndListSt, mold). + subroutine freeBndryList(mold) +c + use amr_module + implicit none + + integer nborCount, mold,nextSpot, i, nextnext + +c :::::::::::::::::::::::::::: freeBndryList ::::::::::::::::::::::::: +c free the linked list of intersecting "boundary" grids for grid 'mold' +c that is no longer active +c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + nborCount = node(bndListNum,mold) ! count for this grid + nextSpot = node(bndListSt,mold) ! first index of this grids nbors + do i = 1, nborCount + nextnext = bndList(nextSpot,nextfree) + call putnod_bnd(nextSpot) + nextSpot = nextnext + end do + + return + end diff --git a/src/3d/init_nodes.f90 b/src/3d/init_nodes.f90 new file mode 100644 index 000000000..b8fad4839 --- /dev/null +++ b/src/3d/init_nodes.f90 @@ -0,0 +1,38 @@ +! ============================================================================ +! Program: AMRClaw +! File: init_nodes.f90 +! Created: 2018-06-15 +! Author: Marsha Berger +! ============================================================================ +! Description: Initialization of rnode and node storage +! ============================================================================ + + +subroutine init_nodes() + + use amr_module + implicit none + + maxgr = 10000 + if (.not.allocated(rnode)) then ! new way, use allocatable arrays, not pointers + allocate(rnode(rsize,maxgr)) + print *, "rnode allocated..." + else + print *, "rnode already allocated!" + endif + if (.not.allocated(node)) then ! new way, use allocatable arrays, not pointers + allocate(node(nsize,maxgr)) + print *, "node allocated..." + else + print *, "rnode already allocated!" + endif + if (.not.allocated(listOfGrids)) then ! include all vars whose size depends on maxgr + allocate(listOfGrids(maxgr)) + print *, "listOfGrids allocated..." + else + print *, "listOfGrids already allocated!" + endif + + +end subroutine init_nodes + diff --git a/src/3d/nodget.f b/src/3d/nodget.f index 34d9de3ec..e5ed12d11 100644 --- a/src/3d/nodget.f +++ b/src/3d/nodget.f @@ -1,10 +1,17 @@ c c ------------------------------------------------------------ c +!> Get first free node of the linked list kept in node +!! array. adjust pointers accordingly. +!! +!! This function is used to create a new grid descriptor, like +!! mptr = new grid_class in c++. +!! integer function nodget() c use amr_module implicit double precision (a-h,o-z) + integer maxgrIncrement/10000/ c c ::::::::::::::::: NODGET ::::::::::::::::::::::::::::::::::::; @@ -18,13 +25,22 @@ integer function nodget() 100 format(' out of nodal space - allowed ',i8,' grids') do level = 1, lfine write(*,101) level,numgrids(level) - 101 format(" level ",i4," has ",i6,'grids') + 101 format(" level ",i4," has ",i7,' grids') end do write(*,*)" Could need twice as many grids as on any given" write(*,*)" level if regridding/birecting" - stop -c -c update pointers +! stop +! new way gets more storage and continues + istatus = 1 + call resize_nodes(maxgr + maxgrIncrement, istatus) + if (istatus > 0) then + write(*,102) maxgr,maxgrIncrement +102 format(' unable to increase nodal space by ',i8,' grids',/, + . 'resize failed') + stop + endif +c +c allocate next node, update pointers c 10 nodget = ndfree ndfree = node(nextfree,ndfree) @@ -72,9 +88,10 @@ integer function nodget_bnd() 101 format(" There are ",i8," total grids", i10," bndry nbors", . " average num/grid ",f10.3) - stop + !stop + call resize_bndryList() c -c ## adjust pointers +c ## adjust pointers for next bndry entry c 10 nodget_bnd = ndfree_bnd ndfree_bnd = bndList(ndfree_bnd,nextfree) @@ -89,6 +106,16 @@ integer function nodget_bnd() c c ----------------------------------------------------------------- c +!> Make (modify) array of grid numbers, listOfGrids, (after sorting them +!! in the linked list so they are in decreasing order of workload, +!! done in arrangeGrid()). +!! +!! This is done every time there is regridding, initial gridding, +!! or restarting. Most often finest level is regridded, so +!! put it last in array. **lbase** is the level that didnt change, so +!! only redo from lbase+1 to lfine. +!! \param[in] lbase all levels from **lbase**+1 to the finest get +!! modifed in the array, listOfGrids subroutine makeGridList(lbase) c use amr_module @@ -123,118 +150,5 @@ subroutine makeGridList(lbase) listStart(lev+1) = levSt + numgrids(lev) end do - return - end -c -c ----------------------------------------------------------------- -c - subroutine initBndryList() - - use amr_module - implicit none - - integer i -c -c need to manage the boundary List too -c - do i = 1, bndListSize - bndList(i,nextfree) = i+1 - end do - - bndList(bndListSize,nextfree) = null - ndfree_bnd = 1 - - end -c -c ----------------------------------------------------------------- -c - subroutine makeBndryList(level) -c - use amr_module - implicit none - - integer level, n, levSt, k, nborCount - integer nodget_bnd, nextSpot, prevNbor, msrc, mptr - integer imin, imax, jmin, jmax, kmin, kmax - integer imlo, imhi, jmlo, jmhi, kmlo, kmhi - integer ixlo, ixhi, jxlo, jxhi, kxlo, kxhi - -c :::::::::::::::::::::::::::: makeBndryList ::::::::::::::::::::::::: -c preprocess each grid to have linked list of other grids at -c same level that supply ghost cells. -c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - -c traverse linked list into array. list already sorted by arrangegrids - levSt = listStart(level) - do n = 1, numgrids(level) - mptr = listOfGrids(levSt+n-1) - imin = node(ndilo,mptr) - nghost ! ghost cells included since - imax = node(ndihi,mptr) + nghost ! this is what you want to fill - jmin = node(ndjlo,mptr) - nghost ! may also use for filval stuff, - jmax = node(ndjhi,mptr) + nghost ! change nghost to mbuff, etc - kmin = node(ndklo,mptr) - nghost ! may also use for filval stuff, - kmax = node(ndkhi,mptr) + nghost ! change nghost to mbuff, etc - nborCount = 0 - - do k = 1, numgrids(level) ! loop over all other grids once to find touching ones - if (k .eq. n) cycle ! dont count yourself as source grid - msrc = listOfgrids(levSt+k-1) - - ! Check if grid mptr and patch intersect - imlo = node(ndilo, msrc) - jmlo = node(ndjlo, msrc) - kmlo = node(ndklo, msrc) - imhi = node(ndihi, msrc) - jmhi = node(ndjhi, msrc) - kmhi = node(ndkhi, msrc) - - ixlo = max(imlo,imin) - ixhi = min(imhi,imax) - jxlo = max(jmlo,jmin) - jxhi = min(jmhi,jmax) - kxlo = max(kmlo,kmin) - kxhi = min(kmhi,kmax) - - if (ixlo .le. ixhi .and. jxlo .le. jxhi - . .and. kxlo .le. kxhi) then ! put on bnd list for mptr - nborCount = nborCount + 1 - nextSpot = nodget_bnd() - bndList(nextSpot,gridNbor) = msrc - ! get spot in bnd list. insert next grid at front to avoid traversing - bndList(nextSpot,nextfree) = node(bndListSt,mptr) - node(bndListSt,mptr) = nextSpot - endif - - end do - -! save final count - node(bndListNum,mptr) = nborcount - end do - - return - end -c -c ----------------------------------------------------------------- -c - subroutine freeBndryList(mold) -c - use amr_module - implicit none - - integer nborCount, mold,nextSpot, i, nextnext - -c :::::::::::::::::::::::::::: freeBndryList ::::::::::::::::::::::::: -c free the linked list of intersecting "boundary" grids for grid 'mold' -c that is no longer active -c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - - nborCount = node(bndListNum,mold) ! count for this grid - nextSpot = node(bndListSt,mold) ! first index of this grids nbors - do i = 1, nborCount - nextnext = bndList(nextSpot,nextfree) - call putnod_bnd(nextSpot) - nextSpot = nextnext - end do - return end diff --git a/src/3d/resize_bndryList.f b/src/3d/resize_bndryList.f new file mode 100644 index 000000000..94e9f7b09 --- /dev/null +++ b/src/3d/resize_bndryList.f @@ -0,0 +1,29 @@ + + subroutine resize_bndryList() + + use amr_module + implicit double precision (a-h,o-z) + integer, allocatable, target, dimension(:,:) :: new_bndList + + ! get new bndry space + new_bndListSize = 1.5*bndListSize + print *,"Expanding size of boundary list from ",bndListSize, + . " to ",new_bndListSize + allocate(new_bndList(new_bndListSize,2),STAT=istatus) + if (istatus > 0) then + write(*,*)" could not get new bndry list space" + stop + endif + new_bndList(1:bndListSize,1:2) = bndList + call move_alloc(new_bndList,bndList) + + ! thread new bndry list space + do i = bndListSize+1, new_bndListSize + bndList(i,nextfree) = i+1 + end do + + bndList(new_bndListSize,nextfree) = null + ndfree_bnd = bndListSize+1 + bndListSize = new_bndListSize + + end diff --git a/src/3d/resize_nodes.f90 b/src/3d/resize_nodes.f90 new file mode 100644 index 000000000..ca6576100 --- /dev/null +++ b/src/3d/resize_nodes.f90 @@ -0,0 +1,65 @@ +! ============================================================================ +! Program: AMRClaw +! File: resize_nodes.f90 +! Author: Marsha Berger +! ============================================================================ +! Description: Resize the node space if have run out +! ============================================================================ + +subroutine resize_nodes(new_size,status) + + use amr_module + implicit none + + integer, intent(out) :: status + integer, intent(in) :: new_size + + integer :: i + + real(kind=8), allocatable, target, dimension(:,:) :: new_rnode + integer, allocatable, target, dimension(:,:) :: new_node + integer, allocatable, target, dimension(:) :: new_listOfGrids + + print *, "Expanding maximum number of grids from ", maxgr," to ", new_size + + ! first for rnode + allocate(new_rnode(rsize,new_size),STAT=status) + if (status > 0) then + return + endif + new_rnode(1:rsize,1:maxgr) = rnode ! new way, use allocatable, not pointer + + call move_alloc(new_rnode,rnode) + + ! next for node + allocate(new_node(nsize,new_size),STAT=status) + if (status > 0) then + return + endif + new_node(1:nsize,1:maxgr) = node ! new way, use allocatable, not pointer + + call move_alloc(new_node,node) + + !! need to rethread new space to be able to use it when new grids requested + do i = maxgr+1, new_size + node(nextfree,i) = i+1 + end do + ! reset last one to null + node(nextfree, new_size) = null + + ! next for listOfGrids + allocate(new_listOfGrids(new_size),STAT=status) + if (status > 0) then + return + endif + new_listOfGrids(1:maxgr) = listOfGrids ! new way, use allocatable, not pointer + + call move_alloc(new_listOfGrids,listOfGrids) + + ! reset maxgr and next free node, to continue + ndfree = maxgr + 1 + maxgr = new_size + + return + +end subroutine resize_nodes diff --git a/src/3d/restrt.f b/src/3d/restrt.f index 6810fb924..9f18f19e1 100644 --- a/src/3d/restrt.f +++ b/src/3d/restrt.f @@ -39,10 +39,12 @@ subroutine restrt(nsteps,time,nvar,naux) open(rstunit,file=trim(rstfile),status='old',form='unformatted') rewind rstunit - read(rstunit) lenmax,lendim,isize + !! new version has flexible node size, so need to read current size maxgr + read(rstunit) lenmax,lendim,isize,maxgr c # need to allocate for dynamic memory: call restrt_alloc(isize) + call restrt_nodes(maxgr) read(rstunit) (alloc(i),i=1,lendim) read(rstunit) hxposs,hyposs,hzposs,possk,icheck @@ -207,8 +209,8 @@ subroutine restrt(nsteps,time,nvar,naux) 15 if (mptr .eq. 0) go to 25 mitot = node(ndihi,mptr)-node(ndilo,mptr)+1+2*nghost mjtot = node(ndjhi,mptr)-node(ndjlo,mptr)+1+2*nghost - mktot = node(ndkhi,mptr)-node(ndklo,mptr)+1+2*nghost - node(store2,mptr) = igetsp(mitot*mjtot*mktot*nvar) + mktot = node(ndkhi,mptr)-node(ndklo,mptr)+1+2*nghost + node(store2,mptr) = igetsp(mitot*mjtot*mktot*nvar) mptr = node(levelptr,mptr) go to 15 25 continue @@ -246,5 +248,6 @@ subroutine restrt(nsteps,time,nvar,naux) call makeBndryList(level) end do c + call initTimers() return end diff --git a/src/3d/restrt_nodes.f90 b/src/3d/restrt_nodes.f90 new file mode 100644 index 000000000..8bf3779e5 --- /dev/null +++ b/src/3d/restrt_nodes.f90 @@ -0,0 +1,45 @@ +! ============================================================================ +! Program: AMRClaw +! File: restrt_nodes.f90 +! Created: 2018-06-15 +! Author: Marsha Berger +! ============================================================================ +! Description: Initialization of nodal storage for restart +! ============================================================================ + + +subroutine restrt_nodes(isize) + + use amr_module + implicit none + integer :: isize + + maxgr = isize + + if (.not.allocated(rnode)) then ! allocatable nodal arrays now + write(6,*)"allocating ",maxgr," -sized rnode array" + allocate(rnode(rsize,maxgr)) + print *, "rnode storage allocated of size ",maxgr," at restart" + else + print *, "rnode storage already allocated!" + endif + + if (.not.allocated(node)) then ! allocatable nodal arrays now + write(6,*)"allocating ",maxgr," -sized node array" + allocate(node(nsize,maxgr)) + print *, "node storage allocated of size ",maxgr," at restart" + else + print *, "node storage already allocated!" + endif + + if (.not.allocated(listOfGrids)) then ! allocatable nodal arrays now + write(6,*)"allocating ",maxgr," -sized listOfGrids " + allocate(listOfGrids(maxgr)) + print *, "listOfGrids allocated of size ",maxgr," at restart" + else + print *, "listOfGrids already allocated!" + endif + +end subroutine restrt_nodes + + diff --git a/src/3d/stst1.f b/src/3d/stst1.f index 3a4bbd146..87600d095 100644 --- a/src/3d/stst1.f +++ b/src/3d/stst1.f @@ -21,6 +21,11 @@ subroutine stst1 c ::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::::::::: c ndfree = 1 + + !! now rnode and node are allocatable to allow resizing + call init_nodes() + + ! node space allocated, now thread it do 10 i = 1, maxgr node(nextfree,i) = i+1 10 continue