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