2D AMRCLAW
intfil.f90
Go to the documentation of this file.
1 !
2 ! ::::::::::::::::::::::: INTFIL ::::::::::::::::::::::::::::::::;
3 ! INTFIL: interpolates values for a patch at the specified level and
4 ! location, using values from grids at LEVEL and coarser, if nec.
5 !
6 ! take the intersection of a grid patch with corners at ilo,ihi,jlo,jhi
7 ! and all grids mptr at LEVEL. If there is a non-null intersection
8 ! copy the solution vaues from mptr (at TIME) into VAL array.
9 ! assumes patch at same levels do straight copy, not skipping
10 ! every intrat or doing any interpolation here,
11 ! assume called in correct order of levels, so that when copying
12 ! is ok to overwrite.
13 !
14 ! N.B.: there are no dummy points around patch, since
15 ! this is not an official "grid" but a "patch".
16 !
17 ! used array marks when point filled. filpatch checks if points left over
18 ! after intersections at specified level.
19 ! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;
20 !
21 !
95 
96 
97 
98 subroutine intfil(val,mi,mj,time,flaguse,nrowst,ncolst,ilo,ihi,jlo,jhi,level,nvar,naux,msrc)
99 
102  use amr_module, only: rnode, ndilo, ndihi, ndjlo, ndjhi, nextfree
103  use amr_module, only: bndlistnum, bndlistst
105 
106  implicit none
107 
108  ! Input
109  integer, intent(in) :: mi, mj, nrowst, ncolst, ilo, ihi, jlo, jhi, level, nvar, naux,msrc
110  real(kind=8), intent(in) :: time
111 
112  ! In/Out
113  integer(kind=1), intent(in out) :: flaguse(ilo:ihi, jlo:jhi)
114  real(kind=8), intent(in out) :: val(nvar,mi,mj)
115 
116  ! Locals
117  integer :: imlo, jmlo, imhi, jmhi, nx, ny, mitot, mjtot
118  integer :: ixlo, ixhi, jxlo, jxhi, locold, locnew, nextSpot
119  integer :: icount, bndNum, bndLoc, levSt
120  integer :: ivar, i, j, mptr, mstart, loc, numg
121  real(kind=8) :: dt, alphac, alphai
122  logical :: t_interpolate
123 
124  integer :: patch_rect(4)
125 
126  real(kind=8), parameter :: t_epsilon = 1.0d-4
127 
128  ! Formats for error statements
129  character(len=*), parameter :: missing_error = &
130  "(' time wanted ',e15.7,' not available from grid ',i4,'level',i4)"
131  character(len=*), parameter :: time_error = &
132  "(' trying to interpolate from previous time values ',/," // &
133  "' for a patch with corners ilo,ihi,jlo,jhi:'" // &
134  ",/,2x,4i10,/," // &
135  "' from source grid ',i4,' at level ',i4,/," // &
136  "' time wanted ',e24.16,' source time is ',e24.16,/," // &
137  "' alphai, t_epsilon ',2e24.16)"
138 
139  patch_rect = [ilo,ihi,jlo,jhi]
140 
141  ! Note that we need a non-dimensionalized t epspatch_rect(1)n as it was a problem
142  ! in tsunami tests ( instead of t_epsilon = dt / 10.d0 )
143 
144  ! Time step at this level
145  dt = possk(level)
146 
147  ! Initialize the flagging where we set things
148  flaguse = 0
149 
150  ! Loop through all grids at this level, initialize to first
151 ! mptr = lstart(level)
152 ! do while (mptr /= 0)
153  if (msrc .eq. -1) then
154  numg = numgrids(level)
155  levst = liststart(level)
156  else
157  bndloc = node(bndlistst,msrc) ! index of first grid in bndList
158  bndnum = node(bndlistnum,msrc)
159  nextspot = node(bndlistst, msrc) ! initialize
160  numg = bndnum
161  endif
162 
163  do icount = 1, numg
164 
165  if (msrc .eq. -1) then
166  mptr = listofgrids(levst+icount-1)
167  else
168  mptr = bndlist(nextspot,gridnbor)
169  endif
170 
171  ! Check if grid mptr and patch intersect
172  imlo = node(ndilo, mptr)
173  jmlo = node(ndjlo, mptr)
174  imhi = node(ndihi, mptr)
175  jmhi = node(ndjhi, mptr)
176 
177  nx = node(ndihi,mptr) - node(ndilo,mptr) + 1
178  ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1
179 
180  mitot = nx + 2 * nghost
181  mjtot = ny + 2 * nghost
182 
183  ixlo = max(imlo,patch_rect(1))
184  ixhi = min(imhi,patch_rect(2))
185  jxlo = max(jmlo,patch_rect(3))
186  jxhi = min(jmhi,patch_rect(4))
187 
188  ! Check to see if grid and patch interesect, if not continue to next
189  ! grid in the list
190  if (ixlo <= ixhi .and. jxlo <= jxhi) then
191 
192  ! grids intersect. figure out what time to use.
193  ! alphai = 1 for new time; 0 for old time
194  alphac = (rnode(timemult,mptr) - time)/dt
195  alphai = 1.d0 - alphac
196 
197  if ((alphai < -t_epsilon) .or. (alphai > 1.d0 + t_epsilon)) then
198  write(outunit,missing_error) time, mptr, level
199  print missing_error, time, mptr, level
200  write(outunit,'(A,E24.16)') 'Line 80', dt
201  write(outunit,time_error) patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
202  print time_error, patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
203  call outtre(mstart,.false.,nvar,naux)
204  stop
205  endif
206 
207  ! Check if we should interpolate in time
208  t_interpolate = .false.
209  if (abs(alphai - 1.d0) < t_epsilon) then
210  ! need no interpolation
211  loc = node(store1,mptr)
212  else if (dabs(alphai) .lt. t_epsilon) then
213  loc = node(store2,mptr)
214  if (level == mxnest) then
215  write(outunit,'(A,E24.16)') 'Line 95', dt
216  write(outunit,time_error) patch_rect,mptr,level,time, &
217  rnode(timemult,mptr),alphai,t_epsilon
218  write(*,time_error) patch_rect,mptr,level,time, &
219  rnode(timemult,mptr),alphai,t_epsilon
220  stop
221  endif
222  else
223  locold = node(store2,mptr)
224  locnew = node(store1,mptr)
225  t_interpolate = .true.
226 
227  ! If we are at the maximum level nesting, abort
228  if (level == mxnest) then
229  write(outunit,'(A,E24.16)') 'Line 107',dt
230  write(outunit,time_error) patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
231  print time_error, patch_rect,mptr,level,time,rnode(timemult,mptr),alphai,t_epsilon
232  stop
233  endif
234  endif
235 
236  ! Actual interpolation
237  if (.not. t_interpolate) then
238  ! No time interp. copy the solution values
239  do ivar = 1, nvar
240  do j = jxlo, jxhi
241  do i = ixlo, ixhi
242  val(ivar,i-patch_rect(1)+nrowst,j-jlo+ncolst) = &
243  alloc(iadd(ivar,i-imlo+nghost+1,j-jmlo+nghost+1))
244  flaguse(i,j) = 1
245  end do
246  end do
247  end do
248  else
249  ! Linear interpolation in time
250  do ivar = 1, nvar
251  do j = jxlo, jxhi
252  do i = ixlo, ixhi
253  val(ivar,i-patch_rect(1)+nrowst,j-jlo+ncolst) = &
254  alloc(iadnew(ivar,i-imlo+nghost+1,j-jmlo+nghost+1))*alphai + &
255  alloc(iadold(ivar,i-imlo+nghost+1,j-jmlo+nghost+1))*alphac
256  flaguse(i,j) = 1
257  end do
258  end do
259  end do
260  endif
261 
262  endif
263 
264  ! Get next grid
265 ! mptr = node(levelptr, mptr)
266  if (msrc .ne. -1) then
267  nextspot = bndlist(nextspot,nextfree)
268  endif
269 
270  end do
271 
272  ! Set used array points which intersect domain boundary to be equal to 1,
273  ! they will be set elsewhere, namely boundary conditions and periodic
274  ! domains
275  if (jhi >= jregsz(level)) then
276  flaguse(patch_rect(1):ihi, max(jregsz(level),jlo):jhi) = 1
277  endif
278 
279  if (jlo < 0) then
280  flaguse(patch_rect(1):ihi, jlo:min(-1,ncolst + jhi - jlo )) = 1
281  endif
282 
283  if (patch_rect(1) < 0) then
284  flaguse(patch_rect(1):min(-1,nrowst + ihi - patch_rect(1)), jlo:jhi) = 1
285  endif
286 
287  if (ihi >= iregsz(level)) then
288  flaguse(max(iregsz(level),patch_rect(1)):ihi, jlo:jhi) = 1
289  endif
290 
291 contains
292 
293  integer pure function iadd(ivar,i,j)
294  implicit none
295  integer, intent(in) :: ivar, i, j
296  iadd = loc + ivar-1 + nvar*((j-1)*mitot+i-1)
297  end function iadd
298 
299  integer pure function iadnew(ivar,i,j)
300  implicit none
301  integer, intent(in) :: ivar, i, j
302  iadnew = locnew + ivar-1 + nvar*((j-1)*mitot+i-1)
303  end function iadnew
304 
305  integer pure function iadold(ivar,i,j)
306  implicit none
307  integer, intent(in) :: ivar, i, j
308  iadold = locold + ivar-1 + nvar*((j-1)*mitot+i-1)
309  end function iadold
310 
311 end subroutine intfil
integer, parameter timemult
current simulation time on this grid
Definition: amr_module.f90:151
integer, parameter bndlistnum
number of grids (on the same level) that border this grid
Definition: amr_module.f90:139
integer, dimension(maxlv) iregsz
Definition: amr_module.f90:198
integer, parameter ndihi
global i index of right border of this grid
Definition: amr_module.f90:111
integer, dimension(maxlv) jregsz
Definition: amr_module.f90:198
integer, dimension(nsize, maxgr) node
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
integer, dimension(0:maxlv+1) liststart
Definition: amr_module.f90:189
integer, dimension(maxlv) numgrids
Definition: amr_module.f90:198
real(kind=8), dimension(rsize, maxgr) rnode
Definition: amr_module.f90:193
integer pure function iadnew(ivar, i, j)
Definition: intfil.f90:300
integer, dimension(maxgr) listofgrids
Definition: amr_module.f90:189
integer pure function iadd(ivar, i, j)
Definition: intfil.f90:294
integer, parameter gridnbor
Definition: amr_module.f90:158
integer, parameter bndlistst
pointer (actually it&#39;s an index in the bndList array) to the first node of a linked list...
Definition: amr_module.f90:136
integer, parameter ndilo
global i index of left border of this grid
Definition: amr_module.f90:108
integer, parameter ndjlo
global j index of lower border of this grid
Definition: amr_module.f90:114
integer pure function iadold(ivar, i, j)
Definition: intfil.f90:306
integer, parameter store1
pointer to the address of memory storing the first copy of solution data on this grid, usually for storing new solution
Definition: amr_module.f90:101
integer, dimension(maxlv) lstart
Definition: amr_module.f90:198
integer, parameter store2
pointer to the address of memory storing the second copy of solution data on this grid...
Definition: amr_module.f90:105
integer, dimension(bndlistsize, 2) bndlist
Definition: amr_module.f90:191
integer, parameter outunit
Definition: amr_module.f90:290
real(kind=8), dimension(maxlv) possk
Definition: amr_module.f90:193
integer, parameter nextfree
Definition: amr_module.f90:154
integer, parameter ndjhi
global j index of upper border of this grid
Definition: amr_module.f90:117
subroutine outtre(mlev, outgrd, nvar, naux)
Output a subtree of the grids.
Definition: outtre.f:11
integer, parameter levelptr
node number (index) of next grid on the same level
Definition: amr_module.f90:35
integer mxnest
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
real(kind=8), dimension(:), allocatable alloc
Definition: amr_module.f90:218