2D AMRCLAW
flux2_dimSplit.f
Go to the documentation of this file.
1 c
2 c
3 c =====================================================
4  subroutine flux2_dimsplit(ixy,maxm,meqn,maux,mbc,mx,
5  & q1d,dtdx1d,aux2,
6  & faddm,faddp,cfl1d,wave,s,
7  & cqxx,rpn2)
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 In dimensionally split version no need for extra arrays amdq,apdq
52 c Use faddm and faddp instead. BUT REMEMBER to negate apdq
53 c since it was previously SUBTRACTED from fadpp.
54 c
55 c
56  use amr_module
57  implicit double precision (a-h,o-z)
58  external rpn2
59  dimension q1d(meqn,1-mbc:maxm+mbc)
60  dimension cqxx(meqn,1-mbc:maxm+mbc)
61  dimension faddm(meqn,1-mbc:maxm+mbc)
62  dimension faddp(meqn,1-mbc:maxm+mbc)
63  dimension dtdx1d(1-mbc:maxm+mbc)
64  dimension aux2(maux,1-mbc:maxm+mbc)
65 c
66  dimension s(mwaves, 1-mbc:maxm+mbc)
67  dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
68 c
69  logical limit
70 c
71  limit = .false.
72  do 5 mw=1,mwaves
73  if (mthlim(mw) .gt. 0) limit = .true.
74  5 continue
75 c
76 c # initialize flux increments:
77 c -----------------------------
78 c
79  do 10 i = 1-mbc, mx+mbc
80  do 20 m=1,meqn
81  faddm(m,i) = 0.d0
82  faddp(m,i) = 0.d0
83  20 continue
84  10 continue
85 c
86 c
87 c # solve Riemann problem at each interface and compute Godunov updates
88 c ---------------------------------------------------------------------
89 c
90  call rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,q1d,q1d,
91  & aux2,aux2,wave,s,faddm,faddp)
92 c
93  faddp = -faddp !! NEGATED, replaces faddp(m,i) = faddp(m,i) - apdq(m,i) in original flux2.
94 c
95 c # dont delete this code yet, in case there is an issue with loop indices
96 ! # Set fadd for the donor-cell upwind method (Godunov)
97 !-- do 40 i=1,mx+1
98 !-- do 40 m=1,meqn
99 !-- faddp(m,i) = faddp(m,i) - apdq(m,i)
100 !-- faddm(m,i) = faddm(m,i) + amdq(m,i)
101 !-- 40 continue
102 c
103 c # compute maximum wave speed for checking Courant number:
104  cfl1d = 0.d0
105  do 50 mw=1,mwaves
106  do 50 i=1,mx+1
107 c # if s>0 use dtdx1d(i) to compute CFL,
108 c # if s<0 use dtdx1d(i-1) to compute CFL:
109  cfl1d = dmax1(cfl1d, dtdx1d(i)*s(mw,i),
110  & -dtdx1d(i-1)*s(mw,i))
111  50 continue
112 c
113  if (method(2).eq.1) go to 130
114 c
115 c # modify F fluxes for second order q_{xx} correction terms:
116 c -----------------------------------------------------------
117 c
118 c # apply limiter to waves:
119  if (limit) call limiter(maxm,meqn,mwaves,mbc,mx,wave,s,mthlim)
120 c
121  do 120 i = 1, mx+1
122 c
123 c # For correction terms below, need average of dtdx in cell
124 c # i-1 and i. Compute these and overwrite dtdx1d:
125 c
126 c # modified in Version 4.3 to use average only in cqxx, not transverse
127  dtdxave = 0.5d0 * (dtdx1d(i-1) + dtdx1d(i))
128 
129 c
130 c # second order corrections:
131 
132  do 120 m=1,meqn
133  cqxx(m,i) = 0.d0
134  do 119 mw=1,mwaves
135 c
136  if (use_fwaves) then
137  abs_sign = dsign(1.d0,s(mw,i))
138  else
139  abs_sign = dabs(s(mw,i))
140  endif
141 
142  cqxx(m,i) = cqxx(m,i) + abs_sign
143  & * (1.d0 - dabs(s(mw,i))*dtdxave) * wave(m,mw,i)
144 c
145  119 continue
146  faddm(m,i) = faddm(m,i) + 0.5d0 * cqxx(m,i)
147  faddp(m,i) = faddp(m,i) + 0.5d0 * cqxx(m,i)
148  120 continue
149 c
150 c
151  130 continue
152 c
153 
154  return
155  end
flux2_dimsplit
subroutine flux2_dimsplit(ixy, maxm, meqn, maux, mbc, mx, q1d, dtdx1d, aux2, faddm, faddp, cfl1d, wave, s, cqxx, rpn2)
Definition: flux2_dimSplit.f:4
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
Generated on Thu Mar 9 2017 15:15:52 for 2D AMRCLAW by   doxygen 1.8.6