diff --git a/olaflow/config.ini b/olaflow/config.ini index 7917bff..2bac493 100644 --- a/olaflow/config.ini +++ b/olaflow/config.ini @@ -15,3 +15,8 @@ out=out_of t0 = 13900 tf = 14300 x0 = -150 + +[post] +pickle = out_post/model.pickle +x = -50 +z = 5 diff --git a/olaflow/processing/olaflow.py b/olaflow/processing/olaflow.py index 48d9165..2de1bc9 100644 --- a/olaflow/processing/olaflow.py +++ b/olaflow/processing/olaflow.py @@ -1,14 +1,55 @@ +import pathlib import re from fluidfoam import readof +import numpy as np class OFModel: def __init__(self, root): self._root = root + self._fields = {} def read_mesh(self): self._x, self._y, self._z = readof.readmesh(str(self._root)) + self._n = self._x.size + + def read_time(self): + _dirs = np.fromiter( + map( + lambda f: f.name, + filter( + lambda f: re.match(r"^[0-9]+(\.[0-9]+)?$", f.name), + self._root.glob("*"), + ), + ), + dtype="U10", + ) + _t = _dirs.astype(np.half) + _sort = np.argsort(_t) + + self._t_dirs = _dirs[_sort] + self._t = _t[_sort] + return self.t + + def read_field(self, field, t): + if not self._root.joinpath(t, field).exists(): + return np.empty((self._n)) + return readof.readfield(self._root, time_name=t, name=field) + + def read_field_all(self, field): + _shape = ( + (self.t.size, self._n) + if readof.typefield(self._root, time_name=self._t_dirs[-1], name=field) + == "scalar" + else (self.t.size, 3, self._n) + ) + + _field = np.empty(_shape, dtype=np.single) + for _f, _dir in zip(_field, self._t_dirs): + _f[:] = self.read_field(field, _dir) + self.fields[field] = _field + return _field def write_field(self, field, values): with open(self._root.joinpath("0", field), "r") as aw_file: @@ -51,3 +92,26 @@ class OFModel: @property def z(self): return self._z + + @property + def coords(self): + return np.stack((self._x, self._y, self._z), axis=1) + + @property + def t(self): + return self._t + + @property + def fields(self): + return self._fields + + @property + def X(self): + return self._X + + @property + def Z(self): + return self._Z + + def FIELD(self, field): + return np.where(self._C > 0, field[:, C], np.nan) diff --git a/olaflow/processing/pickle.py b/olaflow/processing/pickle.py new file mode 100644 index 0000000..1dcf79f --- /dev/null +++ b/olaflow/processing/pickle.py @@ -0,0 +1,39 @@ +import argparse +import configparser +import logging +import pathlib +import pickle + +import matplotlib.pyplot as plt +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("-c", "--config", default="config.ini") +parser.add_argument("-o", "--output", type=pathlib.Path) +args = parser.parse_args() + +logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) +log = logging.getLogger("ola_post") + +log.info("Starting sws -> olaFlow converter") +config = configparser.ConfigParser() +config.read(args.config) +out = pathlib.Path(config.get("post", "pickle")) +out.parent.mkdir(parents=True, exist_ok=True) + +olaflow_root = args.output +model = OFModel(olaflow_root) +model.read_mesh() +model.read_time() +model.read_field_all("alpha.water") +model.read_field_all("porosity") +model.read_field_all("p") +model.read_field_all("p_rgh") +model.read_field_all("U") + +with out.open("wb") as f: + pickle.dump(model, f) diff --git a/olaflow/processing/post.py b/olaflow/processing/post.py new file mode 100644 index 0000000..277b7b7 --- /dev/null +++ b/olaflow/processing/post.py @@ -0,0 +1,80 @@ +import argparse +import configparser +import logging +import pathlib +import pickle + +import matplotlib.pyplot as plt +import matplotlib.animation as animation +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("-c", "--config", default="config.ini") +args = parser.parse_args() + +logging.basicConfig(level=max((10, 20 - 10 * args.verbose))) +log = logging.getLogger("ola_post") + +log.info("Starting sws -> olaFlow converter") +config = configparser.ConfigParser() +config.read(args.config) +out = pathlib.Path(config.get("post", "pickle")) +out.parent.mkdir(parents=True, exist_ok=True) + +with out.open("rb") as f: + model = pickle.load(f) + +x0 = config.getfloat("post", "x") +z0 = config.getfloat("post", "z") +i0 = np.argmin(np.abs((model.x - x0) + 1j * (model.z - z0))) + +X, Z = np.meshgrid(np.unique(model.x), np.unique(model.z)) + +C = np.where( + (model.x[:, None, None].astype(np.single) == X[None, :, :].astype(np.single)) + & (model.z[:, None, None].astype(np.single) == Z[None, :, :].astype(np.single)) +) + +P = np.full((model.t.size, *X.shape), np.nan) +P[:, C[1], C[2]] = model.fields["porosity"][:, C[0]] +AW = np.full((model.t.size, *X.shape), np.nan) +AW[:, C[1], C[2]] = model.fields["alpha.water"][:, C[0]] + +fig, ax = plt.subplots() +tit = ax.text( + 0.5, + 0.95, + f"t={model.t[0]}s", + horizontalalignment="center", + verticalalignment="top", + transform=ax.transAxes, +) +aw_m = ax.pcolormesh(X, Z, AW[0], vmin=0, vmax=1, cmap="Blues", zorder=1) +ax.pcolormesh( + X, + Z, + P[1], + vmin=0, + vmax=1, + cmap="Greys_r", + alpha=np.nan_to_num(1 - P[1])/2, + zorder=1.1, +) +ax.axhline(4.5, ls="-.", lw=1, c="k", alpha=0.2, zorder=1.2) + + +def anim(i): + tit.set_text(f"t={i[0]}s") + aw_m.set_array(i[1]) + return (aw_m,) + + +fig.colorbar(aw_m) +ax.set(xlabel="x (m)", ylabel="z (m)", aspect="equal", facecolor="#bebebe") +ax.grid(c="k", alpha=0.2) +ani = animation.FuncAnimation(fig, anim, frames=zip(model.t, AW), interval=1 / 25) +plt.show()