2D AMRCLAW
Functions/Subroutines
flux2.f File Reference

Go to the source code of this file.

Functions/Subroutines

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. More...
 

Function/Subroutine Documentation

◆ flux2()

subroutine flux2 (   ixy,
  maxm,
  meqn,
  maux,
  mbc,
  mx,
dimension(meqn,1-mbc:maxm+mbc)  q1d,
dimension(1-mbc:maxm+mbc)  dtdx1d,
dimension(maux,1-mbc:maxm+mbc)  aux1,
dimension(maux,1-mbc:maxm+mbc)  aux2,
dimension(maux,1-mbc:maxm+mbc)  aux3,
dimension(meqn,1-mbc:maxm+mbc)  faddm,
dimension(meqn,1-mbc:maxm+mbc)  faddp,
dimension(meqn,1-mbc:maxm+mbc, 2)  gaddm,
dimension(meqn,1-mbc:maxm+mbc, 2)  gaddp,
  cfl1d,
dimension(meqn, mwaves, 1-mbc:maxm+mbc)  wave,
dimension(mwaves, 1-mbc:maxm+mbc)  s,
dimension(meqn,1-mbc:maxm+mbc)  amdq,
dimension(meqn,1-mbc:maxm+mbc)  apdq,
dimension(meqn,1-mbc:maxm+mbc)  cqxx,
dimension(meqn,1-mbc:maxm+mbc)  bmasdq,
dimension(meqn,1-mbc:maxm+mbc)  bpasdq,
external  rpn2,
external  rpt2 
)

Solve Riemann problems based on 1D slice of the 2D grid.

Return fluxes contributed by this slice.

If ixy = 1, it's a slice in x. If ixy = 2, it's a slice in y.

Suppose it's a slice in x, the subroutine call normal Riemann solver to get fluctuations at edges(i-1/2,j) and store the fluctuations in flux form: faddm, faddp. It also call transverse Riemann soler to get up-goning and down-going transverse waves and store them in gaddm and gaddp.

gaddm is value to be added to gm, which is g flux on the lower side of a cell edge. gaddp is value to be added to gp, which is g flux on the upper side of a cell edge. Here gaddm and gaddp have the same value since g flux on upper side and lower side of a cell edge should be the same (it's not like apdq and amdq at the edge). gaddm(*,i,.,1) or gaddp(*,i,.,1) modifies G below cell i; gaddm(*,i,.,2) or gaddp(*,i,.,2) modifies G above cell i;

Definition at line 29 of file flux2.f.

References limiter(), amr_module::method, amr_module::mthlim, amr_module::mwaves, rpn2(), and amr_module::use_fwaves.

Referenced by step2().

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
amr_module::method
integer, dimension(7) method
Definition: amr_module.f90:253
amr_module::mthlim
integer, dimension(:), allocatable mthlim
Definition: amr_module.f90:254
amr_module::use_fwaves
logical use_fwaves
Definition: amr_module.f90:257
limiter
subroutine limiter(maxm, meqn, mwaves, mbc, mx, wave, s, mthlim)
Definition: inlinelimiter.f:5
rpn2
subroutine rpn2(ixy, maxm, meqn, mwaves, maux, mbc, mx, ql, qr, auxl, auxr, wave, s, amdq, apdq)
Definition: rpn2from1d.f:7
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