From b2c4227f83ee23e5bacc924b3301b855a2699077 Mon Sep 17 00:00:00 2001 From: "Edgar P. Burkhart" Date: Wed, 6 Jul 2022 07:55:46 +0200 Subject: [PATCH] Olaflow processing for article --- olaflow/processing/flow_velocity.py | 137 ++++++++++++++++++++++++++++ olaflow/processing/snap.py | 79 ++++++++++++++++ 2 files changed, 216 insertions(+) create mode 100644 olaflow/processing/flow_velocity.py create mode 100644 olaflow/processing/snap.py diff --git a/olaflow/processing/flow_velocity.py b/olaflow/processing/flow_velocity.py new file mode 100644 index 0000000..cfd3133 --- /dev/null +++ b/olaflow/processing/flow_velocity.py @@ -0,0 +1,137 @@ +import argparse +import gzip +from itertools import starmap +import logging +from multiprocessing import pool +import pathlib +import pickle +import sys + +from cycler import cycler +import matplotlib.pyplot as plt +from matplotlib.ticker import MultipleLocator +import numpy as np +from scipy import interpolate + +from .olaflow import OFModel + +parser = argparse.ArgumentParser(description="Post-process olaflow results") +parser.add_argument("-v", "--verbose", action="count", default=0) +parser.add_argument( + "-o", + "--output", + type=pathlib.Path, + help="Output directory for pickled data", + required=True, +) +parser.add_argument( + "-f", + "--func", + type=str, + help="Post-process function to compare", + default="graphUniform", + choices=("graphUniform", "graphUniform2"), +) + +args = parser.parse_args() + +logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) +log = logging.getLogger("ola_post") + +log.info("Plotting comparison of model output") + + +def get_pickle(out): + with ( + path.open("rb") + if (path := out.joinpath("pickle")).exists() + else gzip.open(path.with_suffix(".gz"), "rb") + ) as f: + return pickle.load(f) + + +model = get_pickle(args.output) + +figsize = 15 / 2.54, 6 / 2.54 + +fig, ax_ = plt.subplots( + 2, + figsize=figsize, + dpi=200, + constrained_layout=True, +) + +_ax = ax_[0] +v = np.nanmax(np.abs(np.where( + model.post_fields[args.func]["alpha.water"] > 0.5, + #np.linalg.norm(model.post_fields[args.func]["U"], axis=2), + model.post_fields[args.func]["U"][..., 0], + np.nan, +))) +v150 = np.nanmax(np.abs(np.where( + (model.post_fields[args.func]["alpha.water"] > 0.5) & (model.t[:, None] > 170) & (model.t[:, None] < 200), + #np.linalg.norm(model.post_fields[args.func]["U"], axis=2), + model.post_fields[args.func]["U"][..., 0], + np.nan, +))) +_data = model.post_fields[args.func]["U"][..., 0].T +#_c = _ax.contourf( +# model.t, +# model.post_fields[args.func]["x_U"], +# _data, +# cmap="PiYG", +# #levels=[-15, -10, -5, -2, -1, 0, 1, 2, 5, 10, 15], +# vmin=-np.nanmax(np.abs(_data)), +# vmax=np.nanmax(np.abs(_data)), +# extend="both", +#) +_c = _ax.imshow( + _data[::-1], + cmap="PiYG", + alpha=np.clip(model.post_fields[args.func]["alpha.water"], 0, 1).T[::-1], + extent=( + model.t.min(), + model.t.max(), + model.post_fields[args.func]["x_U"].min(), + model.post_fields[args.func]["x_U"].max(), + ), + vmin=-v150, + vmax=v150, + aspect="auto", +) +_ax.set(xlim=(0, 400)) +_ax.set(facecolor="k") +_ax.xaxis.set_minor_locator(MultipleLocator(10)) +_ax.yaxis.set_minor_locator(MultipleLocator(1)) +_ax.set(ylabel="z (m)") +_ax.axes.set_xticklabels([]) +fig.colorbar(_c, label=f"U (m/s)", ax=_ax) +log.info(f"Vitesse max: {v}m/s") +log.info(f"Vitesse max [170,200]: {v150}m/s") +log.info(f"Écart: {abs(np.nanmax(_data)-17.7)/17.7:%}") + + +x = model.post_fields[args.func]["x_U"] +i0 = np.argmin(np.abs(x - 5)) +_data = model.post_fields[args.func]["U"][..., i0, 0] +_alpha = model.post_fields[args.func]["alpha.water"][..., i0] + +ax = ax_[1] +ax.fill_between(model.t, np.where(_alpha > 0.5, _data, 0), lw=1, color="#898989", edgecolor="k") +#ax.autoscale(True, "x", True) +ax.set(xlim=(0, 400)) +ax.set(xlabel="t (s)", ylabel="U (m/s)") +ax.grid(c="k", alpha=.2) +ax.xaxis.set_minor_locator(MultipleLocator(10)) +ax.yaxis.set_minor_locator(MultipleLocator(2)) + +fig.savefig( + args.output.joinpath( + f"U_{args.func}.pdf" + ) +) +fig.savefig( + args.output.joinpath( + f"U_{args.func}.jpg" + ) +) diff --git a/olaflow/processing/snap.py b/olaflow/processing/snap.py new file mode 100644 index 0000000..5628eed --- /dev/null +++ b/olaflow/processing/snap.py @@ -0,0 +1,79 @@ +import argparse +import gzip +import logging +import multiprocessing as mp +import pathlib +import pickle + +import matplotlib.pyplot as plt +import matplotlib.animation as animation +from matplotlib.gridspec import GridSpec +from matplotlib.ticker import MultipleLocator +import numpy as np +from scipy import interpolate + +from .olaflow import OFModel + +parser = argparse.ArgumentParser(description="Post-process olaflow results") +parser.add_argument("-v", "--verbose", action="count", default=0) +parser.add_argument( + "-o", + "--output", + type=pathlib.Path, + help="Output directory for pickled data", + required=True, +) +parser.add_argument( + "-m", + "--max", + help="Only compute maximum rather than animation", + action="store_true", +) +parser.add_argument( + "-i", + "--initial", + help="Only compute initial domain", + action="store_true", +) +args = parser.parse_args() + +logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) +log = logging.getLogger("ola_post") + +log.info("Animating olaFlow output") +out = args.output +out.mkdir(parents=True, exist_ok=True) + +with ( + path.open("rb") + if (path := out.joinpath("pickle")).exists() + else gzip.open(path.with_suffix(".gz"), "rb") +) as f: + model = pickle.load(f) + +x0, idx0 = np.unique(model.x.astype(np.half), return_inverse=True) +z0, idz0 = np.unique(model.z.astype(np.half), return_inverse=True) + +ix0 = np.argsort(x0) +iz0 = np.argsort(z0)[::-1] + +X, Z = np.meshgrid(x0, z0) + +P = np.full((model.t.size, *X.shape), np.nan) +P[:, iz0[idz0], ix0[idx0]] = model.fields["porosity"] + +AW = np.full((model.t.size, *X.shape), np.nan) +AW[:, iz0[idz0], ix0[idx0]] = model.fields["alpha.water"] + +#U = np.full((model.t.size, *X.shape), np.nan) +#U[:, iz0[idz0], ix0[idx0]] = np.linalg.norm(model.fields["U"], axis=1) + +i0 = np.argmin(np.abs(model.t[:, None] - np.asarray((102, 118, 144.5, 176.5))[None, :]), axis=0) + +fig, ax_ = plt.subplots( + 2, 2, figsize=(15 / 2.54, 4 / 2.54), dpi=200, constrained_layout=True +) +for ax, i in zip(ax_.flatten(), i0): + ax.imshow(AW[i], cmap="Blues", vmin=0, vmax=1) + +fig.savefig(out.joinpath("snap.pdf"))