2D AMRCLAW
flux2.f
Go to the documentation of this file.
1 c
2 c
24 c =====================================================
25  subroutine flux2(ixy,maxm,meqn,maux,mbc,mx,
26  & q1d,dtdx1d,aux1,aux2,aux3,
27  & faddm,faddp,gaddm,gaddp,cfl1d,wave,s,
28  & amdq,apdq,cqxx,bmasdq,bpasdq,rpn2,rpt2)
29 c =====================================================
30 c
31 c # clawpack routine ... modified for AMRCLAW
32 c
33 c # Compute the modification to fluxes f and g that are generated by
34 c # all interfaces along a 1D slice of the 2D grid.
35 c # ixy = 1 if it is a slice in x
36 c # 2 if it is a slice in y
37 c # This value is passed into the Riemann solvers. The flux modifications
38 c # go into the arrays fadd and gadd. The notation is written assuming
39 c # we are solving along a 1D slice in the x-direction.
40 c
41 c # fadd(*,i,.) modifies F to the left of cell i
42 c # gadd(*,i,.,1) modifies G below cell i
43 c # gadd(*,i,.,2) modifies G above cell i
44 c
45 c # The method used is specified by method(2:3):
46 c
47 c method(2) = 1 if only first order increment waves are to be used.
48 c = 2 if second order correction terms are to be added, with
49 c a flux limiter as specified by mthlim.
50 c
51 c method(3) = 0 if no transverse propagation is to be applied.
52 c Increment and perhaps correction waves are propagated
53 c normal to the interface.
54 c = 1 if transverse propagation of increment waves
55 c (but not correction waves, if any) is to be applied.
56 c = 2 if transverse propagation of correction waves is also
57 c to be included.
58 c
59 c Note that if mcapa>0 then the capa array comes into the second
60 c order correction terms, and is already included in dtdx1d:
61 c If ixy = 1 then
62 c dtdx1d(i) = dt/dx if mcapa= 0
63 c = dt/(dx*aux(mcapa,i,jcom)) if mcapa = 1
64 c If ixy = 2 then
65 c dtdx1d(j) = dt/dy if mcapa = 0
66 c = dt/(dy*aux(mcapa,icom,j)) if mcapa = 1
67 c
68 c Notation:
69 c The jump in q (q1d(i,:)-q1d(i-1,:)) is split by rpn2 into
70 c amdq = the left-going flux difference A^- Delta q
71 c apdq = the right-going flux difference A^+ Delta q
72 c Each of these is split by rpt2 into
73 c bmasdq = the down-going transverse flux difference B^- A^* Delta q
74 c bpasdq = the up-going transverse flux difference B^+ A^* Delta q
75 c where A^* represents either A^- or A^+.
76 c
77 c
78  use amr_module
79  implicit double precision (a-h,o-z)
80  external rpn2, rpt2
81  dimension q1d(meqn,1-mbc:maxm+mbc)
82  dimension amdq(meqn,1-mbc:maxm+mbc)
83  dimension apdq(meqn,1-mbc:maxm+mbc)
84  dimension bmasdq(meqn,1-mbc:maxm+mbc)
85  dimension bpasdq(meqn,1-mbc:maxm+mbc)
86  dimension cqxx(meqn,1-mbc:maxm+mbc)
87  dimension faddm(meqn,1-mbc:maxm+mbc)
88  dimension faddp(meqn,1-mbc:maxm+mbc)
89  dimension gaddm(meqn,1-mbc:maxm+mbc, 2)
90  dimension gaddp(meqn,1-mbc:maxm+mbc, 2)
91  dimension dtdx1d(1-mbc:maxm+mbc)
92  dimension aux1(maux,1-mbc:maxm+mbc)
93  dimension aux2(maux,1-mbc:maxm+mbc)
94  dimension aux3(maux,1-mbc:maxm+mbc)
95 c
96  dimension s(mwaves, 1-mbc:maxm+mbc)
97  dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
98 c
99  logical limit
100  common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
101 c
102  limit = .false.
103  do 5 mw=1,mwaves
104  if (mthlim(mw) .gt. 0) limit = .true.
105  5 continue
106 c
107 c # initialize flux increments:
108 c -----------------------------
109 c
110  do 10 i = 1-mbc, mx+mbc
111  do 20 m=1,meqn
112  faddm(m,i) = 0.d0
113  faddp(m,i) = 0.d0
114  gaddm(m,i,1) = 0.d0
115  gaddp(m,i,1) = 0.d0
116  gaddm(m,i,2) = 0.d0
117  gaddp(m,i,2) = 0.d0
118  20 continue
119  10 continue
120 c
121 c
122 c # solve Riemann problem at each interface and compute Godunov updates
123 c ---------------------------------------------------------------------
124 c
125  call rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,q1d,q1d,
126  & aux2,aux2,wave,s,amdq,apdq)
127 c
128 c # Set fadd for the donor-cell upwind method (Godunov)
129  do 40 i=1,mx+1
130  do 40 m=1,meqn
131  faddp(m,i) = faddp(m,i) - apdq(m,i)
132  faddm(m,i) = faddm(m,i) + amdq(m,i)
133  40 continue
134 c
135 c # compute maximum wave speed for checking Courant number:
136  cfl1d = 0.d0
137  do 50 mw=1,mwaves
138  do 50 i=1,mx+1
139 c # if s>0 use dtdx1d(i) to compute CFL,
140 c # if s<0 use dtdx1d(i-1) to compute CFL:
141  cfl1d = dmax1(cfl1d, dtdx1d(i)*s(mw,i),
142  & -dtdx1d(i-1)*s(mw,i))
143  50 continue
144 c
145  if (method(2).eq.1) go to 130
146 c
147 c # modify F fluxes for second order q_{xx} correction terms:
148 c -----------------------------------------------------------
149 c
150 c # apply limiter to waves:
151  if (limit) call limiter(maxm,meqn,mwaves,mbc,mx,wave,s,mthlim)
152 c
153  do 120 i = 1, mx+1
154 c
155 c # For correction terms below, need average of dtdx in cell
156 c # i-1 and i. Compute these and overwrite dtdx1d:
157 c
158 c # modified in Version 4.3 to use average only in cqxx, not transverse
159  dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
160 
161 c
162 c # second order corrections:
163 
164  do 120 m=1,meqn
165  cqxx(m,i) = 0.d0
166  do 119 mw=1,mwaves
167 c
168  if (use_fwaves) then
169  abs_sign = dsign(1.d0,s(mw,i))
170  else
171  abs_sign = dabs(s(mw,i))
172  endif
173 
174  cqxx(m,i) = cqxx(m,i) + abs_sign
175  & * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
176 c
177  119 continue
178  faddm(m,i) = faddm(m,i) + 0.5d0 * cqxx(m,i)
179  faddp(m,i) = faddp(m,i) + 0.5d0 * cqxx(m,i)
180  120 continue
181 c
182 c
183  130 continue
184 c
185  if (method(3).eq.0) go to 999 !# no transverse propagation
186 c
187  if (method(2).gt.1 .and. method(3).eq.2) then
188 c # incorporate cqxx into amdq and apdq so that it is split also.
189  do 150 i = 1, mx+1
190  do 150 m=1,meqn
191  amdq(m,i) = amdq(m,i) + cqxx(m,i)
192  apdq(m,i) = apdq(m,i) - cqxx(m,i)
193  150 continue
194  endif
195 c
196 c
197 c # modify G fluxes for transverse propagation
198 c --------------------------------------------
199 c
200 c
201 c # split the left-going flux difference into down-going and up-going:
202  call rpt2(ixy,1,maxm,meqn,mwaves,maux,mbc,mx,
203  & q1d,q1d,aux1,aux2,aux3,
204  & amdq,bmasdq,bpasdq)
205 c
206 c # modify flux below and above by B^- A^- Delta q and B^+ A^- Delta q:
207  do 160 i = 1, mx+1
208  do 160 m=1,meqn
209  gupdate = 0.5d0*dtdx1d(i-1) * bmasdq(m,i)
210  gaddm(m,i-1,1) = gaddm(m,i-1,1) - gupdate
211  gaddp(m,i-1,1) = gaddp(m,i-1,1) - gupdate
212 c
213  gupdate = 0.5d0*dtdx1d(i-1) * bpasdq(m,i)
214  gaddm(m,i-1,2) = gaddm(m,i-1,2) - gupdate
215  gaddp(m,i-1,2) = gaddp(m,i-1,2) - gupdate
216  160 continue
217 c
218 c # split the right-going flux difference into down-going and up-going:
219  call rpt2(ixy,2,maxm,meqn,mwaves,maux,mbc,mx,
220  & q1d,q1d,aux1,aux2,aux3,
221  & apdq,bmasdq,bpasdq)
222 c
223 c # modify flux below and above by B^- A^+ Delta q and B^+ A^+ Delta q:
224  do 180 i = 1, mx+1
225  do 180 m=1,meqn
226  gupdate = 0.5d0*dtdx1d(i) * bmasdq(m,i)
227  gaddm(m,i,1) = gaddm(m,i,1) - gupdate
228  gaddp(m,i,1) = gaddp(m,i,1) - gupdate
229 c
230  gupdate = 0.5d0*dtdx1d(i) * bpasdq(m,i)
231  gaddm(m,i,2) = gaddm(m,i,2) - gupdate
232  gaddp(m,i,2) = gaddp(m,i,2) - gupdate
233  180 continue
234 c
235  999 continue
236  return
237  end
integer, dimension(7) method
Definition: amr_module.f90:253
integer, dimension(:), allocatable mthlim
Definition: amr_module.f90:254
logical use_fwaves
Definition: amr_module.f90:257
subroutine limiter(maxm, meqn, mwaves, mbc, mx, wave, s, mthlim)
Definition: inlinelimiter.f:5
subroutine rpn2(ixy, maxm, meqn, mwaves, maux, mbc, mx, ql, qr, auxl, auxr, wave, s, amdq, apdq)
Definition: rpn2from1d.f:7
subroutine flux2(ixy, maxm, meqn, maux, mbc, mx, q1d, dtdx1d, aux1, aux2, aux3, faddm, faddp, gaddm, gaddp, cfl1d, wave, s, amdq, apdq, cqxx, bmasdq, bpasdq, rpn2, rpt2)
Solve Riemann problems based on 1D slice of the 2D grid.
Definition: flux2.f:29
The module contains the definition of a "node descriptor" as well as other global variables used duri...
Definition: amr_module.f90:21
integer mwaves
Definition: amr_module.f90:253