2D AMRCLAW
filval_orig.f90
Go to the documentation of this file.
1 !
2 ! :::::::::::::::::::::::::::::: FILVAL ::::::::::::::::::::::::::
3 !
4 ! create and fill coarser (level-1) patch with one extra coarse cell all
5 ! around, plus the ghost cells . will interpolate from this patch to grid mptr
6 ! without needing special boundary code.
7 !
8 ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
9 !
10 !
11 ! ------------------------------------------------------------------
12 !
13 subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, &
14  mjc, xleft, xright, ybot, ytop, nvar, mptr, ilo, ihi, &
15  jlo, jhi, aux, naux)
16 
20 
21  !for setaux timing
22  use amr_module, only: timesetaux, timesetauxcpu
23 
24  implicit none
25 
26  ! Input
27  integer, intent(in) :: mitot, mjtot, level, mic, mjc, nvar, mptr, ilo, ihi
28  integer, intent(in) :: jlo, jhi, naux
29  real(kind=8), intent(in) :: dx, dy, time, xleft, xright, ybot, ytop
30 
31  ! Output
32  real(kind=8), intent(in out) :: val(nvar,mitot,mjtot), aux(naux,mitot,mjtot)
33 
34  ! Local storage
35  integer :: refinement_ratio_x, refinement_ratio_y, iclo, jclo, ichi, jchi, ng
36  integer :: ivar, i, j, ico, jco, ifine, jfine, nx, ny
37  real(kind=8) :: valc(nvar,mic,mjc), auxc(naux,mic,mjc)
38  real(kind=8) :: dx_coarse, dy_coarse, xl, xr, yb, yt, area
39  real(kind=8) :: s1m, s1p, slopex, slopey, xoff, yoff
40  real(kind=8) :: fliparray((mitot+mjtot)*nghost*(nvar+naux))
41  real(kind=8) :: setflags(mitot,mjtot),maxauxdif,aux2(naux,mitot,mjtot)
42  integer :: mjb
43 
44  !for timing
45  integer :: clock_start, clock_finish, clock_rate
46  real(kind=8) cpu_start, cpu_finish
47 
48 
49  ! External function definitions
50  real(kind=8) :: get_max_speed
51 
52  refinement_ratio_x = intratx(level-1)
53  refinement_ratio_y = intraty(level-1)
54  dx_coarse = dx * refinement_ratio_x
55  dy_coarse = dy * refinement_ratio_y
56  xl = xleft - dx_coarse
57  xr = xright + dx_coarse
58  yb = ybot - dy_coarse
59  yt = ytop + dy_coarse
60 
61  ! if topo not yet final then aux is set outside filval (in gfixup)
62  ! and so aux has real data already, (ie dont overwrite here)
63 
64  ! set integer indices for coarser patch enlarged by 1 cell
65  ! (can stick out of domain). proper nesting will insure this one
66  ! call is sufficient.
67  iclo = ilo / refinement_ratio_x - 1
68  jclo = jlo / refinement_ratio_y - 1
69  ichi = (ihi + 1) / refinement_ratio_x - 1 + 1
70  jchi = (jhi + 1) / refinement_ratio_y - 1 + 1
71  ng = 0
72 
73  if (naux == 0) then
74  if (xperdom .or. yperdom .or. spheredom) then
75  call preintcopy(valc,mic,mjc,nvar,iclo,ichi,jclo,jchi,level-1,fliparray)
76  else
77  call intcopy(valc,mic,mjc,nvar,iclo,ichi,jclo,jchi,level-1,1,1)
78  endif
79  else
80  ! intersect grids and copy all (soln and aux)
81  auxc(1,:,:) = needs_to_be_set
82  if (xperdom .or. yperdom .or. spheredom) then
83  call preicall(valc,auxc,mic,mjc,nvar,naux,iclo,ichi,jclo,jchi, &
84  level-1,fliparray)
85  else
86  call icall(valc,auxc,mic,mjc,nvar,naux,iclo,ichi,jclo,jchi,level-1,1,1)
87  endif
88 !!$ do i = 1, mic
89 !!$ do j = 1, mjc
90 !!$ if (auxc(1,:,:) == NEEDS_TO_BE_SET) then
91 !!$ write(*,*)" *** coarsenened new fine grid not completely set
92 !from previously" &
93 !!$ "existing coarse grids ***"
94 !!$ stop
95 !!$ endif
96 !!$ end do
97 !!$ end do
98  ! no ghost cells on coarse enlarged patch. set any remaining
99  ! vals. should only be bcs that stick out of domain.
100  call system_clock(clock_start,clock_rate)
101  call cpu_time(cpu_start)
102  call setaux(ng,mic,mjc,xl,yb,dx_coarse,dy_coarse,naux,auxc)
103  call system_clock(clock_finish,clock_rate)
104  call cpu_time(cpu_finish)
105  timesetaux=timesetaux+clock_finish-clock_start
106  timesetauxcpu=timesetauxcpu+cpu_finish-cpu_start
107  endif
108  call bc2amr(valc,auxc,mic,mjc,nvar,naux,dx_coarse,dy_coarse,level-1,time,xl,xr,yb, &
109  yt,xlower,ylower,xupper,yupper,xperdom,yperdom,spheredom)
110 
111 
112 ! NOTE change in order of code. Since the interp from coarse to fine needs the
113 ! aux
114 ! arrays set already, the fine copy is done first, to set up the aux
115 ! arrays.
116 ! we can do this since we have the flag array to test where to overwrite.
117 
118 ! SO this is no longer overwriting but setting for the first time.
119 ! overwrite interpolated values with fine grid values, if available.
120  nx = mitot - 2*nghost
121  ny = mjtot - 2*nghost
122 
123  if (naux .gt. 0) then
124 ! ## NEEDS_TO_BE_SET is signal that aux array not set.
125 ! ## after calling icall to copy aux from other grids
126 ! ## any remaining NEEDS_TO_BE_SET signals will be set in setaux.
127 ! ## it also signals where soln was copied, so it wont be
128 ! ## overwritten with coarse grid interpolation
129  aux(1,:,:) = needs_to_be_set ! indicates fine cells not yet set.
130 
131  if (xperdom .or. yperdom .or. spheredom) then
132  call preicall(val,aux,mitot,mjtot,nvar,naux,ilo-nghost,ihi+nghost, &
133  jlo-nghost,jhi+nghost,level,1,1,fliparray)
134  else
135  call icall(val,aux,mitot,mjtot,nvar,naux,ilo-nghost,ihi+nghost, &
136  jlo-nghost,jhi+nghost,level,1,1)
137  endif
138  setflags = aux(1,:,:) ! save since will overwrite in setaux when setting all aux vals
139  ! need this so we know where to use coarse grid to set fine solution
140  ! w/o overwriting
141  ! set remaining aux vals not set by copying from prev existing grids
142  call system_clock(clock_start,clock_rate)
143  call cpu_time(cpu_start)
144  call setaux(nghost,nx,ny,xleft,ybot,dx,dy,naux,aux)
145  call system_clock(clock_finish,clock_rate)
146  call cpu_time(cpu_finish)
147  timesetaux=timesetaux+clock_finish-clock_start
148  timesetauxcpu=timesetauxcpu+cpu_finish-cpu_start
149  else ! either no aux exists, or cant reuse yet
150  ! so only call intcopy (which copies soln) and not icall.
151  ! in this case flag q(1,:) to NEEDS_TO_BE_SET flag so wont be
152  ! overwritten
153  ! by coarse grid interp. this is needed due to reversing order of
154  ! work - first copy from fine grids, then interpolate from coarse grids
155  val(1,:,:) = needs_to_be_set
156  if (xperdom .or. yperdom .or. spheredom) then
157  call preintcopy(val,mitot,mjtot,nvar,ilo-nghost,ihi+nghost, &
158  jlo-nghost,jhi+nghost,level,fliparray)
159  else
160  call intcopy(val,mitot,mjtot,nvar,ilo-nghost,ihi+nghost, &
161  jlo-nghost,jhi+nghost,level,1,1)
162  endif
163  setflags = val(1,:,:) ! remaining flags signals need to set
164  endif
165 
166 
167  ! Prepare slopes - use min-mod limiters
168 
169  do j=2, mjc-1
170  do i=2, mic-1
171  do ivar = 1, nvar
172 
173  s1p = valc(ivar,i+1,j) - valc(ivar,i,j)
174  s1m = valc(ivar,i,j) - valc(ivar,i-1,j)
175  slopex = min(abs(s1p), abs(s1m)) &
176  * sign(1.d0,valc(ivar,i+1,j) - valc(ivar,i-1,j))
177  ! if there's a sign change, set slope to 0.
178  if ( s1m*s1p <= 0.d0) slopex = 0.d0
179 
180  s1p = valc(ivar,i,j+1) - valc(ivar,i,j)
181  s1m = valc(ivar,i,j) - valc(ivar,i,j-1)
182  slopey = min(abs(s1p), abs(s1m)) &
183  * sign(1.0d0, valc(ivar,i,j+1) - valc(ivar,i,j-1))
184  if ( s1m*s1p <= 0.d0) slopey = 0.d0
185 
186  ! Interpolate from coarse cells to fine grid to find depth
187  do jco = 1,refinement_ratio_y
188  yoff = (real(jco,kind=8) - 0.5d0) / refinement_ratio_y - 0.5d0
189  jfine = (j-2) * refinement_ratio_y + nghost + jco
190 
191  do ico = 1,refinement_ratio_x
192  xoff = (real(ico,kind=8) - 0.5d0) / refinement_ratio_x - 0.5d0
193  ifine = (i-2) * refinement_ratio_x + nghost + ico
194 
195  if (setflags(ifine,jfine) .eq. needs_to_be_set) then
196  val(ivar,ifine,jfine) = valc(ivar,i,j) + xoff*slopex + yoff*slopey
197  endif
198 
199  end do
200  end do
201 
202  enddo !end of ivar loop
203  enddo !end of coarse i loop
204  enddo !end of coarse j loop
205 
206  ! adjust to conserve kappa*q, but only where coarse grid was interpolated
207  ! so now need to pass setflags to this subr.
208  if (mcapa .ne. 0) then
209  call fixcapaq(val,aux,mitot,mjtot,valc,auxc,mic,mjc,nvar,naux,level-1,setflags)
210  endif
211 
212 !!$! CHECK BY CALLING SETAUX AND SETTING ALL, THEN DIFFING
213 !!$ mjb = 0
214 !!$ if (naux .gt. 0 .and. mjb .eq. 1) then
215 !!$ aux2(1,:,:) = NEEDS_TO_BE_SET ! indicates fine cells not yet set
216 !!$ call setaux(nghost,nx,ny,xleft,ybot,dx,dy,naux,aux2)
217 !!$ maxauxdif = 1.d-13
218 !!$ do i = 1, mitot
219 !!$ do j = 1, mjtot
220 !!$ if (abs(aux(1,i,j)-aux2(1,i,j)) .gt. maxauxdif) then
221 !!$ maxauxdif = abs(aux(1,i,j)-aux2(1,i,j))
222 !!$ write(*,444)i,j,aux(1,i,j),aux2(1,i,j),maxauxdif
223 !!$444 format("i,j = ",2i4," auxs ",2e15.7," maxauxdif ",e12.5)
224 !!$ endif
225 !!$ end do
226 !!$ end do
227 !!$ if (maxauxdif .gt. 2.d-13) then
228 !!$ write(*,*)" maxauxdif = ",maxauxdif," with mitot,mjtot
229 !",mitot,mjtot, &
230 !!$ " on grid ",mptr," level ",level
231 !!$ endif
232 !!$ endif
233 
234 end subroutine filval
235 
236 ! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
237 
238 subroutine dumpaux(aux,naux,mitot,mjtot)
239  implicit none
240  real(kind=8) :: aux(naux,mitot,mjtot)
241  integer :: naux,mitot,mjtot,i,j,iaux
242 
243  do j = 1, mjtot
244  do i = 1, mitot
245  write(*,444) i,j,(aux(iaux,i,j),iaux=1,naux)
246  444 format(2i4,5e12.5)
247  end do
248  end do
249 
250 end subroutine dumpaux
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
subroutine filval(val, mitot, mjtot, dx, dy, level, time, mic, mjc, xleft, xright, ybot, ytop, nvar, mptr, ilo, ihi, jlo, jhi, aux, naux)
Fill grid mptr on level level by copying values from OLD level level grids if available, otherwise by interpolating values from coarser grids.
Definition: filval.f90:6
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
subroutine icall(val, aux, nrow, ncol, nvar, naux, ilo, ihi, jlo, jhi, level, iputst, jputst)
For a rectangle defined on level level and bound by ilo, ihi, jlo, jhi, find intersecting grids at th...
Definition: icall.f:10
subroutine intcopy(val, mitot, mjtot, nvar, ilo, ihi, jlo, jhi, level, iputst, jputst)
For a rectangle that is on level level, described by ilo, ihi, jlo, jhi and made up by mitot mjtot c...
Definition: intcopy.f:14
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
subroutine fixcapaq(val, aux, mitot, mjtot, valc, auxc, mic, mjc, nvar, naux, levc, setflags)
Definition: fixcapaq.f:6
logical spheredom
Definition: amr_module.f90:230
subroutine preintcopy(val, mitot, mjtot, nvar, ilo, ihi, jlo, jhi, level, fliparray)
Definition: preintcopy.f:6
integer mcapa
Definition: amr_module.f90:253
integer, dimension(maxlv) intraty
Definition: amr_module.f90:198
subroutine dumpaux(aux, naux, mitot, mjtot)
Definition: filval.f90:219
subroutine preicall(val, aux, nrow, ncol, nvar, naux, ilo, ihi, jlo, jhi, level, fliparray)
Definition: preicall.f:6
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
real(kind=8), dimension(:), allocatable alloc
Definition: amr_module.f90:218