2D AMRCLAW
Functions/Subroutines
filpatch.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

recursive subroutine filrecur (level, nvar, valbig, aux, naux, t, mitot, mjtot, nrowst, ncolst, ilo, ihi, jlo, jhi, patchOnly, msrc)
 Fill a region (patch) described by: More...
 
integer pure function coarse_index (n, i, j)
 
logical pure function sticksout (iplo, iphi, jplo, jphi)
 

Function/Subroutine Documentation

◆ coarse_index()

integer pure function filrecur::coarse_index ( integer, intent(in)  n,
integer, intent(in)  i,
integer, intent(in)  j 
)

Definition at line 343 of file filpatch.f90.

Referenced by coarse_index(), and filrecur().

343  implicit none
344  integer, intent(in) :: n, i, j
345  coarse_index = n + nvar*(i-1)+nvar*mitot_coarse*(j-1)
integer pure function coarse_index(n, i, j)
Definition: filpatch.f90:343
Here is the caller graph for this function:

◆ filrecur()

recursive subroutine filrecur ( integer, intent(in)  level,
integer, intent(in)  nvar,
real(kind=8), dimension(nvar,mitot,mjtot), intent(inout)  valbig,
real(kind=8), dimension(naux,mitot,mjtot), intent(inout)  aux,
integer, intent(in)  naux,
real(kind=8), intent(in)  t,
integer, intent(in)  mitot,
integer, intent(in)  mjtot,
integer, intent(in)  nrowst,
integer, intent(in)  ncolst,
integer, intent(in)  ilo,
integer, intent(in)  ihi,
integer, intent(in)  jlo,
integer, intent(in)  jhi,
logical  patchOnly,
integer, intent(in)  msrc 
)

Fill a region (patch) described by:

  • global indices of its lower left corner: (ilo, jlo).
  • global indices of its upper right corner: (ihi, jhi).

The patch is on grid msrc and starts from row nrowst and column ncolst of the grid.

It first fills the patch with values obtainable from level level grids. if any left unfilled, then enlarge remaining rectangle of unfilled values by 1 (for later linear interp), and recusively obtain the remaining values from coarser levels.

## Algorithm

procedure filrecur(level, patch_fine, val_fine, firstcall)
fill patch_fine as much as possible by copy values from other level "level" grids
values on patch_fine is stored in array val_fine after the filling
if this is first call to filrecur() (firstcall == true)
don't fill cells outside the computational domain (it will be handled elsewhere)
else
this is a sencond or higher-level (recursive) call to filrecur()
fill cells outside the computational domain by calling bc2amr since they are used for interpolation for finer cells above them
end if
if not all cells in patch_fine is filled (set)
find the smallest rectangular patch on level "level-1", patch_coarse, that contains all unset fine cells in patch_fine
filrecur(level-1, patch_coarse, val_coarse, firstCall=false)
for each unset fine cells in patch_fine
interpolate from val_coarse to fill all unset cells on patch_fine
end for
Parameters
levelAMR level the patch is on
nvarnumber of equations for the system
valbigdata array for solution $q $ (cover the whole grid msrc)
auxdata array for auxiliary variables
nauxnumber of auxiliary variables
tfill the patch with value at time t
mitotnumber of cells in i direction on grid msrc
mjtotnumber of cells in j direction on grid msrc
nrowstlocal i index (relative to lower-left corner of grid msrc) of the left-most cell of the patch to be filled
ncolstlocal j index (relative to lower-left corner of grid msrc) of the lower-most cell of the patch to be filled
iloglobal i index of left-most cell of this patch
ihiglobal i index of right-most cell of this patch
jloglobal j index of lower-most cell of this patch
jhiglobal j index of upper-most cell of this patch
patchOnlyThis is 1) false if this is a first level call to filrecur() and won't set bounday cells outside domain; 2) true if this is a second or higher level (recursive) call to filrecur() and will call bc2amr() to get values for cells outside domain (these cells are only for interpolation to fill cells in first level call to filrecur()
msrcindex of the grid

Definition at line 73 of file filpatch.f90.

References bc2amr(), coarse_index(), amr_module::hxposs, amr_module::hyposs, intfil(), amr_module::intratx, amr_module::intraty, amr_module::iregsz, amr_module::jregsz, amr_module::needs_to_be_set, amr_module::nghost, amr_module::outunit, prefilrecur(), setaux(), amr_module::spheredom, sticksout(), trimbd(), amr_module::xlower, amr_module::xperdom, amr_module::xupper, amr_module::ylower, amr_module::yperdom, and amr_module::yupper.

Referenced by bound(), filrecur(), and prefilrecur().

73 
77 
78  implicit none
79 
80  ! Input
81  integer, intent(in) :: level, nvar, naux, mitot, mjtot, nrowst, ncolst
82  integer, intent(in) :: ilo,ihi,jlo,jhi, msrc
83  real(kind=8), intent(in) :: t
84  logical :: patchonly
85 
86  ! Output
87  real(kind=8), intent(in out) :: valbig(nvar,mitot,mjtot)
88  real(kind=8), intent(in out) :: aux(naux,mitot,mjtot)
89 
90  ! Local storage
91  integer :: i_fine, j_fine, i_coarse, j_coarse, n, k
92  integer :: iplo, iphi, jplo, jphi
93  integer :: mitot_patch, mjtot_patch, mitot_coarse, mjtot_coarse, nghost_patch, lencrse
94  integer :: refinement_ratio_x, refinement_ratio_y
95  integer :: unset_indices(4)
96  real(kind=8) :: dx_fine, dy_fine, dx_coarse, dy_coarse
97  real(kind=8) :: xlow_coarse,ylow_coarse, xlow_fine, ylow_fine, xhi_fine,yhi_fine
98  real(kind=8) :: xcent_fine, xcent_coarse, ycent_fine, ycent_coarse,ratiox,ratioy,floor
99 
100  !for timing
101  integer :: clock_start, clock_finish, clock_rate
102  real(kind=8) :: cpu_start, cpu_finish
103 
104  ! Interpolation variables
105  real(kind=8) :: eta1, eta2, valp10, valm10, valc, valp01, valm01, dupc, dumc
106  real(kind=8) :: ducc, du, fu, dvpc, dvmc, dvcc, dv, fv, valint, uslope, vslope
107 
108  ! Cell set tracking
109  logical :: set
110  integer(kind=1) :: flaguse(ihi-ilo+1, jhi-jlo+1)
111 
112  ! Scratch storage
113  ! use stack-based scratch arrays instead of alloc, since dont really
114  ! need to save beyond these routines, and to allow dynamic memory resizing
115  !
116  ! use 1d scratch arrays that are potentially the same size as
117  ! current grid, since may not coarsen.
118  ! need to make it 1d instead of 2 and do own indexing, since
119  ! when pass it in to subroutines they treat it as having dierent
120  ! dimensions than the max size need to allocate here
121 
122  !-- dimension valcrse((ihi-ilo+2)*(jhi-jlo+2)*nvar) ! NB this is a 1D array
123  !-- dimension auxcrse((ihi-ilo+2)*(jhi-jlo+2)*naux) ! the +2 is to expand on coarse grid to enclose fine
124  ! ### turns out you need 3 rows, forget offset of 1 plus one on each side
125  ! the +3 is to expand on coarse grid to enclose fine
126  real(kind=8) :: valcrse((ihi-ilo+3) * (jhi-jlo+3) * nvar)
127  real(kind=8) :: auxcrse((ihi-ilo+3) * (jhi-jlo+3) * naux)
128 
129  ! We begin by filling values for grids at level level. If all values can be
130  ! filled in this way, we return;
131  mitot_patch = ihi-ilo + 1 ! nrowp
132  mjtot_patch = jhi-jlo + 1
133 
134  dx_fine = hxposs(level)
135  dy_fine = hyposs(level)
136 
137  ! Coordinates of edges of patch (xlp,xrp,ybp,ytp)
138  xlow_fine = xlower + ilo * dx_fine
139  ylow_fine = ylower + jlo * dy_fine
140 
141  ! Fill in the patch as much as possible using values at this level
142  ! note that if only a patch, msrc = -1, otherwise a real grid and intfil
143  ! uses its boundary list
144  ! msrc either -1 (for a patch) or the real grid number
145  call intfil(valbig,mitot,mjtot,t,flaguse,nrowst,ncolst, ilo, &
146  ihi,jlo,jhi,level,nvar,naux,msrc)
147 
148 
149 
150  ! Trimbd returns set = true if all of the entries are filled (=1.).
151  ! set = false, otherwise. If set = true, then no other levels are
152  ! are required to interpolate, and we return.
153  !
154  ! Note that the used array is filled entirely in intfil, i.e. the
155  ! marking done there also takes into account the points filled by
156  ! the boundary conditions. bc2amr will be called later, after all 4
157  ! boundary pieces filled.
158  call trimbd(flaguse,mitot_patch,mjtot_patch,set,unset_indices)
159  ! il,ir,jb,jt = unset_indices(4)
160 
161  ! If set is .true. then all cells have been set and we can skip to setting
162  ! the remaining boundary cells. If it is .false. we need to interpolate
163  ! some values from coarser levels, possibly calling this routine
164  ! recursively.
165  if (.not.set) then
166 
167  ! Error check
168  if (level == 1) then
169  write(outunit,*)" error in filrecur - level 1 not set"
170  write(outunit,'("start at row: ",i4," col ",i4)') nrowst,ncolst
171  print *," error in filrecur - level 1 not set"
172  print *," should not need more recursion "
173  print *," to set patch boundaries"
174  print '("start at row: ",i4," col ",i4)', nrowst,ncolst
175  stop
176  endif
177 
178  ! We begin by initializing the level level arrays, so that we can use
179  ! purely recursive formulation for interpolating.
180  dx_coarse = hxposs(level - 1)
181  dy_coarse = hyposs(level - 1)
182 
183  ! Adjust unset_indices to account for the patch
184  ! isl, isr, jsb, jst
185  unset_indices(1) = unset_indices(1) + ilo - 1
186  unset_indices(2) = unset_indices(2) + ilo - 1
187  unset_indices(3) = unset_indices(3) + jlo - 1
188  unset_indices(4) = unset_indices(4) + jlo - 1
189 
190  ! Coarsened geometry
191  refinement_ratio_x = intratx(level - 1)
192  refinement_ratio_y = intraty(level - 1)
193 
194  ! New patch rectangle (after we have partially filled it in) but in the
195  ! coarse patches [iplo,iphi,jplo,jphi]
196 
197  ! islo = unset_indices(1)
198  ! ishi = unset_indices(2)
199  ! islo ishi
200  ! |----|----|----|----|----|----|----|----|
201  ! If following operation is conducted, it will lead to:
202  !
203  ! iplo iphi
204  ! |-------------------|-------------------|-------------------|-------------------|
205  iplo = (unset_indices(1) - refinement_ratio_x + nghost * refinement_ratio_x) &
206  / refinement_ratio_x - nghost
207  iphi = (unset_indices(2) + refinement_ratio_x) / refinement_ratio_x
208  jplo = (unset_indices(3) - refinement_ratio_y + nghost * refinement_ratio_y) &
209  / refinement_ratio_y - nghost
210  jphi = (unset_indices(4) + refinement_ratio_y) / refinement_ratio_y
211 
212  xlow_coarse = xlower + iplo * dx_coarse
213  ylow_coarse = ylower + jplo * dy_coarse
214 
215  ! Coarse grid number of spatial points (nrowc,ncolc)
216  mitot_coarse = iphi - iplo + 1
217  mjtot_coarse = jphi - jplo + 1
218 
219  ! Check to make sure we created big enough scratch arrays
220  if (mitot_coarse > ihi - ilo + 3 .or. &
221  mjtot_coarse > jhi - jlo + 3) then
222 
223  print *," did not make big enough work space in filrecur "
224  print *," need coarse space with mitot_coarse,mjtot_coarse ",mitot_coarse,mjtot_coarse
225  print *," made space for ilo,ihi,jlo,jhi ",ilo,ihi,jlo,jhi
226  stop
227  endif
228 
229  ! Set the aux array values for the coarse grid, this could be done
230  ! instead in intfil using possibly already available bathy data from the
231  ! grids
232  if (naux > 0) then
233  nghost_patch = 0
234  lencrse = (ihi-ilo+3)*(jhi-jlo+3)*naux ! set 1 component, not all naux
235  do k = 1, lencrse, naux
236  auxcrse(k) = needs_to_be_set ! new system checks initialization before setting aux vals
237  end do
238  call setaux(nghost_patch, mitot_coarse,mjtot_coarse, &
239  xlow_coarse, ylow_coarse, &
240  dx_coarse,dy_coarse,naux,auxcrse)
241  endif
242 
243  ! Fill in the edges of the coarse grid. for recursive calls, patch indices and
244  ! 'coarse grid' indices are the same (no actual coarse grid here, so cant use mptr
245  ! must pass indices. patchOnly argument is thus true
246  if ((xperdom .or. (yperdom .or. spheredom)) .and. sticksout(iplo,iphi,jplo,jphi)) then
247  call prefilrecur(level-1,nvar,valcrse,auxcrse,naux,t,mitot_coarse,mjtot_coarse,1,1, &
248  iplo,iphi,jplo,jphi,iplo,iphi,jplo,jphi,.true.)
249  else
250  call filrecur(level-1,nvar,valcrse,auxcrse,naux,t,mitot_coarse,mjtot_coarse,1,1, &
251  iplo,iphi,jplo,jphi,.true.,-1) ! when going to coarser patch, no source grid (for now at least)
252  endif
253 
254  ! convert to real for use below
255  ratiox = real(refinement_ratio_x,kind=8)
256  ratioy = real(refinement_ratio_y,kind=8)
257 
258  ! fine below refers to level "level"
259  ! coarse below refers to level "level-1"
260  do j_fine = 1,mjtot_patch
261  !j_coarse = 2 + (j_fine - (unset_indices(3) - jlo) - 1) / refinement_ratio_y
262  !eta2 = (-0.5d0 + real(mod(j_fine - 1, refinement_ratio_y),kind=8)) &
263  ! / real(refinement_ratio_y,kind=8)
264 
265  ! we are now processing the j_fine^{th} column of the fine patch
266  ! this is contained in the j_coarse^{th} column of the coarse patch below it
267  ! note that these two are not global indices but local indices
268  j_coarse =floor((j_fine+jlo-1)/ratioy) - jplo + 1
269  ycent_coarse = ylow_coarse + (j_coarse-.5d0)*dy_coarse
270  ycent_fine = ylower + (j_fine-1+jlo + .5d0)*dy_fine
271  eta2 = (ycent_fine-ycent_coarse)/dy_coarse
272  if (abs(eta2) .gt. .5) then
273  write(*,*)" filpatch y indices wrong: eta2 = ",eta2
274  endif
275 
276  do i_fine = 1,mitot_patch
277  !i_coarse = 2 + (i_fine - (unset_indices(1) - ilo) - 1) / refinement_ratio_x
278  !eta1 = (-0.5d0 + real(mod(i_fine - 1, refinement_ratio_x),kind=8)) &
279  ! / real(refinement_ratio_x,kind=8)
280 
281  ! new coarse indexing
282  !i_coarse =(i_fine+ilo-1)/refinement_ratio_x - iplo + 1
283  i_coarse = floor((i_fine+ilo-1)/ratiox) - iplo + 1
284  xcent_coarse = xlow_coarse + (i_coarse-.5d0)*dx_coarse
285  xcent_fine = xlower + (i_fine-1+ilo + .5d0)*dx_fine
286  eta1 = (xcent_fine-xcent_coarse)/dx_coarse
287  if (abs(eta1) .gt. .5) then
288  write(*,*)" filpatch x indices wrong: eta1 = ",eta1
289  endif
290 
291 
292  if (flaguse(i_fine,j_fine) == 0) then
293 
294  do n=1,nvar
295  ! write(*,*)n,i_coarse+1,j_coarse,coarse_index(n,i_coarse + 1,j_coarse)
296  ! QUESTION: why do we interpolate in this way?
297  valp10 = valcrse(coarse_index(n,i_coarse + 1,j_coarse))
298  valm10 = valcrse(coarse_index(n,i_coarse - 1,j_coarse))
299  valc = valcrse(coarse_index(n,i_coarse ,j_coarse))
300  valp01 = valcrse(coarse_index(n,i_coarse ,j_coarse + 1))
301  valm01 = valcrse(coarse_index(n,i_coarse ,j_coarse - 1))
302 
303  dupc = valp10 - valc
304  dumc = valc - valm10
305  ducc = valp10 - valm10
306  du = min(abs(dupc), abs(dumc))
307  du = min(2.d0 * du, 0.5d0 * abs(ducc))
308  fu = max(0.d0, sign(1.d0, dupc * dumc))
309  uslope = du*sign(1.d0,ducc)*fu ! not really - should divide by h
310 
311  dvpc = valp01 - valc
312  dvmc = valc - valm01
313  dvcc = valp01 - valm01
314  dv = min(abs(dvpc), abs(dvmc))
315  dv = min(2.d0 * dv, 0.5d0 * abs(dvcc))
316  fv = max(0.d0,sign(1.d0, dvpc * dvmc))
317  vslope = dv*sign(1.d0,dvcc)*fv ! but faster to put in eta above
318 
319  valint = valc + eta1 * uslope + eta2 * vslope
320 
321  valbig(n,i_fine+nrowst-1,j_fine+ncolst-1) = valint
322 
323  end do
324 
325  endif
326  end do
327  end do
328  end if
329 
330  ! set bcs, whether or not recursive calls needed. set any part of patch that stuck out
331  xhi_fine = xlower + (ihi + 1) * dx_fine
332  yhi_fine = ylower + (jhi + 1) * dy_fine
333  ! only call if a small coarser recursive patch
334  ! otherwise whole grid bcs done from bound
335  if (patchonly) then
336  call bc2amr(valbig,aux,mitot,mjtot,nvar,naux,dx_fine,dy_fine,level,t, &
337  xlow_fine,xhi_fine,ylow_fine,yhi_fine)
338  endif
339 
340 contains
341 
342  integer pure function coarse_index(n,i,j)
343  implicit none
344  integer, intent(in) :: n, i, j
345  coarse_index = n + nvar*(i-1)+nvar*mitot_coarse*(j-1)
346  end function coarse_index
347 
348  logical pure function sticksout(iplo,iphi,jplo,jphi)
349  implicit none
350  integer, intent(in) :: iplo,iphi,jplo,jphi
351  sticksout = (iplo < 0 .or. jplo < 0 .or. &
352  iphi >= iregsz(level - 1) .or. jphi >= jregsz(level - 1))
353  end function sticksout
354 
subroutine bc2amr(val, aux, nrow, ncol, meqn, naux,
Take a grid patch with mesh widths hx,hy, of dimensions nrow by ncol, and set the values of any piece...
Definition: bc2amr.f:88
logical pure function sticksout(iplo, iphi, jplo, jphi)
Definition: filpatch.f90:349
recursive subroutine prefilrecur(level, nvar, valbig, auxbig, naux, time, mitot, mjtot, nrowst, ncolst, ilo, ihi, jlo, jhi, iglo, ighi, jglo, jghi, patchOnly)
For periodic boundary conditions more work needed to fill the piece of the boundary.
Definition: prefilp.f90:20
subroutine trimbd(used, nrow, ncol, set, unset_rect)
Examine the setting status of a patch.
Definition: trimbd.f90:33
integer, dimension(maxlv) iregsz
Definition: amr_module.f90:198
real(kind=8), dimension(maxlv) hyposs
Definition: amr_module.f90:193
real(kind=8), dimension(maxlv) hxposs
Definition: amr_module.f90:193
integer, dimension(maxlv) jregsz
Definition: amr_module.f90:198
subroutine intfil(val, mi, mj, time, flaguse, nrowst, ncolst, ilo, ihi, jlo, jhi, level, nvar, naux, msrc)
Fill values for a patch at the specified level and location, using values from grids at level level O...
Definition: intfil.f90:99
real(kind=8) xupper
Definition: amr_module.f90:231
real(kind=8) xlower
Definition: amr_module.f90:231
logical yperdom
Definition: amr_module.f90:230
real(kind=8) yupper
Definition: amr_module.f90:231
real(kind=8) ylower
Definition: amr_module.f90:231
real(kind=8), parameter needs_to_be_set
Definition: amr_module.f90:168
logical spheredom
Definition: amr_module.f90:230
recursive subroutine filrecur(level, nvar, valbig, aux, naux, t, mitot, mjtot, nrowst, ncolst, ilo, ihi, jlo, jhi, patchOnly, msrc)
Fill a region (patch) described by:
Definition: filpatch.f90:73
integer, dimension(maxlv) intraty
Definition: amr_module.f90:198
integer, parameter outunit
Definition: amr_module.f90:290
logical xperdom
Definition: amr_module.f90:230
integer, dimension(maxlv) intratx
Definition: amr_module.f90:198
integer nghost
Definition: amr_module.f90:232
The module contains the definition of a "node descriptor" as well as other global variables used duri...
Definition: amr_module.f90:21
subroutine setaux(mbc, mx, my, xlower, ylower, dx, dy, maux, aux)
Definition: setaux.f90:2
integer pure function coarse_index(n, i, j)
Definition: filpatch.f90:343
Here is the call graph for this function:
Here is the caller graph for this function:

◆ sticksout()

logical pure function filrecur::sticksout ( integer, intent(in)  iplo,
integer, intent(in)  iphi,
integer, intent(in)  jplo,
integer, intent(in)  jphi 
)

Definition at line 349 of file filpatch.f90.

References amr_module::iregsz, and amr_module::jregsz.

Referenced by filrecur(), and sticksout().

349  implicit none
350  integer, intent(in) :: iplo,iphi,jplo,jphi
351  sticksout = (iplo < 0 .or. jplo < 0 .or. &
352  iphi >= iregsz(level - 1) .or. jphi >= jregsz(level - 1))
logical pure function sticksout(iplo, iphi, jplo, jphi)
Definition: filpatch.f90:349
integer, dimension(maxlv) iregsz
Definition: amr_module.f90:198
integer, dimension(maxlv) jregsz
Definition: amr_module.f90:198
Here is the caller graph for this function: