1
Fork 0
internship/olaflow/processing/diff.py

185 lines
5.8 KiB
Python
Raw Normal View History

2022-05-03 11:53:58 +02:00
import argparse
import gzip
from itertools import starmap
import logging
from multiprocessing import pool
import pathlib
import pickle
2022-05-09 11:23:09 +02:00
import sys
2022-05-03 11:53:58 +02:00
from cycler import cycler
import matplotlib.pyplot as plt
2022-06-24 16:50:38 +02:00
from matplotlib.ticker import MultipleLocator
2022-05-03 11:53:58 +02:00
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",
action="append",
type=pathlib.Path,
help="Post-processing directory",
required=True,
)
2022-05-06 11:14:05 +02:00
parser.add_argument(
"-t",
"--timestep",
type=float,
help="Time-step to compare",
)
parser.add_argument(
"-f",
"--func",
type=str,
help="Post-process function to compare",
default="graphUniform",
2022-05-09 11:23:09 +02:00
choices=("graphUniform", "graphUniform2"),
)
parser.add_argument(
"-y",
"--field",
type=str,
help="Field to compare",
default="alpha.water",
choices=("alpha.water", "U"),
2022-05-06 11:14:05 +02:00
)
2022-05-03 11:53:58 +02:00
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)
models = list(map(get_pickle, args.output))
2022-06-24 16:50:38 +02:00
figsize = 15 / 2.54, 4 / 2.54 * len(models)
2022-05-09 11:23:09 +02:00
fig, ax_ = plt.subplots(
2022-05-09 11:23:09 +02:00
len(models),
2022-06-24 16:50:38 +02:00
figsize=figsize,
dpi=200,
2022-05-09 11:23:09 +02:00
constrained_layout=True,
squeeze=False,
)
ax = ax_[:, 0]
2022-05-09 11:23:09 +02:00
2022-05-06 11:14:05 +02:00
if args.timestep is None:
2022-05-09 11:23:09 +02:00
match args.field:
case "alpha.water":
for i, (_ax, _model) in enumerate(zip(ax, models)):
_ax.contour(
_model.t,
_model.post_fields[args.func][f"x_{args.field}"],
_model.post_fields[args.func][args.field].T,
(0.5,),
colors="k",
)
case "U":
for i, (_ax, _model) in enumerate(zip(ax, models)):
2022-06-24 16:50:38 +02:00
v = np.nanmax(np.abs(np.where(
_model.post_fields[args.func]["alpha.water"] > 0.5,
#np.linalg.norm(_model.post_fields[args.func][args.field], axis=2),
_model.post_fields[args.func][args.field][..., 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][args.field], axis=2),
_model.post_fields[args.func][args.field][..., 0],
np.nan,
)))
_data = _model.post_fields[args.func][args.field][..., 0].T
#_c = _ax.contourf(
# _model.t,
# _model.post_fields[args.func][f"x_{args.field}"],
# _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",
#)
2022-05-09 11:23:09 +02:00
_c = _ax.imshow(
2022-06-24 16:50:38 +02:00
_data[::-1],
cmap="PiYG",
alpha=np.clip(_model.post_fields[args.func]["alpha.water"], 0, 1).T[::-1],
2022-05-09 11:23:09 +02:00
extent=(
_model.t.min(),
_model.t.max(),
_model.post_fields[args.func][f"x_{args.field}"].min(),
_model.post_fields[args.func][f"x_{args.field}"].max(),
),
2022-06-24 16:50:38 +02:00
vmin=-v150,
vmax=v150,
aspect="auto",
2022-05-09 11:23:09 +02:00
)
2022-06-24 16:50:38 +02:00
_ax.set(xlim=(100, 300))
_ax.set(facecolor="k")
_ax.xaxis.set_minor_locator(MultipleLocator(5))
_ax.yaxis.set_minor_locator(MultipleLocator(1))
2022-05-09 11:23:09 +02:00
fig.colorbar(_c, label=f"{args.field} (m/s)", ax=_ax)
2022-06-24 16:50:38 +02:00
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:%}")
2022-05-09 11:23:09 +02:00
case _:
log.error(f"Cannot plot field {args.field} from {args.func}")
sys.exit(1)
2022-05-06 11:14:05 +02:00
2022-05-09 11:23:09 +02:00
for i, (_ax, _model) in enumerate(zip(ax, models)):
2022-06-24 16:50:38 +02:00
_ax.set(xlabel="t (s)", ylabel="z (m)")
if len(models) > 1:
_ax.set(title=f"Case {i}")
#_ax.grid(color="#797979", alpha=0.5)
2022-05-06 11:14:05 +02:00
fig.savefig(
args.output[0].joinpath(
2022-05-09 11:23:09 +02:00
f"diff_{args.func}_{args.field}_{'_'.join([o.name for o in args.output])}.pdf"
2022-05-06 11:14:05 +02:00
)
2022-05-03 11:53:58 +02:00
)
2022-06-24 16:50:38 +02:00
fig.savefig(
args.output[0].joinpath(
f"diff_{args.func}_{args.field}_{'_'.join([o.name for o in args.output])}.jpg"
)
)
2022-05-06 11:14:05 +02:00
else:
2022-05-09 11:23:09 +02:00
match args.field:
case "alpha.water":
for i, (_ax, _model) in enumerate(zip(ax, models)):
_ax.tricontour(
_model.x,
_model.z,
2022-06-24 16:50:38 +02:00
_model.fields[args.field][np.where(_model.t == args.timestep)[0]][
0
],
2022-05-09 11:23:09 +02:00
levels=(0.5,),
colors="k",
)
case _:
log.error(f"Cannot plot field {args.field} from {args.func} at timestep")
sys.exit(1)
2022-05-03 11:53:58 +02:00
2022-05-09 11:23:09 +02:00
for i, (_ax, _model) in enumerate(zip(ax, models)):
2022-05-06 11:14:05 +02:00
_ax.set(xlabel="x (m)", ylabel="z (m)", title=f"Case {i}")
_ax.grid()
2022-05-03 11:53:58 +02:00
2022-05-06 11:14:05 +02:00
fig.savefig(
args.output[0].joinpath(
f"diff_t{args.timestep}_{'_'.join([o.name for o in args.output])}.pdf"
)
)