2D AMRCLAW
Functions/Subroutines
qad.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine qad (valbig, mitot, mjtot, nvar, svdflx, qc1d, lenbc, lratiox, lratioy, hx, hy, maux, aux, auxc1d, delt, mptr)
 For each coarse-fine interface, a Riemann problem between an inner ghost cell value on the fine grid and cell value in the adjacent coarse cell must be solved and added to corresponding location in node(ffluxptr, mptr) for conservative fix later. More...
 

Function/Subroutine Documentation

◆ qad()

subroutine qad ( dimension(nvar,mitot,mjtot)  valbig,
  mitot,
  mjtot,
  nvar,
dimension(nvar,lenbc)  svdflx,
dimension(nvar,lenbc)  qc1d,
  lenbc,
  lratiox,
  lratioy,
  hx,
  hy,
  maux,
dimension(maux,mitot,mjtot)  aux,
dimension(maux,lenbc)  auxc1d,
  delt,
  mptr 
)

For each coarse-fine interface, a Riemann problem between an inner ghost cell value on the fine grid and cell value in the adjacent coarse cell must be solved and added to corresponding location in node(ffluxptr, mptr) for conservative fix later.

Definition at line 12 of file qad.f.

References amr_module::auxtype, amr_module::dbugunit, amr_module::max1d, amr_module::maxaux, amr_module::method, amr_module::mwaves, amr_module::nestlevel, amr_module::nghost, amr_module::node, amr_module::rnode, rpn2(), src1d(), and amr_module::timemult.

Referenced by par_advanc().

12 
13  use amr_module
14  implicit double precision (a-h, o-z)
15 
16 
17  logical qprint
18 
19  dimension valbig(nvar,mitot,mjtot)
20  dimension qc1d(nvar,lenbc)
21  dimension svdflx(nvar,lenbc)
22  dimension aux(maux,mitot,mjtot)
23  dimension auxc1d(maux,lenbc)
24 
25 c
26 c ::::::::::::::::::::::::::: QAD ::::::::::::::::::::::::::::::::::
27 c are added in to coarse grid value, as a conservation fixup.
28 c Done each fine grid time step. If source terms are present, the
29 c coarse grid value is advanced by source terms each fine time step too.
30 
31 c No change needed in this sub. for spherical mapping: correctly
32 c mapped vals already in bcs on this fine grid and coarse saved
33 c vals also properly prepared
34 c
35 c Side 1 is the left side of the fine grid patch. Then go around clockwise.
36 c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
37 c
38 c # local storage
39 c # note that dimension here are bigger than dimensions used
40 c # in rp2, but shouldn't matter since wave is not used in qad
41 c # and for other arrays it is only the last parameter that is wrong
42 c # ok as long as meqn, mwaves < maxvar
43 
44  parameter(max1dp1 = max1d+1)
45  dimension ql(nvar,max1dp1), qr(nvar,max1dp1)
46  dimension wave(nvar,mwaves,max1dp1), s(mwaves,max1dp1)
47  dimension amdq(nvar,max1dp1), apdq(nvar,max1dp1)
48  dimension auxl(maxaux*max1dp1), auxr(maxaux*max1dp1)
49 c
50 c WARNING: auxl,auxr dimensioned at max possible, but used as if
51 c they were dimensioned as the real maux by max1dp1. Would be better
52 c of course to dimension by maux by max1dp1 but this wont work if maux=0
53 c So need to access using your own indexing into auxl,auxr.
54  iaddaux(iaux,i) = iaux + maux*(i-1)
55 
56  data qprint/.false./
57 c
58 c aux is auxiliary array with user parameters needed in Riemann solvers
59 c on fine grid corresponding to valbig
60 c auxc1d is coarse grid stuff from around boundary, same format as qc1d
61 c auxl, auxr are work arrays needed to pass stuff to rpn2
62 c maux is the number of aux variables, which may be zero.
63 c
64 
65  tgrid = rnode(timemult, mptr)
66  if (qprint)
67  . write(dbugunit,*)" working on grid ",mptr," time ",tgrid
68  nc = mjtot-2*nghost
69  nr = mitot-2*nghost
70  level = node(nestlevel, mptr)
71  index = 0
72 
73 c
74 c--------
75 c side 1
76 c--------
77 c
78  do 10 j = nghost+1, mjtot-nghost
79  if (maux.gt.0) then
80  do 5 ma = 1,maux
81  if (auxtype(ma).eq."xleft") then
82 c # Assuming velocity at left-face, this fix
83 c # preserves conservation in incompressible flow:
84  auxl(iaddaux(ma,j-nghost+1)) = aux(ma,nghost+1,j)
85  else
86 c # Normal case -- we set the aux arrays
87 c # from the cell corresponding to q
88  auxl(iaddaux(ma,j-nghost+1)) = aux(ma,nghost,j)
89  endif
90  5 continue
91  endif
92  do 10 ivar = 1, nvar
93  ql(ivar,j-nghost+1) = valbig(ivar,nghost,j)
94  10 continue
95 
96  lind = 0
97  ncrse = (mjtot-2*nghost)/lratioy
98  do 20 jc = 1, ncrse
99  index = index + 1
100  do 25 l = 1, lratioy
101  lind = lind + 1
102  if (maux.gt.0) then
103  do 24 ma=1,maux
104  auxr(iaddaux(ma,lind)) = auxc1d(ma,index)
105  24 continue
106  endif
107  do 25 ivar = 1, nvar
108  25 qr(ivar,lind) = qc1d(ivar,index)
109  20 continue
110 
111  if (qprint) then
112  write(dbugunit,*) 'side 1, ql and qr:'
113  do i=2,nc
114  write(dbugunit,4101) i,qr(1,i-1),ql(1,i)
115  enddo
116  4101 format(i3,4e16.6)
117  if (maux .gt. 0) then
118  write(dbugunit,*) 'side 1, auxr:'
119  do i=2,nc
120  write(dbugunit,4101) i,(auxr(iaddaux(ma,i-1)),ma=1,maux)
121  enddo
122  write(dbugunit,*) 'side 1, auxl:'
123  do i=2,nc
124  write(dbugunit,4101) i,(auxl(iaddaux(ma,i)),ma=1,maux)
125  enddo
126  endif
127  endif
128 
129  call rpn2(1,max1dp1-2*nghost,nvar,mwaves,maux,nghost,
130  . nc+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq)
131 c
132 c we have the wave. for side 1 add into sdflxm
133 c
134  influx = 0
135  do 30 j = 1, nc/lratioy
136  influx = influx + 1
137  jfine = (j-1)*lratioy
138  do 40 ivar = 1, nvar
139  do 50 l = 1, lratioy
140  svdflx(ivar,influx) = svdflx(ivar,influx)
141  . + amdq(ivar,jfine+l+1) * hy * delt
142  . + apdq(ivar,jfine+l+1) * hy * delt
143  50 continue
144  40 continue
145  30 continue
146 
147 c--------
148 c side 2
149 c--------
150 c
151  if (mjtot .eq. 2*nghost+1) then
152 c # a single row of interior cells only happens when using the
153 c # 2d amrclaw code to do a 1d problem with refinement.
154 c # (feature added in Version 4.3)
155 c # skip over sides 2 and 4 in this case
156  go to 299
157  endif
158 
159  do 210 i = nghost+1, mitot-nghost
160  if (maux.gt.0) then
161  do 205 ma = 1,maux
162  auxr(iaddaux(ma,i-nghost)) = aux(ma,i,mjtot-nghost+1)
163  205 continue
164  endif
165  do 210 ivar = 1, nvar
166  qr(ivar,i-nghost) = valbig(ivar,i,mjtot-nghost+1)
167  210 continue
168 
169  lind = 0
170  ncrse = (mitot-2*nghost)/lratiox
171  do 220 ic = 1, ncrse
172  index = index + 1
173  do 225 l = 1, lratiox
174  lind = lind + 1
175  if (maux.gt.0) then
176  do 224 ma=1,maux
177  if (auxtype(ma).eq."yleft") then
178 c # Assuming velocity at bottom-face, this fix
179 c # preserves conservation in incompressible flow:
180  ifine = (ic-1)*lratiox + nghost + l
181  auxl(iaddaux(ma,lind+1)) = aux(ma,ifine,mjtot-nghost+1)
182  else
183  auxl(iaddaux(ma,lind+1)) = auxc1d(ma,index)
184  endif
185  224 continue
186  endif
187  do 225 ivar = 1, nvar
188  225 ql(ivar,lind+1) = qc1d(ivar,index)
189  220 continue
190 
191  if (qprint) then
192  write(dbugunit,*) 'side 2, ql and qr:'
193  do i=1,nr
194  write(dbugunit,4101) i,ql(1,i+1),qr(1,i)
195  enddo
196  if (maux .gt. 0) then
197  write(dbugunit,*) 'side 2, auxr:'
198  do i = 1, nr
199  write(dbugunit,4101) i, (auxr(iaddaux(ma,i)),ma=1,maux)
200  enddo
201  write(dbugunit,*) 'side 2, auxl:'
202  do i = 1, nr
203  write(dbugunit,4101) i, (auxl(iaddaux(ma,i)),ma=1,maux)
204  enddo
205  endif
206  endif
207  call rpn2(2,max1dp1-2*nghost,nvar,mwaves,maux,nghost,
208  . nr+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq)
209 c
210 c we have the wave. for side 2. add into sdflxp
211 c
212  do 230 i = 1, nr/lratiox
213  influx = influx + 1
214  ifine = (i-1)*lratiox
215  do 240 ivar = 1, nvar
216  do 250 l = 1, lratiox
217  svdflx(ivar,influx) = svdflx(ivar,influx)
218  . - amdq(ivar,ifine+l+1) * hx * delt
219  . - apdq(ivar,ifine+l+1) * hx * delt
220  250 continue
221  240 continue
222  230 continue
223 
224  299 continue
225 
226 c--------
227 c side 3
228 c--------
229 c
230  do 310 j = nghost+1, mjtot-nghost
231  if (maux.gt.0) then
232  do 305 ma = 1,maux
233  auxr(iaddaux(ma,j-nghost)) = aux(ma,mitot-nghost+1,j)
234  305 continue
235  endif
236  do 310 ivar = 1, nvar
237  qr(ivar,j-nghost) = valbig(ivar,mitot-nghost+1,j)
238  310 continue
239 
240  lind = 0
241  ncrse = (mjtot-2*nghost)/lratioy
242  do 320 jc = 1, ncrse
243  index = index + 1
244  do 325 l = 1, lratioy
245  lind = lind + 1
246  if (maux.gt.0) then
247  do 324 ma=1,maux
248  if (auxtype(ma).eq."xleft") then
249 c # Assuming velocity at left-face, this fix
250 c # preserves conservation in incompressible flow:
251  jfine = (jc-1)*lratioy + nghost + l
252  auxl(iaddaux(ma,lind+1)) = aux(ma,mitot-nghost+1,jfine)
253  else
254  auxl(iaddaux(ma,lind+1)) = auxc1d(ma,index)
255  endif
256  324 continue
257  endif
258  do 325 ivar = 1, nvar
259  325 ql(ivar,lind+1) = qc1d(ivar,index)
260  320 continue
261 
262  if (qprint) then
263  write(dbugunit,*) 'side 3, ql and qr:'
264  do i=1,nc
265  write(dbugunit,4101) i,ql(1,i),qr(1,i)
266  enddo
267  endif
268  call rpn2(1,max1dp1-2*nghost,nvar,mwaves,maux,nghost,
269  . nc+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq)
270 c
271 c we have the wave. for side 3 add into sdflxp
272 C
273  do 330 j = 1, nc/lratioy
274  influx = influx + 1
275  jfine = (j-1)*lratioy
276  do 340 ivar = 1, nvar
277  do 350 l = 1, lratioy
278  svdflx(ivar,influx) = svdflx(ivar,influx)
279  . - amdq(ivar,jfine+l+1) * hy * delt
280  . - apdq(ivar,jfine+l+1) * hy * delt
281  350 continue
282  340 continue
283  330 continue
284 
285 c--------
286 c side 4
287 c--------
288 c
289  if (mjtot .eq. 2*nghost+1) then
290 c # a single row of interior cells only happens when using the
291 c # 2d amrclaw code to do a 1d problem with refinement.
292 c # (feature added in Version 4.3)
293 c # skip over sides 2 and 4 in this case
294  go to 499
295  endif
296 c
297  do 410 i = nghost+1, mitot-nghost
298  if (maux.gt.0) then
299  do 405 ma = 1,maux
300  if (auxtype(ma).eq."yleft") then
301 c # Assuming velocity at bottom-face, this fix
302 c # preserves conservation in incompressible flow:
303  auxl(iaddaux(ma,i-nghost+1)) = aux(ma,i,nghost+1)
304  else
305  auxl(iaddaux(ma,i-nghost+1)) = aux(ma,i,nghost)
306  endif
307  405 continue
308  endif
309  do 410 ivar = 1, nvar
310  ql(ivar,i-nghost+1) = valbig(ivar,i,nghost)
311  410 continue
312 
313  lind = 0
314  ncrse = (mitot-2*nghost)/lratiox
315  do 420 ic = 1, ncrse
316  index = index + 1
317  do 425 l = 1, lratiox
318  lind = lind + 1
319  if (maux.gt.0) then
320  do 424 ma=1,maux
321  auxr(iaddaux(ma,lind)) = auxc1d(ma,index)
322  424 continue
323  endif
324  do 425 ivar = 1, nvar
325  425 qr(ivar,lind) = qc1d(ivar,index)
326  420 continue
327 
328  if (qprint) then
329  write(dbugunit,*) 'side 4, ql and qr:'
330  do i=1,nr
331  write(dbugunit,4101) i, ql(1,i),qr(1,i)
332  enddo
333  endif
334  call rpn2(2,max1dp1-2*nghost,nvar,mwaves,maux,nghost,
335  . nr+1-2*nghost,ql,qr,auxl,auxr,wave,s,amdq,apdq)
336 c
337 c we have the wave. for side 4. add into sdflxm
338 c
339  do 430 i = 1, nr/lratiox
340  influx = influx + 1
341  ifine = (i-1)*lratiox
342  do 440 ivar = 1, nvar
343  do 450 l = 1, lratiox
344  svdflx(ivar,influx) = svdflx(ivar,influx)
345  . + amdq(ivar,ifine+l+1) * hx * delt
346  . + apdq(ivar,ifine+l+1) * hx * delt
347  450 continue
348  440 continue
349  430 continue
350 
351  499 continue
352 
353 c # for source terms:
354  if (method(5) .ne. 0) then ! should I test here if index=0 and all skipped?
355  call src1d(nvar,nghost,lenbc,qc1d,maux,auxc1d,tgrid,delt)
356 c # how can this be right - where is the integrated src term used?
357  endif
358 
359  return
amr_module::dbugunit
integer, parameter dbugunit
Definition: amr_module.f90:293
amr_module::method
integer, dimension(7) method
Definition: amr_module.f90:253
amr_module::timemult
integer, parameter timemult
current simulation time on this grid
Definition: amr_module.f90:151
src1d
subroutine src1d(meqn, mbc, mx1d, q1d, maux, aux1d, t, dt)
Definition: src1d.f90:2
amr_module::max1d
integer, parameter max1d
Definition: amr_module.f90:181
amr_module::node
integer, dimension(nsize, maxgr) node
Definition: amr_module.f90:198
amr_module::nestlevel
integer, parameter nestlevel
AMR level of the grid.
Definition: amr_module.f90:44
amr_module::rnode
real(kind=8), dimension(rsize, maxgr) rnode
Definition: amr_module.f90:193
amr_module::maxaux
integer, parameter maxaux
Definition: amr_module.f90:184
rpn2
subroutine rpn2(ixy, maxm, meqn, mwaves, maux, mbc, mx, ql, qr, auxl, auxr, wave, s, amdq, apdq)
Definition: rpn2from1d.f:7
amr_module::auxtype
character(len=10), dimension(:), allocatable auxtype
Definition: amr_module.f90:252
amr_module::nghost
integer nghost
Definition: amr_module.f90:232
amr_module
The module contains the definition of a "node descriptor" as well as other global variables used duri...
Definition: amr_module.f90:21
amr_module::mwaves
integer mwaves
Definition: amr_module.f90:253
Here is the call graph for this function:
Here is the caller graph for this function:
Generated on Sun May 7 2017 10:59:40 for 2D AMRCLAW by   doxygen 1.8.14