Skip to content

Commit

Permalink
Merge pull request #222 from mjberger/flexMax_3d
Browse files Browse the repository at this point in the history 
Flex max 3d
  • Loading branch information
@mandli
mandli authored Jun 17, 2018
2 parents 8329ac1 + 1df7f72 commit fc882fd
Show file tree
Hide file tree
Showing 13 changed files with 453 additions and 242 deletions.
5 changes: 5 additions & 0 deletions src/3d/Makefile.amr_3d
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 \


27 changes: 20 additions & 7 deletions src/3d/amr_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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), &
Expand Down
4 changes: 3 additions & 1 deletion src/3d/check.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/3d/domain.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
179 changes: 68 additions & 111 deletions src/3d/gauges_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

!
! -------------------------------------------------------------------------
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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...
Expand All @@ -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) &
Expand All @@ -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) &
Expand All @@ -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 =============================================

Expand Down
Loading

0 comments on commit fc882fd

Please sign in to comment.