2D AMRCLAW
filpatch_orig.f90
Go to the documentation of this file.
1 ! :::::::::::::::::::::::::::: FILPATCH :::::::::::::::::::::::::;
2 !
3 ! fill the portion of valbig from rows nrowst
4 ! and cols ncolst
5 ! the patch can also be described by the corners (xlp,ybp) by (xrp,ytp).
6 ! vals are needed at time t, and level level,
7 !
8 ! first fill with values obtainable from the level level
9 ! grids. if any left unfilled, then enlarge remaining rectangle of
10 ! unfilled values by 1 (for later linear interp), and recusively
11 ! obtain the remaining values from coarser levels.
12 !
13 ! :::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::;
14 !
15 recursive subroutine filrecur(level,nvar,valbig,aux,naux,t,mx,my, &
16  nrowst,ncolst,ilo,ihi,jlo,jhi)
17 
21 
22  !for setaux timing
23  use amr_module, only: timesetaux, timesetauxcpu
24 
25  implicit none
26 
27  ! Input
28  integer, intent(in) :: level, nvar, naux, mx, my, nrowst, ncolst
29  integer, intent(in) :: ilo,ihi,jlo,jhi
30 ! integer, intent(in) :: fill_indices(4)
31  real(kind=8), intent(in) :: t
32 
33  ! Output
34  real(kind=8), intent(in out) :: valbig(nvar,mx,my)
35  real(kind=8), intent(in out) :: aux(naux,mx,my)
36 
37  ! Local storage
38  integer :: i_fine, j_fine, i_coarse, j_coarse, n, k
39  integer :: iplo, iphi, jplo, jphi
40  integer :: mx_patch, my_patch, mx_coarse, my_coarse, nghost_patch, lencrse
41  integer :: refinement_ratio_x, refinement_ratio_y
42  integer :: unset_indices(4)
43  real(kind=8) :: dx_fine, dy_fine, dx_coarse, dy_coarse
44  real(kind=8) :: xlow_coarse,ylow_coarse, xlow_fine, ylow_fine, xhi_fine,yhi_fine
45 
46  !for setaux timing
47  integer :: clock_start, clock_finish, clock_rate
48  real(kind=8) :: cpu_start, cpu_finish
49 
50  ! Interpolation variables
51  real(kind=8) :: eta1, eta2, valp10, valm10, valc, valp01, valm01, dupc, dumc
52  real(kind=8) :: ducc, du, fu, dvpc, dvmc, dvcc, dv, fv, valint
53 
54  ! Cell set tracking
55  logical :: set
56  integer(kind=1) :: flaguse(ihi-ilo+1, jhi-jlo+1)
57 
58  ! Scratch storage
59  ! use stack-based scratch arrays instead of alloc, since dont really
60  ! need to save beyond these routines, and to allow dynamic memory resizing
61  !
62  ! use 1d scratch arrays that are potentially the same size as
63  ! current grid, since may not coarsen.
64  ! need to make it 1d instead of 2 and do own indexing, since
65  ! when pass it in to subroutines they treat it as having dierent
66  ! dimensions than the max size need to allocate here
67 
68  !-- dimension valcrse((ihi-ilo+2)*(jhi-jlo+2)*nvar) ! NB this is a 1D
69  !array
70  !-- dimension auxcrse((ihi-ilo+2)*(jhi-jlo+2)*naux) ! the +2 is to
71  !expand on coarse grid to enclose fine
72  ! ### turns out you need 3 rows, forget offset of 1 plus one on each side
73  ! the +3 is to expand on coarse grid to enclose fine
74 !!$ real(kind=8) :: valcrse((fill_indices(2) - fill_indices(1) + 3) &
75 !!$ * (fill_indices(4) - fill_indices(3) + 3) *nvar)
76 !!$ real(kind=8) :: auxcrse((fill_indices(2) - fill_indices(1) + 3) &
77 !!$ * (fill_indices(4) - fill_indices(3) + 3) *naux)
78  real(kind=8) :: valcrse((ihi-ilo+3) * (jhi-jlo+3) * nvar)
79  real(kind=8) :: auxcrse((ihi-ilo+3) * (jhi-jlo+3) * naux)
80  ! We begin by filling values for grids at level level. If all values can be
81  ! filled in this way, we return;
82 !!$ mx_patch = fill_indices(2) - fill_indices(1) + 1 ! nrowp
83 !!$ my_patch = fill_indices(4) - fill_indices(3) + 1 ! ncolp
84  mx_patch = ihi-ilo + 1 ! nrowp
85  my_patch = jhi-jlo + 1
86 
87  dx_fine = hxposs(level)
88  dy_fine = hyposs(level)
89 
90  ! Coordinates of edges of patch (xlp,xrp,ybp,ytp)
91 !!$ fill_rect = [xlower + fill_indices(1) * dx_fine, &
92 !!$ xlower + (fill_indices(2) + 1) * dx_fine, &
93 !!$ ylower + fill_indices(3) * dy_fine, &
94 !!$ ylower + (fill_indices(4) + 1) * dy_fine]
95  xlow_fine = xlower + ilo * dx_fine
96  ylow_fine = ylower + jlo * dy_fine
97 
98  ! Fill in the patch as much as possible using values at this level
99 !!$ call intfil(valbig,mx,my,t,flaguse,nrowst,ncolst,fill_indices(1), &
100 !!$
101 !fill_indices(2),fill_indices(3),fill_indices(4),level,nvar,naux)
102  call intfil(valbig,mx,my,t,flaguse,nrowst,ncolst, ilo, &
103  ihi,jlo,jhi,level,nvar,naux)
104 
105 
106  ! Trimbd returns set = true if all of the entries are filled (=1.).
107  ! set = false, otherwise. If set = true, then no other levels are
108  ! are required to interpolate, and we return.
109  !
110  ! Note that the used array is filled entirely in intfil, i.e. the
111  ! marking done there also takes into account the points filled by
112  ! the boundary conditions. bc2amr will be called later, after all 4
113  ! boundary pieces filled.
114  call trimbd(flaguse,mx_patch,my_patch,set,unset_indices) ! il,ir,jb,jt = unset_indices(4)
115 
116  ! If set is .true. then all cells have been set and we can skip to setting
117  ! the remaining boundary cells. If it is .false. we need to interpolate
118  ! some values from coarser levels, possibly calling this routine
119  ! recursively.
120  if (.not.set) then
121 
122  ! Error check
123  if (level == 1) then
124  write(outunit,*)" error in filrecur - level 1 not set"
125  write(outunit,'("start at row: ",i4," col ",i4)') nrowst,ncolst
126  print *," error in filrecur - level 1 not set"
127  print *," should not need more recursion "
128  print *," to set patch boundaries"
129  print '("start at row: ",i4," col ",i4)', nrowst,ncolst
130  stop
131  endif
132 
133  ! We begin by initializing the level level arrays, so that we can use
134  ! purely recursive formulation for interpolating.
135  dx_coarse = hxposs(level - 1)
136  dy_coarse = hyposs(level - 1)
137 
138  ! Adjust unset_indices to account for the patch
139  ! isl, isr, jsb, jst
140  unset_indices(1) = unset_indices(1) + ilo - 1
141  unset_indices(2) = unset_indices(2) + ilo - 1
142  unset_indices(3) = unset_indices(3) + jlo - 1
143  unset_indices(4) = unset_indices(4) + jlo - 1
144 
145  ! Coarsened geometry
146  refinement_ratio_x = intratx(level - 1)
147  refinement_ratio_y = intraty(level - 1)
148 
149  ! New patch rectangle (after we have partially filled it in) but in the
150  ! coarse patches [iplo,iphi,jplo,jphi]
151 
152  iplo = (unset_indices(1) - refinement_ratio_x + nghost * refinement_ratio_x) &
153  / refinement_ratio_x - nghost
154  iphi = (unset_indices(2) + refinement_ratio_x) / refinement_ratio_x
155  jplo = (unset_indices(3) - refinement_ratio_y + nghost * refinement_ratio_y) &
156  / refinement_ratio_y - nghost
157  jphi = (unset_indices(4) + refinement_ratio_y) / refinement_ratio_y
158 
159  xlow_coarse = xlower + iplo * dx_coarse
160  ylow_coarse = ylower + jplo * dy_coarse
161 
162  ! Coarse grid number of spatial points (nrowc,ncolc)
163  mx_coarse = iphi - iplo + 1
164  my_coarse = jphi - jplo + 1
165 
166  ! Check to make sure we created big enough scratch arrays
167  if (mx_coarse > ihi - ilo + 3 .or. &
168  my_coarse > jhi - jlo + 3) then
169 
170  print *," did not make big enough work space in filrecur "
171  print *," need coarse space with mx_coarse,my_coarse ",mx_coarse,my_coarse
172  print *," made space for ilo,ihi,jlo,jhi ",ilo,ihi,jlo,jhi
173  stop
174  endif
175 
176  ! Set the aux array values for the coarse grid, this could be done
177  ! instead in intfil using possibly already available bathy data from the
178  ! grids
179  if (naux > 0) then
180  nghost_patch = 0
181  lencrse = (ihi-ilo+3)*(jhi-jlo+3)*naux ! set 1 component, not all naux
182  do k = 1, lencrse, naux
183  auxcrse(k) = needs_to_be_set ! new system checks initialization before setting aux vals
184  end do
185  call system_clock(clock_start, clock_rate)
186  call cpu_time(cpu_start)
187  call setaux(nghost_patch, mx_coarse,my_coarse, &
188  xlow_coarse, ylow_coarse, &
189  dx_coarse,dy_coarse,naux,auxcrse)
190  call system_clock(clock_finish, clock_rate)
191  call cpu_time(cpu_finish)
192  timesetaux = timesetaux + clock_finish - clock_start
193  timesetauxcpu = timesetaux + cpu_finish - cpu_start
194  endif
195 
196  ! Fill in the edges of the coarse grid
197  if ((xperdom .or. (yperdom .or. spheredom)) .and. sticksout(iplo,iphi,jplo,jphi)) then
198  call prefilrecur(level-1,nvar,valcrse,auxcrse,naux,t,mx_coarse,my_coarse,1,1,iplo,iphi,jplo,jphi)
199  else
200  call filrecur(level-1,nvar,valcrse,auxcrse,naux,t,mx_coarse,my_coarse,1,1,iplo,iphi,jplo,jphi)
201  endif
202 
203  do i_fine = 1,mx_patch
204  i_coarse = 2 + (i_fine - (unset_indices(1) - ilo) - 1) / refinement_ratio_x
205  eta1 = (-0.5d0 + real(mod(i_fine - 1, refinement_ratio_x),kind=8)) &
206  / real(refinement_ratio_x,kind=8)
207 
208  do j_fine = 1,my_patch
209  j_coarse = 2 + (j_fine - (unset_indices(3) - jlo) - 1) / refinement_ratio_y
210  eta2 = (-0.5d0 + real(mod(j_fine - 1, refinement_ratio_y),kind=8)) &
211  / real(refinement_ratio_y,kind=8)
212 
213  if (flaguse(i_fine,j_fine) == 0) then
214 
215  do n=1,nvar
216 
217  valp10 = valcrse(coarse_index(n,i_coarse + 1,j_coarse))
218  valm10 = valcrse(coarse_index(n,i_coarse - 1,j_coarse))
219  valc = valcrse(coarse_index(n,i_coarse ,j_coarse))
220  valp01 = valcrse(coarse_index(n,i_coarse ,j_coarse + 1))
221  valm01 = valcrse(coarse_index(n,i_coarse ,j_coarse - 1))
222 
223  dupc = valp10 - valc
224  dumc = valc - valm10
225  ducc = valp10 - valm10
226  du = min(abs(dupc), abs(dumc))
227  du = min(2.d0 * du, 0.5d0 * abs(ducc))
228  fu = max(0.d0, sign(1.d0, dupc * dumc))
229 
230  dvpc = valp01 - valc
231  dvmc = valc - valm01
232  dvcc = valp01 - valm01
233  dv = min(abs(dvpc), abs(dvmc))
234  dv = min(2.d0 * dv, 0.5d0 * abs(dvcc))
235  fv = max(0.d0,sign(1.d0, dvpc * dvmc))
236 
237  valint = valc + eta1 * du * sign(1.d0, ducc) * fu &
238  + eta2 * dv * sign(1.d0, dvcc) * fv
239 
240 
241  valbig(n,i_fine+nrowst-1,j_fine+ncolst-1) = valint
242 
243  end do
244 
245  endif
246  end do
247  end do
248  end if
249 
250  ! set bcs, whether or not recursive calls needed. set any part of patch
251  ! that stuck out
252  xhi_fine = xlower + (ihi + 1) * dx_fine
253  yhi_fine = ylower + (jhi + 1) * dy_fine
254  call bc2amr(valbig,aux,mx,my,nvar,naux,dx_fine,dy_fine,level,t,xlow_fine,xhi_fine, &
255  ylow_fine,yhi_fine,xlower,ylower,xupper,yupper,xperdom,yperdom,spheredom)
256 
257 contains
258 
259  integer pure function coarse_index(n,i,j)
260  implicit none
261  integer, intent(in) :: n, i, j
262  coarse_index = n + nvar*(i-1)+nvar*mx_coarse*(j-1)
263  end function coarse_index
264 
265  logical pure function sticksout(iplo,iphi,jplo,jphi)
266  implicit none
267  integer, intent(in) :: iplo,iphi,jplo,jphi
268  sticksout = (iplo < 0 .or. jplo < 0 .or. &
269  iphi >= iregsz(level - 1) .or. jphi >= jregsz(level - 1))
270  end function sticksout
271 
272 end subroutine filrecur
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