2D AMRCLAW
stepgrid.f
Go to the documentation of this file.
1 c
2 c -------------------------------------------------------------
3 c
4  subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
5  & nvar,xlow,ylow,time,mptr,maux,aux)
6 c
7 c
8 c ::::::::::::::::::: STEPGRID ::::::::::::::::::::::::::::::::::::
9 c take a time step on a single grid. overwrite solution array q.
10 c A modified version of the clawpack routine step2 is used.
11 c
12 c return fluxes in fm,fp and gm,gp.
13 c patch has room for ghost cells (mbc of them) around the grid.
14 c everything is the enlarged size (mitot by mjtot).
15 c
16 c mbc = number of ghost cells (= lwidth)
17 c mptr = grid number (for debugging)
18 c xlow,ylow = lower left corner of enlarged grid (including ghost cells).
19 c dt = incoming time step
20 c dx,dy = mesh widths for this grid
21 c dtnew = return suggested new time step for this grid's soln.
22 c
23 c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
24 
25  use amr_module
26  implicit double precision (a-h,o-z)
27  external rpn2,rpt2
28 
29  common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
30 
31  parameter(msize=max1d+4)
32  parameter(mwork=msize*(maxvar*maxvar + 13*maxvar + 3*maxaux +2))
33 
34  dimension q(nvar,mitot,mjtot)
35  dimension fp(nvar,mitot,mjtot),gp(nvar,mitot,mjtot)
36  dimension fm(nvar,mitot,mjtot),gm(nvar,mitot,mjtot)
37  dimension aux(maux,mitot,mjtot)
38 c dimension work(mwork)
39 
40  logical debug, dump
41  data debug/.false./, dump/.false./
42 
43 c
44 c # set tcom = time. This is in the common block comxyt that could
45 c # be included in the Riemann solver, for example, if t is explicitly
46 c # needed there.
47 
48  tcom = time
49 
50  if (dump) then
51  write(outunit,*) "dumping grid ",mptr," at time ",time
52  do i = 1, mitot
53  do j = 1, mjtot
54  write(outunit,545) i,j,(q(ivar,i,j),ivar=1,nvar)
55 c . ,(aux(ivar,i,j),ivar=1,maux)
56  545 format(2i4,5e15.7)
57  end do
58  end do
59  endif
60 c
61  meqn = nvar
62  mx = mitot - 2*mbc
63  my = mjtot - 2*mbc
64  maxm = max(mx,my) !# size for 1d scratch array
65  mbig = maxm
66  xlowmbc = xlow + mbc*dx
67  ylowmbc = ylow + mbc*dy
68 
69 c # method(2:7) and mthlim
70 c # are set in the amr2ez file (read by amr)
71 c
72  method(1) = 0
73 c
74 c
75 c # fluxes initialized in step2
76 c
77 C mwork0 = (maxm+2*mbc)*(12*meqn + mwaves + meqn*mwaves + 2)
78 c
79 C if (mwork .lt. mwork0) then
80 C write(outunit,*) 'CLAW2 ERROR... mwork must be increased to ',
81 C & mwork0
82 C write(* ,*) 'CLAW2 ERROR... mwork must be increased to ',
83 C & mwork0
84 C stop
85 C endif
86 c
87 c
88 c # partition work array into pieces needed for local storage in
89 c # step2 routine. Find starting index of each piece:
90 c
91 C i0faddm = 1
92 C i0faddp = i0faddm + (maxm+2*mbc)*meqn
93 C i0gaddm = i0faddp + (maxm+2*mbc)*meqn
94 C i0gaddp = i0gaddm + 2*(maxm+2*mbc)*meqn
95 C i0q1d = i0gaddp + 2*(maxm+2*mbc)*meqn
96 C i0dtdx1 = i0q1d + (maxm+2*mbc)*meqn
97 C i0dtdy1 = i0dtdx1 + (maxm+2*mbc)
98 C i0aux1 = i0dtdy1 + (maxm+2*mbc)
99 C i0aux2 = i0aux1 + (maxm+2*mbc)*maux
100 C i0aux3 = i0aux2 + (maxm+2*mbc)*maux
101 c
102 c
103 C i0next = i0aux3 + (maxm+2*mbc)*maux !# next free space
104 C mused = i0next - 1 !# space already used
105 C mwork1 = mwork - mused !# remaining space (passed to step2)
106 
107 c
108 c
109  call b4step2(mbc,mx,my,nvar,q,
110  & xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux)
111 c
112 c
113 c # take one step on the conservation law:
114 c
115  call step2(mbig,nvar,maux,
116  & mbc,mx,my,
117  & q,aux,dx,dy,dt,cflgrid,
118  & fm,fp,gm,gp,rpn2,rpt2)
119 c
120 c
121  write(outunit,1001) mptr, node(nestlevel,mptr),cflgrid
122  1001 format(' Courant # of grid', i4,
123  & ' on level', i3, ' is ', e10.3)
124 c
125 
126 !$OMP CRITICAL (cflm)
127 
128  cfl_level = dmax1(cfl_level,cflgrid)
129 
130 !$OMP END CRITICAL (cflm)
131 
132 c
133 c # update q
134  dtdx = dt/dx
135  dtdy = dt/dy
136  do 50 j=mbc+1,mjtot-mbc
137  do 50 i=mbc+1,mitot-mbc
138  do 50 m=1,nvar
139  if (mcapa.eq.0) then
140 c
141 c # no capa array. Standard flux differencing:
142 
143  q(m,i,j) = q(m,i,j)
144  & - dtdx * (fm(m,i+1,j) - fp(m,i,j))
145  & - dtdy * (gm(m,i,j+1) - gp(m,i,j))
146  else
147 c # with capa array.
148  q(m,i,j) = q(m,i,j)
149  & - (dtdx * (fm(m,i+1,j) - fp(m,i,j))
150  & + dtdy * (gm(m,i,j+1) - gp(m,i,j))) / aux(mcapa,i,j)
151  endif
152 
153  50 continue
154 c
155 c
156  if (method(5).eq.1) then
157 c # with source term: use Godunov splitting
158  call src2(nvar,mbc,mx,my,xlowmbc,ylowmbc,dx,dy,
159  & q,maux,aux,time,dt)
160  endif
161 c
162 c
163 c
164 c # output fluxes for debugging purposes:
165  if (debug) then
166  write(dbugunit,*)" fluxes for grid ",mptr
167 c do 830 j = mbc+1, mjtot-1
168  do 830 i = mbc+1, mitot-1
169  do 830 j = mbc+1, mjtot-1
170  write(dbugunit,831) i,j,fm(1,i,j),fp(1,i,j),
171  . gm(1,i,j),gp(1,i,j)
172  do 830 m = 2, meqn
173  write(dbugunit,832) fm(m,i,j),fp(m,i,j),
174  . gm(m,i,j),gp(m,i,j)
175  831 format(2i4,4d16.6)
176  832 format(8x,4d16.6)
177  830 continue
178  endif
179 
180 c
181 c
182 c For variable time stepping, use max speed seen on this grid to
183 c choose the allowable new time step dtnew. This will later be
184 c compared to values seen on other grids.
185 c
186  if (cflgrid .gt. 0.d0) then
187  dtnew = dt*cfl/cflgrid
188  else
189 c # velocities are all zero on this grid so there's no
190 c # time step restriction coming from this grid.
191  dtnew = rinfinity
192  endif
193 
194 c # give a warning if Courant number too large...
195 c
196  if (cflgrid .gt. cflv1) then
197  write(*,810) cflgrid
198  write(outunit,810) cflgrid, cflv1
199  810 format('*** WARNING *** Courant number =', d12.4,
200  & ' is larger than input cfl_max = ', d12.4)
201  endif
202 c
203  if (dump) then
204  write(outunit,*) "dumping grid ",mptr," after stepgrid"
205  do i = 1, mitot
206  do j = 1, mjtot
207  write(outunit,545) i,j,(q(ivar,i,j),ivar=1,nvar)
208  end do
209  end do
210  endif
211  return
212  end
213 
214 
b4step2
subroutine b4step2(mbc, mx, my, meqn, q, xlower, ylower, dx, dy, t, dt, maux, aux)
Definition: b4step2.f90:2
src2
subroutine src2(meqn, mbc, mx, my, xlower, ylower, dx, dy, q, maux, aux, t, dt)
Definition: src2.f90:1
amr_module
Definition: amr_module.f90:1
stepgrid
subroutine stepgrid(q, fm, fp, gm, gp, mitot, mjtot, mbc, dt, dtnew, dx, dy, nvar, xlow, ylow, time, mptr, maux, aux)
Definition: stepgrid.f:4
rpn2
subroutine rpn2(ixy, maxm, meqn, mwaves, maux, mbc, mx, ql, qr, auxl, auxr, wave, s, amdq, apdq)
Definition: rpn2from1d.f:5
step2
subroutine step2(maxm, meqn, maux, mbc, mx, my, qold, aux, dx, dy, dt, cflgrid, fm, fp, gm, gp, rpn2, rpt2)
Definition: step2.f90:1
Generated on Thu Mar 9 2017 15:15:52 for 2D AMRCLAW by   doxygen 1.8.6