pyvista scripts for compositional field visualizations

cookbooks/composition-reaction/doc/plotting_scripts/plot_figure.py

@@ -0,0 +1,130 @@
+import numpy as np
+import pyvista as pv
+import matplotlib.pyplot as plt
+from cmcrameri import cm
+
+# inputs/config
+solution_file = "composition-reaction/output-composition-reaction/solution/solution-00000.pvtu"
+output_png = "composition-reaction/doc/0.png"
+
+c1_name = "C_1"
+c2_name = "C_2"
+temp_name = "T"
+
+# domain
+x0 = 0.0
+x1 = 2.0
+y0 = 0.0
+y1 = 1.0
+
+# grid resolution
+nx = 1600
+n_levels = 160
+
+# turns long sampled data into a 2d grid
+def reshape_data(sampled, name):
+    return np.asarray(sampled[name]).reshape((nx, ny), order="F").T
+
+# load solution and trim to the region of interest
+mesh = pv.read(solution_file)
+mesh = mesh.clip_box(bounds=[x0, x1, y0, y1, 0, 0], invert=False)
+
+width = x1 - x0
+height = y1 - y0
+ny = int(nx * height / width)  # keep pixels square
+
+# sample onto a regular grid
+grid = pv.ImageData()
+grid.dimensions = (nx, ny, 1)
+grid.origin = (x0, y0, 0.0)
+grid.spacing = (width / (nx - 1), height / (ny - 1), 1.0)
+sampled = grid.sample(mesh)
+
+x = np.linspace(x0, x1, nx)
+y = np.linspace(y0, y1, ny)
+X, Y = np.meshgrid(x, y)
+
+# fields
+T = reshape_data(sampled, temp_name)
+C1 = reshape_data(sampled, c1_name)
+C2 = reshape_data(sampled, c2_name)
+valid = reshape_data(sampled, "vtkValidPointMask").astype(bool)
+
+# remove bad points
+T = np.where(valid, T, np.nan)
+C1 = np.where(valid, C1, np.nan)
+C2 = np.where(valid, C2, np.nan)
+
+# plotting
+fig_w = 14
+fig_h = fig_w * height / width
+fig, ax = plt.subplots(figsize=(fig_w, fig_h), dpi=220)
+
+# temperature field
+levels = np.linspace(np.nanmin(T), np.nanmax(T), n_levels)
+cf = ax.contourf(
+    X,
+    Y,
+    T,
+    levels=levels,
+    cmap=cm.vik,
+    zorder=1,
+)
+
+# black line: material from the bottom
+ax.contour(
+    X,
+    Y,
+    C1,
+    levels=[0.5],
+    colors="black",
+    linewidths=1.3,
+    zorder=2,
+)
+
+# white line: material made by reaction
+ax.contour(
+    X,
+    Y,
+    C2,
+    levels=[0.5],
+    colors="white",
+    linewidths=1.1,
+    zorder=3,
+)
+
+ax.set_xlim(x0, x1)
+ax.set_ylim(y0, y1)
+ax.set_aspect("equal")
+
+# axis ticks and labels
+x_ticks = np.linspace(x0, x1, 5)
+y_ticks = np.linspace(y0, y1, 5)
+ax.set_xticks(x_ticks)
+ax.set_yticks(y_ticks)
+ax.set_xticklabels([f"{v:g}" for v in x_ticks], fontsize=9)
+ax.set_yticklabels([f"{v:g}" for v in y_ticks], fontsize=9)
+ax.set_xlabel("X", fontsize=11)
+ax.set_ylabel("Y", fontsize=11)
+ax.tick_params(direction="out", length=4, width=0.8)
+
+# leave room for scalar bar
+fig.subplots_adjust(left=0.06, right=0.90, bottom=0.10, top=0.99)
+
+# scalar bar for temperature
+ticks = np.linspace(np.nanmin(T), np.nanmax(T), 6)
+cax = fig.add_axes([0.92, 0.25, 0.015, 0.50])
+cb = fig.colorbar(cf, cax=cax, ticks=ticks)
+cb.ax.set_title("T", fontsize=9, pad=6)
+cb.ax.set_yticklabels([f"{v:g}" for v in ticks])
+cb.ax.tick_params(labelsize=8)
+
+# save final figure
+plt.savefig(
+    output_png,
+    dpi=300,
+    bbox_inches="tight",
+    pad_inches=0.01,
+    facecolor="white"
+)
+plt.close(fig)

cookbooks/composition_active/doc/plotting_scripts/plot_composition.py

@@ -0,0 +1,98 @@
+import numpy as np
+import pyvista as pv
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+from matplotlib.colors import ListedColormap
+from cmcrameri import cm
+
+# inputs/config
+solution_file = "composition_active/output-composition-active/solution/solution-00000.pvtu"
+output_png    = "composition_active/doc/visit0007.png"
+comp_name     = "C_1"
+
+# domain bounds (model coordinates)
+x0, x1 = 0.0, 2.0
+y0, y1 = 0.0, 1.0
+
+# sampling grid resolution and plot tuning
+nx           = 1600
+n_arrows_x   = 70
+quiver_scale = 25
+
+def reshape_data(sampled, name):
+    # pyvista returns flattened Fortran-order arrays; reshape back to (ny, nx) image layout
+    return np.asarray(sampled[name]).reshape((nx, ny), order="F").T
+
+# build a custom colormap: blue->cyan->green->yellow->red across 0..1
+lower = cm.lapaz(np.linspace(0.15, 1.0, 128))
+upper = cm.lajolla(np.linspace(0.0, 1.0, 128))
+active_cmap = ListedColormap(np.vstack([lower, upper]))
+
+# load solution and trim to the region of interest
+mesh = pv.read(solution_file)
+mesh = mesh.clip_box(bounds=[x0, x1, y0, y1, 0, 0], invert=False)
+
+width, height = x1 - x0, y1 - y0
+ny = int(nx * height / width)  # keep pixels square
+
+# resample the unstructured mesh onto a regular image grid
+grid = pv.ImageData(
+    dimensions=(nx, ny, 1),
+    origin=(x0, y0, 0.0),
+    spacing=(width / (nx-1), height / (ny-1), 1.0),
+)
+sampled = grid.sample(mesh)
+
+# coordinate arrays for plotting/quiver
+x = np.linspace(x0, x1, nx)
+y = np.linspace(y0, y1, ny)
+X, Y = np.meshgrid(x, y)
+
+# extract composition field, validity mask, and velocity components
+C1    = reshape_data(sampled, comp_name)
+valid = reshape_data(sampled, "vtkValidPointMask").astype(bool)
+vel   = np.asarray(sampled["velocity"]).reshape((nx, ny, 3), order="F").transpose(1, 0, 2)
+
+# mask out points outside the original mesh
+C1 = np.where(valid, C1, np.nan)
+U  = np.where(valid, vel[:, :, 0], np.nan)
+V  = np.where(valid, vel[:, :, 1], np.nan)
+
+# imshow renders every pixel directly
+norm = mcolors.Normalize(vmin=0.0, vmax=1.0) # normalize to [0,1]
+rgba = active_cmap(norm(np.clip(C1, 0.0, 1.0)))
+rgba[~valid, 3] = 0.0  # transparent outside domain
+
+# plotting
+fig, ax = plt.subplots(figsize=(16, 16 * height / width), dpi=220)
+ax.set_facecolor("white")
+
+# composition field as an image
+ax.imshow(rgba, origin="lower", extent=[x0, x1, y0, y1],
+          interpolation="bilinear", zorder=1)
+
+# velocity arrows, subsampled on a regular grid so they're not too dense
+step = max(1, nx // n_arrows_x)
+ax.quiver(
+    x[::step], y[::step, None],
+    U[::step, ::step], V[::step, ::step],
+    color="black", angles="xy", scale_units="xy",
+    scale=quiver_scale, width=0.0010, pivot="mid", zorder=2)
+
+# axes formatting
+ax.set(xlim=(x0, x1), ylim=(y0, y1), aspect="equal")
+ax.set_xlabel("X", fontsize=11)
+ax.set_ylabel("Y", fontsize=11)
+ax.tick_params(direction="out", length=4, width=0.8, labelsize=9)
+fig.subplots_adjust(left=0.06, right=0.91, bottom=0.10, top=0.99)
+
+# colorbar on the right side
+cax = fig.add_axes([0.925, 0.25, 0.015, 0.50])
+cb  = fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=active_cmap),
+                    cax=cax, ticks=np.linspace(0, 1, 6))
+cb.ax.set_title("C_1", fontsize=9, pad=6)
+cb.ax.tick_params(labelsize=8)
+
+# save final figure
+plt.savefig(output_png, dpi=300, bbox_inches="tight", pad_inches=0.01, facecolor="white")
+plt.close(fig)

cookbooks/composition_active/doc/plotting_scripts/plot_temperature.py

@@ -0,0 +1,92 @@
+import numpy as np
+import pyvista as pv
+import matplotlib.pyplot as plt
+import matplotlib.colors as mcolors
+from cmcrameri import cm
+
+# inputs/config
+solution_file = "composition_active/output-composition-active/solution/solution-00079.pvtu"
+output_png    = "composition_active/doc/visit0003.png"
+
+# sampling grid resolution and plot tuning
+nx         = 1200
+n_arrows_x = 70
+threshold  = 0.5
+
+def to_2d(sampled, name):
+    # pyvista returns flattened Fortran-order arrays; reshape back to (ny, nx) image layout
+    return np.asarray(sampled[name]).reshape((nx, ny), order="F").T
+
+# load solution and get domain bounds
+mesh = pv.read(solution_file)
+xmin, xmax, ymin, ymax = mesh.bounds[:4]
+width, height = xmax - xmin, ymax - ymin
+ny = int(nx * height / width)  # keep pixels square
+
+# resample the unstructured mesh onto a regular image grid
+grid = pv.ImageData(
+    dimensions=(nx, ny, 1),
+    origin=(xmin, ymin, 0.0),
+    spacing=(width / (nx-1), height / (ny-1), 1.0),
+)
+sampled = grid.sample(mesh)
+
+# mask out points outside the original mesh
+valid = to_2d(sampled, "vtkValidPointMask").astype(bool)
+
+def get_field(name):
+    return np.where(valid, to_2d(sampled, name), np.nan)
+
+# extract fields needed for the plot
+T   = get_field("T")
+C1  = get_field("C_1")
+vel = np.asarray(sampled["velocity"]).reshape((nx, ny, 3), order="F").transpose(1, 0, 2)
+U   = np.where(valid, vel[:, :, 0], np.nan)
+V   = np.where(valid, vel[:, :, 1], np.nan)
+speed = np.sqrt(U**2 + V**2)
+
+# plotting
+fig, ax = plt.subplots(figsize=(16, 16 * height / width), dpi=220)
+extent = [xmin, xmax, ymin, ymax]
+
+# temperature background
+t_img = ax.imshow(T, origin="lower", extent=extent,
+                   cmap=cm.vik, norm=mcolors.Normalize(vmin=np.nanmin(T), vmax=np.nanmax(T)),
+                   interpolation="bilinear", zorder=1)
+
+# velocity arrows colored by speed
+step     = max(1, nx // n_arrows_x)
+spd_norm = mcolors.Normalize(vmin=np.nanmin(speed), vmax=np.nanmax(speed))
+ax.quiver(
+    np.linspace(xmin, xmax, nx)[::step],
+    np.linspace(ymin, ymax, ny)[::step, None],
+    U[::step, ::step], V[::step, ::step], speed[::step, ::step],
+    cmap=cm.batlow, norm=spd_norm,
+    angles="xy", scale_units="xy", scale=32,
+    width=0.0012, pivot="mid", zorder=3)
+
+# overlay C_1 = 0.5 contour as the compositional boundary
+ax.contour(C1, levels=[threshold], colors="black", linewidths=1.5,
+           extent=extent, origin="lower", zorder=4)
+
+# axes formatting
+ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax), aspect="equal")
+ax.set_xlabel("X", fontsize=11)
+ax.set_ylabel("Y", fontsize=11)
+ax.tick_params(direction="out", length=4, width=0.8, labelsize=9)
+fig.subplots_adjust(left=0.06, right=0.86, bottom=0.08, top=0.99)
+
+# temperature colorbar (top right)
+cb1 = fig.colorbar(t_img, cax=fig.add_axes([0.885, 0.56, 0.018, 0.28]))
+cb1.ax.set_title("T", fontsize=9, pad=6)
+cb1.ax.tick_params(labelsize=8)
+
+# speed colorbar (bottom right)
+cb2 = fig.colorbar(plt.cm.ScalarMappable(norm=spd_norm, cmap=cm.batlow),
+                    cax=fig.add_axes([0.885, 0.20, 0.018, 0.28]))
+cb2.ax.set_title("speed", fontsize=9, pad=6)
+cb2.ax.tick_params(labelsize=8)
+
+# save final figure
+plt.savefig(output_png, dpi=300, bbox_inches="tight", pad_inches=0.01, facecolor="white")
+plt.close(fig)