Skip to content

Commit 697aed8

Browse files
zingalezingale
committed
revert the last change
1 parent 0d2a13f commit 697aed8

5 files changed

Lines changed: 3189 additions & 0 deletions

File tree

‎advection_fv4/interface.py‎

Lines changed: 302 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,302 @@
1+
import numpy as np
2+
from numba import njit
3+
4+
5+
@njit(cache=True)
6+
def states(a, qx, qy, ng, idir):
7+
r"""
8+
Predict the cell-centered state to the edges in one-dimension using the
9+
reconstructed, limited slopes. We use a fourth-order Godunov method.
10+
11+
Our convention here is that:
12+
13+
``al[i,j]`` will be :math:`al_{i-1/2,j}`,
14+
15+
``al[i+1,j]`` will be :math:`al_{i+1/2,j}`.
16+
17+
Parameters
18+
----------
19+
a : ndarray
20+
The cell-centered state.
21+
qx, qy : int
22+
The dimensions of `a`.
23+
ng : int
24+
The number of ghost cells
25+
idir : int
26+
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
27+
28+
Returns
29+
-------
30+
out : ndarray, ndarray
31+
The state predicted to the left and right edges.
32+
"""
33+
34+
al = np.zeros((qx, qy))
35+
ar = np.zeros((qx, qy))
36+
37+
a_int = np.zeros((qx, qy))
38+
dafm = np.zeros((qx, qy))
39+
dafp = np.zeros((qx, qy))
40+
d2af = np.zeros((qx, qy))
41+
d2ac = np.zeros((qx, qy))
42+
d3a = np.zeros((qx, qy))
43+
44+
C2 = 1.25
45+
C3 = 0.1
46+
47+
nx = qx - 2 * ng
48+
ny = qy - 2 * ng
49+
ilo = ng
50+
ihi = ng + nx
51+
jlo = ng
52+
jhi = ng + ny
53+
54+
# we need interface values on all faces of the domain
55+
if (idir == 1):
56+
57+
for i in range(ilo - 2, ihi + 3):
58+
for j in range(jlo - 1, jhi + 1):
59+
60+
# interpolate to the edges
61+
a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \
62+
(1.0 / 12.0) * (a[i - 2, j] + a[i + 1, j])
63+
64+
al[i, j] = a_int[i, j]
65+
ar[i, j] = a_int[i, j]
66+
67+
for i in range(ilo - 2, ihi + 3):
68+
for j in range(jlo - 1, jhi + 1):
69+
# these live on cell-centers
70+
dafm[i, j] = a[i, j] - a_int[i, j]
71+
dafp[i, j] = a_int[i + 1, j] - a[i, j]
72+
73+
# these live on cell-centers
74+
d2af[i, j] = 6.0 * (a_int[i, j] - 2.0 *
75+
a[i, j] + a_int[i + 1, j])
76+
77+
for i in range(ilo - 3, ihi + 3):
78+
for j in range(jlo - 1, jhi + 1):
79+
d2ac[i, j] = a[i - 1, j] - 2.0 * a[i, j] + a[i + 1, j]
80+
81+
for i in range(ilo - 2, ihi + 3):
82+
for j in range(jlo - 1, jhi + 1):
83+
# this lives on the interface
84+
d3a[i, j] = d2ac[i, j] - d2ac[i - 1, j]
85+
86+
# this is a look over cell centers, affecting
87+
# i-1/2,R and i+1/2,L
88+
for i in range(ilo - 1, ihi + 1):
89+
for j in range(jlo - 1, jhi + 1):
90+
91+
# limit? MC Eq. 24 and 25
92+
if (dafm[i, j] * dafp[i, j] <= 0.0 or
93+
(a[i, j] - a[i - 2, j]) * (a[i + 2, j] - a[i, j]) <= 0.0):
94+
95+
# we are at an extrema
96+
97+
s = np.copysign(1.0, d2ac[i, j])
98+
if (s == np.copysign(1.0, d2ac[i - 1, j]) and s == np.copysign(1.0, d2ac[i + 1, j]) and
99+
s == np.copysign(1.0, d2af[i, j])):
100+
# MC Eq. 26
101+
d2a_lim = s * min(abs(d2af[i, j]), C2 * abs(d2ac[i - 1, j]),
102+
C2 * abs(d2ac[i, j]), C2 * abs(d2ac[i + 1, j]))
103+
else:
104+
d2a_lim = 0.0
105+
106+
if (abs(d2af[i, j]) <= 1.e-12 * max(abs(a[i - 2, j]), abs(a[i - 1, j]),
107+
abs(a[i, j]), abs(a[i + 1, j]), abs(a[i + 2, j]))):
108+
rho = 0.0
109+
else:
110+
# MC Eq. 27
111+
rho = d2a_lim / d2af[i, j]
112+
113+
if (rho < 1.0 - 1.e-12):
114+
# we may need to limit -- these quantities are at cell-centers
115+
d3a_min = min(d3a[i - 1, j], d3a[i, j],
116+
d3a[i + 1, j], d3a[i + 2, j])
117+
d3a_max = max(d3a[i - 1, j], d3a[i, j],
118+
d3a[i + 1, j], d3a[i + 2, j])
119+
120+
if (C3 * max(abs(d3a_min), abs(d3a_max)) <= (d3a_max - d3a_min)):
121+
# limit
122+
if (dafm[i, j] * dafp[i, j] < 0.0):
123+
# Eqs. 29, 30
124+
ar[i, j] = a[i, j] - rho * \
125+
dafm[i, j] # note: typo in Eq 29
126+
al[i + 1, j] = a[i, j] + rho * dafp[i, j]
127+
elif (abs(dafm[i, j]) >= 2.0 * abs(dafp[i, j])):
128+
# Eq. 31
129+
ar[i, j] = a[i, j] - 2.0 * \
130+
(1.0 - rho) * dafp[i, j] - rho * dafm[i, j]
131+
elif (abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j])):
132+
# Eq. 32
133+
al[i + 1, j] = a[i, j] + 2.0 * \
134+
(1.0 - rho) * dafm[i, j] + rho * dafp[i, j]
135+
136+
else:
137+
# if Eqs. 24 or 25 didn't hold we still may need to limit
138+
if (abs(dafm[i, j]) >= 2.0 * abs(dafp[i, j])):
139+
ar[i, j] = a[i, j] - 2.0 * dafp[i, j]
140+
141+
if (abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j])):
142+
al[i + 1, j] = a[i, j] + 2.0 * dafm[i, j]
143+
144+
elif (idir == 2):
145+
146+
for i in range(ilo - 1, ihi + 1):
147+
for j in range(jlo - 2, jhi + 3):
148+
149+
# interpolate to the edges
150+
a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \
151+
(1.0 / 12.0) * (a[i, j - 2] + a[i, j + 1])
152+
153+
al[i, j] = a_int[i, j]
154+
ar[i, j] = a_int[i, j]
155+
156+
for i in range(ilo - 1, ihi + 1):
157+
for j in range(jlo - 2, jhi + 3):
158+
# these live on cell-centers
159+
dafm[i, j] = a[i, j] - a_int[i, j]
160+
dafp[i, j] = a_int[i, j + 1] - a[i, j]
161+
162+
# these live on cell-centers
163+
d2af[i, j] = 6.0 * (a_int[i, j] - 2.0 *
164+
a[i, j] + a_int[i, j + 1])
165+
166+
for i in range(ilo - 1, ihi + 1):
167+
for j in range(jlo - 3, jhi + 3):
168+
d2ac[i, j] = a[i, j - 1] - 2.0 * a[i, j] + a[i, j + 1]
169+
170+
for i in range(ilo - 1, ihi + 1):
171+
for j in range(jlo - 2, jhi + 2):
172+
# this lives on the interface
173+
d3a[i, j] = d2ac[i, j] - d2ac[i, j - 1]
174+
175+
# this is a look over cell centers, affecting
176+
# j-1/2,R and j+1/2,L
177+
for i in range(ilo - 1, ihi + 1):
178+
for j in range(jlo - 1, jhi + 1):
179+
180+
# limit? MC Eq. 24 and 25
181+
if (dafm[i, j] * dafp[i, j] <= 0.0 or
182+
(a[i, j] - a[i, j - 2]) * (a[i, j + 2] - a[i, j]) <= 0.0):
183+
184+
# we are at an extrema
185+
186+
s = np.copysign(1.0, d2ac[i, j])
187+
if (s == np.copysign(1.0, d2ac[i, j - 1]) and s == np.copysign(1.0, d2ac[i, j + 1]) and
188+
s == np.copysign(1.0, d2af[i, j])):
189+
# MC Eq. 26
190+
d2a_lim = s * min(abs(d2af[i, j]), C2 * abs(d2ac[i, j - 1]),
191+
C2 * abs(d2ac[i, j]), C2 * abs(d2ac[i, j + 1]))
192+
else:
193+
d2a_lim = 0.0
194+
195+
if (abs(d2af[i, j]) <= 1.e-12 * max(abs(a[i, j - 2]), abs(a[i, j - 1]),
196+
abs(a[i, j]), abs(a[i, j + 1]), abs(a[i, j + 2]))):
197+
rho = 0.0
198+
else:
199+
# MC Eq. 27
200+
rho = d2a_lim / d2af[i, j]
201+
202+
if (rho < 1.0 - 1.e-12):
203+
# we may need to limit -- these quantities are at cell-centers
204+
d3a_min = min(d3a[i, j - 1], d3a[i, j],
205+
d3a[i, j + 1], d3a[i, j + 2])
206+
d3a_max = max(d3a[i, j - 1], d3a[i, j],
207+
d3a[i, j + 1], d3a[i, j + 2])
208+
209+
if (C3 * max(abs(d3a_min), abs(d3a_max)) <= (d3a_max - d3a_min)):
210+
# limit
211+
if (dafm[i, j] * dafp[i, j] < 0.0):
212+
# Eqs. 29, 30
213+
ar[i, j] = a[i, j] - rho * \
214+
dafm[i, j] # note: typo in Eq 29
215+
al[i, j + 1] = a[i, j] + rho * dafp[i, j]
216+
elif (abs(dafm[i, j]) >= 2.0 * abs(dafp[i, j])):
217+
# Eq. 31
218+
ar[i, j] = a[i, j] - 2.0 * \
219+
(1.0 - rho) * dafp[i, j] - rho * dafm[i, j]
220+
elif (abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j])):
221+
# Eq. 32
222+
al[i, j + 1] = a[i, j] + 2.0 * \
223+
(1.0 - rho) * dafm[i, j] + rho * dafp[i, j]
224+
225+
else:
226+
# if Eqs. 24 or 25 didn't hold we still may need to limit
227+
if (abs(dafm[i, j]) >= 2.0 * abs(dafp[i, j])):
228+
ar[i, j] = a[i, j] - 2.0 * dafp[i, j]
229+
230+
if (abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j])):
231+
al[i, j + 1] = a[i, j] + 2.0 * dafm[i, j]
232+
233+
return al, ar
234+
235+
236+
@njit(cache=True)
237+
def states_nolimit(a, qx, qy, ng, idir):
238+
r"""
239+
Predict the cell-centered state to the edges in one-dimension using the
240+
reconstructed slopes (and without slope limiting). We use a fourth-order
241+
Godunov method.
242+
243+
Our convention here is that:
244+
245+
``al[i,j]`` will be :math:`al_{i-1/2,j}`,
246+
247+
``al[i+1,j]`` will be :math:`al_{i+1/2,j}`.
248+
249+
Parameters
250+
----------
251+
a : ndarray
252+
The cell-centered state.
253+
qx, qy : int
254+
The dimensions of `a`.
255+
ng : int
256+
The number of ghost cells
257+
idir : int
258+
Are we predicting to the edges in the x-direction (1) or y-direction (2)?
259+
260+
Returns
261+
-------
262+
out : ndarray, ndarray
263+
The state predicted to the left and right edges.
264+
"""
265+
266+
a_int = np.zeros((qx, qy))
267+
al = np.zeros((qx, qy))
268+
ar = np.zeros((qx, qy))
269+
270+
nx = qx - 2 * ng
271+
ny = qy - 2 * ng
272+
ilo = ng
273+
ihi = ng + nx
274+
jlo = ng
275+
jhi = ng + ny
276+
277+
# we need interface values on all faces of the domain
278+
if (idir == 1):
279+
280+
for i in range(ilo - 2, ihi + 3):
281+
for j in range(jlo - 1, jhi + 1):
282+
283+
# interpolate to the edges
284+
a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \
285+
(1.0 / 12.0) * (a[i - 2, j] + a[i + 1, j])
286+
287+
al[i, j] = a_int[i, j]
288+
ar[i, j] = a_int[i, j]
289+
290+
elif (idir == 2):
291+
292+
for i in range(ilo - 1, ihi + 1):
293+
for j in range(jlo - 2, jhi + 3):
294+
295+
# interpolate to the edges
296+
a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \
297+
(1.0 / 12.0) * (a[i, j - 2] + a[i, j + 1])
298+
299+
al[i, j] = a_int[i, j]
300+
ar[i, j] = a_int[i, j]
301+
302+
return al, ar

0 commit comments

Comments
 (0)