From 903285056a2411715e69c1b7cb323e9ed7b4f549 Mon Sep 17 00:00:00 2001
From: "Edgar P. Burkhart" <git@edgarpierre.fr>
Date: Wed, 20 Apr 2022 13:54:22 +0200
Subject: [PATCH] Velocity zoomed animation

---
 olaflow/processing/animateU.py | 104 +++++++++++++++++++++++++++++++++
 1 file changed, 104 insertions(+)
 create mode 100644 olaflow/processing/animateU.py

diff --git a/olaflow/processing/animateU.py b/olaflow/processing/animateU.py
new file mode 100644
index 0000000..18c7298
--- /dev/null
+++ b/olaflow/processing/animateU.py
@@ -0,0 +1,104 @@
+import argparse
+import configparser
+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
+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("Animating olaFlow output")
+config = configparser.ConfigParser()
+config.read(args.config)
+out = pathlib.Path(config.get("post", "out"))
+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 = config.getfloat("post", "x")
+z0 = config.getfloat("post", "z")
+
+flt = np.where((model.x >= -60) & (model.x <= -20) & (model.z >= 0) & (model.z <= 10))[
+    0
+]
+x0, idx0 = np.unique(model.x[flt].astype(np.half), return_inverse=True)
+z0, idz0 = np.unique(model.z[flt].astype(np.half), return_inverse=True)
+
+X, Z = np.meshgrid(x0, z0)
+
+P = np.full((model.t.size, *X.shape), np.nan)
+P[:, idz0, idx0] = model.fields["porosity"][:, flt]
+
+AW = np.full((model.t.size, *X.shape), np.nan)
+AW[:, idz0, idx0] = model.fields["alpha.water"][:, flt]
+
+watl = z0[np.argmax((AW > 0.5)[:, ::-1, :], axis=1)]
+
+U = np.full((model.t.size, 2, *X.shape), np.nan)
+UU = np.full((model.t.size, *X.shape), np.nan)
+U[..., idz0, idx0] = model.fields["U"][..., flt][:, (0, 2)]
+UU[..., idz0, idx0] = np.linalg.norm(model.fields["U"][..., flt], axis=1)
+
+figU = plt.figure(figsize=(19.2, 10.8), dpi=100)
+gsU = GridSpec(2, 1, figure=figU, height_ratios=[1, 0.05])
+axU = figU.add_subplot(gsU[0])
+caxu1 = figU.add_subplot(gsU[1])
+# caxu2 = figU.add_subplot(gsU[2])
+
+alp = np.clip(np.nan_to_num(AW), 0, 1)
+axU.pcolormesh(X, Z, P[1], vmin=0, vmax=1, cmap="Greys_r")
+u_m = axU.quiver(
+    X,
+    Z,
+    *U[0],
+    UU[0],
+    alpha=alp[0],
+    cmap="spring",
+    clim=(0, np.nanquantile(UU, 0.99)),
+)
+# (wat_p,) = axU.plot(x0, watl[0])
+
+axU.set(xlabel="x (m)", ylabel="z (m)", aspect="equal", facecolor="#bebebe")
+axU.grid(c="k", alpha=0.2)
+titU = axU.text(
+    0.5,
+    0.95,
+    f"t={model.t[0]}s",
+    horizontalalignment="center",
+    verticalalignment="top",
+    transform=axU.transAxes,
+)
+
+figU.colorbar(u_m, label=r"$U$", cax=caxu1, shrink=0.6, orientation="horizontal")
+
+
+def animU(i):
+    titU.set_text(f"t={model.t[i]}s")
+    u_m.set_UVC(*U[i], UU[i])
+    u_m.set_alpha(alp[i])
+    # wat_p.set_ydata(watl[i])
+
+
+aniU = animation.FuncAnimation(figU, animU, frames=model.t.size, interval=1 / 24)
+
+aniU.save(out.joinpath("animUzoom.mp4"), fps=24)