temp plotter
This commit is contained in:
Magnus Ulimoen
parent
7f0b5b8430
commit
c1895a6c33
|
|
@ -0,0 +1,264 @@
|
|||
#! /bin/env python3
|
||||
import matplotlib as mpl
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
import os
|
||||
|
||||
from argparse import ArgumentParser
|
||||
|
||||
|
||||
def gridlines(obj, x, y):
|
||||
for j in range(1, x.shape[0] - 1):
|
||||
obj.plot(x[j, :], y[j, :], color="#7f7f7f", linewidth=0.1, alpha=0.3)
|
||||
for j in range(1, x.shape[1] - 1):
|
||||
obj.plot(x[:, j], y[:, j], color="#7f7f7f", linewidth=0.1, alpha=0.3)
|
||||
|
||||
obj.plot(x[0, :], y[0, :], color="#7f7f7f", linewidth=0.2)
|
||||
obj.plot(x[-1, :], y[-1, :], color="#7f7f7f", linewidth=0.2)
|
||||
obj.plot(x[:, 0], y[:, 0], color="#7f7f7f", linewidth=0.2)
|
||||
obj.plot(x[:, -1], y[:, -1], color="#7f7f7f", linewidth=0.2)
|
||||
|
||||
|
||||
def plot_all(grids, error: bool, save: bool, filename="figure.png"):
|
||||
sym_cmap = plt.get_cmap("PiYG") # Symmetric around zero
|
||||
if error:
|
||||
e_cmap = sym_cmap
|
||||
else:
|
||||
e_cmap = plt.get_cmap("Greys")
|
||||
|
||||
f, axarr = plt.subplots(2, 2)
|
||||
|
||||
min_rho = min(np.min(g["rho"]) for g in grids)
|
||||
max_rho = max(np.max(g["rho"]) for g in grids)
|
||||
if error:
|
||||
r = 1.2 * max(abs(min_rho), abs(max_rho))
|
||||
rho_levels = np.linspace(-r, r, 34)
|
||||
else:
|
||||
r = 1.2 * max(abs(min_rho - 1), abs(max_rho - 1))
|
||||
rho_levels = np.linspace(1 - r, 1 + r, 34)
|
||||
|
||||
min_rhou = min(np.min(g["rhou"]) for g in grids)
|
||||
max_rhou = max(np.max(g["rhov"]) for g in grids)
|
||||
if error:
|
||||
r = 1.2 * max(abs(min_rhou), abs(max_rhou))
|
||||
rhou_levels = np.linspace(-r, r, 20)
|
||||
else:
|
||||
r = 1.2 * max(abs(min_rhou - 1), abs(max_rhou - 1))
|
||||
rhou_levels = np.linspace(1 - r, 1 + r, 20)
|
||||
|
||||
min_rhov = min(np.min(g["rhov"]) for g in grids)
|
||||
max_rhov = max(np.max(g["rhov"]) for g in grids)
|
||||
r = 1.2 * max(abs(min_rhov), abs(max_rhov))
|
||||
rhov_levels = np.linspace(-r, r, 20)
|
||||
|
||||
min_e = min(np.min(g["e"]) for g in grids)
|
||||
max_e = max(np.max(g["e"]) for g in grids)
|
||||
if error:
|
||||
r = max(abs(min_e), abs(max_e))
|
||||
e_levels = np.linspace(-r, r, 20)
|
||||
else:
|
||||
e_levels = np.linspace(min_e, max_e)
|
||||
|
||||
for g in grids:
|
||||
x = g["x"]
|
||||
y = g["y"]
|
||||
axarr[0, 0].contourf(x, y, g["rho"], cmap=sym_cmap, levels=rho_levels)
|
||||
gridlines(axarr[0, 0], x, y)
|
||||
|
||||
axarr[0, 1].contourf(x, y, g["rhou"], cmap=sym_cmap, levels=rhou_levels)
|
||||
gridlines(axarr[0, 1], x, y)
|
||||
|
||||
axarr[1, 0].contourf(x, y, g["rhov"], cmap=sym_cmap, levels=rhov_levels)
|
||||
gridlines(axarr[1, 0], x, y)
|
||||
|
||||
axarr[1, 1].contourf(x, y, g["e"], cmap=e_cmap, levels=e_levels)
|
||||
gridlines(axarr[1, 1], x, y)
|
||||
|
||||
axarr[0, 0].set_title(r"$\rho$")
|
||||
axarr[0, 0].set_xlabel("x")
|
||||
axarr[0, 0].set_ylabel("y")
|
||||
norm = mpl.colors.Normalize(vmin=rho_levels[0], vmax=rho_levels[-1])
|
||||
sm = plt.cm.ScalarMappable(cmap=sym_cmap, norm=norm)
|
||||
sm.set_array([])
|
||||
plt.colorbar(sm, ax=axarr[0, 0])
|
||||
|
||||
axarr[0, 1].set_title(r"$\rho u$")
|
||||
axarr[0, 1].set_xlabel("x")
|
||||
axarr[0, 1].set_ylabel("y")
|
||||
norm = mpl.colors.Normalize(vmin=rhou_levels[0], vmax=rhou_levels[-1])
|
||||
sm = plt.cm.ScalarMappable(cmap=sym_cmap, norm=norm)
|
||||
sm.set_array([])
|
||||
plt.colorbar(sm, ax=axarr[0, 1])
|
||||
|
||||
axarr[1, 0].set_title(r"$\rho v$")
|
||||
axarr[1, 0].set_xlabel("x")
|
||||
axarr[1, 0].set_ylabel("y")
|
||||
norm = mpl.colors.Normalize(vmin=rhov_levels[0], vmax=rhov_levels[-1])
|
||||
sm = plt.cm.ScalarMappable(cmap=sym_cmap, norm=norm)
|
||||
sm.set_array([])
|
||||
plt.colorbar(sm, ax=axarr[1, 0])
|
||||
|
||||
axarr[1, 1].set_title(r"$e$")
|
||||
axarr[1, 1].set_xlabel("x")
|
||||
axarr[1, 1].set_ylabel("y")
|
||||
norm = mpl.colors.Normalize(vmin=e_levels[0], vmax=e_levels[-1])
|
||||
sm = plt.cm.ScalarMappable(cmap=e_cmap, norm=norm)
|
||||
sm.set_array([])
|
||||
plt.colorbar(sm, ax=axarr[1, 1])
|
||||
|
||||
if save:
|
||||
plt.savefig(filename, bbox_inches="tight", dpi=600)
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
def plot_total_error(grids, save: bool, filename="figure.png"):
|
||||
cmap = plt.get_cmap("Greys")
|
||||
|
||||
total_err = [
|
||||
np.abs(g["rho"]) + np.abs(g["rhou"]) + np.abs(g["rhov"]) + np.abs(g["e"])
|
||||
for g in grids
|
||||
]
|
||||
|
||||
r = max(np.max(err) for err in total_err)
|
||||
|
||||
levels = np.linspace(0, r, 30)
|
||||
|
||||
for g, err in zip(grids, total_err):
|
||||
x = g["x"]
|
||||
y = g["y"]
|
||||
|
||||
plt.contourf(x, y, err, cmap=cmap, levels=levels)
|
||||
gridlines(plt, x, y)
|
||||
|
||||
plt.title("Total error")
|
||||
norm = mpl.colors.Normalize(vmin=levels[0], vmax=levels[-1])
|
||||
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
|
||||
sm.set_array([])
|
||||
plt.colorbar(sm)
|
||||
|
||||
plt.xlabel("x")
|
||||
plt.ylabel("y")
|
||||
|
||||
if save:
|
||||
plt.savefig(args.output, bbox_inches="tight", dpi=600)
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
def plot_pressure(grids, save: bool, filename="figure.png"):
|
||||
cmap = plt.get_cmap("RdGy")
|
||||
gamma = 1.4 # Assumption might be wrong
|
||||
Mach = 0.5
|
||||
|
||||
p = [
|
||||
(gamma - 1) * (g["e"] - (g["rhou"] ** 2 + g["rhov"] ** 2) / (2 * g["rho"]))
|
||||
for g in grids
|
||||
]
|
||||
|
||||
flat_p = np.array([])
|
||||
for p_ in p:
|
||||
flat_p = np.append(flat_p, p_)
|
||||
|
||||
max_p = np.max(flat_p)
|
||||
min_p = np.min(flat_p)
|
||||
|
||||
p_inf = 1 / (gamma * Mach ** 2)
|
||||
|
||||
r = max(max_p - p_inf, p_inf - min_p)
|
||||
|
||||
levels = np.linspace(p_inf - r, p_inf + r, 30)
|
||||
|
||||
for g, p_ in zip(grids, p):
|
||||
x = g["x"]
|
||||
y = g["y"]
|
||||
|
||||
plt.contourf(x, y, p_, cmap=cmap, levels=levels)
|
||||
gridlines(plt, x, y)
|
||||
|
||||
plt.title("Pressure")
|
||||
norm = mpl.colors.Normalize(vmin=levels[0], vmax=levels[-1])
|
||||
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
|
||||
sm.set_array([])
|
||||
plt.colorbar(sm)
|
||||
|
||||
plt.xlabel("x")
|
||||
plt.ylabel("y")
|
||||
|
||||
if save:
|
||||
plt.savefig(filename, bbox_inches="tight", dpi=600)
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
def read_from_file(filename):
|
||||
grids = []
|
||||
|
||||
with open(filename, "rb") as f:
|
||||
ngrids = int(np.fromfile(f, dtype=np.uint32, count=1))
|
||||
for i in range(ngrids):
|
||||
(neta, nxi) = np.fromfile(f, dtype=np.uint32, count=2)
|
||||
|
||||
x = np.fromfile(f, dtype=np.double, count=neta * nxi)
|
||||
x = x.reshape((neta, nxi))
|
||||
|
||||
y = np.fromfile(f, dtype=np.double, count=neta * nxi)
|
||||
y = y.reshape((neta, nxi))
|
||||
|
||||
rho = np.fromfile(f, dtype=np.double, count=neta * nxi)
|
||||
rho = rho.reshape((neta, nxi))
|
||||
|
||||
rhou = np.fromfile(f, dtype=np.double, count=neta * nxi)
|
||||
rhou = rhou.reshape((neta, nxi))
|
||||
|
||||
rhov = np.fromfile(f, dtype=np.double, count=neta * nxi)
|
||||
rhov = rhov.reshape((neta, nxi))
|
||||
|
||||
e = np.fromfile(f, dtype=np.double, count=neta * nxi)
|
||||
e = e.reshape((neta, nxi))
|
||||
|
||||
grids.append(
|
||||
{"x": x, "y": y, "rho": rho, "rhou": rhou, "rhov": rhov, "e": e}
|
||||
)
|
||||
|
||||
return grids
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
parser = ArgumentParser(description="Plot a solution from the eulersolver")
|
||||
parser.add_argument("filename", metavar="filename", type=str)
|
||||
parser.add_argument(
|
||||
"-e",
|
||||
help="Scale is centered around zero (implies -a)",
|
||||
action="store_true",
|
||||
dest="error",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-te", help="Plots total error", action="store_true", dest="total_error"
|
||||
)
|
||||
parser.add_argument("-s", help="Save figure", action="store_true", dest="save")
|
||||
parser.add_argument(
|
||||
"-o",
|
||||
help="Output of saved figure",
|
||||
type=str,
|
||||
default="figure.png",
|
||||
dest="output",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-a", help="Show all four variables", action="store_true", dest="all"
|
||||
)
|
||||
|
||||
args = parser.parse_args()
|
||||
filename = args.filename
|
||||
if not os.path.isfile(filename):
|
||||
filename = "solution{:03}.bin".format(int(filename))
|
||||
|
||||
grids = read_from_file(filename)
|
||||
|
||||
if args.all or args.error:
|
||||
plot_all(grids, args.error, args.save, args.output)
|
||||
elif args.total_error:
|
||||
plot_total_error(grids, args.save, args.output)
|
||||
else:
|
||||
plot_pressure(grids, args.save, args.output)
|
||||
|
|
@ -118,6 +118,7 @@ impl<T: operators::UpwindOperator> System<T> {
|
|||
*y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
|
||||
});
|
||||
}
|
||||
std::mem::swap(&mut self.fnext, &mut self.fnow);
|
||||
return;
|
||||
}
|
||||
_ => {
|
||||
|
|
@ -127,29 +128,29 @@ impl<T: operators::UpwindOperator> System<T> {
|
|||
|
||||
let bt = vec![
|
||||
euler::BoundaryTerms {
|
||||
north: self.fnext[1].north(),
|
||||
south: self.fnext[1].south(),
|
||||
north: self.fnext[1].south(),
|
||||
south: self.fnext[1].north(),
|
||||
east: self.fnext[3].west(),
|
||||
west: self.fnext[3].east(),
|
||||
},
|
||||
euler::BoundaryTerms {
|
||||
north: self.fnext[0].north(),
|
||||
south: self.fnext[0].south(),
|
||||
north: self.fnext[0].south(),
|
||||
south: self.fnext[0].north(),
|
||||
east: self.fnext[2].west(),
|
||||
west: self.fnext[2].east(),
|
||||
},
|
||||
euler::BoundaryTerms {
|
||||
north: self.fnext[3].north(),
|
||||
south: self.fnext[3].south(),
|
||||
east: self.fnext[0].west(),
|
||||
west: self.fnext[0].east(),
|
||||
},
|
||||
euler::BoundaryTerms {
|
||||
north: self.fnext[2].north(),
|
||||
south: self.fnext[2].south(),
|
||||
north: self.fnext[3].south(),
|
||||
south: self.fnext[3].north(),
|
||||
east: self.fnext[1].west(),
|
||||
west: self.fnext[1].east(),
|
||||
},
|
||||
euler::BoundaryTerms {
|
||||
north: self.fnext[2].south(),
|
||||
south: self.fnext[2].north(),
|
||||
east: self.fnext[0].west(),
|
||||
west: self.fnext[0].east(),
|
||||
},
|
||||
];
|
||||
|
||||
for ((((prev, fut), grid), wb), bt) in self
|
||||
|
|
@ -166,7 +167,7 @@ impl<T: operators::UpwindOperator> System<T> {
|
|||
}
|
||||
}
|
||||
|
||||
fn mesh(x: (f64, f64, usize), y: (f64, f64, usize)) -> (Array2<f64>, Array2<f64>) {
|
||||
fn mesh(y: (f64, f64, usize), x: (f64, f64, usize)) -> (Array2<f64>, Array2<f64>) {
|
||||
let arrx = Array1::linspace(x.0, x.1, x.2);
|
||||
let arry = Array1::linspace(y.0, y.1, y.2);
|
||||
|
||||
|
|
@ -174,25 +175,26 @@ fn mesh(x: (f64, f64, usize), y: (f64, f64, usize)) -> (Array2<f64>, Array2<f64>
|
|||
let mut gy = arry.broadcast((x.2, y.2)).unwrap();
|
||||
gy.swap_axes(0, 1);
|
||||
|
||||
(gx.into_owned(), gy.into_owned())
|
||||
(gy.into_owned(), gx.into_owned())
|
||||
}
|
||||
|
||||
fn main() {
|
||||
let n = 20;
|
||||
let nx = 32;
|
||||
let ny = 31;
|
||||
|
||||
let mut grids = Vec::with_capacity(4);
|
||||
|
||||
let (x0, y0) = mesh((-5.0, 0.0, n), (0.0, 5.0, n));
|
||||
grids.push(grid::Grid::<operators::Upwind4>::new(x0, y0).unwrap());
|
||||
let (y0, x0) = mesh((0.0, 5.0, ny), (-5.0, 0.0, nx));
|
||||
grids.push(grid::Grid::<operators::Upwind9>::new(x0, y0).unwrap());
|
||||
|
||||
let (x1, y1) = mesh((-5.0, 0.0, n), (-5.0, 0.0, n));
|
||||
grids.push(grid::Grid::<operators::Upwind4>::new(x1, y1).unwrap());
|
||||
let (y1, x1) = mesh((-5.0, 0.0, ny), (-5.0, 0.0, nx));
|
||||
grids.push(grid::Grid::<operators::Upwind9>::new(x1, y1).unwrap());
|
||||
|
||||
let (x2, y2) = mesh((0.0, 5.0, n), (-5.0, 0.0, n));
|
||||
grids.push(grid::Grid::<operators::Upwind4>::new(x2, y2).unwrap());
|
||||
let (y2, x2) = mesh((-5.0, 0.0, ny), (0.0, 5.0, nx));
|
||||
grids.push(grid::Grid::<operators::Upwind9>::new(x2, y2).unwrap());
|
||||
|
||||
let (x3, y3) = mesh((0.0, 5.0, n), (0.0, 5.0, n));
|
||||
grids.push(grid::Grid::<operators::Upwind4>::new(x3, y3).unwrap());
|
||||
let (y3, x3) = mesh((0.0, 5.0, ny), (0.0, 5.0, nx));
|
||||
grids.push(grid::Grid::<operators::Upwind9>::new(x3, y3).unwrap());
|
||||
|
||||
let mut sys = System::new(grids);
|
||||
sys.vortex(
|
||||
|
|
@ -205,32 +207,42 @@ fn main() {
|
|||
eps: 1.0,
|
||||
},
|
||||
);
|
||||
sys.advance(0.05);
|
||||
let t = 0.2;
|
||||
let dt = 0.1 / (std::cmp::max(nx, ny) as Float);
|
||||
for _ in 0..(t / dt) as _ {
|
||||
sys.advance(dt);
|
||||
}
|
||||
|
||||
/*
|
||||
let bt0 = euler::BoundaryTerms {
|
||||
north: sys1.field().north(),
|
||||
south: sys1.field().south(),
|
||||
east: sys3.field().west(),
|
||||
west: sys3.field().east(),
|
||||
};
|
||||
let bt1 = euler::BoundaryTerms {
|
||||
north: sys0.field().north(),
|
||||
south: sys0.field().south(),
|
||||
east: sys2.field().west(),
|
||||
west: sys2.field().east(),
|
||||
};
|
||||
let bt2 = euler::BoundaryTerms {
|
||||
north: sys3.field().north(),
|
||||
south: sys3.field().south(),
|
||||
east: sys0.field().west(),
|
||||
west: sys0.field().east(),
|
||||
};
|
||||
let bt3 = euler::BoundaryTerms {
|
||||
north: sys2.field().north(),
|
||||
south: sys2.field().south(),
|
||||
east: sys1.field().west(),
|
||||
west: sys1.field().east(),
|
||||
};
|
||||
*/
|
||||
println!("{}", sys.fnow[0].e()[(0, 0)]);
|
||||
dump_to_file(&sys);
|
||||
}
|
||||
|
||||
fn dump_to_file<T: sbp::operators::UpwindOperator>(sys: &System<T>) {
|
||||
use std::io::prelude::*;
|
||||
let file = std::fs::File::create("output").unwrap();
|
||||
let mut file = std::io::BufWriter::new(file);
|
||||
let ngrids = sys.grids.len();
|
||||
file.write_all(&(ngrids as u32).to_le_bytes()).unwrap();
|
||||
for (grid, s) in sys.grids.iter().zip(&sys.fnow) {
|
||||
file.write_all(&(grid.ny() as u32).to_le_bytes()).unwrap();
|
||||
file.write_all(&(grid.nx() as u32).to_le_bytes()).unwrap();
|
||||
for x in grid.x().as_slice().unwrap() {
|
||||
file.write_all(&(x.to_le_bytes())).unwrap();
|
||||
}
|
||||
for y in grid.y().as_slice().unwrap() {
|
||||
file.write_all(&(y.to_le_bytes())).unwrap();
|
||||
}
|
||||
for rho in s.rho().as_slice().unwrap() {
|
||||
file.write_all(&(rho.to_le_bytes())).unwrap();
|
||||
}
|
||||
for rhou in s.rhou().as_slice().unwrap() {
|
||||
file.write_all(&(rhou.to_le_bytes())).unwrap();
|
||||
}
|
||||
for rhov in s.rhov().as_slice().unwrap() {
|
||||
file.write_all(&(rhov.to_le_bytes())).unwrap();
|
||||
}
|
||||
for e in s.e().as_slice().unwrap() {
|
||||
file.write_all(&(e.to_le_bytes())).unwrap();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue