-
Notifications
You must be signed in to change notification settings - Fork 0
/
plots.py
182 lines (161 loc) · 5.77 KB
/
plots.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
# Make plots
import matplotlib.pyplot as plt
from argparse import Namespace
import numpy as np
import os
from common import *
from filters import *
# plt.style.use('fivethirtyeight')
plt.style.use('./presentation.mplstyle')
def filterplots():
x = np.linspace(-2, 2, 100)
fig, ax = plt.subplots(2,2, figsize=(14,8))
ax[0][0].plot(x, cesaro(x), label='Cesaro')
ax[0][1].plot(x, raisedcos(x), label="Raised cosine")
ax[1][0].plot(x, lanczos(x), label="Lanczos")
ax[1][1].plot(x, exponential(x, 2), label="Exponential2")
for i in range(2):
for j in range(2):
ax[i][j].legend()
ax[i][j].grid()
fig.savefig("filters.svg")
# Plot all the ggb polys
def ggbplots():
import ggb
x = np.linspace(-1, 1, 2048)
fig, ax = plt.subplots(figsize=(14,8))
l = 3
ax.set_ylim((-20, 20))
for n in range(5):
ax.plot(x, ggb.C(x, n, l), label=f"$k={n}$")
ax.legend()
fig.savefig('ggbs.svg')
# Plot the resolution of the given FFT
def plot_resolution(c, ax, **kwargs):
k = fftshift(freqs(len(c)))
ax.semilogy(k, np.abs(fftshift(c)), **kwargs)
return
def plot_error_fft(uh, exact, ax, **kwargs):
x = exact[0]
if len(uh) < len(x):
u = ifft_at(x, uh)
e = np.abs(u - exact[1])
ax[1].plot(x, e, **kwargs)
else:
print("WARNING: Interpolating exact not supported yet.\nSupply exact solution on finer grid.")
return
def plot_error(c, exact, ax, **kwargs):
print("print_error not implemented yet.")
# Fourier interpolate the IFFT
def smoothplot(v, ax,nn=2048, **plotargs):
n = len(v)
w = pad(v, (nn - n)//2)
dy = (xmax - xmin) / n
dz = (xmax - xmin) / nn
y = cgrid(n)
z = np.linspace(y[0], y[-1]+dy-dz, nn, endpoint=True)
ax.plot(z, ifft(w).real, **plotargs)
return z, ifft(w)
def smooth_and_error(solax, errax, v, exact, nn=2048, **plotargs):
if np.any(exact == None):
print("Exact not found.")
return smoothplot(v, solax, **plotargs)
plot_and_error(solax, errax, exact[0], ifft_at(exact[0], v).real, exact, **plotargs)
def plot_and_error(solax, errax, x, u, exact, **plotargs):
solax.plot(x, u, **plotargs)
if np.any(exact == None):
print("Exact not found.")
return
if len(x) > len(exact[0]):
ei = np.interp(x, exact[0], exact[1])
errax.semilogy(x, np.abs(u-ei), **plotargs)
elif (len(x) < len(exact[0])):
ui = np.interp(exact[0], x, u)
errax.semilogy(exact[0], np.abs(ui-exact[1]), **plotargs)
else:
errax.semilogy(x, np.abs(u-exact[1]), **plotargs)
errax.set_ylim((1e-6, 1e1))
def convergence_plot(exactfile, filenames, saveas, **kwargs):
for fn in filenames:
if not os.path.isfile(fn):
print(fn + " not found or is invalid.")
exit(1)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1])
fig.tight_layout()
ex = np.loadtxt(exact).transpose()
ax[0].plot(ex[0] * 2 * np.pi, ex[1], 'ko', markersize=0.1, label="Exact")
x = np.arange(xmin, xmax, len(ex[0]))
for fn in filenames:
ku = np.loadtxt(fn)
k = ku[0]
u = ku[1]
N = len(k)
ax[0].plot(x, ifft(u[-1]).real, label=str(N))
plot_resolution(u[-1], ax[1], {'linewidth':0.5, 'markersize':0.5})
fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1])
x = np.linspace(xmin, xmax, 1000)
ax[0].plot(x, initial_condition(x), linewidth=0.1, color='k', label="init")
ax[0].legend()
fig.savefig(saveas)
def solplot(sols, args, plotname):
return
nn = 2048
fig, ax = plt.subplots(nrows=1, ncols=2, figsize = (14, 8), width_ratios=[3,1])
# Initial
x = np.linspace(xmin, xmax, 1000)
ax[0].plot(x, initial_condition(x), linewidth=1, color='k', label="Init")
# All the modes
for (t, uf, v) in sols:
n = len(uf)
label = str(n) + f", t={np.round(t, 3)}"
plots.smoothplot(uf, ax[0], label=label, linewidth=1)
plots.plot_resolution(uf, ax[1], linewidth=0.5, markersize=0.5)
if args.ggb:
label += ",ggb"
if sigma != filters.no_filter:
label += "," + args.filter
ax[0].plot(cgrid(n), v, label=label)
elif sigma != filters.no_filter:
label += "," + args.filter
plots.smoothplot(uf, ax[0])
if args.show_markers:
ax[0].plot(cgrid(n), ifft(uf), "+", color= "red", markersize=5)
# Exact
if args.exact != None:
exact = np.loadtxt(args.exact).transpose()
ax[0].plot(exact[0] * 2 * np.pi, exact[1], 'ko', markersize=0.8, label="Exact")
ax[0].grid(visible=True)
ax[0].legend()
fig.tight_layout()
fig.savefig(plotname)
plt.close()
print("Saved plot to "+plotname)
# convergence_plot('a', ['16.txt', '32.txt', '64.txt'])
def aliasing_plots():
N = 4
n = 1
x = np.arange(0, 2*np.pi, 2 * np.pi / N)
xx = np.linspace(0, 2*np.pi, 1000)
y1 = np.sin(n*xx)
y2 = np.sin((n+2*N)*xx)
y3 = np.sin(n*x)
fig, ax = plt.subplots(figsize=(14,8))
ax.plot(xx, y1, label=f"$\sin({n}x)$")
ax.plot(xx, y2, label=f"$\sin({n+2*N}x)$")
ax.plot(x, y3, ls="None", marker='o', markersize=25, alpha=0.5)
ax.legend()
fig.tight_layout()
fig.savefig('aliasing.svg')
plt.close()
if __name__ == "__main__":
aliasing_plots()
ggbplots()
filterplots()
'''
initial_condition = ic.sin
args = Namespace(L=3, N=[16, 64], Tf=1.0, cfl=0.1, exact=None, filter='no_filter', ggb=False, ic='sin', pde='burgers', show_markers=False)
sols, prefix, _ = main.run(args)
# prefix = 'burgers-1.0-sin-no_filter-gFalse'
# sols = [np.loadtxt(file) for file in [prefix+'-16.txt', prefix+'-64']]
solplot(sol, args, initial_condition, prefix+'.svg')
'''