Olaflow output pickling and animating

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

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)

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)

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()