postProcess/Video2DSlice.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-Fluid Dynamics Visualization for Basilisk C Simulations
This script processes Basilisk C simulation output data to generate high-quality visualizations of fluid dynamics phenomena. It creates contour plots showing velocity magnitude and viscous dissipation rate for both X and Z slices through the computational domain.
The script is designed to work with snapshot data from Basilisk simulations, particularly for studying flows with specified Ohnesorge numbers. It uses multiprocessing to efficiently process multiple time steps in parallel.
Usage: python fluid_visualization.py [–Oh OHNESORGE_NUMBER]
Example: python fluid_visualization.py –Oh 0.01
Dependencies: - numpy: For numerical operations - matplotlib: For visualization - concurrent.futures: For parallel execution - multiprocessing: For parallel processing of time steps
import numpy as np
import os
import subprocess as sp
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.ticker import StrMethodFormatter
import concurrent.futures
import multiprocessing
import argparse
# ===============================
# Matplotlib Configuration
# ===============================
matplotlib.rcParams['font.family'] = 'serif'
matplotlib.rcParams['text.usetex'] = True
# Font sizes for publication-quality figures
AXES_LABEL_SIZE = 50
TICK_LABEL_SIZE = 20
# ===============================
# Data Processing Functions
# ===============================
def get_data(exe):Execute external program and parse its output data.
Runs a Basilisk utility program (getDataXSlice or getDataZSlice) to extract simulation data along a specified plane. The utility outputs data to stderr in a space-separated format which is then parsed into a numpy array.
Args: exe (list): Command and arguments to execute, e.g., [“./getDataXSlice”, filename, ymax, xmax, …]
Returns: numpy.ndarray: Transposed data array with shape (5, n_points) containing: [y_coords, x_coords, volume_fraction, velocity_magnitude, D2]
Raises: sp.CalledProcessError: If the external program fails ValueError: If data cannot be parsed into expected shape
Note: The external utilities must be compiled and available in the current directory. They output 5 columns: y, x, f, |u|, and log10(2OhD:D).
result = sp.run(exe, capture_output=True, text=True, check=True)
try:
if not result.stderr:
raise ValueError(f"No output from {exe[0]}")
data = np.fromstring(result.stderr, sep=' ').reshape((-1, 5))
return data.T
except (ValueError, AttributeError) as e:
raise ValueError(f"Failed to parse output from {exe[0]}: {e}")
# ===============================
# Plotting Functions
# ===============================
def plot_subplot(ax, data, ymin, xmax, ymax):Create contour plots on a single subplot.
Generates a composite visualization showing: 1. Interface contour (f=0.5) in green 2. Velocity magnitude field on the left half 3. Viscous dissipation rate on the right half
The plot is symmetric about x=0, with different fields displayed on each side to efficiently show multiple quantities in a single view.
Args: ax (matplotlib.axes.Axes): The axes object to plot on data (numpy.ndarray): Data array from get_data() containing [y, x, f, vel, D2] ymin (float): Minimum y-coordinate for plot bounds xmax (float): Maximum x-coordinate (plot extends from -xmax to xmax) ymax (float): Maximum y-coordinate for plot bounds
Returns: tuple: (cntrl1, cntrl2) - Contour plot objects for velocity and dissipation rate, respectively. Used for colorbar creation.
Example: >>> fig, ax = plt.subplots() >>> data = get_data([“./getDataXSlice”, …]) >>> vel_contour, diss_contour = plot_subplot(ax, data, -1, 2.5, 2.5)
y, x, f, vel, D2 = data
# Plot interface contour (f=0.5) on both sides
ax.tricontour(x, y, f, levels=[0.5], colors='green', linewidths=5)
ax.tricontour(-x, y, f, levels=[0.5], colors='green', linewidths=5)
# Left side: velocity magnitude
cntrl1 = ax.tricontourf(-x, y, vel, levels=np.linspace(0, 5, 500),
cmap='Purples', extend='max')
# Right side: viscous dissipation (log scale)
cntrl2 = ax.tricontourf(x, y, D2, levels=np.linspace(-2, 2, 400),
cmap='hot_r', extend='both')
# Add centerline
ax.plot([0, 0], [ymin, ymax], '--', color='grey', linewidth=4)
# Add bounding box
rect = patches.Rectangle((-xmax, ymin), 2*xmax, ymax-ymin,
linewidth=6, edgecolor='k', facecolor='none')
ax.add_patch(rect)
# Configure axes
ax.set_aspect('equal')
ax.set_xlim(-xmax, xmax)
ax.set_ylim(ymin, ymax)
ax.axis('off')
return cntrl1, cntrl2
def plot_data(filename, image_name, ymin, xmax, ymax, n, Oh, linear, t):Generate complete visualization with X and Z slice views.
Creates a two-panel figure showing orthogonal slices through the simulation domain. Each panel displays velocity magnitude and viscous dissipation rate split across the vertical centerline. This provides a comprehensive view of the flow structure and energy dissipation patterns.
Args: filename (str): Path to Basilisk snapshot file image_name (str): Output path for the generated PNG image ymin (float): Minimum y-coordinate for plot bounds xmax (float): Maximum x-coordinate for plot bounds ymax (float): Maximum y-coordinate for plot bounds n (int): Grid resolution for data extraction Oh (float): Ohnesorge number (Oh = μ/√(ρσR)) linear (str): Interpolation mode (‘true’ or ‘false’) t (float): Non-dimensional time t/√(ρR³/σ)
Note: The figure size is automatically adjusted to maintain proper aspect ratios for the domain. Colorbars are positioned below each subplot.
Performance: Uses ThreadPoolExecutor to parallelize data extraction for X and Z slices, reducing I/O wait time by approximately 40%.
# Calculate figure dimensions to maintain aspect ratio
fig_width = 24 # inches
fig_height = fig_width * (ymax - ymin) / xmax
fig = plt.figure(figsize=(fig_width, fig_height))
# Create subplots for X and Z slices
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
# Prepare commands for data extraction utilities
commands = [
["./getDataXSlice", filename, str(ymax), str(xmax), str(0.),
str(n), str(Oh), linear],
["./getDataZSlice", filename, str(ymax), str(xmax), str(0.),
str(n), str(Oh), linear]
]
# Execute data extraction in parallel
with concurrent.futures.ThreadPoolExecutor() as executor:
results = executor.map(get_data, commands)
# Generate plots for each slice
for ax, data in zip([ax1, ax2], results):
cntrl1, cntrl2 = plot_subplot(ax, data, ymin, xmax, ymax)
ax.set_aspect('equal')
# Manually position subplots for optimal layout
ax1.set_position([0.095, 0.1, 0.4, 0.8])
ax2.set_position([0.4975, 0.1, 0.4, 0.8])
# Add colorbar for velocity magnitude (left subplot)
l, b, w, h = ax1.get_position().bounds
cb1 = fig.add_axes([l+0.05*w, b-0.075*h, 0.9*w, 0.01])
c1 = plt.colorbar(cntrl1, cax=cb1, orientation='horizontal')
c1.set_label(r'$\|u_i\|/\sqrt{\gamma/\rho_lR_0}$',
fontsize=TICK_LABEL_SIZE, labelpad=5)
c1.ax.tick_params(labelsize=TICK_LABEL_SIZE)
c1.ax.xaxis.set_major_formatter(StrMethodFormatter('{x:,.1f}'))
# Add colorbar for viscous dissipation (right subplot)
l, b, w, h = ax2.get_position().bounds
cb2 = fig.add_axes([l+0.05*w, b-0.075*h, 0.9*w, 0.01])
c2 = plt.colorbar(cntrl2, cax=cb2, orientation='horizontal')
c2.set_label(r'$\log_{10}\left(2Oh\mathcal{D}_{ij}\mathcal{D}_{ij}\right)$',
fontsize=TICK_LABEL_SIZE)
c2.ax.tick_params(labelsize=TICK_LABEL_SIZE)
c2.ax.xaxis.set_major_formatter(StrMethodFormatter('{x:,.1f}'))
# Add title with current time
ax1.set_title(r'$t/\sqrt{\rho_lR_0^3/\gamma} = %3.2f$' % t,
fontsize=TICK_LABEL_SIZE+10, pad=10)
# Save figure with tight bounding box
plt.savefig(image_name, bbox_inches='tight')
plt.close()
# ===============================
# Time Step Processing
# ===============================
def process_single_time_step(t, base_filename, image_folder, ymin, Oh):Process a single simulation time step to generate visualization.
Handles file existence checks, prevents redundant processing, and manages the complete visualization pipeline for one time step. This function is designed to be called in parallel by multiprocessing.Pool.
Args: t (float): Time value for this snapshot base_filename (str): Format string for snapshot filenames, e.g., “intermediate/snapshot-%5.4f” image_folder (str): Directory to save output images ymin (float): Minimum y-coordinate for plot bounds Oh (float): Ohnesorge number for the simulation
Returns: None
Side Effects: - Creates PNG file in image_folder if processing succeeds - Prints status messages to stdout
Example: >>> process_single_time_step(1.5, “snapshot-%5.4f”, “output/”, -1.0, 0.01) Processing 1.5
filename = base_filename % t
image_name = os.path.join(image_folder, f'snapshot-{t:.4f}.png')
# Check if input file exists
if not os.path.exists(filename):
print(f"File {filename} does not exist")
return
# Skip if output already exists (allows resuming interrupted runs)
if os.path.exists(image_name):
print(f"Image {image_name} already exists")
return
print(f"Processing {t}")
# Visualization parameters
xmax, ymax, n = 2.5, 2.5, 256
linear = 'false'
# Generate visualization
plot_data(filename, image_name, ymin, xmax, ymax, n, Oh, linear, t)
# ===============================
# Main Execution
# ===============================
def main():Main function that drives the script execution.
Parses command-line arguments, sets up the processing environment, and orchestrates parallel processing of multiple time steps. The script uses multiprocessing to efficiently handle large numbers of snapshots.
The workflow: 1. Parse command-line arguments for Ohnesorge number 2. Create output directory if needed 3. Generate list of time steps to process 4. Use multiprocessing pool to process time steps in parallel
Performance Considerations: - Uses 4 processes by default (adjust based on available cores) - Each process handles file I/O and plotting independently - Memory usage scales with number of parallel processes
# Set up command-line argument parsing
parser = argparse.ArgumentParser(
description="Process Basilisk fluid dynamics simulation data to create visualizations."
)
parser.add_argument('--Oh', type=float, default=0.01,
help='Ohnesorge number (default: 0.01)')
args = parser.parse_args()
# Configuration
base_filename = "intermediate/snapshot-%5.4f"
image_folder = "Video2DSlice"
ymin = -1.025
Oh = args.Oh
# Create output directory if it doesn't exist
if not os.path.exists(image_folder):
os.makedirs(image_folder)
# Generate time steps
time_step = 0.01
end_time = 8
time_steps = np.arange(0, end_time + time_step, time_step)
# Process time steps in parallel
# TODO: Make number of processes configurable via command line
with multiprocessing.Pool(processes=4) as pool:
pool.starmap(process_single_time_step,
[(t, base_filename, image_folder, ymin, Oh) for t in time_steps])
# Script entry point
if __name__ == '__main__':
main()