Skip to content

Commit

Permalink
apply formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
@gardner48
gardner48 committed Jun 30, 2024
1 parent e737583 commit f604770
Show file tree
Hide file tree
Showing 40 changed files with 3,345 additions and 2,086 deletions.
83 changes: 42 additions & 41 deletions benchmarks/advection_reaction_3D/scripts/compare_error.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,8 @@
import glob
import sys
import matplotlib
matplotlib.use('Agg')

matplotlib.use("Agg")
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import pandas as pd
Expand All @@ -25,75 +26,75 @@
# load pickled data
def load_data(file):
data = np.load(file)
m = data['mesh']
t = data['t']
u = data['u']
v = data['v']
w = data['w']
m = data["mesh"]
t = data["t"]
u = data["u"]
v = data["v"]
w = data["w"]

hx = m[0,1] - m[0,0]
hy = m[1,1] - m[1,0]
hz = m[2,1] - m[2,0]
hx = m[0, 1] - m[0, 0]
hy = m[1, 1] - m[1, 0]
hz = m[2, 1] - m[2, 0]

return { 'm': m, 'h': (hx,hy,hz), 't': t, 'u': u, 'v': v, 'w': w }
return {"m": m, "h": (hx, hy, hz), "t": t, "u": u, "v": v, "w": w}


# grid function norm
def norm_3Dgrid(h, x, q=1):
hx,hy,hz = h
hx, hy, hz = h
s = np.shape(x)
return (hx*hy*hz*np.sum(np.abs(x)**q, axis=(1,2,3)))**(1./q)
return (hx * hy * hz * np.sum(np.abs(x) ** q, axis=(1, 2, 3))) ** (1.0 / q)


# load data files
np111 = load_data('np-111/output-with-h-8.33e-02.npz')
np211 = load_data('np-211/output-with-h-8.33e-02.npz')
np311 = load_data('np-311/output-with-h-8.33e-02.npz')
np131 = load_data('np-131/output-with-h-8.33e-02.npz')
np113 = load_data('np-113/output-with-h-8.33e-02.npz')
np911 = load_data('np-911/output-with-h-8.33e-02.npz')
np111 = load_data("np-111/output-with-h-8.33e-02.npz")
np211 = load_data("np-211/output-with-h-8.33e-02.npz")
np311 = load_data("np-311/output-with-h-8.33e-02.npz")
np131 = load_data("np-131/output-with-h-8.33e-02.npz")
np113 = load_data("np-113/output-with-h-8.33e-02.npz")
np911 = load_data("np-911/output-with-h-8.33e-02.npz")
# np133 = load_data('np-133/output-with-h-8.33e-02.npz')
np313 = load_data('np-313/output-with-h-8.33e-02.npz')
np331 = load_data('np-331/output-with-h-8.33e-02.npz')
np333 = load_data('np-333/output-with-h-8.33e-02.npz')
np313 = load_data("np-313/output-with-h-8.33e-02.npz")
np331 = load_data("np-331/output-with-h-8.33e-02.npz")
np333 = load_data("np-333/output-with-h-8.33e-02.npz")
# np666 = load_data('np-666/output-with-h-8.33e-02.npz')

for component in ['u', 'v', 'w']:
for component in ["u", "v", "w"]:
# Reference solution
ref = np111[component]

# Now compute E(h) = ||U(h) - \bar{U}(h)|| using the grid-function norm
E_np211 = norm_3Dgrid(np211['h'], np211[component] - ref)
E_np311 = norm_3Dgrid(np311['h'], np311[component] - ref)
E_np131 = norm_3Dgrid(np131['h'], np131[component] - ref)
E_np113 = norm_3Dgrid(np113['h'], np113[component] - ref)
E_np911 = norm_3Dgrid(np911['h'], np911[component] - ref)
E_np211 = norm_3Dgrid(np211["h"], np211[component] - ref)
E_np311 = norm_3Dgrid(np311["h"], np311[component] - ref)
E_np131 = norm_3Dgrid(np131["h"], np131[component] - ref)
E_np113 = norm_3Dgrid(np113["h"], np113[component] - ref)
E_np911 = norm_3Dgrid(np911["h"], np911[component] - ref)
# E_np133 = norm_3Dgrid(np133['h'], np133[component] - ref)
E_np313 = norm_3Dgrid(np313['h'], np313[component] - ref)
E_np331 = norm_3Dgrid(np331['h'], np331[component] - ref)
E_np333 = norm_3Dgrid(np333['h'], np333[component] - ref)
E_np313 = norm_3Dgrid(np313["h"], np313[component] - ref)
E_np331 = norm_3Dgrid(np331["h"], np331[component] - ref)
E_np333 = norm_3Dgrid(np333["h"], np333[component] - ref)
# E_np666 = norm_3Dgrid(np666['h'], np666[component] - ref)

# Plot error across time
X, Y = np.meshgrid(np111['m'][0,:], np111['t'])
X, Y = np.meshgrid(np111["m"][0, :], np111["t"])
# fig = plt.figure()
# ax = plt.subplot(311, projection='3d')
# ax.plot_surface(X, Y, np.abs(np911[component][:,:,0,0] - ref[:,:,0,0]))
# ax = plt.subplot(312, projection='3d')
# ax.plot_surface(X, Y, np.abs(np911[component][:,0,:,0] - ref[:,0,:,0]))
# ax = plt.subplot(313, projection='3d')
# ax.plot_surface(X, Y, np.abs(np911[component][:,0,0,:] - ref[:,0,0,:]))
plt.plot(np111['t'], E_np211)
plt.plot(np111['t'], E_np131)
plt.plot(np111['t'], E_np113)
plt.plot(np111['t'], E_np911)
plt.plot(np111["t"], E_np211)
plt.plot(np111["t"], E_np131)
plt.plot(np111["t"], E_np113)
plt.plot(np111["t"], E_np911)
# plt.plot(np111['t'], E_np133)
plt.plot(np111['t'], E_np313)
plt.plot(np111['t'], E_np331)
plt.plot(np111['t'], E_np333)
plt.plot(np111["t"], E_np313)
plt.plot(np111["t"], E_np331)
plt.plot(np111["t"], E_np333)
# plt.plot(np111['t'], E_np666)
# plt.legend(['2 1 1', '3 1 1', '1 3 3', '3 1 3', '3 3 1', '3 3 3', '6 6 6'])
# plt.legend(['3 1 1', '1 3 1', '1 1 3', '9 1 1', '1 3 3', '3 1 3', '3 3 1'])
plt.ylabel('||E(hx,hy,hz)||')
plt.xlabel('time')
plt.savefig('compare-error-plot-%s.png' % component)
plt.ylabel("||E(hx,hy,hz)||")
plt.xlabel("time")
plt.savefig("compare-error-plot-%s.png" % component)
77 changes: 40 additions & 37 deletions benchmarks/advection_reaction_3D/scripts/compute_error.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,8 @@
import glob
import sys
import matplotlib
matplotlib.use('Agg')

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
Expand All @@ -24,65 +25,67 @@
# load pickled data
def load_data(file):
data = np.load(file)
m = data['mesh']
t = data['t']
u = data['u']
v = data['v']
w = data['w']
m = data["mesh"]
t = data["t"]
u = data["u"]
v = data["v"]
w = data["w"]

hx = m[0,1] - m[0,0]
hy = m[1,1] - m[1,0]
hz = m[2,1] - m[2,0]
hx = m[0, 1] - m[0, 0]
hy = m[1, 1] - m[1, 0]
hz = m[2, 1] - m[2, 0]

return { 'm': m, 'h': (hx,hy,hz), 't': t, 'u': u, 'v': v, 'w': w }
return {"m": m, "h": (hx, hy, hz), "t": t, "u": u, "v": v, "w": w}


# grid function norm
def norm_3Dgrid(h, x, q=1):
hx,hy,hz = h
return (hx*hy*hz*np.sum(np.abs(x)**q, axis=(1,2,3)))**(1/q)
hx, hy, hz = h
return (hx * hy * hz * np.sum(np.abs(x) ** q, axis=(1, 2, 3))) ** (1 / q)


# computer order of accuracy p
def calc_order(h1, Eh1, h2, Eh2):
return np.log( Eh1/Eh2 ) / np.log( np.prod(h1)/np.prod(h2) )
return np.log(Eh1 / Eh2) / np.log(np.prod(h1) / np.prod(h2))


# load data files
h_over_8 = load_data('middle-h/output-with-h-1.04e-02.npz')
h_over_4 = load_data('large-h/output-with-h-2.08e-02.npz')
h_over_8 = load_data("middle-h/output-with-h-1.04e-02.npz")
h_over_4 = load_data("large-h/output-with-h-2.08e-02.npz")
# h_over_2 = load_data('larger-h/output-with-h-4.16e-02.npz')
h_over_1 = load_data('largest-h/output-with-h-8.33e-02.npz')
h_over_1 = load_data("largest-h/output-with-h-8.33e-02.npz")

for component in ['u', 'v', 'w']:
for component in ["u", "v", "w"]:
# Restrict reference to the coarsest grid
ref = h_over_8[component][:,::8,::8,::8]
ref = h_over_8[component][:, ::8, ::8, ::8]

# Now compute E(h) = ||U(h) - \bar{U}(h)|| using the grid-function norm
Eh_over_4 = norm_3Dgrid(h_over_4['h'], h_over_4[component][:,::4,::4,::4] - ref)
Eh_over_1 = norm_3Dgrid(h_over_1['h'], h_over_1[component][:,:,:,:] - ref)
Eh_over_4 = norm_3Dgrid(h_over_4["h"], h_over_4[component][:, ::4, ::4, ::4] - ref)
Eh_over_1 = norm_3Dgrid(h_over_1["h"], h_over_1[component][:, :, :, :] - ref)

# Compute order p as in O(h^p)
p = calc_order(h_over_1['h'], Eh_over_1, h_over_4['h'], Eh_over_4)
print('min p for %s component: %.4f' % (component, np.min(p)))
p = calc_order(h_over_1["h"], Eh_over_1, h_over_4["h"], Eh_over_4)
print("min p for %s component: %.4f" % (component, np.min(p)))

# Plot error across time
plt.figure()
plt.plot(h_over_8['t'], Eh_over_4, 'r-')
plt.plot(h_over_8['t'], Eh_over_1, 'b-')
plt.ylabel('||E(hx,hy,hz)||')
plt.xlabel('time')
plt.savefig('error-in-time-plot-%s.png' % component)
plt.plot(h_over_8["t"], Eh_over_4, "r-")
plt.plot(h_over_8["t"], Eh_over_1, "b-")
plt.ylabel("||E(hx,hy,hz)||")
plt.xlabel("time")
plt.savefig("error-in-time-plot-%s.png" % component)

# Plot error norm with respect to h
plt.figure()
x = np.array([np.prod(h_over_4['h']), np.prod(h_over_1['h'])])
plt.plot(x, x, 'k-')
plt.plot(x, x**2, 'k-')
plt.plot(x, [np.linalg.norm(Eh_over_4, np.Inf), np.linalg.norm(Eh_over_1, np.Inf)], 'r-')
plt.legend(['1st order', '2nd order', 'actual'])
plt.ylabel('|| ||E(hx,hy,hz)|| ||_inf')
plt.xlabel('hx * hy * hz')
plt.yscale('log')
plt.xscale('log')
plt.savefig('error-plot-%s.png' % component)
x = np.array([np.prod(h_over_4["h"]), np.prod(h_over_1["h"])])
plt.plot(x, x, "k-")
plt.plot(x, x**2, "k-")
plt.plot(
x, [np.linalg.norm(Eh_over_4, np.Inf), np.linalg.norm(Eh_over_1, np.Inf)], "r-"
)
plt.legend(["1st order", "2nd order", "actual"])
plt.ylabel("|| ||E(hx,hy,hz)|| ||_inf")
plt.xlabel("hx * hy * hz")
plt.yscale("log")
plt.xscale("log")
plt.savefig("error-plot-%s.png" % component)
Loading

0 comments on commit f604770

Please sign in to comment.