From c1895a6c339d359e686dbfb3bcfb30555d7daa00 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Wed, 1 Apr 2020 22:37:01 +0200 Subject: [PATCH] temp plotter --- sbp/eulerplot | 264 ++++++++++++++++++++++++++++++++++++++ sbp/examples/multigrid.rs | 112 ++++++++-------- 2 files changed, 326 insertions(+), 50 deletions(-) create mode 100755 sbp/eulerplot diff --git a/sbp/eulerplot b/sbp/eulerplot new file mode 100755 index 0000000..a3ceac0 --- /dev/null +++ b/sbp/eulerplot @@ -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) diff --git a/sbp/examples/multigrid.rs b/sbp/examples/multigrid.rs index b0b0806..929af0d 100644 --- a/sbp/examples/multigrid.rs +++ b/sbp/examples/multigrid.rs @@ -118,6 +118,7 @@ impl System { *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 System { 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 System { } } -fn mesh(x: (f64, f64, usize), y: (f64, f64, usize)) -> (Array2, Array2) { +fn mesh(y: (f64, f64, usize), x: (f64, f64, usize)) -> (Array2, Array2) { 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, Array2 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::::new(x0, y0).unwrap()); + let (y0, x0) = mesh((0.0, 5.0, ny), (-5.0, 0.0, nx)); + grids.push(grid::Grid::::new(x0, y0).unwrap()); - let (x1, y1) = mesh((-5.0, 0.0, n), (-5.0, 0.0, n)); - grids.push(grid::Grid::::new(x1, y1).unwrap()); + let (y1, x1) = mesh((-5.0, 0.0, ny), (-5.0, 0.0, nx)); + grids.push(grid::Grid::::new(x1, y1).unwrap()); - let (x2, y2) = mesh((0.0, 5.0, n), (-5.0, 0.0, n)); - grids.push(grid::Grid::::new(x2, y2).unwrap()); + let (y2, x2) = mesh((-5.0, 0.0, ny), (0.0, 5.0, nx)); + grids.push(grid::Grid::::new(x2, y2).unwrap()); - let (x3, y3) = mesh((0.0, 5.0, n), (0.0, 5.0, n)); - grids.push(grid::Grid::::new(x3, y3).unwrap()); + let (y3, x3) = mesh((0.0, 5.0, ny), (0.0, 5.0, nx)); + grids.push(grid::Grid::::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(sys: &System) { + 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(); + } + } }