From a6b0ad179111e2706d1622bf91a5789eb8e64e8d Mon Sep 17 00:00:00 2001 From: Philip Korsager Nickel Date: Tue, 9 Dec 2025 13:35:00 +0100 Subject: [PATCH 1/3] finally it seems to be working! --- conf/experiment/benchmarking/timings.yaml | 4 +- .../validation/convergence/spectral.yaml | 4 +- conf/experiment/validation/ghia/spectral.yaml | 4 +- scripts/plot_lic.py | 461 ++++++++++++++ scripts/test_pyvista_plots.py | 67 ++ src/shared/plotting/ldc/__init__.py | 6 +- src/shared/plotting/ldc/fields.py | 169 ----- src/shared/plotting/ldc/orchestrator.py | 34 +- src/shared/plotting/ldc/pyvista_fields.py | 589 ++++++++++++++++++ src/shared/plotting/ldc/validation.py | 3 +- src/solvers/__init__.py | 8 +- src/solvers/datastructures.py | 2 +- src/solvers/spectral/__init__.py | 4 +- src/solvers/spectral/fmg.py | 45 -- .../spectral/operators/transfer_operators.py | 1 - src/solvers/spectral/vmg.py | 459 -------------- 16 files changed, 1151 insertions(+), 709 deletions(-) create mode 100644 scripts/plot_lic.py create mode 100644 scripts/test_pyvista_plots.py create mode 100644 src/shared/plotting/ldc/pyvista_fields.py delete mode 100644 src/solvers/spectral/fmg.py delete mode 100644 src/solvers/spectral/vmg.py diff --git a/conf/experiment/benchmarking/timings.yaml b/conf/experiment/benchmarking/timings.yaml index 253f7da..bc7988d 100644 --- a/conf/experiment/benchmarking/timings.yaml +++ b/conf/experiment/benchmarking/timings.yaml @@ -21,6 +21,4 @@ hydra: # Sweep over spectral solver variants # sg: Single Grid (no multigrid) # fsg: Full Single Grid - # vmg: V-cycle MultiGrid - # fmg: Full MultiGrid - solver: spectral/sg,spectral/fsg,spectral/vmg,spectral/fmg + solver: spectral/sg,spectral/fsg diff --git a/conf/experiment/validation/convergence/spectral.yaml b/conf/experiment/validation/convergence/spectral.yaml index a907def..0ac085f 100644 --- a/conf/experiment/validation/convergence/spectral.yaml +++ b/conf/experiment/validation/convergence/spectral.yaml @@ -1,6 +1,6 @@ # @package _global_ # Ghia validation - Spectral solvers -# Sweeps over SG, FSG, VMG, FMG solvers at N=15 +# Sweeps over SG, FSG solvers defaults: - override /solver: spectral/sg @@ -15,6 +15,6 @@ hydra: sweeper: params: # Sweep over spectral solver types - solver: spectral/fsg #,spectral/fsg + solver: spectral/fsg Re: 1000 N: 80, 90, 100, 110, 120, 130, 140 diff --git a/conf/experiment/validation/ghia/spectral.yaml b/conf/experiment/validation/ghia/spectral.yaml index ed25fca..cfec235 100644 --- a/conf/experiment/validation/ghia/spectral.yaml +++ b/conf/experiment/validation/ghia/spectral.yaml @@ -1,6 +1,6 @@ # @package _global_ # Ghia validation - Spectral solvers -# Sweeps over SG, FSG, VMG, FMG solvers at N=15 +# Sweeps over SG, FSG solvers defaults: - override /solver: spectral/sg @@ -15,6 +15,6 @@ hydra: sweeper: params: # Sweep over spectral solver types - solver: spectral/fsg #,spectral/vmg #, spectral/vmg + solver: spectral/fsg Re: 400, 1000 N: 12, 14, 16, 18, 20, 24, 28, 32, 36, 40, 46, 50, 56, 60, 64, 68, 74, 82, 86, 90 diff --git a/scripts/plot_lic.py b/scripts/plot_lic.py new file mode 100644 index 0000000..e2d18b5 --- /dev/null +++ b/scripts/plot_lic.py @@ -0,0 +1,461 @@ +"""Line Integral Convolution (LIC) visualization for lid-driven cavity flow. + +LIC convolves vector field orientation with local texture over a kernel length, +preserving all details of the velocity field geometry at pixel scale. + +Usage: + python scripts/plot_lic.py # Uses bundled Re=1000, N=40 solution + python scripts/plot_lic.py --npz path/to/solution.npz + python scripts/plot_lic.py --vts path/to/solution.vts + python scripts/plot_lic.py --run-id + +References: + - Cabral & Leedom (1993): "Imaging Vector Fields Using Line Integral Convolution" + - SciPy Cookbook LIC implementation +""" + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import Normalize +from scipy.ndimage import map_coordinates +from scipy.interpolate import RectBivariateSpline +from pathlib import Path +import sys + +# Default bundled solution (Re=1000, N=40, spectral FSG solver) +DEFAULT_SOLUTION = Path(__file__).parent.parent / "data" / "ldc_solution_Re1000_N40.npz" + + +def generate_lic_fast(u: np.ndarray, v: np.ndarray, length: int = 30) -> np.ndarray: + """Generate Line Integral Convolution texture from velocity field. + + Uses vectorized operations with scipy's map_coordinates for speed. + + Parameters + ---------- + u : np.ndarray + x-component of velocity, shape (ny, nx) + v : np.ndarray + y-component of velocity, shape (ny, nx) + length : int + Kernel length for integration (number of steps in each direction) + + Returns + ------- + np.ndarray + LIC texture, shape (ny, nx) + """ + ny, nx = u.shape + + # Generate white noise texture + np.random.seed(42) # For reproducibility + noise = np.random.rand(ny, nx).astype(np.float32) + + # Normalize velocity to unit vectors (avoid division by zero) + magnitude = np.sqrt(u**2 + v**2) + magnitude = np.maximum(magnitude, 1e-10) + u_norm = (u / magnitude).astype(np.float32) + v_norm = (v / magnitude).astype(np.float32) + + # Create coordinate grids + y_coords, x_coords = np.mgrid[0:ny, 0:nx].astype(np.float32) + + # Accumulator for LIC + lic = np.zeros((ny, nx), dtype=np.float32) + weights = np.zeros((ny, nx), dtype=np.float32) + + # Forward and backward integration + for direction in [1, -1]: + x = x_coords.copy() + y = y_coords.copy() + + for step in range(length): + # Clip to valid range + x_clipped = np.clip(x, 0, nx - 1.001) + y_clipped = np.clip(y, 0, ny - 1.001) + + # Sample noise at current positions + coords = np.array([y_clipped.ravel(), x_clipped.ravel()]) + sampled = map_coordinates(noise, coords, order=1, mode='constant', cval=0) + sampled = sampled.reshape(ny, nx) + + # Check which points are still in bounds + in_bounds = (x >= 0) & (x < nx) & (y >= 0) & (y < ny) + + # Accumulate + lic += sampled * in_bounds + weights += in_bounds.astype(np.float32) + + # Get velocity at current positions + coords_int = np.array([y_clipped.ravel(), x_clipped.ravel()]) + dx = map_coordinates(u_norm, coords_int, order=1, mode='constant', cval=0).reshape(ny, nx) + dy = map_coordinates(v_norm, coords_int, order=1, mode='constant', cval=0).reshape(ny, nx) + + # Step in direction + x += direction * dx * 0.5 + y += direction * dy * 0.5 + + # Normalize + weights = np.maximum(weights, 1) + lic = lic / weights + + return lic + + +def load_solution_from_npz(npz_path: Path) -> tuple: + """Load solution from NPZ file. + + Parameters + ---------- + npz_path : Path + Path to solution.npz file + + Returns + ------- + tuple + (x_nodes, y_nodes, U_2d, V_2d, Re, N) + """ + data = np.load(npz_path) + return ( + data["x_nodes"], + data["y_nodes"], + data["U"], + data["V"], + float(data["Re"]), + int(data["N"]), + ) + + +def load_solution_from_vts(vts_path: Path) -> dict: + """Load solution from VTS file using pyvista. + + Parameters + ---------- + vts_path : Path + Path to solution.vts file + + Returns + ------- + dict + Dictionary with x, y, u, v, p arrays + """ + import pyvista as pv + + grid = pv.read(str(vts_path)) + points = grid.points + + return { + "x": points[:, 0], + "y": points[:, 1], + "u": grid["u"], + "v": grid["v"], + "p": grid["pressure"], + } + + +def spectral_interpolate_to_fine_grid( + x_nodes: np.ndarray, + y_nodes: np.ndarray, + U_2d: np.ndarray, + V_2d: np.ndarray, + n_fine: int = 512, +) -> tuple: + """Interpolate solution from spectral nodes to fine uniform grid. + + Uses the fact that spectral solutions can be evaluated at arbitrary + points via polynomial interpolation (spectrally accurate). + + Parameters + ---------- + x_nodes : np.ndarray + 1D array of x node positions (Chebyshev/Legendre nodes) + y_nodes : np.ndarray + 1D array of y node positions + U_2d : np.ndarray + u velocity on spectral grid, shape (ny, nx) + V_2d : np.ndarray + v velocity on spectral grid, shape (ny, nx) + n_fine : int + Number of points in fine grid + + Returns + ------- + tuple + (x_fine, y_fine, U_fine, V_fine) on uniform grid + """ + # Create fine uniform grid + x_fine = np.linspace(x_nodes.min(), x_nodes.max(), n_fine) + y_fine = np.linspace(y_nodes.min(), y_nodes.max(), n_fine) + + # Use RectBivariateSpline for smooth interpolation + # This is spectrally accurate when source points are spectral nodes + U_spline = RectBivariateSpline(y_nodes, x_nodes, U_2d, kx=3, ky=3) + V_spline = RectBivariateSpline(y_nodes, x_nodes, V_2d, kx=3, ky=3) + + U_fine = U_spline(y_fine, x_fine) + V_fine = V_spline(y_fine, x_fine) + + return x_fine, y_fine, U_fine, V_fine + + +def restructure_to_2d(fields: dict) -> tuple: + """Convert flattened fields to 2D arrays. + + Parameters + ---------- + fields : dict + Dictionary with x, y, u, v arrays (flattened) + + Returns + ------- + tuple + (x_unique, y_unique, U_2d, V_2d) + """ + x, y = fields["x"], fields["y"] + u, v = fields["u"], fields["v"] + + x_unique = np.sort(np.unique(x)) + y_unique = np.sort(np.unique(y)) + nx, ny = len(x_unique), len(y_unique) + + # Build index maps + x_to_idx = {val: i for i, val in enumerate(x_unique)} + y_to_idx = {val: i for i, val in enumerate(y_unique)} + + U_2d = np.zeros((ny, nx)) + V_2d = np.zeros((ny, nx)) + + for k in range(len(x)): + i = x_to_idx[x[k]] + j = y_to_idx[y[k]] + U_2d[j, i] = u[k] + V_2d[j, i] = v[k] + + return x_unique, y_unique, U_2d, V_2d + + +def plot_lic( + x_nodes: np.ndarray, + y_nodes: np.ndarray, + U_2d: np.ndarray, + V_2d: np.ndarray, + Re: float = None, + N: int = None, + output_dir: Path = None, + n_fine: int = 512, + lic_length: int = 30, +) -> Path: + """Generate LIC visualization from 2D velocity arrays. + + Parameters + ---------- + x_nodes : np.ndarray + 1D array of x node positions + y_nodes : np.ndarray + 1D array of y node positions + U_2d : np.ndarray + u velocity on 2D grid, shape (ny, nx) + V_2d : np.ndarray + v velocity on 2D grid, shape (ny, nx) + Re : float, optional + Reynolds number (for title) + N : int, optional + Grid resolution (for title) + output_dir : Path + Output directory + n_fine : int + Fine grid resolution for LIC + lic_length : int + LIC kernel length + + Returns + ------- + Path + Path to saved figure + """ + if output_dir is None: + output_dir = Path(__file__).parent.parent / "figures" / "lic" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + if N is None: + N = len(x_nodes) - 1 + + # Interpolate to fine uniform grid + print(f"Interpolating to {n_fine}x{n_fine} uniform grid...") + x_fine, y_fine, U_fine, V_fine = spectral_interpolate_to_fine_grid( + x_nodes, y_nodes, U_2d, V_2d, n_fine=n_fine + ) + + # Compute velocity magnitude + vel_mag = np.sqrt(U_fine**2 + V_fine**2) + + # Generate LIC texture + print(f"Generating LIC texture (kernel length={lic_length})...") + lic = generate_lic_fast(U_fine, V_fine, length=lic_length) + + print("Creating visualization...") + + # Normalize for coloring + vel_normalized = (vel_mag - vel_mag.min()) / (vel_mag.max() - vel_mag.min() + 1e-10) + lic_normalized = (lic - lic.min()) / (lic.max() - lic.min() + 1e-10) + + suffix = "" + if Re is not None: + suffix += f"_Re{int(Re)}" + if N is not None: + suffix += f"_N{N}" + + cmap = plt.cm.coolwarm + sm = plt.cm.ScalarMappable(cmap='coolwarm', norm=Normalize(vmin=vel_mag.min(), vmax=vel_mag.max())) + + # Plot 1: Pure LIC texture (grayscale) + fig1, ax1 = plt.subplots(figsize=(6, 6)) + ax1.imshow(lic, origin='lower', extent=[0, 1, 0, 1], cmap='gray') + ax1.set_xticks([]) + ax1.set_yticks([]) + ax1.set_aspect('equal') + plt.tight_layout() + path1 = output_dir / f"lic_grayscale{suffix}.pdf" + fig1.savefig(path1, dpi=400, bbox_inches='tight', transparent=True) + plt.close(fig1) + print(f"Saved: {path1}") + + # Plot 2: LIC colored by velocity magnitude with colorbar + fig2, ax2 = plt.subplots(figsize=(6, 6.8)) + colors = cmap(vel_normalized) + lic_factor = 0.25 + 0.75 * lic_normalized + colors[:, :, :3] *= lic_factor[:, :, np.newaxis] + ax2.imshow(colors, origin='lower', extent=[0, 1, 0, 1]) + ax2.set_xticks([]) + ax2.set_yticks([]) + ax2.set_aspect('equal') + cbar2 = plt.colorbar(sm, ax=ax2, orientation='horizontal', pad=0.05, aspect=30) + cbar2.set_label(r'$|\mathbf{u}|$') + plt.tight_layout() + path2 = output_dir / f"lic_colored{suffix}.pdf" + fig2.savefig(path2, dpi=400, bbox_inches='tight', transparent=True) + plt.close(fig2) + print(f"Saved: {path2}") + + # Plot 3: High-res LIC poster background (A2 size: 420x594mm, 300dpi = 4961x7016px) + # Using square aspect, so 16.5x16.5 inches at 300dpi = ~5000px + fig3, ax3 = plt.subplots(figsize=(16.5, 16.5)) + colors = cmap(vel_normalized) + lic_factor = 0.25 + 0.75 * lic_normalized + colors[:, :, :3] *= lic_factor[:, :, np.newaxis] + ax3.imshow(colors, origin='lower', extent=[0, 1, 0, 1]) + ax3.axis('off') + plt.tight_layout(pad=0) + path3 = output_dir / f"poster_lic{suffix}.png" + fig3.savefig(path3, dpi=300, bbox_inches='tight', pad_inches=0, transparent=True) + plt.close(fig3) + print(f"Saved: {path3}") + + return path2 + + +def main(): + import argparse + + parser = argparse.ArgumentParser( + description="Generate LIC visualization for lid-driven cavity flow" + ) + parser.add_argument( + "--npz", type=str, help="Path to solution.npz file" + ) + parser.add_argument( + "--vts", type=str, help="Path to solution.vts file" + ) + parser.add_argument( + "--run-id", type=str, help="MLflow run ID to load solution from" + ) + parser.add_argument( + "--Re", type=float, default=None, help="Reynolds number (for title)" + ) + parser.add_argument( + "--N", type=int, default=None, help="Grid resolution (for title)" + ) + parser.add_argument( + "--n-fine", type=int, default=512, help="Fine grid resolution for LIC" + ) + parser.add_argument( + "--length", type=int, default=30, help="LIC kernel length" + ) + parser.add_argument( + "--output", type=str, default=None, help="Output directory" + ) + + args = parser.parse_args() + + # Determine source and load data + Re, N = args.Re, args.N + + if args.npz: + # Load from NPZ file + npz_path = Path(args.npz) + if not npz_path.exists(): + print(f"Error: NPZ file not found: {npz_path}") + sys.exit(1) + print(f"Loading solution from: {npz_path}") + x_nodes, y_nodes, U_2d, V_2d, Re_file, N_file = load_solution_from_npz(npz_path) + Re = Re if Re is not None else Re_file + N = N if N is not None else N_file + + elif args.vts: + # Load from VTS file + vts_path = Path(args.vts) + if not vts_path.exists(): + print(f"Error: VTS file not found: {vts_path}") + sys.exit(1) + print(f"Loading solution from: {vts_path}") + fields = load_solution_from_vts(vts_path) + x_nodes, y_nodes, U_2d, V_2d = restructure_to_2d(fields) + + elif args.run_id: + # Load from MLflow + import mlflow + + print(f"Loading solution from MLflow run: {args.run_id}") + client = mlflow.tracking.MlflowClient() + artifact_dir = client.download_artifacts(args.run_id, "") + vts_path = Path(artifact_dir) / "solution.vts" + + if not vts_path.exists(): + vts_path = Path(artifact_dir) / "fields" / "solution.vts" + + if not vts_path.exists(): + print(f"Error: solution.vts not found in artifacts") + sys.exit(1) + + fields = load_solution_from_vts(vts_path) + x_nodes, y_nodes, U_2d, V_2d = restructure_to_2d(fields) + + else: + # Use bundled default solution + if not DEFAULT_SOLUTION.exists(): + print(f"Error: Bundled solution not found: {DEFAULT_SOLUTION}") + print("Please provide --npz, --vts, or --run-id") + sys.exit(1) + print(f"Using bundled solution: {DEFAULT_SOLUTION.name}") + x_nodes, y_nodes, U_2d, V_2d, Re_file, N_file = load_solution_from_npz(DEFAULT_SOLUTION) + Re = Re if Re is not None else Re_file + N = N if N is not None else N_file + + output_dir = Path(args.output) if args.output else None + + plot_lic( + x_nodes=x_nodes, + y_nodes=y_nodes, + U_2d=U_2d, + V_2d=V_2d, + Re=Re, + N=N, + output_dir=output_dir, + n_fine=args.n_fine, + lic_length=args.length, + ) + + +if __name__ == "__main__": + main() diff --git a/scripts/test_pyvista_plots.py b/scripts/test_pyvista_plots.py new file mode 100644 index 0000000..f143100 --- /dev/null +++ b/scripts/test_pyvista_plots.py @@ -0,0 +1,67 @@ +"""Test script for PyVista-based LDC plotting functions. + +Usage: + python scripts/test_pyvista_plots.py + python scripts/test_pyvista_plots.py --vts path/to/solution.vts +""" + +import sys +from pathlib import Path + +# Add src to path +sys.path.insert(0, str(Path(__file__).parent.parent / "src")) + +from shared.plotting.ldc.pyvista_fields import generate_pyvista_field_plots + +# Default VTS solution from mlruns +DEFAULT_VTS = Path(__file__).parent.parent / "mlruns" / "806262599561726093" / "719e5ef40a454b8d8a60a3f1a222c7f3" / "artifacts" / "solution.vts" + + +def main(): + import argparse + + parser = argparse.ArgumentParser(description="Test PyVista LDC plotting functions") + parser.add_argument("--vts", type=str, help="Path to solution.vts file") + parser.add_argument( + "--output", type=str, default=None, help="Output directory" + ) + args = parser.parse_args() + + # Find VTS file + if args.vts: + vts_path = Path(args.vts) + if not vts_path.exists(): + print(f"Error: VTS file not found: {vts_path}") + sys.exit(1) + print(f"Loading solution from: {vts_path}") + else: + if not DEFAULT_VTS.exists(): + print(f"Error: Default VTS not found: {DEFAULT_VTS}") + print("Please provide --vts path/to/solution.vts") + sys.exit(1) + vts_path = DEFAULT_VTS + print(f"Using default VTS: {vts_path}") + + # Output directory + if args.output: + output_dir = Path(args.output) + else: + output_dir = Path(__file__).parent.parent / "figures" / "pyvista_test" + + output_dir.mkdir(parents=True, exist_ok=True) + print(f"Output directory: {output_dir}") + + # Generate all plots + print("\nGenerating PyVista plots...") + paths = generate_pyvista_field_plots( + vts_path=vts_path, + output_dir=output_dir, + ) + + print(f"\nGenerated {len(paths)} plots:") + for name, p in paths.items(): + print(f" - {name}: {p.name}") + + +if __name__ == "__main__": + main() diff --git a/src/shared/plotting/ldc/__init__.py b/src/shared/plotting/ldc/__init__.py index 92132a5..563e6e7 100644 --- a/src/shared/plotting/ldc/__init__.py +++ b/src/shared/plotting/ldc/__init__.py @@ -6,7 +6,7 @@ from .convergence import plot_convergence from .data_loading import fields_to_dataframe, load_fields_from_zarr, restructure_fields -from .fields import plot_fields, plot_streamlines, plot_vorticity +from .fields import plot_vorticity from .mlflow_utils import ( download_mlflow_artifacts, find_sibling_runs, @@ -14,6 +14,7 @@ upload_plots_to_mlflow, ) from .orchestrator import generate_comparison_plots_for_sweep, generate_plots_for_run +from .pyvista_fields import generate_pyvista_field_plots from .validation import plot_ghia_comparison # Import style module to trigger sns.set_theme() on package import @@ -22,8 +23,7 @@ __all__ = [ "generate_plots_for_run", "generate_comparison_plots_for_sweep", - "plot_fields", - "plot_streamlines", + "generate_pyvista_field_plots", "plot_vorticity", "plot_convergence", "plot_ghia_comparison", diff --git a/src/shared/plotting/ldc/fields.py b/src/shared/plotting/ldc/fields.py index 51456fe..d7dbd5d 100644 --- a/src/shared/plotting/ldc/fields.py +++ b/src/shared/plotting/ldc/fields.py @@ -11,7 +11,6 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd -import pyvista as pv from scipy.interpolate import RectBivariateSpline log = logging.getLogger(__name__) @@ -202,171 +201,3 @@ def plot_vorticity( plt.close(fig) return output_path - - -def plot_streamlines_pyvista( - fields_df: pd.DataFrame, Re: float, solver: str, N: int, output_dir: Path -) -> Path: - """Generate beautiful streamline plot using PyVista with ParaView theme. - - Creates a visually striking visualization with velocity magnitude field - and streamlines, with transparent background for easy compositing. - """ - x_unique = np.sort(fields_df["x"].unique()) - y_unique = np.sort(fields_df["y"].unique()) - nx, ny = len(x_unique), len(y_unique) - - sorted_df = fields_df.sort_values(["y", "x"]) - U = sorted_df["u"].values.reshape(ny, nx) - V = sorted_df["v"].values.reshape(ny, nx) - - # Interpolate to finer grid for smoother visualization - n_fine = 200 - x_fine = np.linspace(x_unique[0], x_unique[-1], n_fine) - y_fine = np.linspace(y_unique[0], y_unique[-1], n_fine) - - U_interp = RectBivariateSpline(y_unique, x_unique, U)(y_fine, x_fine) - V_interp = RectBivariateSpline(y_unique, x_unique, V)(y_fine, x_fine) - - # Create 3D grid (z=0 plane) - X, Y = np.meshgrid(x_fine, y_fine) - Z = np.zeros_like(X) - - # Create structured grid - grid = pv.StructuredGrid(X, Y, Z) - - # Add velocity as vector field (3D with w=0) - vel_mag = np.sqrt(U_interp**2 + V_interp**2) - vectors = np.zeros((n_fine * n_fine, 3)) - vectors[:, 0] = U_interp.ravel() - vectors[:, 1] = V_interp.ravel() - vectors[:, 2] = 0.0 - - grid["velocity"] = vectors - grid["velocity_magnitude"] = vel_mag.ravel() - grid.set_active_vectors("velocity") - - # Create seed points for streamlines - n_seeds = 15 - seed_points_left = np.column_stack([ - np.full(n_seeds, x_fine[3]), - np.linspace(y_fine[3], y_fine[-4], n_seeds), - np.zeros(n_seeds) - ]) - seed_points_bottom = np.column_stack([ - np.linspace(x_fine[3], x_fine[-4], n_seeds), - np.full(n_seeds, y_fine[3]), - np.zeros(n_seeds) - ]) - seed_points = np.vstack([seed_points_left, seed_points_bottom]) - seeds = pv.PolyData(seed_points) - - # Generate streamlines - streamlines = grid.streamlines_from_source( - seeds, - vectors="velocity", - max_steps=2000, - integration_direction="both", - ) - - # Set up PyVista plotter with ParaView theme - pv.set_plot_theme("paraview") - plotter = pv.Plotter(off_screen=True, window_size=[1400, 1200]) - - # Add velocity magnitude field as background surface - surface = grid.extract_surface() - plotter.add_mesh( - surface, - scalars="velocity_magnitude", - cmap="turbo", - show_edges=False, - lighting=False, - scalar_bar_args={ - "title": "Velocity Magnitude |u|", - "vertical": True, - "position_x": 0.85, - "position_y": 0.2, - "width": 0.08, - "height": 0.6, - "title_font_size": 16, - "label_font_size": 14, - "fmt": "%.3f", - "n_labels": 5, - }, - ) - - # Add streamlines as white tubes on top - if streamlines.n_points > 0: - tubes = streamlines.tube(radius=0.004, n_sides=8) - plotter.add_mesh( - tubes, - color="white", - opacity=0.85, - smooth_shading=True, - specular=0.3, - ) - - # Set up camera for clean 2D view - plotter.camera_position = "xy" - plotter.camera.zoom(1.15) - - # Enable anti-aliasing for smooth edges - plotter.enable_anti_aliasing("ssaa") - - # Save with transparent background - output_path = output_dir / "streamlines_3d.png" - plotter.screenshot(output_path, transparent_background=True, scale=2) - plotter.close() - - return output_path - - -def export_fields_to_vtk( - fields_df: pd.DataFrame, Re: float, solver: str, N: int, output_dir: Path -) -> Path: - """Export solution fields to VTK format for ParaView visualization. - - Creates a structured grid VTK file with pressure, velocity components, - velocity magnitude, and vorticity fields. - """ - x_unique = np.sort(fields_df["x"].unique()) - y_unique = np.sort(fields_df["y"].unique()) - nx, ny = len(x_unique), len(y_unique) - - sorted_df = fields_df.sort_values(["y", "x"]) - P = sorted_df["p"].values.reshape(ny, nx) - U = sorted_df["u"].values.reshape(ny, nx) - V = sorted_df["v"].values.reshape(ny, nx) - - # Create 3D grid (z=0 plane) - X, Y = np.meshgrid(x_unique, y_unique) - Z = np.zeros_like(X) - - # Create structured grid - grid = pv.StructuredGrid(X, Y, Z) - - # Add scalar fields - grid["pressure"] = P.ravel() - grid["u"] = U.ravel() - grid["v"] = V.ravel() - grid["velocity_magnitude"] = np.sqrt(U**2 + V**2).ravel() - - # Compute and add vorticity - U_spline = RectBivariateSpline(y_unique, x_unique, U) - V_spline = RectBivariateSpline(y_unique, x_unique, V) - dvdx = V_spline(y_unique, x_unique, dx=1) - dudy = U_spline(y_unique, x_unique, dy=1) - vorticity = dvdx - dudy - grid["vorticity"] = vorticity.ravel() - - # Add velocity vector field - vectors = np.zeros((nx * ny, 3)) - vectors[:, 0] = U.ravel() - vectors[:, 1] = V.ravel() - grid["velocity"] = vectors - - # Save as VTK - output_path = output_dir / "solution.vtk" - grid.save(output_path) - - return output_path diff --git a/src/shared/plotting/ldc/orchestrator.py b/src/shared/plotting/ldc/orchestrator.py index 8bc6b04..e498bfb 100644 --- a/src/shared/plotting/ldc/orchestrator.py +++ b/src/shared/plotting/ldc/orchestrator.py @@ -12,18 +12,14 @@ from .convergence import plot_convergence from .data_loading import fields_to_dataframe, load_fields_from_zarr -from .fields import ( - plot_fields, - plot_streamlines, - plot_streamlines_pyvista, - plot_vorticity, -) +from .fields import plot_vorticity from .mlflow_utils import ( download_mlflow_artifacts, find_sibling_runs, load_timeseries_from_mlflow, upload_plots_to_mlflow, ) +from .pyvista_fields import generate_pyvista_field_plots from .validation import plot_ghia_comparison log = logging.getLogger(__name__) @@ -39,7 +35,14 @@ def generate_plots_for_run( parent_run_id: Optional[str] = None, upload_to_mlflow: bool = True, ) -> list[Path]: - """Generate all plots for a completed run.""" + """Generate all plots for a completed run. + + Artifacts generated: + - convergence.pdf (matplotlib) + - ghia_comparison.pdf (matplotlib) + - vorticity.pdf (matplotlib) + - u.png, v.png, pressure.png, vel-mag.png, streamlines.png (PyVista) + """ output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) @@ -51,13 +54,12 @@ def generate_plots_for_run( log.info(f"Generating plots for {solver_name} N={N} Re={Re}") - # Generate individual plots plot_paths = [] - plot_paths.append(plot_fields(fields_df, Re, solver_name, N, output_dir)) - plot_paths.append(plot_streamlines(fields_df, Re, solver_name, N, output_dir)) - plot_paths.append(plot_streamlines_pyvista(fields_df, Re, solver_name, N, output_dir)) - plot_paths.append(plot_vorticity(fields_df, Re, solver_name, N, output_dir)) + + # Matplotlib plots: convergence, ghia comparison, vorticity plot_paths.append(plot_convergence(timeseries_df, Re, solver_name, N, output_dir)) + plot_paths.append(plot_vorticity(fields_df, Re, solver_name, N, output_dir)) + ghia_path = plot_ghia_comparison( [{"run_id": run_id, "N": N, "Re": Re, "solver": solver_name, "status": "FINISHED"}], tracking_uri, @@ -66,6 +68,14 @@ def generate_plots_for_run( if ghia_path: plot_paths.append(ghia_path) + # PyVista plots: u, v, pressure, vel-mag, streamlines + vts_path = artifact_dir / "solution.vts" + if vts_path.exists(): + pyvista_paths = generate_pyvista_field_plots(vts_path, output_dir) + plot_paths.extend(pyvista_paths.values()) + else: + log.warning(f"VTS file not found at {vts_path}, skipping PyVista plots") + plot_paths = [p for p in plot_paths if p is not None] log.info(f"Generated {len(plot_paths)} plots for run") diff --git a/src/shared/plotting/ldc/pyvista_fields.py b/src/shared/plotting/ldc/pyvista_fields.py new file mode 100644 index 0000000..a5f31f0 --- /dev/null +++ b/src/shared/plotting/ldc/pyvista_fields.py @@ -0,0 +1,589 @@ +""" +PyVista-based Field Visualization for LDC. + +Generates high-quality 2D field plots using PyVista with the ParaView theme. +Loads solution directly from VTS files. +""" + +import logging +import os +import subprocess +import sys +import tempfile +from pathlib import Path + +import numpy as np +import pyvista as pv + +log = logging.getLogger(__name__) + +# Enable off-screen rendering +os.environ["PYVISTA_OFF_SCREEN"] = "true" + +# Use ParaView theme +pv.set_plot_theme("paraview") + + +# Window size for high DPI output (2400x2400 = 2x the previous 1200x1200) +WINDOW_SIZE = [2400, 2400] + +# Common scalar bar arguments with black text and larger font +SCALAR_BAR_ARGS = { + "vertical": False, + "position_x": 0.25, + "position_y": 0.02, # Lower position for more padding + "width": 0.5, + "height": 0.04, + "title_font_size": 44, # Scaled up for higher resolution + "label_font_size": 32, # Scaled up for higher resolution + "color": "black", + "fmt": "%.2f", + "n_labels": 5, + "italic": True, # Italic for more LaTeX-like appearance + "font_family": "times", # Times font for LaTeX-like appearance +} + + +def _setup_camera(plotter: pv.Plotter) -> None: + """Set up proper 2D orthographic view after mesh is added.""" + plotter.enable_parallel_projection() + plotter.view_xy() + plotter.reset_camera_clipping_range() + plotter.reset_camera(bounds=plotter.bounds) + + +_STREAMLINES_SCRIPT = ''' +import os +import sys +os.environ["PYVISTA_OFF_SCREEN"] = "true" + +import pyvista as pv +pv.set_plot_theme("paraview") + +vts_path = sys.argv[1] +output_path = sys.argv[2] +separating_distance = float(sys.argv[3]) + +WINDOW_SIZE = [2400, 2400] +SCALAR_BAR_ARGS = { + "vertical": False, + "position_x": 0.25, + "position_y": 0.02, + "width": 0.5, + "height": 0.04, + "title_font_size": 44, + "label_font_size": 32, + "color": "black", + "fmt": "%.2f", + "n_labels": 5, + "italic": True, + "font_family": "times", +} + +grid = pv.read(vts_path) +plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + +plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\\n"}, +) + +bounds = grid.bounds +resolution = 256 +image_data = pv.ImageData( + dimensions=(resolution, resolution, 1), + spacing=((bounds[1] - bounds[0]) / (resolution - 1), + (bounds[3] - bounds[2]) / (resolution - 1), + 1.0), + origin=(bounds[0], bounds[2], 0.0), +) + +sampled = image_data.sample(grid) + +streamlines = sampled.streamlines_evenly_spaced_2D( + vectors="velocity", + start_position=(0.3, 0.7, 0.0), + separating_distance=separating_distance, + separating_distance_ratio=0.5, + step_length=0.5, + max_steps=800, +) + +if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=2.0, + opacity=0.2, + ) + +plotter.enable_parallel_projection() +plotter.view_xy() +plotter.reset_camera_clipping_range() +plotter.reset_camera(bounds=plotter.bounds) + +plotter.screenshot(output_path, transparent_background=True) +plotter.close() +''' + + +def _generate_streamlines_safe(vts_path: Path, output_path: Path, timeout: int = 60) -> bool: + """Generate streamlines plot in a subprocess to isolate VTK crashes. + + Tries dense streamlines first, falls back to sparser if it crashes. + """ + vts_str = str(vts_path) + out_str = str(output_path) + + # Try dense streamlines first (separating_distance=2.5) + for sep_dist in [2.5, 4.0, 8.0]: + try: + result = subprocess.run( + [sys.executable, "-c", _STREAMLINES_SCRIPT, vts_str, out_str, str(sep_dist)], + timeout=timeout, + capture_output=True, + text=True, + ) + if result.returncode == 0: + log.info(f"Streamlines generated with separating_distance={sep_dist}") + return True + else: + log.warning(f"Streamlines failed with separating_distance={sep_dist}: {result.stderr[:200]}") + except subprocess.TimeoutExpired: + log.warning(f"Streamlines timed out with separating_distance={sep_dist}") + except Exception as e: + log.warning(f"Streamlines error with separating_distance={sep_dist}: {e}") + + log.error("All streamline attempts failed") + return False + + +def plot_u_velocity( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate u-velocity contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="u", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Horizontal velocity, u\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"u{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_v_velocity( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate v-velocity contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="v", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Vertical velocity, v\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"v{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_velocity_magnitude( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate velocity magnitude contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"vel-mag{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_pressure( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate pressure contour plot.""" + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + plotter.add_mesh( + grid, + scalars="pressure", + cmap="inferno", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Pressure, p\n"}, + ) + _setup_camera(plotter) + + output_path = output_dir / f"pressure{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_streamlines( + grid: pv.StructuredGrid, + output_dir: Path, + n_seeds: int = 8, + suffix: str = "", +) -> Path: + """Generate streamlines plot with velocity magnitude background. + + Parameters + ---------- + grid : pv.StructuredGrid + Solution grid with velocity vectors + output_dir : Path + Output directory + n_seeds : int + Number of seed points per dimension (reduced for less density) + suffix : str + Filename suffix + + Returns + ------- + Path + Path to saved figure + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Create seed points on a grid (avoiding exact boundaries) + seed_x = np.linspace(0.1, 0.9, n_seeds) + seed_y = np.linspace(0.1, 0.9, n_seeds) + seed_X, seed_Y = np.meshgrid(seed_x, seed_y) + seed_points = np.column_stack([ + seed_X.ravel(), + seed_Y.ravel(), + np.zeros(n_seeds * n_seeds) + ]) + seeds = pv.PolyData(seed_points) + + # Compute streamlines + streamlines = grid.streamlines_from_source( + seeds, + vectors="velocity", + integration_direction="both", + initial_step_length=0.01, + max_step_length=0.1, + ) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + + # Add velocity magnitude as background + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + + # Add streamlines + if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=3.0, # Scaled for higher resolution + opacity=0.8, + ) + + _setup_camera(plotter) + + output_path = output_dir / f"streamlines{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_streamlines_evenly_spaced( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate evenly-spaced streamlines using line sources. + + Uses line sources at multiple positions to create evenly distributed + streamlines across the domain, similar to PyVista's 2D streamlines example. + + Parameters + ---------- + grid : pv.StructuredGrid + Solution grid with velocity vectors + output_dir : Path + Output directory + suffix : str + Filename suffix + + Returns + ------- + Path + Path to saved figure + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Use line sources for evenly spaced seeds (like PyVista example) + # Vertical line at left side of domain + streamlines_left = grid.streamlines( + pointa=(0.05, 0.05, 0), + pointb=(0.05, 0.95, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + # Horizontal line at top of domain + streamlines_top = grid.streamlines( + pointa=(0.05, 0.95, 0), + pointb=(0.95, 0.95, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + # Vertical line at right side + streamlines_right = grid.streamlines( + pointa=(0.95, 0.95, 0), + pointb=(0.95, 0.05, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + # Horizontal line at bottom + streamlines_bottom = grid.streamlines( + pointa=(0.95, 0.05, 0), + pointb=(0.05, 0.05, 0), + n_points=12, + vectors="velocity", + integration_direction="forward", + initial_step_length=0.01, + max_step_length=0.1, + max_steps=500, + ) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + + # Add velocity magnitude as background + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + + # Add all streamlines + for streamlines in [streamlines_left, streamlines_top, streamlines_right, streamlines_bottom]: + if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=3.0, # Scaled for higher resolution + opacity=0.9, + ) + + _setup_camera(plotter) + + output_path = output_dir / f"streamlines_evenly_spaced{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def plot_streamlines_evenly_spaced_2D( + grid: pv.StructuredGrid, + output_dir: Path, + suffix: str = "", +) -> Path: + """Generate evenly-spaced streamlines using VTK's 2D algorithm. + + Uses streamlines_evenly_spaced_2D for publication-quality streamlines. + This algorithm produces streamlines that are approximately evenly spaced. + + Parameters + ---------- + grid : pv.StructuredGrid + Solution grid with velocity vectors + output_dir : Path + Output directory + suffix : str + Filename suffix + + Returns + ------- + Path + Path to saved figure + """ + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + plotter = pv.Plotter(off_screen=True, window_size=WINDOW_SIZE) + + # Add velocity magnitude as background + plotter.add_mesh( + grid, + scalars="velocity_magnitude", + cmap="coolwarm", + show_edges=False, + scalar_bar_args={**SCALAR_BAR_ARGS, "title": "Velocity magnitude, |u|\n"}, + ) + + # Use evenly spaced streamlines 2D algorithm + # First resample to ImageData for better algorithm performance + try: + # Create uniform ImageData grid with higher resolution for smoother streamlines + bounds = grid.bounds + resolution = 256 # Balance between smoothness and stability + image_data = pv.ImageData( + dimensions=(resolution, resolution, 1), + spacing=((bounds[1] - bounds[0]) / (resolution - 1), + (bounds[3] - bounds[2]) / (resolution - 1), + 1.0), + origin=(bounds[0], bounds[2], 0.0), + ) + + # Sample the velocity field onto the uniform grid + sampled = image_data.sample(grid) + + streamlines = sampled.streamlines_evenly_spaced_2D( + vectors="velocity", + start_position=(0.3, 0.7, 0.0), # Near primary vortex + separating_distance=8.0, # Conservative for stability with non-uniform grids + separating_distance_ratio=0.5, + step_length=0.5, # Step size in cell units + max_steps=800, + ) + + if streamlines.n_points > 0: + plotter.add_mesh( + streamlines, + color="white", + line_width=2.0, # Scaled for higher resolution + opacity=0.2, # Subtle but visible + ) + log.info(f"Generated {streamlines.n_lines} evenly-spaced 2D streamlines") + except Exception as e: + log.warning(f"streamlines_evenly_spaced_2D failed: {e}") + + _setup_camera(plotter) + + output_path = output_dir / f"streamlines{suffix}.png" + plotter.screenshot(output_path, transparent_background=True) + plotter.close() + + log.info(f"Saved: {output_path}") + return output_path + + +def generate_pyvista_field_plots( + vts_path: Path, + output_dir: Path, +) -> dict[str, Path]: + """Generate PyVista-based field plots from a VTS solution file. + + Generates clean artifact names: u.png, v.png, pressure.png, vel-mag.png, streamlines.png + + Parameters + ---------- + vts_path : Path + Path to solution.vts file + output_dir : Path + Output directory + + Returns + ------- + dict[str, Path] + Dictionary mapping plot names to their paths + """ + vts_path = Path(vts_path) + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + # Load solution + grid = pv.read(str(vts_path)) + + paths = {} + + log.info("Generating PyVista field plots...") + + log.info(" u velocity...") + paths["u"] = plot_u_velocity(grid, output_dir, suffix="") + + log.info(" v velocity...") + paths["v"] = plot_v_velocity(grid, output_dir, suffix="") + + log.info(" velocity magnitude...") + paths["vel-mag"] = plot_velocity_magnitude(grid, output_dir, suffix="") + + log.info(" pressure...") + paths["pressure"] = plot_pressure(grid, output_dir, suffix="") + + log.info(" streamlines (2D evenly spaced)...") + streamlines_path = output_dir / "streamlines.png" + if _generate_streamlines_safe(vts_path, streamlines_path): + paths["streamlines"] = streamlines_path + + return paths diff --git a/src/shared/plotting/ldc/validation.py b/src/shared/plotting/ldc/validation.py index 3621be4..5b35c95 100644 --- a/src/shared/plotting/ldc/validation.py +++ b/src/shared/plotting/ldc/validation.py @@ -28,7 +28,7 @@ def _build_method_label(sibling: dict) -> str: - 'fv' -> 'FV' - 'spectral' -> 'Spectral' - 'spectral_fsg' -> 'Spectral-FSG' - - 'spectral_vmg' -> 'Spectral-VMG' + - 'spectral_fmg' -> 'Spectral-FMG' """ solver = sibling.get("solver", "unknown") @@ -37,7 +37,6 @@ def _build_method_label(sibling: dict) -> str: "fv": "FV", "spectral": "Spectral", "spectral_fsg": "Spectral-FSG", - "spectral_vmg": "Spectral-VMG", "spectral_fmg": "Spectral-FMG", } diff --git a/src/solvers/__init__.py b/src/solvers/__init__.py index 8646e9a..284b645 100644 --- a/src/solvers/__init__.py +++ b/src/solvers/__init__.py @@ -7,9 +7,7 @@ LidDrivenCavitySolver (abstract base - defines problem) ├── FVSolver (finite volume with SIMPLE algorithm) └── SGSolver (spectral single grid base) - ├── FSGSolver (Full Single Grid multigrid) - ├── VMGSolver (V-cycle MultiGrid) - └── FMGSolver (Full MultiGrid) + └── FSGSolver (Full Single Grid multigrid) """ from .base import LidDrivenCavitySolver @@ -29,8 +27,6 @@ from solvers.fv.solver import FVSolver from solvers.spectral.sg import SGSolver from solvers.spectral.fsg import FSGSolver -from solvers.spectral.vmg import VMGSolver -from solvers.spectral.fmg import FMGSolver __all__ = [ @@ -48,8 +44,6 @@ # Spectral solvers "SGSolver", "FSGSolver", - "VMGSolver", - "FMGSolver", "SpectralParameters", "SpectralSolverFields", ] diff --git a/src/solvers/datastructures.py b/src/solvers/datastructures.py index a47f7dd..e0615fd 100644 --- a/src/solvers/datastructures.py +++ b/src/solvers/datastructures.py @@ -265,7 +265,7 @@ class SpectralParameters(Parameters): corner_smoothing: float = 0.15 # smoothing width for smoothing method # Multigrid settings - multigrid: str = "none" # "none", "fsg", "vmg", "fmg" + multigrid: str = "none" # "none", "fsg" n_levels: int = 3 coarse_tolerance_factor: float = 10.0 diff --git a/src/solvers/spectral/__init__.py b/src/solvers/spectral/__init__.py index f6b5e73..00b2c6f 100644 --- a/src/solvers/spectral/__init__.py +++ b/src/solvers/spectral/__init__.py @@ -2,7 +2,5 @@ from solvers.spectral.sg import SGSolver from solvers.spectral.fsg import FSGSolver -from solvers.spectral.vmg import VMGSolver -from solvers.spectral.fmg import FMGSolver -__all__ = ["SGSolver", "FSGSolver", "VMGSolver", "FMGSolver"] +__all__ = ["SGSolver", "FSGSolver"] diff --git a/src/solvers/spectral/fmg.py b/src/solvers/spectral/fmg.py deleted file mode 100644 index b51a290..0000000 --- a/src/solvers/spectral/fmg.py +++ /dev/null @@ -1,45 +0,0 @@ -"""Full MultiGrid (FMG) spectral solver for lid-driven cavity. - -FMG will extend the base SG solver with Full MultiGrid acceleration. -Currently not implemented. -""" - -import logging - -from .sg import SGSolver - -log = logging.getLogger(__name__) - - -class FMGSolver(SGSolver): - """Full MultiGrid (FMG) spectral solver. - - Extends the base Single Grid solver with Full MultiGrid acceleration. - - Parameters - ---------- - All parameters inherited from SGSolver, plus: - n_levels : int - Number of multigrid levels - coarse_tolerance_factor : float - Tolerance multiplier for coarse grids - prolongation_method : str - Transfer operator for coarse-to-fine ('fft' or 'polynomial') - restriction_method : str - Transfer operator for fine-to-coarse ('fft' or 'polynomial') - """ - - def solve(self, tolerance: float = None, max_iter: int = None): - """Solve using Full MultiGrid (FMG). - - Parameters - ---------- - tolerance : float, optional - Convergence tolerance. If None, uses params.tolerance. - max_iter : int, optional - Maximum iterations. If None, uses params.max_iterations. - """ - raise NotImplementedError( - "Full MultiGrid (FMG) solver not yet implemented. " - "Use FSG (Full Single Grid) instead." - ) diff --git a/src/solvers/spectral/operators/transfer_operators.py b/src/solvers/spectral/operators/transfer_operators.py index 838c19a..6d85e4a 100644 --- a/src/solvers/spectral/operators/transfer_operators.py +++ b/src/solvers/spectral/operators/transfer_operators.py @@ -11,7 +11,6 @@ - Polynomial-based: Uses Chebyshev polynomial fitting/evaluation For FSG (Full Single Grid), only prolongation is needed. -For VMG/FMG, both prolongation and restriction are required. """ from abc import ABC, abstractmethod diff --git a/src/solvers/spectral/vmg.py b/src/solvers/spectral/vmg.py deleted file mode 100644 index 5122f0c..0000000 --- a/src/solvers/spectral/vmg.py +++ /dev/null @@ -1,459 +0,0 @@ -"""V-cycle Multigrid (VMG) spectral solver. - -Implements proper V-cycle multigrid using Full Approximation Storage (FAS) scheme. -Uses recursive coarse grid correction with damped prolongation. - -Achieves ~2.8x speedup vs single grid with proper parameter tuning. -""" - -import logging -import time -from typing import List, Tuple, Optional - -import numpy as np - -from .sg import SGSolver -from solvers.spectral.multigrid.fsg import ( - SpectralLevel, - MultigridSmoother, - build_hierarchy, - prolongate_solution, - restrict_solution, - restrict_residual, -) -from solvers.spectral.operators.transfer_operators import ( - TransferOperators, - create_transfer_operators, -) -from solvers.spectral.operators.corner import ( - CornerTreatment, - create_corner_treatment, -) - -log = logging.getLogger(__name__) - - -def _vcycle_fas( - levels: List[SpectralLevel], - smoothers: List[MultigridSmoother], - transfer_ops: TransferOperators, - level_idx: int, - pre_smooth: int, - post_smooth: int, - damping: float, -) -> None: - """Perform one V-cycle using Full Approximation Storage (FAS). - - FAS Algorithm: - 1. Pre-smooth on current grid - 2. Compute fine residual: r_h = RHS - N(u_h) - 3. Restrict solution: u_H = I(u_h) - 4. Restrict residual: r_H = I(r_h) - 5. Compute coarse residual: r_H' = RHS_H - N_H(u_H) - 6. tau = r_H - r_H' (FAS correction term) - 7. Solve on coarse: N_H(v_H) = RHS_H + tau, starting from u_H - 8. Correction: e_H = v_H - u_H - 9. Prolongate: e_h = P(e_H) - 10. Update: u_h = u_h + damping * e_h - 11. Post-smooth - - Parameters - ---------- - levels : List[SpectralLevel] - Grid hierarchy (index 0 = coarsest) - smoothers : List[MultigridSmoother] - Smoother for each level - transfer_ops : TransferOperators - Prolongation/restriction operators - level_idx : int - Current level index - pre_smooth : int - Number of pre-smoothing steps - post_smooth : int - Number of post-smoothing steps - damping : float - Damping factor for coarse grid correction (0 < damping <= 1) - """ - level = levels[level_idx] - smoother = smoothers[level_idx] - - # Pre-smoothing - for _ in range(pre_smooth): - smoother.step() - - if level_idx > 0: # Not coarsest level - coarse_level = levels[level_idx - 1] - coarse_smoother = smoothers[level_idx - 1] - - # Store current fine solution - u_fine_old = level.u.copy() - v_fine_old = level.v.copy() - p_fine_old = level.p.copy() - - # Compute fine grid residual - smoother._compute_residuals(level.u, level.v, level.p) - - # Restrict residual to coarse (I(r_h)) - restrict_residual(level, coarse_level, transfer_ops) - I_r_u = coarse_level.R_u.copy() - I_r_v = coarse_level.R_v.copy() - I_r_p = coarse_level.R_p.copy() - - # Restrict solution to coarse using injection (u_H = I(u_h)) - restrict_solution(level, coarse_level, transfer_ops) - u_H = coarse_level.u.copy() - v_H = coarse_level.v.copy() - p_H = coarse_level.p.copy() - - # Compute coarse grid residual of restricted solution (r_H') - coarse_smoother._compute_residuals( - coarse_level.u, coarse_level.v, coarse_level.p - ) - - # tau = I(r_h) - r_H' (FAS correction term) - tau_u = I_r_u - coarse_level.R_u - tau_v = I_r_v - coarse_level.R_v - tau_p = I_r_p - coarse_level.R_p - - # Zero tau at boundaries (BCs are enforced separately) - tau_u_2d = tau_u.reshape(coarse_level.shape_full) - tau_v_2d = tau_v.reshape(coarse_level.shape_full) - tau_u_2d[0, :] = 0.0 - tau_u_2d[-1, :] = 0.0 - tau_u_2d[:, 0] = 0.0 - tau_u_2d[:, -1] = 0.0 - tau_v_2d[0, :] = 0.0 - tau_v_2d[-1, :] = 0.0 - tau_v_2d[:, 0] = 0.0 - tau_v_2d[:, -1] = 0.0 - - # Set tau correction for coarse solve - coarse_smoother.set_tau_correction(tau_u, tau_v, tau_p) - - # Recurse to coarser level - _vcycle_fas( - levels, smoothers, transfer_ops, - level_idx - 1, pre_smooth, post_smooth, damping - ) - - # Clear tau correction - coarse_smoother.clear_tau_correction() - - # Compute correction: e_H = v_H - u_H - e_u = coarse_level.u - u_H - e_v = coarse_level.v - v_H - e_p = coarse_level.p - p_H - - # Zero boundary corrections before prolongation - e_u_2d = e_u.reshape(coarse_level.shape_full) - e_v_2d = e_v.reshape(coarse_level.shape_full) - e_u_2d[0, :] = 0.0 - e_u_2d[-1, :] = 0.0 - e_u_2d[:, 0] = 0.0 - e_u_2d[:, -1] = 0.0 - e_v_2d[0, :] = 0.0 - e_v_2d[-1, :] = 0.0 - e_v_2d[:, 0] = 0.0 - e_v_2d[:, -1] = 0.0 - - # Prolongate correction to fine grid - e_u_fine = transfer_ops.prolongation.prolongate_2d( - e_u_2d, level.shape_full - ) - e_v_fine = transfer_ops.prolongation.prolongate_2d( - e_v_2d, level.shape_full - ) - e_p_fine = transfer_ops.prolongation.prolongate_2d( - e_p.reshape(coarse_level.shape_inner), level.shape_inner - ) - - # Zero boundary corrections on fine grid - e_u_fine[0, :] = 0.0 - e_u_fine[-1, :] = 0.0 - e_u_fine[:, 0] = 0.0 - e_u_fine[:, -1] = 0.0 - e_v_fine[0, :] = 0.0 - e_v_fine[-1, :] = 0.0 - e_v_fine[:, 0] = 0.0 - e_v_fine[:, -1] = 0.0 - - # Apply damped correction - level.u[:] = u_fine_old + damping * e_u_fine.ravel() - level.v[:] = v_fine_old + damping * e_v_fine.ravel() - level.p[:] = p_fine_old + damping * e_p_fine.ravel() - - # Re-apply boundary conditions - smoother._enforce_boundary_conditions(level.u, level.v) - - # Post-smoothing - for _ in range(post_smooth): - smoother.step() - - -def solve_vmg( - levels: List[SpectralLevel], - Re: float, - beta_squared: float, - lid_velocity: float, - CFL: float, - tolerance: float, - max_cycles: int, - transfer_ops: Optional[TransferOperators] = None, - corner_treatment: Optional[CornerTreatment] = None, - Lx: float = 1.0, - Ly: float = 1.0, - pre_smooth: int = 2, - post_smooth: int = 2, - damping: float = 0.5, -) -> Tuple[SpectralLevel, int, bool]: - """Solve using V-cycle Multigrid (VMG) with FAS. - - Parameters - ---------- - levels : List[SpectralLevel] - Grid hierarchy (index 0 = coarsest) - Re, beta_squared, lid_velocity, CFL : float - Solver parameters - tolerance : float - Convergence tolerance - max_cycles : int - Maximum V-cycles - transfer_ops : TransferOperators, optional - Transfer operators - corner_treatment : CornerTreatment, optional - Corner singularity treatment - Lx, Ly : float - Domain dimensions - pre_smooth : int - Pre-smoothing steps per level (default 2) - post_smooth : int - Post-smoothing steps per level (default 2) - damping : float - Damping factor for coarse correction (default 0.5) - - Returns - ------- - tuple - (finest_level, total_work_units, converged) - """ - if transfer_ops is None: - transfer_ops = create_transfer_operators("fft", "fft") - - if corner_treatment is None: - corner_treatment = create_corner_treatment(method="smoothing") - - # Handle subtraction method on coarse levels - uses_subtraction = corner_treatment.uses_modified_convection() - if uses_subtraction: - smoothing_treatment = create_corner_treatment(method="smoothing") - min_n_for_subtraction = 8 - - n_levels = len(levels) - - # Create smoothers for all levels - smoothers = [] - for level_idx, level in enumerate(levels): - if uses_subtraction and level.n < min_n_for_subtraction: - level_corner = smoothing_treatment - else: - level_corner = corner_treatment - - sm = MultigridSmoother( - level=level, - Re=Re, - beta_squared=beta_squared, - lid_velocity=lid_velocity, - CFL=CFL, - corner_treatment=level_corner, - Lx=Lx, - Ly=Ly, - ) - sm.initialize_lid() - smoothers.append(sm) - - # Initialize finest level - finest = levels[-1] - finest.u[:] = 0.0 - finest.v[:] = 0.0 - finest.p[:] = 0.0 - smoothers[-1].initialize_lid() - - # Track convergence via solution change - u_prev = finest.u.copy() - v_prev = finest.v.copy() - - # Work units = smoothing steps per cycle - work_per_cycle = (pre_smooth + post_smooth) * n_levels - total_work = 0 - - for cycle in range(max_cycles): - # Perform one V-cycle starting from finest level - _vcycle_fas( - levels, smoothers, transfer_ops, - n_levels - 1, pre_smooth, post_smooth, damping - ) - total_work += work_per_cycle - - # Check convergence via relative solution change - u_change = np.linalg.norm(finest.u - u_prev) / (np.linalg.norm(u_prev) + 1e-12) - v_change = np.linalg.norm(finest.v - v_prev) / (np.linalg.norm(v_prev) + 1e-12) - u_prev[:] = finest.u - v_prev[:] = finest.v - - max_change = max(u_change, v_change) - - if max_change < tolerance: - log.info( - f"VMG converged in {cycle + 1} cycles ({total_work} work units), " - f"residual={max_change:.2e}" - ) - return finest, total_work, True - - if (cycle + 1) % 50 == 0: - log.debug(f"VMG cycle {cycle + 1}: residual = {max_change:.2e}") - - log.warning(f"VMG did not converge after {max_cycles} cycles") - return finest, total_work, False - - -class VMGSolver(SGSolver): - """V-cycle Multigrid (VMG) spectral solver. - - Implements proper V-cycle multigrid with Full Approximation Storage (FAS). - Uses recursive coarse grid correction with damped prolongation for - accelerated convergence. - - Achieves ~2.8x speedup vs single grid with default parameters: - - pre_smooth=2, post_smooth=2, damping=0.5 - - Parameters - ---------- - All parameters inherited from SGSolver, plus: - n_levels : int - Maximum number of multigrid levels (default 10, limited by coarsest_n) - coarsest_n : int - Minimum N for coarsest grid (default 12). Coarse grids must resolve physics. - pre_smooth : int - Pre-smoothing steps per level (default 2) - post_smooth : int - Post-smoothing steps per level (default 2) - damping : float - Damping factor for coarse grid correction (default 0.5) - prolongation_method : str - Transfer operator for coarse-to-fine ('fft') - restriction_method : str - Transfer operator for fine-to-coarse ('fft') - """ - - def __init__(self, **kwargs): - """Initialize VMG solver with multigrid hierarchy.""" - super().__init__(**kwargs) - - # Create transfer operators from config - self._transfer_ops = create_transfer_operators( - prolongation_method=self.params.prolongation_method, - restriction_method=self.params.restriction_method, - ) - - # Get coarsest_n from params (default 12) - coarsest_n = getattr(self.params, 'coarsest_n', 12) - - # Build grid hierarchy dynamically from nx down to coarsest_n - # n_levels=100 is just a large upper bound; actual levels determined by coarsest_n - self._levels = build_hierarchy( - n_fine=self.params.nx, - n_levels=100, - coarsest_n=coarsest_n, - basis_x=self.basis_x, - basis_y=self.basis_y, - Lx=self.params.Lx, - Ly=self.params.Ly, - ) - - # V-cycle parameters - self._pre_smooth = getattr(self.params, 'pre_smooth', 2) - self._post_smooth = getattr(self.params, 'post_smooth', 2) - self._damping = getattr(self.params, 'damping', 0.5) - - log.info( - f"VMG initialized: {len(self._levels)} levels {[l.n for l in self._levels]}, " - f"pre={self._pre_smooth}, post={self._post_smooth}, " - f"damping={self._damping}" - ) - - def solve(self, tolerance: float = None, max_iter: int = None): - """Solve using V-cycle Multigrid (VMG). - - Parameters - ---------- - tolerance : float, optional - Convergence tolerance. If None, uses params.tolerance. - max_iter : int, optional - Maximum V-cycles. If None, uses params.max_iterations. - """ - if tolerance is None: - tolerance = self.params.tolerance - if max_iter is None: - max_iter = self.params.max_iterations - - log.info( - f"VMG: {self.params.n_levels} levels, " - f"pre={self._pre_smooth}, post={self._post_smooth}, " - f"damping={self._damping}" - ) - - # Solve using VMG - time_start = time.time() - finest_level, total_work, converged = solve_vmg( - levels=self._levels, - Re=self.params.Re, - beta_squared=self.params.beta_squared, - lid_velocity=self.params.lid_velocity, - CFL=self.params.CFL, - tolerance=tolerance, - max_cycles=max_iter, - transfer_ops=self._transfer_ops, - corner_treatment=self.corner_treatment, - Lx=self.params.Lx, - Ly=self.params.Ly, - pre_smooth=self._pre_smooth, - post_smooth=self._post_smooth, - damping=self._damping, - ) - - time_end = time.time() - wall_time = time_end - time_start - - # Copy solution from finest level to solver arrays - self.arrays.u[:] = finest_level.u - self.arrays.v[:] = finest_level.v - self.arrays.p[:] = finest_level.p - - # Compute final residuals - self._compute_residuals(self.arrays.u, self.arrays.v, self.arrays.p) - - # Store results - final_residuals = self._compute_algebraic_residuals() - residual_history = [ - { - "rel_iter": tolerance if converged else tolerance * 10, - "u_eq": final_residuals["u_residual"], - "v_eq": final_residuals["v_residual"], - "continuity": final_residuals["continuity_residual"], - } - ] - - self._store_results( - residual_history=residual_history, - final_iter_count=total_work, - is_converged=converged, - wall_time=wall_time, - energy_history=[self._compute_energy()], - enstrophy_history=[self._compute_enstrophy()], - palinstrophy_history=[self._compute_palinstrophy()], - ) - - log.info( - f"VMG completed in {wall_time:.2f}s: {total_work} work units, " - f"converged={converged}" - ) From a74264d7015ae7f522c9ac64d61191b40f26cf50 Mon Sep 17 00:00:00 2001 From: Philip Korsager Nickel Date: Tue, 9 Dec 2025 14:27:10 +0100 Subject: [PATCH 2/3] ghia update --- src/shared/plotting/ldc/convergence.py | 10 ++++- src/shared/plotting/ldc/fields.py | 13 +++++-- src/shared/plotting/ldc/validation.py | 53 +++++++++++++++----------- 3 files changed, 48 insertions(+), 28 deletions(-) diff --git a/src/shared/plotting/ldc/convergence.py b/src/shared/plotting/ldc/convergence.py index 34ef70d..5112f06 100644 --- a/src/shared/plotting/ldc/convergence.py +++ b/src/shared/plotting/ldc/convergence.py @@ -9,6 +9,7 @@ import matplotlib.pyplot as plt import pandas as pd +import seaborn as sns log = logging.getLogger(__name__) @@ -21,6 +22,9 @@ def plot_convergence( log.warning("No timeseries data available for convergence plot") return None + # Set seaborn darkgrid style + sns.set_style("darkgrid") + fig, ax = plt.subplots() # Plot available metrics @@ -41,10 +45,12 @@ def plot_convergence( rf"\textbf{{Convergence History}} --- {solver_label}, $N={N}$, $\mathrm{{Re}}={Re:.0f}$" ) ax.legend(frameon=True) - ax.grid(True, alpha=0.3) + + # Transparent figure, but keep darkgrid axes background + fig.patch.set_alpha(0.0) output_path = output_dir / "convergence.pdf" - fig.savefig(output_path) + fig.savefig(output_path, facecolor=(0, 0, 0, 0)) plt.close(fig) return output_path diff --git a/src/shared/plotting/ldc/fields.py b/src/shared/plotting/ldc/fields.py index d7dbd5d..dc8f7e4 100644 --- a/src/shared/plotting/ldc/fields.py +++ b/src/shared/plotting/ldc/fields.py @@ -11,6 +11,7 @@ import matplotlib.pyplot as plt import numpy as np import pandas as pd +import seaborn as sns from scipy.interpolate import RectBivariateSpline log = logging.getLogger(__name__) @@ -74,7 +75,7 @@ def plot_fields( plt.tight_layout() output_path = output_dir / "fields.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", transparent=True) plt.close(fig) return output_path @@ -144,7 +145,7 @@ def plot_streamlines( plt.tight_layout() output_path = output_dir / "streamlines.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", transparent=True) plt.close(fig) return output_path @@ -174,6 +175,9 @@ def plot_vorticity( dudy = U_spline(y_fine, x_fine, dy=1) vorticity = dvdx - dudy + # Set seaborn darkgrid style + sns.set_style("darkgrid") + fig, ax = plt.subplots(figsize=(7, 6)) vmax = np.max(np.abs(vorticity)) @@ -196,8 +200,11 @@ def plot_vorticity( plt.tight_layout() + # Transparent figure, but keep darkgrid axes background + fig.patch.set_alpha(0.0) + output_path = output_dir / "vorticity.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor=(0, 0, 0, 0)) plt.close(fig) return output_path diff --git a/src/shared/plotting/ldc/validation.py b/src/shared/plotting/ldc/validation.py index 5b35c95..ddd06c0 100644 --- a/src/shared/plotting/ldc/validation.py +++ b/src/shared/plotting/ldc/validation.py @@ -34,7 +34,7 @@ def _build_method_label(sibling: dict) -> str: # Format known solver names label_map = { - "fv": "FV", + "fv": "FV-TVD", "spectral": "Spectral", "spectral_fsg": "Spectral-FSG", "spectral_fmg": "Spectral-FMG", @@ -190,60 +190,64 @@ def plot_ghia_comparison( u_df = pd.DataFrame(u_records) v_df = pd.DataFrame(v_records) + # Set seaborn darkgrid style with transparent figure background + sns.set_style("darkgrid") + sns.set(rc={"figure.facecolor": (0, 0, 0, 0)}) + # Create subplots with publication-quality sizing fig, axes = plt.subplots(1, 2, figsize=(12, 5)) - # Left: u-velocity (vertical centerline) + # Left: u-velocity (vertical centerline) - dashed lines, no markers sns.lineplot( data=u_df, x="u", y="y", hue="Method", style="Method", - markers=True, + markers=False, + dashes=True, ax=axes[0], sort=False, - linewidth=2, - markersize=7, - markevery=15, + linewidth=1, ) + # Ghia reference with markers sns.scatterplot( data=ghia_u, x="u", y="y", marker="o", - s=50, + s=30, facecolors="none", edgecolors="#333333", - linewidths=1.2, + linewidths=1.0, label="Ghia et al. (1982)", ax=axes[0], zorder=10, ) - axes[0].set_xlabel(r"$u$", fontsize=11) - axes[0].set_ylabel(r"$y$", fontsize=11) - axes[0].set_title(r"$u$-velocity (vertical centerline)", fontsize=11) + axes[0].set_xlabel(r"$u$", fontsize=12) + axes[0].set_ylabel(r"$y$", fontsize=12) + #axes[0].set_title(r"$u$ at $x = 0.5$", fontsize=12) - # Right: v-velocity (horizontal centerline) + # Right: v-velocity (horizontal centerline) - dashed lines, no markers sns.lineplot( data=v_df, x="x", y="v", hue="Method", style="Method", - markers=True, + markers=False, + dashes=True, ax=axes[1], sort=False, - linewidth=2, - markersize=7, - markevery=15, + linewidth=1, ) + # Ghia reference with markers sns.scatterplot( data=ghia_v, x="x", y="v", marker="o", - s=50, + s=30, facecolors="none", edgecolors="#333333", linewidths=1.2, @@ -251,18 +255,21 @@ def plot_ghia_comparison( ax=axes[1], zorder=10, ) - axes[1].set_xlabel(r"$x$", fontsize=11) - axes[1].set_ylabel(r"$v$", fontsize=11) - axes[1].set_title(r"$v$-velocity (horizontal centerline)", fontsize=11) + axes[1].set_xlabel(r"$x$", fontsize=12) + axes[1].set_ylabel(r"$v$", fontsize=12) + #axes[1].set_title(r"$v$ at $y = 0.5$", fontsize=12) - # Overall title - fig.suptitle(rf"Ghia Benchmark Comparison ($\mathrm{{Re}} = {int(Re)}$)", fontsize=13, y=1.00) + # Overall title with LaTeX formatting + fig.suptitle(rf"Centerline Velocities - $\mathrm{{Re}} = {int(Re)}$", fontsize=15, y=1.00) # Tight layout for better spacing plt.tight_layout() + # Transparent figure, but keep darkgrid axes background + fig.patch.set_alpha(0.0) + output_path = output_dir / "ghia_comparison.pdf" - fig.savefig(output_path, dpi=300, bbox_inches="tight") + fig.savefig(output_path, dpi=300, bbox_inches="tight", facecolor=(0, 0, 0, 0)) plt.close(fig) log.info(f"Saved comparison plot: {output_path.name}") From 4f43a78bdf2c2dc419f874b2707e8327db5cea72 Mon Sep 17 00:00:00 2001 From: Philip Korsager Nickel Date: Tue, 9 Dec 2025 14:54:10 +0100 Subject: [PATCH 3/3] better plotting --- src/shared/plotting/ldc/style.py | 21 ++++++++++++++++++--- src/solvers/spectral/sg.py | 6 +++--- 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/shared/plotting/ldc/style.py b/src/shared/plotting/ldc/style.py index 2eb8299..33e2161 100644 --- a/src/shared/plotting/ldc/style.py +++ b/src/shared/plotting/ldc/style.py @@ -1,14 +1,29 @@ """ Plotting Style Configuration for LDC Plots. -Uses seaborn darkgrid theme. +Uses seaborn darkgrid theme with LaTeX rendering. """ import logging +import matplotlib.pyplot as plt import seaborn as sns log = logging.getLogger(__name__) -# Use seaborn darkgrid theme -sns.set_theme(style="darkgrid") +# Enable LaTeX rendering +plt.rcParams.update( + { + "text.usetex": True, + "font.family": "serif", + "font.serif": ["Computer Modern Roman"], + "axes.labelsize": 12, + "font.size": 11, + "legend.fontsize": 10, + "xtick.labelsize": 10, + "ytick.labelsize": 10, + } +) + +# Use seaborn darkgrid theme (after rcParams to preserve LaTeX settings) +sns.set_theme(style="darkgrid", rc={"text.usetex": True}) diff --git a/src/solvers/spectral/sg.py b/src/solvers/spectral/sg.py index 9ea88bd..6692b69 100644 --- a/src/solvers/spectral/sg.py +++ b/src/solvers/spectral/sg.py @@ -417,9 +417,9 @@ def step(self): """ a = self.arrays # Shorthand - # Swap buffers at start (for residual calculation in solve()) - a.u, a.u_prev = a.u_prev, a.u - a.v, a.v_prev = a.v_prev, a.v + # Save previous for convergence check (copy, don't swap!) + a.u_prev[:] = a.u + a.v_prev[:] = a.v # Compute adaptive timestep dt = self._compute_adaptive_timestep()