diff --git a/CHANGELOG.md b/CHANGELOG.md index ef321a3..53ece98 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- Added the **Aerodynamics 3D** Advanced Problem: a lightweight incompressible periodic + Cartesian structured-grid solver with 3D Fourier projection, immersed/penalized sphere, + ellipsoid, box, and finite-span NACA 0012 wing obstacles, cached visual result views, + Cd/Cl/Cs histories, and shared volumetric-history performance guardrails. +- Corrected 3D Brinkman penalization, rejected obstacles crossing the periodic box, normalized + streamline integration, and upgraded 3D flow, streamline, and slice result views to cached, + stable-scale animations with continuous analytical obstacle surfaces. - Reordered Advanced Results tabs and display defaults so the most representative visual views appear first, while diagnostics and conservation views remain secondary. - Polished tooltip input-modality handling, clarified optional ODE event detection, improved narrow-window text wrapping, and refined Standard Results sizing and responsive control layout. - Modernized Advanced result view controls with responsive grouped layouts, accessible multi-item selectors, immediate readonly-selector updates, and explicit updates for editable FPUT inputs. diff --git a/docs/complex-problems.md b/docs/complex-problems.md index 6890e81..4e3a289 100644 --- a/docs/complex-problems.md +++ b/docs/complex-problems.md @@ -1,6 +1,6 @@ # Advanced Problems Guide -`Advanced Problems` is a nine-plugin subsystem where each problem contributes a +`Advanced Problems` is a ten-plugin subsystem where each problem contributes a specialized UI, solver, and result dialog. ## Gravitational N-Body Dynamics @@ -68,6 +68,13 @@ Plugins are registered lazily in `src/complex_problems/problem_registry.py`. | `schrodinger_td` | Schrodinger TD (1D/2D) | Split-operator spectral time-dependent Schrodinger solver with configurable potentials and packet states | density/phase animation, momentum spectrum, expectations, invariant curves | | `antenna_radiation` | Antenna Radiation | Far-field patterns for dipole, loop, patch-like aperture, and uniform linear array models | gain/directivity maps, polar cuts, 3D pattern, field metrics | | `aerodynamics_2d` | Aerodynamics 2D | 2D incompressible flow around obstacles using projection/Stokes-style approximations | speed/vorticity/pressure views, drag/lift curves, streamlines | +| `aerodynamics_3d` | Aerodynamics 3D | Incompressible 3D structured-grid flow in a periodic Cartesian domain with immersed/penalized obstacles; lightweight educational/scientific scope, not industrial CFD | velocity-perturbation/velocity/vorticity vector animation, opt-in cached 3D streamline animation, animated XY/XZ/YZ slices, Cd/Cl/Cs and diagnostics | + +The Aerodynamics 3D v1 configuration requires the complete obstacle to remain inside the +Cartesian periodic box; configurations whose geometry crosses an x, y, or z boundary are +rejected. Results retain only `u`, `v`, `w`, and pressure volumetric histories. Derived slice +histories and streamline polylines are prepared lazily in the Results window and are not added +to the solver result object. | `pipe_flow` | Pipe Flow | Steady Darcy-Weisbach and transient 1D pressure-wave pipe-flow models | pressure/velocity profiles, geometry plots, Reynolds/friction metrics, transient maps | | `gravitational_n_body` | Gravitational N-Body Dynamics | Softened Newtonian multi-body systems and curated orbit studies | orbit animation, COM diagnostics, energy and angular-momentum drift | | `fput_experiment` | Fermi-Pasta-Ulam-Tsingou Experiment | Fixed-end alpha/beta FPUT chains | recurrence, modal energy, strain, Hamiltonian drift, scaling fits | diff --git a/src/complex_problems/aerodynamics_3d/__init__.py b/src/complex_problems/aerodynamics_3d/__init__.py new file mode 100644 index 0000000..91c05fe --- /dev/null +++ b/src/complex_problems/aerodynamics_3d/__init__.py @@ -0,0 +1 @@ +"""Three-dimensional incompressible aerodynamics Advanced Problem.""" diff --git a/src/complex_problems/aerodynamics_3d/model.py b/src/complex_problems/aerodynamics_3d/model.py new file mode 100644 index 0000000..1982347 --- /dev/null +++ b/src/complex_problems/aerodynamics_3d/model.py @@ -0,0 +1,398 @@ +"""Geometry and finite-difference helpers for 3D aerodynamics.""" + +from __future__ import annotations + +import numpy as np + +_OBSTACLE_SHAPES = {"sphere", "ellipsoid", "box", "naca0012_wing"} +_BOUNDARY_ERROR = ( + "Aerodynamics 3D v1 requires the obstacle to remain fully inside the periodic box; " + "move the obstacle inward or reduce its dimensions." +) + + +def build_periodic_domain( + *, nx: int, ny: int, nz: int, lx: float, ly: float, lz: float +) -> tuple[ + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + float, + float, + float, +]: + """Build periodic coordinates and meshes stored in ``(z, y, x)`` order.""" + if min(nx, ny, nz) < 4: + raise ValueError("nx, ny, and nz must be >= 4.") + if min(lx, ly, lz) <= 0: + raise ValueError("lx, ly, and lz must be positive.") + x = np.linspace(0.0, lx, nx, endpoint=False) + y = np.linspace(0.0, ly, ny, endpoint=False) + z = np.linspace(0.0, lz, nz, endpoint=False) + zz, yy, xx = np.meshgrid(z, y, x, indexing="ij") + return x, y, z, xx, yy, zz, lx / nx, ly / ny, lz / nz + + +def rotate_coordinates( + X: np.ndarray, + Y: np.ndarray, + *, + center_x: float, + center_y: float, + angle_deg: float, +) -> tuple[np.ndarray, np.ndarray]: + """Rotate x/y coordinates about the obstacle center around the z axis.""" + angle = np.deg2rad(angle_deg) + cosine = np.cos(angle) + sine = np.sin(angle) + xr = X - center_x + yr = Y - center_y + return cosine * xr + sine * yr, -sine * xr + cosine * yr + + +def _naca_thickness(x_local: np.ndarray, chord: float, thickness_ratio: float) -> np.ndarray: + """Return the half-thickness of a symmetric NACA 00xx section.""" + xc = np.clip(x_local / chord, 0.0, 1.0) + return ( + 5.0 + * thickness_ratio + * chord + * (0.2969 * np.sqrt(xc) - 0.1260 * xc - 0.3516 * xc**2 + 0.2843 * xc**3 - 0.1015 * xc**4) + ) + + +def obstacle_half_extents( + *, + shape: str, + diameter: float = 0.4, + size_x: float = 0.7, + size_y: float = 0.4, + size_z: float = 0.4, + chord: float = 0.8, + span: float = 0.8, + thickness_ratio: float = 0.12, + attack_deg: float = 0.0, +) -> tuple[float, float, float]: + """Return world-aligned half extents for an obstacle geometry.""" + obstacle_shape = shape.strip().lower() + if obstacle_shape not in _OBSTACLE_SHAPES: + raise ValueError(f"Unknown obstacle shape '{shape}'.") + positive = { + "diameter": diameter, + "size_x": size_x, + "size_y": size_y, + "size_z": size_z, + "chord": chord, + "span": span, + "thickness_ratio": thickness_ratio, + } + if any(value <= 0 for value in positive.values()): + raise ValueError("Obstacle dimensions and thickness ratio must be positive.") + + if obstacle_shape == "sphere": + radius = diameter / 2.0 + return radius, radius, radius + + if obstacle_shape == "ellipsoid": + local_x, local_y, local_z = size_x / 2.0, size_y / 2.0, size_z / 2.0 + elif obstacle_shape == "box": + local_x, local_y, local_z = size_x / 2.0, size_y / 2.0, size_z / 2.0 + else: + local_x = chord / 2.0 + sample_x = np.linspace(0.0, chord, 2049) + local_y = float(np.max(_naca_thickness(sample_x, chord, thickness_ratio))) + local_z = span / 2.0 + + angle = np.deg2rad(attack_deg) + cosine, sine = abs(float(np.cos(angle))), abs(float(np.sin(angle))) + return cosine * local_x + sine * local_y, sine * local_x + cosine * local_y, local_z + + +def validate_obstacle_fits_domain( + *, + shape: str, + center_x: float, + center_y: float, + center_z: float, + lx: float, + ly: float, + lz: float, + diameter: float = 0.4, + size_x: float = 0.7, + size_y: float = 0.4, + size_z: float = 0.4, + chord: float = 0.8, + span: float = 0.8, + thickness_ratio: float = 0.12, + attack_deg: float = 0.0, +) -> None: + """Reject obstacles whose complete geometry crosses a Cartesian boundary.""" + if min(lx, ly, lz) <= 0: + raise ValueError("Domain lengths must be positive.") + half_x, half_y, half_z = obstacle_half_extents( + shape=shape, + diameter=diameter, + size_x=size_x, + size_y=size_y, + size_z=size_z, + chord=chord, + span=span, + thickness_ratio=thickness_ratio, + attack_deg=attack_deg, + ) + lower = (center_x - half_x, center_y - half_y, center_z - half_z) + upper = (center_x + half_x, center_y + half_y, center_z + half_z) + bounds = ((lower[0], upper[0], lx), (lower[1], upper[1], ly), (lower[2], upper[2], lz)) + if any(lo < 0.0 or hi > limit for lo, hi, limit in bounds): + raise ValueError(_BOUNDARY_ERROR) + + +def _rotate_surface_xy( + local_x: np.ndarray, + local_y: np.ndarray, + *, + center_x: float, + center_y: float, + attack_deg: float, +) -> tuple[np.ndarray, np.ndarray]: + """Rotate local surface coordinates back into world x/y coordinates.""" + angle = np.deg2rad(attack_deg) + cosine, sine = np.cos(angle), np.sin(angle) + return ( + center_x + cosine * local_x - sine * local_y, + center_y + sine * local_x + cosine * local_y, + ) + + +def obstacle_surface_mesh( + *, + shape: str, + center_x: float, + center_y: float, + center_z: float, + diameter: float = 0.4, + size_x: float = 0.7, + size_y: float = 0.4, + size_z: float = 0.4, + chord: float = 0.8, + span: float = 0.8, + thickness_ratio: float = 0.12, + attack_deg: float = 0.0, + resolution: int = 28, +) -> list[tuple[np.ndarray, np.ndarray, np.ndarray]]: + """Return continuous surface patches using the same geometry as the mask.""" + if resolution < 4: + raise ValueError("resolution must be at least 4.") + obstacle_shape = shape.strip().lower() + if obstacle_shape not in _OBSTACLE_SHAPES: + raise ValueError(f"Unknown obstacle shape '{shape}'.") + + patches: list[tuple[np.ndarray, np.ndarray, np.ndarray]] = [] + if obstacle_shape in {"sphere", "ellipsoid"}: + angle = np.linspace(0.0, 2.0 * np.pi, resolution) + polar = np.linspace(0.0, np.pi, max(12, resolution // 2)) + aa, pp = np.meshgrid(angle, polar) + if obstacle_shape == "sphere": + radii = (diameter / 2.0, diameter / 2.0, diameter / 2.0) + else: + radii = (size_x / 2.0, size_y / 2.0, size_z / 2.0) + local_x = radii[0] * np.sin(pp) * np.cos(aa) + local_y = radii[1] * np.sin(pp) * np.sin(aa) + local_z = center_z + radii[2] * np.cos(pp) + world_x, world_y = _rotate_surface_xy( + local_x, + local_y, + center_x=center_x, + center_y=center_y, + attack_deg=attack_deg, + ) + return [(world_x, world_y, local_z)] + + if obstacle_shape == "box": + half_x, half_y, half_z = size_x / 2.0, size_y / 2.0, size_z / 2.0 + values_x = np.linspace(-half_x, half_x, resolution) + values_y = np.linspace(-half_y, half_y, resolution) + values_z = np.linspace(-half_z, half_z, resolution) + face_xy_x, face_xy_y = np.meshgrid(values_x, values_y) + face_xz_x, face_xz_z = np.meshgrid(values_x, values_z) + face_yz_y, face_yz_z = np.meshgrid(values_y, values_z) + for local_x, local_y, local_z in ( + (face_xy_x, face_xy_y, np.full_like(face_xy_x, -half_z)), + (face_xy_x, face_xy_y, np.full_like(face_xy_x, half_z)), + (face_xz_x, np.full_like(face_xz_x, -half_y), face_xz_z), + (face_xz_x, np.full_like(face_xz_x, half_y), face_xz_z), + (np.full_like(face_yz_y, -half_x), face_yz_y, face_yz_z), + (np.full_like(face_yz_y, half_x), face_yz_y, face_yz_z), + ): + world_x, world_y = _rotate_surface_xy( + local_x, + local_y, + center_x=center_x, + center_y=center_y, + attack_deg=attack_deg, + ) + patches.append((world_x, world_y, center_z + local_z)) + return patches + + local_x_values = np.linspace(0.0, chord, resolution) + local_z_values = np.linspace(-span / 2.0, span / 2.0, max(12, resolution // 2)) + local_x, local_z = np.meshgrid(local_x_values, local_z_values) + local_thickness = _naca_thickness(local_x, chord, thickness_ratio) + for sign in (-1.0, 1.0): + world_x, world_y = _rotate_surface_xy( + local_x - chord / 2.0, + sign * local_thickness, + center_x=center_x, + center_y=center_y, + attack_deg=attack_deg, + ) + patches.append((world_x, world_y, center_z + local_z)) + return patches + + +def build_obstacle_mask( + *, + shape: str, + X: np.ndarray, + Y: np.ndarray, + Z: np.ndarray, + center_x: float, + center_y: float, + center_z: float, + diameter: float = 0.4, + size_x: float = 0.7, + size_y: float = 0.4, + size_z: float = 0.4, + chord: float = 0.8, + span: float = 0.8, + thickness_ratio: float = 0.12, + attack_deg: float = 0.0, +) -> tuple[np.ndarray, float, float]: + """Build an immersed obstacle mask and its reference area and length. + + All dimensions are full dimensions. The returned mask has the same + ``(z, y, x)`` shape as the coordinate arrays. + """ + obstacle_shape = shape.strip().lower() + if obstacle_shape not in _OBSTACLE_SHAPES: + raise ValueError(f"Unknown obstacle shape '{shape}'.") + positive = { + "diameter": diameter, + "size_x": size_x, + "size_y": size_y, + "size_z": size_z, + "chord": chord, + "span": span, + "thickness_ratio": thickness_ratio, + } + if any(value <= 0 for value in positive.values()): + raise ValueError("Obstacle dimensions and thickness ratio must be positive.") + + xp, yp = rotate_coordinates(X, Y, center_x=center_x, center_y=center_y, angle_deg=attack_deg) + zp = Z - center_z + if obstacle_shape == "sphere": + radius = diameter / 2.0 + mask = xp * xp + yp * yp + zp * zp <= radius * radius + reference_area = np.pi * radius * radius + reference_length = diameter + elif obstacle_shape == "ellipsoid": + mask = (xp / (size_x / 2.0)) ** 2 + (yp / (size_y / 2.0)) ** 2 + ( + zp / (size_z / 2.0) + ) ** 2 <= 1.0 + reference_area = np.pi * (size_y / 2.0) * (size_z / 2.0) + reference_length = size_x + elif obstacle_shape == "box": + mask = ( + (np.abs(xp) <= size_x / 2.0) + & (np.abs(yp) <= size_y / 2.0) + & (np.abs(zp) <= size_z / 2.0) + ) + reference_area = size_y * size_z + reference_length = size_x + else: + x_local = xp + chord / 2.0 + mask = ( + (x_local >= 0.0) + & (x_local <= chord) + & (np.abs(yp) <= _naca_thickness(x_local, chord, thickness_ratio)) + & (np.abs(zp) <= span / 2.0) + ) + reference_area = chord * span + reference_length = chord + + return mask.astype(bool), float(reference_area), float(reference_length) + + +def ddx_periodic(field: np.ndarray, dx: float) -> np.ndarray: + """Centered periodic derivative along the last (x) axis.""" + result = np.empty_like(field, dtype=float) + result[..., 1:-1] = (field[..., 2:] - field[..., :-2]) / (2.0 * dx) + result[..., 0] = (field[..., 1] - field[..., -1]) / (2.0 * dx) + result[..., -1] = (field[..., 0] - field[..., -2]) / (2.0 * dx) + return result + + +def ddy_periodic(field: np.ndarray, dy: float) -> np.ndarray: + """Centered periodic derivative along the y axis.""" + result = np.empty_like(field, dtype=float) + result[:, 1:-1, ...] = (field[:, 2:, ...] - field[:, :-2, ...]) / (2.0 * dy) + result[:, 0, ...] = (field[:, 1, ...] - field[:, -1, ...]) / (2.0 * dy) + result[:, -1, ...] = (field[:, 0, ...] - field[:, -2, ...]) / (2.0 * dy) + return result + + +def ddz_periodic(field: np.ndarray, dz: float) -> np.ndarray: + """Centered periodic derivative along the first (z) axis.""" + result = np.empty_like(field, dtype=float) + result[1:-1, ...] = (field[2:, ...] - field[:-2, ...]) / (2.0 * dz) + result[0, ...] = (field[1, ...] - field[-1, ...]) / (2.0 * dz) + result[-1, ...] = (field[0, ...] - field[-2, ...]) / (2.0 * dz) + return result + + +def laplacian_periodic(field: np.ndarray, dx: float, dy: float, dz: float) -> np.ndarray: + """Periodic centered finite-difference Laplacian.""" + return ( + (np.roll(field, -1, axis=-1) - 2.0 * field + np.roll(field, 1, axis=-1)) / dx**2 + + (np.roll(field, -1, axis=-2) - 2.0 * field + np.roll(field, 1, axis=-2)) / dy**2 + + (np.roll(field, -1, axis=-3) - 2.0 * field + np.roll(field, 1, axis=-3)) / dz**2 + ) + + +def divergence_periodic( + u: np.ndarray, v: np.ndarray, w: np.ndarray, dx: float, dy: float, dz: float +) -> np.ndarray: + """Centered finite-difference divergence for diagnostics outside projection.""" + return ddx_periodic(u, dx) + ddy_periodic(v, dy) + ddz_periodic(w, dz) + + +def spectral_divergence_l2( + u: np.ndarray, + v: np.ndarray, + w: np.ndarray, + kx: np.ndarray, + ky: np.ndarray, + kz: np.ndarray, + *, + fluid_mask: np.ndarray | None = None, +) -> float: + """Return the L2 norm of the Fourier divergence.""" + divergence_hat = 1j * (kx * np.fft.fftn(u) + ky * np.fft.fftn(v) + kz * np.fft.fftn(w)) + divergence = np.real(np.fft.ifftn(divergence_hat)) + if fluid_mask is not None and np.any(fluid_mask): + divergence = divergence[fluid_mask] + return float(np.sqrt(np.mean(divergence * divergence))) + + +def vorticity_periodic( + u: np.ndarray, v: np.ndarray, w: np.ndarray, dx: float, dy: float, dz: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return the centered finite-difference vorticity vector.""" + return ( + ddy_periodic(w, dy) - ddz_periodic(v, dz), + ddz_periodic(u, dz) - ddx_periodic(w, dx), + ddx_periodic(v, dx) - ddy_periodic(u, dy), + ) diff --git a/src/complex_problems/aerodynamics_3d/problem.py b/src/complex_problems/aerodynamics_3d/problem.py new file mode 100644 index 0000000..00e34c2 --- /dev/null +++ b/src/complex_problems/aerodynamics_3d/problem.py @@ -0,0 +1,33 @@ +"""Plugin entrypoint for 3D aerodynamics simulations.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from complex_problems.base import ProblemDescriptor + +if TYPE_CHECKING: + from tkinter import Tk, Toplevel + + +class Aerodynamics3DProblem: + """Complex problem plugin for lightweight 3D incompressible flow.""" + + descriptor = ProblemDescriptor( + id="aerodynamics_3d", + name="Aerodynamics 3D", + description=( + "Incompressible 3D structured-grid flow in a periodic Cartesian domain " + "with immersed/penalized obstacles; a lightweight educational/scientific " + "solver, not an industrial CFD replacement." + ), + ) + + def open_dialog(self, parent: "Tk | Toplevel") -> None: + """Open the aerodynamics configuration dialog.""" + from complex_problems.aerodynamics_3d.ui import Aerodynamics3DDialog + + Aerodynamics3DDialog(parent) + + +PROBLEM = Aerodynamics3DProblem() diff --git a/src/complex_problems/aerodynamics_3d/result_dialog.py b/src/complex_problems/aerodynamics_3d/result_dialog.py new file mode 100644 index 0000000..cd9c31a --- /dev/null +++ b/src/complex_problems/aerodynamics_3d/result_dialog.py @@ -0,0 +1,1106 @@ +"""Advanced result views for the 3D aerodynamics solver.""" + +from __future__ import annotations + +import tkinter as tk +from dataclasses import dataclass +from pathlib import Path +from tkinter import filedialog, messagebox, ttk +from typing import TYPE_CHECKING, Any, Callable + +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +from complex_problems.aerodynamics_3d.model import obstacle_surface_mesh, vorticity_periodic +from complex_problems.aerodynamics_3d.solver import Aerodynamics3DResult +from complex_problems.common.result_dialog_ui import ( + AdvancedResultShell, + AdvancedResultSize, + close_embedded_figures, + format_result_summary, + make_view_controls, + reset_embedded_animation, +) +from config import generate_output_basename, get_output_dir +from frontend.plot_embed import embed_animation_plot_in_tk, embed_plot_in_tk +from frontend.theme import get_font +from frontend.window_utils import make_modal +from plotting import ( + create_image_animation_plot, + create_solution_plot, + export_animated_figure_to_mp4, +) +from plotting.animation_metadata import attach_animation_metadata +from utils import get_logger + +if TYPE_CHECKING: + from matplotlib.figure import Figure + +logger = get_logger(__name__) + +_FLOW_DISPLAYS = ("Velocity perturbation", "Velocity vectors", "Vorticity vectors") +_SLICE_FIELDS = ("Vorticity magnitude", "Speed", "Pressure", "u", "v", "w") + + +def derived_speed(u: np.ndarray, v: np.ndarray, w: np.ndarray) -> np.ndarray: + """Compute speed without retaining a separate history.""" + return np.sqrt(u * u + v * v + w * w) + + +def derived_vorticity( + u: np.ndarray, v: np.ndarray, w: np.ndarray, *, dx: float, dy: float, dz: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Compute vorticity components and magnitude on demand.""" + omega = vorticity_periodic(u, v, w, dx, dy, dz) + magnitude = np.sqrt(sum(component * component for component in omega)) + return (*omega, magnitude) + + +def _velocity_interpolator(result: Aerodynamics3DResult, frame: int) -> RegularGridInterpolator: + """Create one vector-valued interpolator for a saved velocity frame.""" + axes = (result.z, result.y, result.x) + values = np.stack((result.u[frame], result.v[frame], result.w[frame]), axis=-1) + return RegularGridInterpolator( + axes, + values, + bounds_error=False, + fill_value=np.nan, + ) + + +def _nearest_uniform_indices(values: np.ndarray, coordinates: np.ndarray) -> np.ndarray: + """Return nearest indexes on a uniformly spaced coordinate axis.""" + if len(values) < 2: + return np.zeros(len(coordinates), dtype=int) + spacing = float(values[1] - values[0]) + indexes = np.rint((coordinates - float(values[0])) / spacing).astype(int) + return np.clip(indexes, 0, len(values) - 1) + + +def _points_are_obstacles(result: Aerodynamics3DResult, points: np.ndarray) -> np.ndarray: + """Check many points against the nearest stored obstacle cells.""" + if len(points) == 0: + return np.zeros(0, dtype=bool) + ix = _nearest_uniform_indices(result.x, points[:, 0]) + iy = _nearest_uniform_indices(result.y, points[:, 1]) + iz = _nearest_uniform_indices(result.z, points[:, 2]) + return result.obstacle_mask[iz, iy, ix] + + +def _batch_tangent( + interpolator: RegularGridInterpolator, + points: np.ndarray, + *, + lower: np.ndarray, + upper: np.ndarray, + speed_epsilon: float, +) -> tuple[np.ndarray, np.ndarray]: + """Interpolate and normalize a batch of RK stages in one SciPy call.""" + tangents = np.zeros_like(points) + in_bounds = np.all((points >= lower) & (points <= upper), axis=1) + if not np.any(in_bounds): + return tangents, in_bounds + values = np.asarray(interpolator(points[in_bounds][:, (2, 1, 0)]), dtype=float) + magnitudes = np.linalg.norm(values, axis=1) + finite = np.all(np.isfinite(values), axis=1) & np.isfinite(magnitudes) + usable = finite & (magnitudes > speed_epsilon) + valid_indexes = np.flatnonzero(in_bounds) + usable_indexes = valid_indexes[usable] + tangents[usable_indexes] = values[usable] / magnitudes[usable, None] + valid = np.zeros(len(points), dtype=bool) + valid[usable_indexes] = True + return tangents, valid + + +def trace_streamlines_3d( + seeds: np.ndarray, + *, + result: Aerodynamics3DResult, + frame: int = 0, + step_size: float | None = None, + max_steps: int = 300, + speed_epsilon: float = 1.0e-8, + _interpolator: RegularGridInterpolator | None = None, +) -> tuple[np.ndarray, ...]: + """Trace many normalized-tangent RK4 streamlines in vectorized batches.""" + if not 0 <= frame < len(result.t): + raise ValueError("frame is outside the saved result history.") + seed_points = np.asarray(seeds, dtype=float) + if seed_points.ndim != 2 or seed_points.shape[1] != 3: + raise ValueError("seeds must have shape (n, 3).") + if len(seed_points) == 0: + return () + if max_steps < 1: + return tuple(point[None, :] for point in seed_points) + step = ( + step_size + or min( + result.x[1] - result.x[0], + result.y[1] - result.y[0], + result.z[1] - result.z[0], + ) + * 0.45 + ) + if step <= 0: + raise ValueError("step_size must be positive.") + + interpolator = _interpolator or _velocity_interpolator(result, frame) + lower = np.array([result.x[0], result.y[0], result.z[0]], dtype=float) + upper = np.array([result.x[-1], result.y[-1], result.z[-1]], dtype=float) + current = seed_points.copy() + histories: list[list[np.ndarray]] = [[point.copy()] for point in seed_points] + active = ~_points_are_obstacles(result, current) + active &= np.all((current >= lower) & (current <= upper), axis=1) + + for _ in range(max_steps): + indexes = np.flatnonzero(active) + if len(indexes) == 0: + break + points = current[indexes] + + k1, valid = _batch_tangent( + interpolator, + points, + lower=lower, + upper=upper, + speed_epsilon=speed_epsilon, + ) + active[indexes[~valid]] = False + indexes, points, k1 = indexes[valid], points[valid], k1[valid] + if len(indexes) == 0: + continue + + k2, valid = _batch_tangent( + interpolator, + points + 0.5 * step * k1, + lower=lower, + upper=upper, + speed_epsilon=speed_epsilon, + ) + active[indexes[~valid]] = False + indexes, points, k1, k2 = indexes[valid], points[valid], k1[valid], k2[valid] + if len(indexes) == 0: + continue + + k3, valid = _batch_tangent( + interpolator, + points + 0.5 * step * k2, + lower=lower, + upper=upper, + speed_epsilon=speed_epsilon, + ) + active[indexes[~valid]] = False + indexes, points, k1, k2, k3 = ( + indexes[valid], + points[valid], + k1[valid], + k2[valid], + k3[valid], + ) + if len(indexes) == 0: + continue + + k4, valid = _batch_tangent( + interpolator, + points + step * k3, + lower=lower, + upper=upper, + speed_epsilon=speed_epsilon, + ) + active[indexes[~valid]] = False + indexes, points, k1, k2, k3, k4 = ( + indexes[valid], + points[valid], + k1[valid], + k2[valid], + k3[valid], + k4[valid], + ) + if len(indexes) == 0: + continue + + next_points = points + step * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0 + valid_next = np.all((next_points >= lower) & (next_points <= upper), axis=1) + if np.any(valid_next): + valid_next &= ~_points_are_obstacles(result, next_points) + active[indexes[~valid_next]] = False + survivor_indexes = indexes[valid_next] + survivor_points = next_points[valid_next] + current[survivor_indexes] = survivor_points + for streamline_index, point in zip(survivor_indexes, survivor_points): + histories[int(streamline_index)].append(point.copy()) + + return tuple(np.asarray(points, dtype=float) for points in histories) + + +def trace_streamline_3d( + seed: tuple[float, float, float], + *, + result: Aerodynamics3DResult, + frame: int = 0, + step_size: float | None = None, + max_steps: int = 300, + speed_epsilon: float = 1.0e-8, + _interpolator: RegularGridInterpolator | None = None, +) -> np.ndarray: + """Trace one streamline using the vectorized RK4 implementation.""" + return trace_streamlines_3d( + np.asarray([seed], dtype=float), + result=result, + frame=frame, + step_size=step_size, + max_steps=max_steps, + speed_epsilon=speed_epsilon, + _interpolator=_interpolator, + )[0] + + +def make_streamline_seeds( + result: Aerodynamics3DResult, density: int = 4 +) -> list[tuple[float, float, float]]: + """Seed a modest y-z grid near the upstream x boundary.""" + count = max(2, min(9, int(density))) + x_seed = float(result.x[0] + 1.5 * (result.x[1] - result.x[0])) + ys = np.linspace(result.y[0], result.y[-1], count + 2)[1:-1] + zs = np.linspace(result.z[0], result.z[-1], count + 2)[1:-1] + return [(x_seed, float(y_value), float(z_value)) for z_value in zs for y_value in ys] + + +@dataclass(frozen=True) +class StreamlineAnimationCache: + """Cached streamline polylines for every saved solver frame.""" + + seed_density: int + lines_by_frame: tuple[tuple[np.ndarray, ...], ...] + + @property + def frames(self) -> tuple[tuple[np.ndarray, ...], ...]: + """Expose frame lines with a concise name for callers and tests.""" + return self.lines_by_frame + + @property + def frame_count(self) -> int: + """Return the deterministic number of cached animation frames.""" + return len(self.lines_by_frame) + + +def build_streamline_cache( + result: Aerodynamics3DResult, density: int = 4 +) -> StreamlineAnimationCache: + """Trace and cache all streamline frames while reusing one interpolator per frame.""" + seeds = make_streamline_seeds(result, density) + frames: list[tuple[np.ndarray, ...]] = [] + seed_points = np.asarray(seeds, dtype=float) + for frame in range(len(result.t)): + interpolator = _velocity_interpolator(result, frame) + frames.append( + trace_streamlines_3d( + seed_points, + result=result, + frame=frame, + _interpolator=interpolator, + ) + ) + return StreamlineAnimationCache( + seed_density=int(density), + lines_by_frame=tuple(frames), + ) + + +def slice_index_count(result: Aerodynamics3DResult, plane: str) -> int: + """Return the valid number of slice indexes for one plane.""" + counts = {"XY": len(result.z), "XZ": len(result.y), "YZ": len(result.x)} + try: + return counts[plane.upper()] + except KeyError as exc: + raise ValueError("plane must be XY, XZ, or YZ") from exc + + +def slice_center_index(result: Aerodynamics3DResult, plane: str) -> int: + """Return the central valid index for one slice plane.""" + return slice_index_count(result, plane) // 2 + + +def _slice_data( + result: Aerodynamics3DResult, data: np.ndarray, *, plane: str, index: int +) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, str]: + """Select one physical-coordinate 2D slice.""" + plane_name = plane.upper() + selected = max(0, min(int(index), slice_index_count(result, plane_name) - 1)) + if plane_name == "XY": + return ( + result.x, + result.y, + data[selected], + float(result.z[selected]), + f"z = {result.z[selected]:.4g}", + ) + if plane_name == "XZ": + return ( + result.x, + result.z, + data[:, selected, :], + float(result.y[selected]), + f"y = {result.y[selected]:.4g}", + ) + return ( + result.y, + result.z, + data[:, :, selected], + float(result.x[selected]), + f"x = {result.x[selected]:.4g}", + ) + + +def slice_field( + result: Aerodynamics3DResult, *, frame: int, plane: str, index: int, field: str +) -> tuple[np.ndarray, np.ndarray, np.ndarray, str]: + """Return one physical-coordinate slice for the selected plane and field.""" + if not 0 <= frame < len(result.t): + raise ValueError("frame is outside the saved result history.") + if field == "Speed": + data = derived_speed(result.u[frame], result.v[frame], result.w[frame]) + elif field == "Vorticity magnitude": + _, _, _, data = derived_vorticity( + result.u[frame], + result.v[frame], + result.w[frame], + dx=float(result.x[1] - result.x[0]), + dy=float(result.y[1] - result.y[0]), + dz=float(result.z[1] - result.z[0]), + ) + elif field in {"Pressure", "u", "v", "w"}: + data = {"Pressure": result.pressure, "u": result.u, "v": result.v, "w": result.w}[field][ + frame + ] + else: + raise ValueError(f"Unknown slice field '{field}'.") + x_axis, y_axis, values, _, coordinate = _slice_data(result, data, plane=plane, index=index) + return x_axis, y_axis, values, coordinate + + +@dataclass(frozen=True) +class SliceAnimationPayload: + """One cached 2D time history prepared for a single slice selector.""" + + t: np.ndarray + x_axis: np.ndarray + y_axis: np.ndarray + frames: np.ndarray + plane: str + index: int + field: str + coordinate: float + coordinate_label: str + + +def prepare_slice_history( + result: Aerodynamics3DResult, *, plane: str, index: int, field: str +) -> SliceAnimationPayload: + """Prepare only the requested 2D history, deriving fields frame by frame.""" + plane_name = plane.upper() + selected = max(0, min(int(index), slice_index_count(result, plane_name) - 1)) + slices: list[np.ndarray] = [] + x_axis: np.ndarray | None = None + y_axis: np.ndarray | None = None + coordinate = 0.0 + coordinate_label = "" + for frame in range(len(result.t)): + x_axis, y_axis, values, coordinate_label = slice_field( + result, frame=frame, plane=plane_name, index=selected, field=field + ) + slices.append(np.asarray(values, dtype=float)) + if plane_name == "XY": + coordinate = float(result.z[selected]) + elif plane_name == "XZ": + coordinate = float(result.y[selected]) + else: + coordinate = float(result.x[selected]) + if x_axis is None or y_axis is None: + raise ValueError("The result history must contain at least one frame.") + return SliceAnimationPayload( + t=result.t, + x_axis=x_axis, + y_axis=y_axis, + frames=np.stack(slices), + plane=plane_name, + index=selected, + field=field, + coordinate=coordinate, + coordinate_label=coordinate_label, + ) + + +@dataclass(frozen=True) +class _FlowPayload: + """Cached data used by both the embedded flow animation and MP4 export.""" + + result: Aerodynamics3DResult + display: str + arrow_density: int + magnitude_min: float + magnitude_max: float + + +@dataclass(frozen=True) +class _StreamlinePayload: + """Cached streamlines shared by the embedded view and MP4 export.""" + + result: Aerodynamics3DResult + cache: StreamlineAnimationCache + + +def _obstacle_kwargs(result: Aerodynamics3DResult) -> dict[str, Any]: + """Read continuous obstacle parameters from solver metadata.""" + metadata = result.metadata + return { + "shape": metadata.get("obstacle_shape", "sphere"), + "center_x": float(metadata.get("obstacle_center_x", result.x[len(result.x) // 2])), + "center_y": float(metadata.get("obstacle_center_y", result.y[len(result.y) // 2])), + "center_z": float(metadata.get("obstacle_center_z", result.z[len(result.z) // 2])), + "diameter": float(metadata.get("obstacle_diameter", 0.4)), + "size_x": float(metadata.get("obstacle_size_x", 0.7)), + "size_y": float(metadata.get("obstacle_size_y", 0.4)), + "size_z": float(metadata.get("obstacle_size_z", 0.4)), + "chord": float(metadata.get("obstacle_chord", 0.8)), + "span": float(metadata.get("obstacle_span", 0.8)), + "thickness_ratio": float(metadata.get("obstacle_thickness_ratio", 0.12)), + "attack_deg": float(metadata.get("obstacle_attack_deg", 0.0)), + } + + +def _draw_obstacle_surface(axis: Any, result: Aerodynamics3DResult) -> None: + """Draw the analytical continuous obstacle surface.""" + for surface_x, surface_y, surface_z in obstacle_surface_mesh( + **_obstacle_kwargs(result), resolution=28 + ): + axis.plot_surface( + surface_x, + surface_y, + surface_z, + color="#65717f", + alpha=0.68, + linewidth=0.0, + antialiased=True, + shade=True, + ) + + +def _decimated_indices( + shape: tuple[int, int, int], target: int = 900 +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return a regular subsampling whose product is close to the arrow target.""" + stride = max(1, int(np.ceil((np.prod(shape) / target) ** (1.0 / 3.0)))) + indices = np.indices(shape) + return ( + indices[0, ::stride, ::stride, ::stride], + indices[1, ::stride, ::stride, ::stride], + indices[2, ::stride, ::stride, ::stride], + ) + + +def flow_vector_components( + result: Aerodynamics3DResult, *, frame: int, display: str +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Return the selected visualization-only vector field for one frame.""" + if display == "Velocity perturbation": + u_inf = float(result.metadata.get("u_inf", 0.0)) + return result.u[frame] - u_inf, result.v[frame], result.w[frame] + if display == "Velocity vectors": + return result.u[frame], result.v[frame], result.w[frame] + if display == "Vorticity vectors": + return derived_vorticity( + result.u[frame], + result.v[frame], + result.w[frame], + dx=float(result.x[1] - result.x[0]), + dy=float(result.y[1] - result.y[0]), + dz=float(result.z[1] - result.z[0]), + )[:3] + raise ValueError(f"Unknown flow display '{display}'.") + + +def _arrow_length_bounds(result: Aerodynamics3DResult) -> tuple[float, float]: + """Return fixed visual arrow-length bounds for one result domain.""" + spans = ( + float(result.x[-1] - result.x[0]), + float(result.y[-1] - result.y[0]), + float(result.z[-1] - result.z[0]), + ) + reference = max(min(spans), 1.0e-12) + return 0.02 * reference, 0.08 * reference + + +def _relative_magnitudes( + magnitude: np.ndarray, + magnitude_min: float, + magnitude_max: float, +) -> np.ndarray: + """Normalize magnitudes against one global range shared by all frames.""" + lower = float(magnitude_min) + upper = float(magnitude_max) + if upper <= lower: + return np.full_like(magnitude, 0.5, dtype=float) + return np.clip((magnitude - lower) / (upper - lower), 0.0, 1.0) + + +def _scaled_arrow_lengths( + magnitude: np.ndarray, + magnitude_min: float, + magnitude_max: float, + *, + minimum: float, + maximum: float, + bias: float = 15.0, +) -> np.ndarray: + """Map a global magnitude range logarithmically into fixed visual bounds.""" + relative = _relative_magnitudes(magnitude, magnitude_min, magnitude_max) + normalized = np.log1p(bias * relative) / np.log1p(bias) + return minimum + (maximum - minimum) * normalized + + +def _flow_magnitude_range( + result: Aerodynamics3DResult, + display: str, + arrow_density: int, +) -> tuple[float, float]: + """Return global nonzero min/max for all arrows displayed across all frames.""" + indices = _decimated_indices( + result.obstacle_mask.shape, + target=max(300, arrow_density * 250), + ) + zz, yy, xx = indices + fluid = ~result.obstacle_mask[zz, yy, xx] + minimum = np.inf + maximum = 0.0 + for frame in range(len(result.t)): + components = flow_vector_components(result, frame=frame, display=display) + sampled = tuple(component[zz, yy, xx][fluid] for component in components) + magnitude = derived_speed(*sampled) + usable = magnitude[np.isfinite(magnitude) & (magnitude > 1.0e-12)] + if usable.size: + minimum = min(minimum, float(np.min(usable))) + maximum = max(maximum, float(np.max(usable))) + if not np.isfinite(minimum) or maximum <= 0.0: + return 0.0, 1.0 + return minimum, maximum + + +def _create_flow_figure(payload: _FlowPayload) -> Figure: + """Build an animated 3D vector view from cached result arrays.""" + import matplotlib.pyplot as plt + from matplotlib.colors import LinearSegmentedColormap + + result = payload.result + figure = plt.figure() + axis: Any = figure.add_subplot(111, projection="3d") + indices = _decimated_indices( + result.obstacle_mask.shape, target=max(300, payload.arrow_density * 250) + ) + zz, yy, xx = indices + fluid = ~result.obstacle_mask[zz, yy, xx] + xs, ys, zs = result.x[xx][fluid], result.y[yy][fluid], result.z[zz][fluid] + cmap = LinearSegmentedColormap.from_list( + "aero_vectors", ("#14253d", "#1f4e79", "#2c7fb8", "#66c2a5") + ) + + def draw(frame: int) -> None: + elevation, azimuth = axis.elev, axis.azim + axis.clear() + _draw_obstacle_surface(axis, result) + u, v, w = flow_vector_components(result, frame=frame, display=payload.display) + u, v, w = u[zz, yy, xx][fluid], v[zz, yy, xx][fluid], w[zz, yy, xx][fluid] + magnitude = derived_speed(u, v, w) + valid = np.isfinite(magnitude) & (magnitude > 1.0e-12) + directions = [np.zeros_like(component) for component in (u, v, w)] + for direction, component in zip(directions, (u, v, w)): + np.divide(component, magnitude, out=direction, where=valid) + relative = _relative_magnitudes( + magnitude, + payload.magnitude_min, + payload.magnitude_max, + ) + colors = cmap(0.05 + 0.75 * relative) + minimum_length, maximum_length = _arrow_length_bounds(result) + arrow_lengths = _scaled_arrow_lengths( + magnitude, + payload.magnitude_min, + payload.magnitude_max, + minimum=minimum_length, + maximum=maximum_length, + ) + visible = valid + axis.quiver( + xs[visible], + ys[visible], + zs[visible], + directions[0][visible] * arrow_lengths[visible], + directions[1][visible] * arrow_lengths[visible], + directions[2][visible] * arrow_lengths[visible], + colors=colors[visible], + length=1.0, + normalize=False, + linewidth=0.9, + alpha=0.86, + ) + axis.set_xlim(result.x[0], result.x[-1]) + axis.set_ylim(result.y[0], result.y[-1]) + axis.set_zlim(result.z[0], result.z[-1]) + axis.set_box_aspect( + (result.x[-1] - result.x[0], result.y[-1] - result.y[0], result.z[-1] - result.z[0]) + ) + axis.set_xlabel("x") + axis.set_ylabel("y") + axis.set_zlabel("z") + axis.grid(False) + axis.set_title(f"{payload.display} (t={result.t[frame]:.3g})") + axis.view_init(elev=elevation, azim=azimuth) + + draw(0) + + def update(frame: int) -> None: + draw(max(0, min(frame, len(result.t) - 1))) + figure.canvas.draw_idle() + + return attach_animation_metadata( + figure, + update=update, + n_points=len(result.t), + frame_label="t", + frame_coordinates=result.t, + ) + + +def _create_streamline_figure(payload: _StreamlinePayload) -> Figure: + """Build a 3D animation from already cached streamline polylines.""" + import matplotlib.pyplot as plt + + result = payload.result + cache = payload.cache + figure = plt.figure() + axis: Any = figure.add_subplot(111, projection="3d") + _draw_obstacle_surface(axis, result) + line_artists: list[Any] = [] + + def draw_lines(frame: int) -> None: + line_artists.clear() + for line in cache.lines_by_frame[frame]: + if len(line) > 1: + (artist,) = axis.plot( + line[:, 0], + line[:, 1], + line[:, 2], + color="#1f4e79", + linewidth=1.8, + alpha=0.92, + ) + line_artists.append(artist) + + draw_lines(0) + axis.set_xlim(result.x[0], result.x[-1]) + axis.set_ylim(result.y[0], result.y[-1]) + axis.set_zlim(result.z[0], result.z[-1]) + axis.set_box_aspect( + (result.x[-1] - result.x[0], result.y[-1] - result.y[0], result.z[-1] - result.z[0]) + ) + axis.set_xlabel("x") + axis.set_ylabel("y") + axis.set_zlabel("z") + axis.grid(False) + axis.set_title(f"3D streamlines (t={result.t[0]:.3g})") + figure.tight_layout() + + def update(frame: int) -> None: + index = max(0, min(frame, cache.frame_count - 1)) + for artist in line_artists: + artist.remove() + draw_lines(index) + axis.set_title(f"3D streamlines (t={result.t[index]:.3g})") + figure.canvas.draw_idle() + + return attach_animation_metadata( + figure, + update=update, + n_points=cache.frame_count, + frame_label="t", + frame_coordinates=result.t, + ) + + +def _line_limits(values: np.ndarray) -> tuple[float, float]: + """Return readable limits for a possibly constant diagnostic.""" + lower, upper = float(np.min(values)), float(np.max(values)) + if lower == upper: + margin = max(1.0, abs(lower) * 0.1) + return lower - margin, upper + margin + margin = 0.05 * (upper - lower) + return lower - margin, upper + margin + + +class Aerodynamics3DResultDialog: + """Tabbed visual workspace for a completed 3D aerodynamics solve.""" + + def __init__(self, parent: tk.Tk | tk.Toplevel, *, result: Aerodynamics3DResult) -> None: + self.parent = parent + self._result = result + self._flow_range_cache: dict[tuple[str, int], tuple[float, float]] = {} + self._streamline_cache: dict[int, StreamlineAnimationCache] = {} + self._slice_cache: dict[tuple[str, int, str], SliceAnimationPayload] = {} + self._canvases: dict[str, object | None] = {} + self._stream_initialized = False + self.win = tk.Toplevel(parent) + self.win.title("Aerodynamics 3D Results") + self._build_ui() + self._shell.finish(AdvancedResultSize(1460, 940, 980, 650)) + make_modal(self.win, parent) + + def _on_close(self) -> None: + close_embedded_figures(self, tuple(f"_canvas_{name}" for name in self._canvases)) + self.win.destroy() + + def _build_ui(self) -> None: + magnitudes = self._result.magnitudes + summary = format_result_summary( + ( + ("Reynolds number", f"{magnitudes['reynolds']:.1f}"), + ("Mean tail Cd", f"{magnitudes['mean_cd_tail']:+.3e}"), + ("RMS Cl", f"{magnitudes['rms_cl']:.3e}"), + ("RMS Cs", f"{magnitudes['rms_cs']:.3e}"), + ("Maximum L2 divergence", f"{magnitudes['max_divergence_l2']:.3e}"), + ) + ) + self._shell = AdvancedResultShell( + self.win, title="Aerodynamics 3D Results", summary=summary, close_command=self._on_close + ) + notebook = self._shell.notebook + tabs: list[tuple[str, Callable[[ttk.Frame], None]]] = [ + ("3D Flow", self._build_flow_tab), + ("3D Streamlines", self._build_streamlines_tab), + ("Slices", self._build_slices_tab), + ("Forces", self._build_forces_tab), + ("Diagnostics", self._build_diagnostics_tab), + ] + for label, builder in tabs: + tab = ttk.Frame(notebook) + notebook.add(tab, text=label) + if label == "3D Streamlines": + self._stream_tab = tab + else: + builder(tab) + notebook.bind("<>", self._on_tab_changed) + + def _on_tab_changed(self, _event: object) -> None: + if ( + not self._stream_initialized + and self._shell.notebook.tab("current", "text") == "3D Streamlines" + ): + self._stream_initialized = True + self._build_streamlines_tab(self._stream_tab) + + def _build_flow_tab(self, parent: ttk.Frame) -> None: + controls = make_view_controls(parent) + group = controls.add_group(requested_width=260) + self._flow_display = tk.StringVar(value="Velocity perturbation") + ttk.Label(group, text="Display:", style="Small.TLabel").pack(side=tk.LEFT, padx=(0, 4)) + combo = ttk.Combobox( + group, + textvariable=self._flow_display, + values=_FLOW_DISPLAYS, + state="readonly", + width=22, + font=get_font(), + ) + combo.pack(side=tk.LEFT) + self._flow_density = tk.StringVar(value="4") + group = controls.add_group(requested_width=150) + ttk.Label(group, text="Arrow density:", style="Small.TLabel").pack( + side=tk.LEFT, padx=(0, 4) + ) + density_combo = ttk.Combobox( + group, + textvariable=self._flow_density, + values=("2", "3", "4", "5", "6"), + state="readonly", + width=5, + font=get_font(), + ) + density_combo.pack(side=tk.LEFT) + combo.bind("<>", lambda _event: self._update_flow()) + density_combo.bind("<>", lambda _event: self._update_flow()) + frame = ttk.Frame(parent) + frame.pack(fill=tk.BOTH, expand=True) + self._flow_frame = frame + self._update_flow() + + def _update_flow(self) -> None: + reset_embedded_animation(self._flow_frame, self._canvas("flow")) + display = self._flow_display.get() + density = int(self._flow_density.get()) + range_key = (display, density) + if range_key not in self._flow_range_cache: + self._flow_range_cache[range_key] = _flow_magnitude_range( + self._result, + display, + density, + ) + magnitude_min, magnitude_max = self._flow_range_cache[range_key] + payload = _FlowPayload( + self._result, + display, + density, + magnitude_min, + magnitude_max, + ) + canvas = embed_animation_plot_in_tk( + _create_flow_figure(payload), + self._flow_frame, + on_export_mp4=lambda duration: self._export_flow(payload, duration), + ) + self._set_canvas("flow", canvas) + + def _build_streamlines_tab(self, parent: ttk.Frame) -> None: + controls = make_view_controls(parent) + density_var = tk.StringVar(value="4") + group = controls.add_group(requested_width=205) + ttk.Label(group, text="Seed density:", style="Small.TLabel").pack(side=tk.LEFT, padx=(0, 4)) + density_combo = ttk.Combobox( + group, + textvariable=density_var, + values=("2", "3", "4", "5", "6"), + state="readonly", + width=5, + font=get_font(), + ) + density_combo.pack(side=tk.LEFT) + self._stream_density = density_var + self._stream_build_button = ttk.Button( + controls.add_group(requested_width=240), + text="Build streamline animation", + command=lambda: self._request_streamline_build(int(density_var.get()), ask=False), + ) + self._stream_build_button.pack(side=tk.LEFT) + self._stream_status = ttk.Label(parent, text="", style="Small.TLabel") + self._stream_status.pack(fill=tk.X, padx=8, pady=(0, 4)) + target = ttk.Frame(parent) + target.pack(fill=tk.BOTH, expand=True) + self._stream_frame = target + + def density_changed(_event: object) -> None: + self._request_streamline_build(int(density_var.get()), ask=True) + + density_combo.bind("<>", density_changed) + self._request_streamline_build(4, ask=True) + + def _request_streamline_build(self, density: int, *, ask: bool) -> None: + if density in self._streamline_cache: + self._show_streamline_cache(self._streamline_cache[density]) + return + if ask and not messagebox.askyesno( + "Build streamline animation?", + "DifferentialLab needs to trace the streamline set for all saved frames.\n\n" + "This may take some time. The result will be cached for this Results window.\n\n" + "Build it now?", + parent=self.win, + ): + self._stream_status.configure( + text="Streamline animation not built. Use Build streamline animation when ready." + ) + return + self._build_streamline_cache(density) + + def _build_streamline_cache(self, density: int) -> None: + self._stream_build_button.configure(state=tk.DISABLED) + self._stream_status.configure(text="Tracing streamlines for all saved frames…") + self.win.configure(cursor="watch") + self.win.update_idletasks() + try: + cache = build_streamline_cache(self._result, density) + self._streamline_cache[density] = cache + self._show_streamline_cache(cache) + except Exception as exc: + logger.error("Streamline animation cache failed: %s", exc, exc_info=True) + self._stream_status.configure(text="Streamline animation could not be built.") + messagebox.showerror("Streamline animation", str(exc), parent=self.win) + finally: + self.win.configure(cursor="") + self._stream_build_button.configure(state=tk.NORMAL) + + def _show_streamline_cache(self, cache: StreamlineAnimationCache) -> None: + reset_embedded_animation(self._stream_frame, self._canvas("streamlines")) + self._stream_status.configure( + text=f"Cached {cache.frame_count} frames at seed density {cache.seed_density}." + ) + payload = _StreamlinePayload(self._result, cache) + canvas = embed_animation_plot_in_tk( + _create_streamline_figure(payload), + self._stream_frame, + on_export_mp4=lambda duration: self._export_streamlines(payload, duration), + ) + self._set_canvas("streamlines", canvas) + + def _build_slices_tab(self, parent: ttk.Frame) -> None: + controls = make_view_controls(parent) + self._slice_plane = tk.StringVar(value="XY") + self._slice_field = tk.StringVar(value="Vorticity magnitude") + self._slice_index = tk.StringVar(value=str(slice_center_index(self._result, "XY"))) + group = controls.add_group(requested_width=150) + ttk.Label(group, text="Plane:", style="Small.TLabel").pack(side=tk.LEFT, padx=(0, 4)) + plane_combo = ttk.Combobox( + group, + textvariable=self._slice_plane, + values=("XY", "XZ", "YZ"), + state="readonly", + width=5, + font=get_font(), + ) + plane_combo.pack(side=tk.LEFT) + group = controls.add_group(requested_width=210) + ttk.Label(group, text="Coordinate/index:", style="Small.TLabel").pack( + side=tk.LEFT, padx=(0, 4) + ) + index_combo = ttk.Combobox( + group, + textvariable=self._slice_index, + values=tuple(str(i) for i in range(slice_index_count(self._result, "XY"))), + state="readonly", + width=6, + font=get_font(), + ) + index_combo.pack(side=tk.LEFT) + self._slice_index_combo = index_combo + group = controls.add_group(requested_width=250) + ttk.Label(group, text="Field:", style="Small.TLabel").pack(side=tk.LEFT, padx=(0, 4)) + field_combo = ttk.Combobox( + group, + textvariable=self._slice_field, + values=_SLICE_FIELDS, + state="readonly", + width=18, + font=get_font(), + ) + field_combo.pack(side=tk.LEFT) + self._slice_coordinate = ttk.Label(parent, text="", style="Small.TLabel") + self._slice_coordinate.pack(fill=tk.X, padx=8, pady=(0, 4)) + target = ttk.Frame(parent) + target.pack(fill=tk.BOTH, expand=True) + self._slice_frame = target + + def redraw(_event: object | None = None) -> None: + plane = self._slice_plane.get() + index = int(self._slice_index.get()) + field = self._slice_field.get() + key = (plane, index, field) + payload = self._slice_cache.get(key) + if payload is None: + payload = prepare_slice_history(self._result, plane=plane, index=index, field=field) + self._slice_cache[key] = payload + self._slice_coordinate.configure( + text=f"Physical slice coordinate: {payload.coordinate_label}" + ) + reset_embedded_animation(target, self._canvas("slices")) + figure = create_image_animation_plot( + payload.t, + payload.frames, + title=f"{payload.field} ({payload.coordinate_label})", + xlabel=plane[0].lower(), + ylabel=plane[1].lower(), + x_coordinates=payload.x_axis, + y_coordinates=payload.y_axis, + symmetric_color_range=payload.field in {"u", "v", "w"}, + ) + self._set_canvas("slices", embed_animation_plot_in_tk(figure, target)) + + def plane_changed(_event: object) -> None: + plane = self._slice_plane.get() + index_combo.configure( + values=tuple(str(i) for i in range(slice_index_count(self._result, plane))) + ) + self._slice_index.set(str(slice_center_index(self._result, plane))) + redraw() + + plane_combo.bind("<>", plane_changed) + index_combo.bind("<>", redraw) + field_combo.bind("<>", redraw) + redraw() + + def _build_forces_tab(self, parent: ttk.Frame) -> None: + values = np.vstack( + [self._result.drag_coeff, self._result.lift_coeff, self._result.side_force_coeff] + ) + figure = create_solution_plot( + self._result.t, + values, + title="Aerodynamic force coefficients", + xlabel="t", + ylabel="coefficient", + selected_derivatives=[0, 1, 2], + labels=["Cd", "Cl", "Cs"], + ) + self._set_canvas("forces", embed_plot_in_tk(figure, parent)) + + def _build_diagnostics_tab(self, parent: ttk.Frame) -> None: + import matplotlib.pyplot as plt + + figure, axes = plt.subplots(2, 1, sharex=True) + axes[0].plot(self._result.t, self._result.divergence_l2, label="divergence L2") + axes[0].set_ylabel("L2") + axes[0].legend() + axes[1].plot(self._result.t, self._result.max_speed, label="maximum speed") + axes[1].set_xlabel("t") + axes[1].set_ylabel("speed") + axes[1].legend() + for axis in axes: + axis.grid(True, alpha=0.3) + axis.set_ylim(*_line_limits(axis.lines[0].get_ydata())) + figure.tight_layout() + self._set_canvas("diagnostics", embed_plot_in_tk(figure, parent)) + + def _canvas(self, name: str) -> object | None: + return self._canvases.get(name) + + def _set_canvas(self, name: str, canvas: object) -> None: + self._canvases[name] = canvas + setattr(self, f"_canvas_{name}", canvas) + + def _export_flow(self, payload: _FlowPayload, duration_seconds: float) -> None: + self._export_figure( + _create_flow_figure(payload), duration_seconds, prefix="aerodynamics_3d_flow" + ) + + def _export_streamlines(self, payload: _StreamlinePayload, duration_seconds: float) -> None: + self._export_figure( + _create_streamline_figure(payload), + duration_seconds, + prefix="aerodynamics_3d_streamlines", + ) + + def _export_figure(self, figure: Figure, duration_seconds: float, *, prefix: str) -> None: + default_path = get_output_dir() / f"{generate_output_basename(prefix=prefix)}.mp4" + filepath_str = filedialog.asksaveasfilename( + parent=self.win, + defaultextension=".mp4", + initialfile=default_path.name, + initialdir=str(default_path.parent), + filetypes=[("MP4 video", "*.mp4"), ("All files", "*.*")], + ) + if not filepath_str: + return + try: + export_animated_figure_to_mp4( + figure, Path(filepath_str), duration_seconds=duration_seconds + ) + messagebox.showinfo( + "Animation export saved", + f"Animation was saved to:\n{filepath_str}", + parent=self.win, + ) + except RuntimeError as exc: + logger.warning("MP4 export failed (ffmpeg): %s", exc) + messagebox.showerror( + "Animation was not saved", + str(exc) + "\n\nInstall ffmpeg and ensure it is in your PATH.", + parent=self.win, + ) + except Exception as exc: + logger.error("MP4 export failed: %s", exc, exc_info=True) + messagebox.showerror("Animation was not saved", str(exc), parent=self.win) diff --git a/src/complex_problems/aerodynamics_3d/solver.py b/src/complex_problems/aerodynamics_3d/solver.py new file mode 100644 index 0000000..20dbcf5 --- /dev/null +++ b/src/complex_problems/aerodynamics_3d/solver.py @@ -0,0 +1,349 @@ +"""Lightweight 3D incompressible Navier--Stokes solver on a periodic grid.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Any + +import numpy as np + +from complex_problems.aerodynamics_3d.model import ( + build_obstacle_mask, + build_periodic_domain, + ddx_periodic, + ddy_periodic, + ddz_periodic, + laplacian_periodic, + spectral_divergence_l2, + validate_obstacle_fits_domain, +) +from utils import get_logger + +logger = get_logger(__name__) + +_APPROXIMATIONS = {"nonlinear_ns", "stokes"} +_SHAPES = {"sphere", "ellipsoid", "box", "naca0012_wing"} + + +@dataclass +class Aerodynamics3DResult: + """Core saved histories and compact diagnostics for a 3D solve.""" + + x: np.ndarray + y: np.ndarray + z: np.ndarray + t: np.ndarray + u: np.ndarray + v: np.ndarray + w: np.ndarray + pressure: np.ndarray + obstacle_mask: np.ndarray + drag_coeff: np.ndarray + lift_coeff: np.ndarray + side_force_coeff: np.ndarray + divergence_l2: np.ndarray + max_speed: np.ndarray + metadata: dict[str, Any] = field(default_factory=dict) + magnitudes: dict[str, float] = field(default_factory=dict) + + +def _wave_numbers( + *, nx: int, ny: int, nz: int, dx: float, dy: float, dz: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Return broadcast Fourier wave numbers in ``(z, y, x)`` order.""" + + def frequencies(size: int, spacing: float) -> np.ndarray: + values = 2.0 * np.pi * np.fft.fftfreq(size, d=spacing) + # The Nyquist coefficient has no distinct negative partner. Removing + # it keeps real-valued projected fields Hermitian in every dimension. + if size % 2 == 0: + values[size // 2] = 0.0 + return values + + kx = frequencies(nx, dx)[None, None, :] + ky = frequencies(ny, dy)[None, :, None] + kz = frequencies(nz, dz)[:, None, None] + k2 = kx * kx + ky * ky + kz * kz + return ( + np.broadcast_to(kx, (nz, ny, nx)), + np.broadcast_to(ky, (nz, ny, nx)), + np.broadcast_to(kz, (nz, ny, nx)), + k2, + ) + + +def project_velocity_fft( + u_star: np.ndarray, + v_star: np.ndarray, + w_star: np.ndarray, + *, + dt: float, + kx: np.ndarray, + ky: np.ndarray, + kz: np.ndarray, + k2: np.ndarray, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Apply the 3D Fourier Helmholtz projection and recover projection pressure. + + The pressure convention is ``u_new = u_star - dt * grad(p)``. The zero + Fourier mode is deliberately left unchanged and has zero pressure. + """ + if dt <= 0: + raise ValueError("dt must be positive.") + u_hat = np.fft.fftn(u_star) + v_hat = np.fft.fftn(v_star) + w_hat = np.fft.fftn(w_star) + d = kx * u_hat + ky * v_hat + kz * w_hat + mask = k2 > 0.0 + u_new_hat = u_hat.copy() + v_new_hat = v_hat.copy() + w_new_hat = w_hat.copy() + u_new_hat[mask] -= kx[mask] * d[mask] / k2[mask] + v_new_hat[mask] -= ky[mask] * d[mask] / k2[mask] + w_new_hat[mask] -= kz[mask] * d[mask] / k2[mask] + pressure_hat = np.zeros_like(d, dtype=complex) + pressure_hat[mask] = -1j * d[mask] / (dt * k2[mask]) + pressure_hat[0, 0, 0] = 0.0 + return ( + np.real(np.fft.ifftn(u_new_hat)), + np.real(np.fft.ifftn(v_new_hat)), + np.real(np.fft.ifftn(w_new_hat)), + np.real(np.fft.ifftn(pressure_hat)), + ) + + +def _sample_indices(n_steps: int, sample_every: int) -> list[int]: + """Return deterministic saved-step indexes including the final step.""" + indices = list(range(0, n_steps + 1, sample_every)) + if indices[-1] != n_steps: + indices.append(n_steps) + return indices + + +def solve_aerodynamics_3d( + *, + approximation: str = "nonlinear_ns", + nx: int = 48, + ny: int = 32, + nz: int = 32, + lx: float = 4.0, + ly: float = 2.0, + lz: float = 2.0, + t_max: float = 1.0, + dt: float = 0.002, + sample_every: int = 10, + rho: float = 1.0, + nu: float = 0.01, + u_inf: float = 1.0, + penalization: float = 0.005, + obstacle_shape: str = "sphere", + obstacle_center_x: float = 1.0, + obstacle_center_y: float = 1.0, + obstacle_center_z: float = 1.0, + obstacle_diameter: float = 0.4, + obstacle_size_x: float = 0.7, + obstacle_size_y: float = 0.4, + obstacle_size_z: float = 0.4, + obstacle_chord: float = 0.8, + obstacle_span: float = 0.8, + obstacle_thickness_ratio: float = 0.12, + obstacle_attack_deg: float = 0.0, +) -> Aerodynamics3DResult: + """Solve incompressible flow in a periodic Cartesian box with penalization.""" + if approximation not in _APPROXIMATIONS: + raise ValueError(f"approximation must be one of {sorted(_APPROXIMATIONS)}") + obstacle_shape = obstacle_shape.strip().lower() + if obstacle_shape not in _SHAPES: + raise ValueError(f"obstacle_shape must be one of {sorted(_SHAPES)}") + if t_max <= 0 or dt <= 0: + raise ValueError("t_max and dt must be positive.") + if sample_every < 1: + raise ValueError("sample_every must be >= 1.") + if rho <= 0 or nu <= 0 or u_inf <= 0 or penalization <= 0: + raise ValueError("rho, nu, u_inf, and penalization must be positive.") + + validate_obstacle_fits_domain( + shape=obstacle_shape, + center_x=obstacle_center_x, + center_y=obstacle_center_y, + center_z=obstacle_center_z, + lx=lx, + ly=ly, + lz=lz, + diameter=obstacle_diameter, + size_x=obstacle_size_x, + size_y=obstacle_size_y, + size_z=obstacle_size_z, + chord=obstacle_chord, + span=obstacle_span, + thickness_ratio=obstacle_thickness_ratio, + attack_deg=obstacle_attack_deg, + ) + + x, y, z, X, Y, Z, dx, dy, dz = build_periodic_domain(nx=nx, ny=ny, nz=nz, lx=lx, ly=ly, lz=lz) + obstacle_mask, reference_area, reference_length = build_obstacle_mask( + shape=obstacle_shape, + X=X, + Y=Y, + Z=Z, + center_x=obstacle_center_x, + center_y=obstacle_center_y, + center_z=obstacle_center_z, + diameter=obstacle_diameter, + size_x=obstacle_size_x, + size_y=obstacle_size_y, + size_z=obstacle_size_z, + chord=obstacle_chord, + span=obstacle_span, + thickness_ratio=obstacle_thickness_ratio, + attack_deg=obstacle_attack_deg, + ) + if not np.any(obstacle_mask): + raise ValueError("Obstacle mask is empty; increase its size or move it into the domain.") + + n_steps = int(np.ceil(t_max / dt)) + t_end = n_steps * dt + sample_indices = _sample_indices(n_steps, sample_every) + n_samples = len(sample_indices) + shape = (nz, ny, nx) + u = np.full(shape, u_inf, dtype=float) + v = np.zeros(shape, dtype=float) + w = np.zeros(shape, dtype=float) + pressure = np.zeros(shape, dtype=float) + u[obstacle_mask] = 0.0 + v[obstacle_mask] = 0.0 + w[obstacle_mask] = 0.0 + mask_f = obstacle_mask.astype(float) + fluid_mask = ~obstacle_mask + kx, ky, kz, k2 = _wave_numbers(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz) + u, v, w, pressure = project_velocity_fft(u, v, w, dt=dt, kx=kx, ky=ky, kz=kz, k2=k2) + + u_hist = np.empty((n_samples, *shape), dtype=float) + v_hist = np.empty_like(u_hist) + w_hist = np.empty_like(u_hist) + p_hist = np.empty_like(u_hist) + t_hist = np.empty(n_samples, dtype=float) + cd_hist = np.empty(n_samples, dtype=float) + cl_hist = np.empty(n_samples, dtype=float) + cs_hist = np.empty(n_samples, dtype=float) + div_l2_hist = np.empty(n_samples, dtype=float) + max_speed_hist = np.empty(n_samples, dtype=float) + q_ref = 0.5 * rho * u_inf * u_inf + cell_volume = dx * dy * dz + tau_drive = 0.2 # Matches the 2D plugin's mean-flow relaxation timescale. + sample_pos = 0 + + def save_sample(step: int) -> None: + nonlocal sample_pos + u_hist[sample_pos] = u + v_hist[sample_pos] = v + w_hist[sample_pos] = w + p_hist[sample_pos] = pressure + t_hist[sample_pos] = step * dt + reaction = rho * cell_volume / penalization * mask_f + cd_hist[sample_pos] = float(np.sum(reaction * u) / (q_ref * reference_area)) + cl_hist[sample_pos] = float(np.sum(reaction * v) / (q_ref * reference_area)) + cs_hist[sample_pos] = float(np.sum(reaction * w) / (q_ref * reference_area)) + div_l2_hist[sample_pos] = spectral_divergence_l2(u, v, w, kx, ky, kz) + max_speed_hist[sample_pos] = float(np.max(np.sqrt(u * u + v * v + w * w))) + sample_pos += 1 + + save_sample(0) + logger.info( + "Solving aerodynamics_3d: approx=%s grid=%dx%dx%d t=[0,%g] dt=%g shape=%s", + approximation, + nx, + ny, + nz, + t_end, + dt, + obstacle_shape, + ) + + for step in range(1, n_steps + 1): + ux, uy, uz = ddx_periodic(u, dx), ddy_periodic(u, dy), ddz_periodic(u, dz) + vx, vy, vz = ddx_periodic(v, dx), ddy_periodic(v, dy), ddz_periodic(v, dz) + wx, wy, wz = ddx_periodic(w, dx), ddy_periodic(w, dy), ddz_periodic(w, dz) + lap_u = laplacian_periodic(u, dx, dy, dz) + lap_v = laplacian_periodic(v, dx, dy, dz) + lap_w = laplacian_periodic(w, dx, dy, dz) + if approximation == "stokes": + conv_u = conv_v = conv_w = 0.0 + else: + conv_u = u * ux + v * uy + w * uz + conv_v = u * vx + v * vy + w * vz + conv_w = u * wx + v * wy + w * wz + + mean_u = float(np.mean(u[fluid_mask])) if np.any(fluid_mask) else float(np.mean(u)) + drive_x = (u_inf - mean_u) / tau_drive + pen_u = -mask_f * u / penalization + pen_v = -mask_f * v / penalization + pen_w = -mask_f * w / penalization + u_star = u + dt * (-conv_u + nu * lap_u + pen_u + drive_x) + v_star = v + dt * (-conv_v + nu * lap_v + pen_v) + w_star = w + dt * (-conv_w + nu * lap_w + pen_w) + u, v, w, pressure = project_velocity_fft( + u_star, v_star, w_star, dt=dt, kx=kx, ky=ky, kz=kz, k2=k2 + ) + if sample_pos < n_samples and step == sample_indices[sample_pos]: + save_sample(step) + + tail = max(3, len(cd_hist) // 4) + magnitudes = { + "reynolds": float(u_inf * reference_length / nu), + "mean_cd_tail": float(np.mean(cd_hist[-tail:])), + "rms_cl": float(np.sqrt(np.mean(cl_hist * cl_hist))), + "rms_cs": float(np.sqrt(np.mean(cs_hist * cs_hist))), + "max_divergence_l2": float(np.max(div_l2_hist)), + "max_speed": float(np.max(max_speed_hist)), + } + metadata = { + "approximation": approximation, + "nx": int(nx), + "ny": int(ny), + "nz": int(nz), + "lx": float(lx), + "ly": float(ly), + "lz": float(lz), + "dt": float(dt), + "t_max": float(t_end), + "sample_every": int(sample_every), + "rho": float(rho), + "nu": float(nu), + "u_inf": float(u_inf), + "penalization": float(penalization), + "obstacle_shape": obstacle_shape, + "obstacle_center_x": float(obstacle_center_x), + "obstacle_center_y": float(obstacle_center_y), + "obstacle_center_z": float(obstacle_center_z), + "obstacle_diameter": float(obstacle_diameter), + "obstacle_size_x": float(obstacle_size_x), + "obstacle_size_y": float(obstacle_size_y), + "obstacle_size_z": float(obstacle_size_z), + "obstacle_chord": float(obstacle_chord), + "obstacle_span": float(obstacle_span), + "obstacle_thickness_ratio": float(obstacle_thickness_ratio), + "obstacle_attack_deg": float(obstacle_attack_deg), + "reference_area": float(reference_area), + "reference_length": float(reference_length), + "stored_fields": ("u", "v", "w", "pressure"), + "derived_fields": ("speed", "vorticity_vector", "vorticity_magnitude"), + } + return Aerodynamics3DResult( + x=x, + y=y, + z=z, + t=t_hist, + u=u_hist, + v=v_hist, + w=w_hist, + pressure=p_hist, + obstacle_mask=obstacle_mask, + drag_coeff=cd_hist, + lift_coeff=cl_hist, + side_force_coeff=cs_hist, + divergence_l2=div_l2_hist, + max_speed=max_speed_hist, + metadata=metadata, + magnitudes=magnitudes, + ) diff --git a/src/complex_problems/aerodynamics_3d/ui.py b/src/complex_problems/aerodynamics_3d/ui.py new file mode 100644 index 0000000..795c8a1 --- /dev/null +++ b/src/complex_problems/aerodynamics_3d/ui.py @@ -0,0 +1,259 @@ +"""Configuration dialog for the 3D aerodynamics problem.""" + +from __future__ import annotations + +import tkinter as tk +from tkinter import ttk +from typing import cast + +from complex_problems.aerodynamics_3d.model import validate_obstacle_fits_domain +from complex_problems.aerodynamics_3d.solver import solve_aerodynamics_3d +from complex_problems.common import ( + add_how_to_config_section, + parse_float, + parse_positive_float, + parse_positive_int, +) +from complex_problems.common.dialog_ui import ( + AdvancedDialogShell, + AdvancedDialogSize, + make_labeled_combo, + make_labeled_entry, + make_labeled_spinbox, + run_solver_dialog, +) +from config import get_env_from_schema +from frontend.performance_guard import ( + assess_aerodynamics_3d_request, + confirm_performance_advisory, +) +from frontend.ui_dialogs.tooltip import ToolTip +from frontend.window_utils import make_modal + +_APPROXIMATIONS = ("nonlinear_ns", "stokes") +_SHAPES = ("sphere", "ellipsoid", "box", "naca0012_wing") + + +class Aerodynamics3DDialog: + """Scrollable, progressive-disclosure configuration dialog.""" + + def __init__(self, parent: tk.Tk | tk.Toplevel) -> None: + self.parent = parent + self.win = tk.Toplevel(parent) + self.win.title("Aerodynamics 3D") + self._build_ui() + self._shell.finish(AdvancedDialogSize(980, 800, 700, 540)) + make_modal(self.win, parent) + self._initial_focus.focus_set() + + def _build_ui(self) -> None: + pad = int(get_env_from_schema("UI_PADDING")) + self._shell = AdvancedDialogShell( + self.win, + title="Aerodynamics 3D", + description=( + "Solve incompressible 3D structured-grid flow in a periodic Cartesian domain " + "with an immersed/penalized obstacle. This is a lightweight educational and " + "scientific solver, not an industrial CFD replacement." + ), + pad=pad, + ) + body = self._shell.body + add_how_to_config_section( + body, self._shell.scroll, problem_id="aerodynamics_3d", pad=pad, wraplength=820 + ) + + model = ttk.LabelFrame(body, text="Model", padding=pad) + model.pack(fill=tk.X, pady=(0, pad)) + row = ttk.Frame(model) + row.pack(fill=tk.X, pady=pad // 2) + self._approx_var = tk.StringVar(value="nonlinear_ns") + self._shape_var = tk.StringVar(value="sphere") + self._initial_focus = make_labeled_combo( + row, "Approximation", self._approx_var, _APPROXIMATIONS, width=15 + ) + row = ttk.Frame(model) + row.pack(fill=tk.X, pady=pad // 2) + shape_combo = make_labeled_combo(row, "Obstacle shape", self._shape_var, _SHAPES, width=18) + shape_combo.bind("<>", lambda _event: self._refresh_geometry()) + + domain = ttk.LabelFrame(body, text="Domain and grid", padding=pad) + domain.pack(fill=tk.X, pady=(0, pad)) + row = ttk.Frame(domain) + row.pack(fill=tk.X, pady=pad // 2) + self._nx_var, self._ny_var, self._nz_var = ( + tk.StringVar(value=value) for value in ("48", "32", "32") + ) + make_labeled_spinbox(row, "Nₓ", self._nx_var, from_=4, to=2048, width=7) + make_labeled_spinbox(row, "Nᵧ", self._ny_var, from_=4, to=2048, width=7) + make_labeled_spinbox(row, "N_z", self._nz_var, from_=4, to=2048, width=7) + row = ttk.Frame(domain) + row.pack(fill=tk.X, pady=pad // 2) + self._lx_var, self._ly_var, self._lz_var = ( + tk.StringVar(value=value) for value in ("4.0", "2.0", "2.0") + ) + make_labeled_entry(row, "Lₓ", self._lx_var, width=7) + make_labeled_entry(row, "Lᵧ", self._ly_var, width=7) + make_labeled_entry(row, "L_z", self._lz_var, width=7) + row = ttk.Frame(domain) + row.pack(fill=tk.X, pady=pad // 2) + self._t_max_var, self._dt_var, self._sample_every_var = ( + tk.StringVar(value=value) for value in ("1.0", "0.002", "10") + ) + make_labeled_entry(row, "tₘₐₓ", self._t_max_var, width=8) + make_labeled_entry(row, "Δt", self._dt_var, width=8) + make_labeled_entry(row, "Sample every", self._sample_every_var, width=10) + ToolTip(row, "Lower Δt and/or Sample every values store more volumetric animation frames.") + + physics = ttk.LabelFrame(body, text="Flow properties", padding=pad) + physics.pack(fill=tk.X, pady=(0, pad)) + row = ttk.Frame(physics) + row.pack(fill=tk.X, pady=pad // 2) + self._rho_var, self._nu_var, self._u_inf_var, self._penal_var = ( + tk.StringVar(value=value) for value in ("1.0", "0.01", "1.0", "0.005") + ) + make_labeled_entry(row, "ρ", self._rho_var, width=8) + make_labeled_entry(row, "ν", self._nu_var, width=8) + make_labeled_entry(row, "U∞", self._u_inf_var, width=8) + make_labeled_entry(row, "Penalization", self._penal_var, width=10) + + geometry = ttk.LabelFrame(body, text="Obstacle geometry", padding=pad) + geometry.pack(fill=tk.X) + row = ttk.Frame(geometry) + row.pack(fill=tk.X, pady=pad // 2) + self._center_x_var, self._center_y_var, self._center_z_var = ( + tk.StringVar(value=value) for value in ("1.0", "1.0", "1.0") + ) + make_labeled_entry(row, "Center x", self._center_x_var, width=8) + make_labeled_entry(row, "Center y", self._center_y_var, width=8) + make_labeled_entry(row, "Center z", self._center_z_var, width=8) + self._geometry_body = ttk.Frame(geometry) + self._geometry_body.pack(fill=tk.X) + self._geometry_rows: dict[str, ttk.Frame] = {} + self._make_shape_rows(pad) + self._refresh_geometry() + self._shell.add_footer_button("Close", self.win.destroy) + self._shell.add_footer_button("Solve", self._on_solve, primary=True) + + def _make_shape_rows(self, pad: int) -> None: + for shape in _SHAPES: + frame = ttk.Frame(self._geometry_body) + self._geometry_rows[shape] = frame + if shape == "sphere": + self._diameter_var = tk.StringVar(value="0.4") + make_labeled_entry(frame, "Diameter", self._diameter_var, width=8) + elif shape in ("ellipsoid", "box"): + self._size_x_var = tk.StringVar(value="0.7") + self._size_y_var = tk.StringVar(value="0.4") + self._size_z_var = tk.StringVar(value="0.4") + self._attack_var = tk.StringVar(value="0.0") + make_labeled_entry(frame, "Size x", self._size_x_var, width=8) + make_labeled_entry(frame, "Size y", self._size_y_var, width=8) + make_labeled_entry(frame, "Size z", self._size_z_var, width=8) + make_labeled_entry(frame, "Attack (deg)", self._attack_var, width=10) + else: + self._chord_var = tk.StringVar(value="0.8") + self._span_var = tk.StringVar(value="0.8") + self._thickness_var = tk.StringVar(value="0.12") + self._wing_attack_var = tk.StringVar(value="0.0") + make_labeled_entry(frame, "Chord", self._chord_var, width=8) + make_labeled_entry(frame, "Span", self._span_var, width=8) + make_labeled_entry(frame, "Thickness ratio", self._thickness_var, width=12) + make_labeled_entry(frame, "Attack (deg)", self._wing_attack_var, width=10) + ToolTip(frame, "The finite-span NACA 0012 wing has no sweep, taper, or twist.") + + def _refresh_geometry(self) -> None: + for frame in self._geometry_rows.values(): + frame.pack_forget() + self._geometry_rows[self._shape_var.get()].pack(fill=tk.X, pady=4) + self._shell.refresh() + + def _collect_inputs(self) -> dict[str, object]: + nx = parse_positive_int(self._nx_var.get(), name="Nₓ", min_value=4) + ny = parse_positive_int(self._ny_var.get(), name="Nᵧ", min_value=4) + nz = parse_positive_int(self._nz_var.get(), name="N_z", min_value=4) + values = { + "lx": parse_positive_float(self._lx_var.get(), name="Lₓ"), + "ly": parse_positive_float(self._ly_var.get(), name="Lᵧ"), + "lz": parse_positive_float(self._lz_var.get(), name="L_z"), + "t_max": parse_positive_float(self._t_max_var.get(), name="tₘₐₓ"), + "dt": parse_positive_float(self._dt_var.get(), name="Δt"), + "sample_every": parse_positive_int(self._sample_every_var.get(), name="Sample every"), + "rho": parse_positive_float(self._rho_var.get(), name="ρ"), + "nu": parse_positive_float(self._nu_var.get(), name="ν"), + "u_inf": parse_positive_float(self._u_inf_var.get(), name="U∞"), + "penalization": parse_positive_float(self._penal_var.get(), name="Penalization"), + "obstacle_center_x": parse_float(self._center_x_var.get(), name="Center x"), + "obstacle_center_y": parse_float(self._center_y_var.get(), name="Center y"), + "obstacle_center_z": parse_float(self._center_z_var.get(), name="Center z"), + } + shape = self._shape_var.get() + values.update( + { + "approximation": self._approx_var.get(), + "obstacle_shape": shape, + "nx": nx, + "ny": ny, + "nz": nz, + } + ) + if shape == "sphere": + values["obstacle_diameter"] = parse_positive_float( + self._diameter_var.get(), name="Diameter" + ) + values["obstacle_attack_deg"] = 0.0 + elif shape in ("ellipsoid", "box"): + values["obstacle_size_x"] = parse_positive_float(self._size_x_var.get(), name="Size x") + values["obstacle_size_y"] = parse_positive_float(self._size_y_var.get(), name="Size y") + values["obstacle_size_z"] = parse_positive_float(self._size_z_var.get(), name="Size z") + values["obstacle_attack_deg"] = parse_float(self._attack_var.get(), name="Attack") + else: + values["obstacle_chord"] = parse_positive_float(self._chord_var.get(), name="Chord") + values["obstacle_span"] = parse_positive_float(self._span_var.get(), name="Span") + values["obstacle_thickness_ratio"] = parse_positive_float( + self._thickness_var.get(), name="Thickness ratio" + ) + values["obstacle_attack_deg"] = parse_float(self._wing_attack_var.get(), name="Attack") + validate_obstacle_fits_domain( + shape=cast(str, values["obstacle_shape"]), + center_x=cast(float, values["obstacle_center_x"]), + center_y=cast(float, values["obstacle_center_y"]), + center_z=cast(float, values["obstacle_center_z"]), + lx=cast(float, values["lx"]), + ly=cast(float, values["ly"]), + lz=cast(float, values["lz"]), + diameter=cast(float, values.get("obstacle_diameter", 0.4)), + size_x=cast(float, values.get("obstacle_size_x", 0.7)), + size_y=cast(float, values.get("obstacle_size_y", 0.4)), + size_z=cast(float, values.get("obstacle_size_z", 0.4)), + chord=cast(float, values.get("obstacle_chord", 0.8)), + span=cast(float, values.get("obstacle_span", 0.8)), + thickness_ratio=cast(float, values.get("obstacle_thickness_ratio", 0.12)), + attack_deg=cast(float, values.get("obstacle_attack_deg", 0.0)), + ) + return values + + def _confirm_heavy_request(self, params: dict[str, object], window: tk.Toplevel) -> bool: + advisory = assess_aerodynamics_3d_request( + nx=cast(int, params["nx"]), + ny=cast(int, params["ny"]), + nz=cast(int, params["nz"]), + t_max=cast(float, params["t_max"]), + dt=cast(float, params["dt"]), + sample_every=cast(int, params["sample_every"]), + ) + return confirm_performance_advisory(window, advisory) + + def _on_solve(self) -> None: + from complex_problems.aerodynamics_3d.result_dialog import Aerodynamics3DResultDialog + + run_solver_dialog( + parent=self.parent, + window=self.win, + collect_inputs=self._collect_inputs, + solver=solve_aerodynamics_3d, + message="Solving 3D aerodynamics...", + result_parent=self.parent, + result_dialog_factory=Aerodynamics3DResultDialog, + confirm_run=self._confirm_heavy_request, + ) diff --git a/src/complex_problems/problem_docs.py b/src/complex_problems/problem_docs.py index 236185e..2873fce 100644 --- a/src/complex_problems/problem_docs.py +++ b/src/complex_problems/problem_docs.py @@ -149,6 +149,29 @@ class ProblemDoc: "Streamlines and centerline diagnostic profiles.", ), ), + "aerodynamics_3d": ProblemDoc( + problem_type="Incompressible 3D structured-grid flow around obstacles", + extended_description=( + "Lightweight incompressible flow solver on a periodic Cartesian domain with " + "immersed/Brinkman-style penalization. It is an educational/scientific " + "structured-grid solver, not an industrial CFD replacement." + ), + equation_summary="∂u/∂t + (u·∇)u = −∇p + ν∇²u + fₚₑₙ + f_drive, ∇·u = 0", + config_options_summary=( + "Choose nonlinear Navier–Stokes or the Stokes approximation.", + "Set a periodic Cartesian Nx × Ny × Nz domain, time step, and saved-frame cadence.", + "Configure density, viscosity, mean-flow speed, and penalization strength.", + "Select a sphere, ellipsoid, box, or finite-span NACA 0012 wing with " + "shape-specific dimensions.", + ), + visualizations_summary=( + "Animated 3D velocity or vorticity vectors with an obstacle boundary.", + "RK4-traced 3D streamlines from a selected saved frame.", + "Responsive XY, XZ, and YZ slices for speed, vorticity, pressure, or " + "velocity components.", + "Cd/Cl/Cs force histories and factual divergence/speed diagnostics.", + ), + ), "pipe_flow": ProblemDoc( problem_type="1D steady/transient hydraulic modeling", extended_description=( diff --git a/src/complex_problems/problem_registry.py b/src/complex_problems/problem_registry.py index 81fdd47..8f6ecca 100644 --- a/src/complex_problems/problem_registry.py +++ b/src/complex_problems/problem_registry.py @@ -96,6 +96,7 @@ def __len__(self) -> int: ProblemRegistration(module_path="complex_problems.schrodinger_td.problem"), ProblemRegistration(module_path="complex_problems.antenna_radiation.problem"), ProblemRegistration(module_path="complex_problems.aerodynamics_2d.problem"), + ProblemRegistration(module_path="complex_problems.aerodynamics_3d.problem"), ProblemRegistration(module_path="complex_problems.pipe_flow.problem"), ProblemRegistration(module_path="complex_problems.gravitational_n_body.problem"), ProblemRegistration(module_path="complex_problems.fput_experiment.problem"), diff --git a/src/frontend/performance_guard.py b/src/frontend/performance_guard.py index 88cb5a8..7db86d3 100644 --- a/src/frontend/performance_guard.py +++ b/src/frontend/performance_guard.py @@ -316,6 +316,31 @@ def assess_aerodynamics_request( ) +def assess_aerodynamics_3d_request( + *, + nx: int, + ny: int, + nz: int, + t_max: float, + dt: float, + sample_every: int, +) -> PerformanceAdvisory | None: + """Assess retained 3D velocity/pressure history before GUI execution.""" + if t_max <= 0 or dt <= 0 or sample_every < 1: + return None + n_steps = int(math.ceil(t_max / dt)) + n_samples = len(range(0, n_steps + 1, sample_every)) + if n_steps % sample_every != 0: + n_samples += 1 + return assess_time_history_request( + label="3D aerodynamics", + frames=n_samples, + points_per_frame=max(1, nx * ny * nz), + array_count=4, + bytes_per_value=8, + ) + + def assess_fput_request( *, n_particles: int, diff --git a/tests/frontend/test_performance_guard.py b/tests/frontend/test_performance_guard.py index c2081e5..69cddc4 100644 --- a/tests/frontend/test_performance_guard.py +++ b/tests/frontend/test_performance_guard.py @@ -6,6 +6,7 @@ from frontend.performance_guard import ( PerformanceAdvisory, + assess_aerodynamics_3d_request, assess_membrane_request, assess_parameters_dialog_request, assess_pde_3d_request, @@ -108,6 +109,23 @@ def test_small_membrane_history_has_no_advisory() -> None: ) +def test_default_aerodynamics_3d_history_has_no_advisory() -> None: + assert ( + assess_aerodynamics_3d_request(nx=48, ny=32, nz=32, t_max=1.0, dt=0.002, sample_every=10) + is None + ) + + +def test_huge_aerodynamics_3d_history_requires_confirmation() -> None: + advisory = assess_aerodynamics_3d_request( + nx=100, ny=100, nz=100, t_max=2.0, dt=0.002, sample_every=2 + ) + assert advisory is not None + assert advisory.severity == "confirm" + assert advisory.title == "Large 3D aerodynamics history" + assert "4 array(s)" in advisory.message + + def test_confirm_performance_advisory_uses_warning_and_confirmation_dialogs() -> None: warn = PerformanceAdvisory( severity="warn", diff --git a/tests/test_advanced_problem_contract.py b/tests/test_advanced_problem_contract.py index 9122ab1..892ffc5 100644 --- a/tests/test_advanced_problem_contract.py +++ b/tests/test_advanced_problem_contract.py @@ -14,6 +14,7 @@ import pytest from complex_problems.aerodynamics_2d.solver import solve_aerodynamics_2d +from complex_problems.aerodynamics_3d.solver import solve_aerodynamics_3d from complex_problems.antenna_radiation.solver import solve_antenna_radiation from complex_problems.coupled_oscillators.solver import solve_coupled_oscillators from complex_problems.fput_experiment.solver import solve_fput @@ -136,6 +137,27 @@ def _solve_membrane() -> object: obstacle_size_y=0.2, ), ), + _SmokeCase( + "aerodynamics_3d", + lambda: solve_aerodynamics_3d( + approximation="stokes", + nx=8, + ny=8, + nz=8, + lx=2.0, + ly=2.0, + lz=2.0, + t_max=0.004, + dt=0.002, + sample_every=2, + rho=1.0, + nu=0.03, + u_inf=0.8, + penalization=0.01, + obstacle_shape="sphere", + obstacle_diameter=0.5, + ), + ), _SmokeCase( "pipe_flow", lambda: solve_pipe_flow( @@ -208,6 +230,7 @@ def loaded_plugin_modules() -> set[str]: "schrodinger_td", "antenna_radiation", "aerodynamics_2d", + "aerodynamics_3d", "pipe_flow", "gravitational_n_body", "fput_experiment", @@ -258,10 +281,10 @@ def test_registered_plugin_contract_is_lazy_and_valid_in_isolated_interpreter() def test_registered_plugin_descriptors_are_unique() -> None: - """All nine registered plugin descriptors have distinct identifiers.""" + """All ten registered plugin descriptors have distinct identifiers.""" descriptors = ProblemRegistry(_REGISTRATIONS).get_descriptors() - assert len(descriptors) == len(_REGISTRATIONS) == 9 + assert len(descriptors) == len(_REGISTRATIONS) == 10 assert set(descriptors) == {case.problem_id for case in _SMOKE_CASES} diff --git a/tests/test_aerodynamics_3d.py b/tests/test_aerodynamics_3d.py new file mode 100644 index 0000000..cce062c --- /dev/null +++ b/tests/test_aerodynamics_3d.py @@ -0,0 +1,195 @@ +"""Focused numerical tests for the 3D aerodynamics Advanced Problem.""" + +from __future__ import annotations + +import numpy as np + +from complex_problems.aerodynamics_3d import model +from complex_problems.aerodynamics_3d.solver import ( + _wave_numbers, + project_velocity_fft, + solve_aerodynamics_3d, +) + + +def _grid() -> tuple[np.ndarray, np.ndarray, np.ndarray]: + x = np.linspace(0.0, 4.0, 16, endpoint=False) + y = np.linspace(0.0, 2.0, 12, endpoint=False) + z = np.linspace(0.0, 2.0, 10, endpoint=False) + zz, yy, xx = np.meshgrid(z, y, x, indexing="ij") + return xx, yy, zz + + +def test_supported_obstacle_masks_and_reference_scales_are_nonempty() -> None: + X, Y, Z = _grid() + common = dict(X=X, Y=Y, Z=Z, center_x=1.5, center_y=1.0, center_z=1.0) + cases = ( + {"shape": "sphere", "diameter": 0.5}, + {"shape": "ellipsoid", "size_x": 0.8, "size_y": 0.5, "size_z": 0.4}, + {"shape": "box", "size_x": 0.8, "size_y": 0.5, "size_z": 0.4}, + {"shape": "naca0012_wing", "chord": 0.8, "span": 0.8}, + ) + for case in cases: + mask, area, length = model.build_obstacle_mask(**common, **case) + assert mask.shape == X.shape + assert np.any(mask) + assert np.isfinite(area) and area > 0 + assert np.isfinite(length) and length > 0 + + +def test_attack_angle_changes_a_non_axisymmetric_mask() -> None: + X, Y, Z = _grid() + kwargs = dict( + shape="box", + X=X, + Y=Y, + Z=Z, + center_x=1.5, + center_y=1.0, + center_z=1.0, + size_x=1.0, + size_y=0.4, + size_z=0.4, + ) + mask_zero, _, _ = model.build_obstacle_mask(**kwargs, attack_deg=0.0) + mask_rotated, _, _ = model.build_obstacle_mask(**kwargs, attack_deg=27.0) + assert np.any(mask_zero) + assert np.any(mask_rotated) + assert not np.array_equal(mask_zero, mask_rotated) + + +def test_obstacle_boundary_validation_accepts_center_and_rejects_crossing() -> None: + model.validate_obstacle_fits_domain( + shape="sphere", + center_x=2.0, + center_y=1.0, + center_z=1.0, + lx=4.0, + ly=2.0, + lz=2.0, + diameter=0.4, + ) + with np.testing.assert_raises_regex(ValueError, "fully inside the periodic box"): + model.validate_obstacle_fits_domain( + shape="sphere", + center_x=0.1, + center_y=1.0, + center_z=1.0, + lx=4.0, + ly=2.0, + lz=2.0, + diameter=0.4, + ) + + +def test_rotated_box_boundary_validation_uses_world_aligned_extent() -> None: + with np.testing.assert_raises_regex(ValueError, "fully inside the periodic box"): + model.validate_obstacle_fits_domain( + shape="box", + center_x=0.45, + center_y=1.0, + center_z=1.0, + lx=4.0, + ly=2.0, + lz=2.0, + size_x=1.0, + size_y=0.4, + size_z=0.4, + attack_deg=45.0, + ) + + +def test_fft_projection_removes_nontrivial_spectral_divergence() -> None: + nx, ny, nz = 12, 10, 8 + dx, dy, dz = 4.0 / nx, 2.0 / ny, 2.0 / nz + kx, ky, kz, k2 = _wave_numbers(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz) + rng = np.random.default_rng(42) + u_star, v_star, w_star = rng.normal(size=(3, nz, ny, nx)) + before = np.sqrt( + np.mean( + np.real( + np.fft.ifftn( + 1j + * ( + kx * np.fft.fftn(u_star) + + ky * np.fft.fftn(v_star) + + kz * np.fft.fftn(w_star) + ) + ) + ) + ** 2 + ) + ) + u, v, w, pressure = project_velocity_fft( + u_star, v_star, w_star, dt=0.01, kx=kx, ky=ky, kz=kz, k2=k2 + ) + after = model.spectral_divergence_l2(u, v, w, kx, ky, kz) + assert before > 1.0e-3 + assert after < 1.0e-12 + assert np.all(np.isfinite(pressure)) + np.testing.assert_allclose(np.mean(u), np.mean(u_star)) + np.testing.assert_allclose(np.mean(v), np.mean(v_star)) + np.testing.assert_allclose(np.mean(w), np.mean(w_star)) + + +def test_short_stokes_and_nonlinear_runs_store_only_core_volumes() -> None: + for approximation in ("stokes", "nonlinear_ns"): + result = solve_aerodynamics_3d( + approximation=approximation, + nx=8, + ny=8, + nz=8, + lx=2.0, + ly=2.0, + lz=2.0, + t_max=0.006, + dt=0.002, + sample_every=2, + nu=0.03, + penalization=0.01, + obstacle_diameter=0.5, + ) + assert result.u.shape == (3, 8, 8, 8) + assert result.v.shape == result.w.shape == result.pressure.shape == result.u.shape + assert result.obstacle_mask.shape == (8, 8, 8) + assert ( + result.drag_coeff.shape + == result.lift_coeff.shape + == result.side_force_coeff.shape + == (3,) + ) + assert result.divergence_l2.shape == result.max_speed.shape == (3,) + assert all( + np.all(np.isfinite(array)) + for array in ( + result.u, + result.v, + result.w, + result.pressure, + result.drag_coeff, + result.divergence_l2, + ) + ) + assert not hasattr(result, "speed") + assert not hasattr(result, "vorticity") + + +def test_penalization_timescale_changes_solved_obstacle_state() -> None: + common = dict( + approximation="stokes", + nx=10, + ny=10, + nz=10, + lx=2.0, + ly=2.0, + lz=2.0, + t_max=0.02, + dt=0.001, + sample_every=20, + nu=0.03, + obstacle_diameter=0.5, + ) + fast = solve_aerodynamics_3d(**common, penalization=0.002) + slow = solve_aerodynamics_3d(**common, penalization=0.02) + obstacle = fast.obstacle_mask + assert np.max(np.abs(fast.u[-1][obstacle] - slow.u[-1][obstacle])) > 1.0e-5 diff --git a/tests/test_aerodynamics_3d_result_dialog_ui.py b/tests/test_aerodynamics_3d_result_dialog_ui.py new file mode 100644 index 0000000..ec8c0a0 --- /dev/null +++ b/tests/test_aerodynamics_3d_result_dialog_ui.py @@ -0,0 +1,171 @@ +"""Pure result-helper tests for 3D aerodynamics views.""" + +from __future__ import annotations + +from dataclasses import replace +from unittest.mock import patch + +import numpy as np + +import complex_problems.aerodynamics_3d.result_dialog as result_dialog +from complex_problems.aerodynamics_3d.result_dialog import ( + _flow_magnitude_range, + _scaled_arrow_lengths, + build_streamline_cache, + prepare_slice_history, + slice_center_index, + slice_field, + slice_index_count, + trace_streamline_3d, +) +from complex_problems.aerodynamics_3d.solver import Aerodynamics3DResult + + +def _uniform_result() -> Aerodynamics3DResult: + x = np.linspace(0.0, 4.0, 16, endpoint=False) + y = np.linspace(0.0, 2.0, 8, endpoint=False) + z = np.linspace(0.0, 2.0, 8, endpoint=False) + shape = (1, len(z), len(y), len(x)) + return Aerodynamics3DResult( + x=x, + y=y, + z=z, + t=np.array([0.0]), + u=np.ones(shape), + v=np.zeros(shape), + w=np.zeros(shape), + pressure=np.zeros(shape), + obstacle_mask=np.zeros(shape[1:], dtype=bool), + drag_coeff=np.zeros(1), + lift_coeff=np.zeros(1), + side_force_coeff=np.zeros(1), + divergence_l2=np.zeros(1), + max_speed=np.ones(1), + ) + + +def test_uniform_positive_x_streamline_is_straight() -> None: + result = _uniform_result() + slow = replace(result, u=result.u * 0.2) + fast = replace(result, u=result.u * 5.0) + line = trace_streamline_3d((0.2, 0.8, 0.8), result=slow, step_size=0.1, max_steps=10) + fast_line = trace_streamline_3d((0.2, 0.8, 0.8), result=fast, step_size=0.1, max_steps=10) + assert len(line) > 3 + np.testing.assert_allclose(line, fast_line, atol=1.0e-12) + np.testing.assert_allclose(line[:, 1], line[0, 1], atol=1.0e-12) + np.testing.assert_allclose(line[:, 2], line[0, 2], atol=1.0e-12) + assert np.all(np.diff(line[:, 0]) > 0) + + +def test_slice_helper_returns_requested_physical_coordinate() -> None: + result = _uniform_result() + x_axis, y_axis, values, coordinate = slice_field( + result, frame=0, plane="XY", index=3, field="Speed" + ) + np.testing.assert_array_equal(x_axis, result.x) + np.testing.assert_array_equal(y_axis, result.y) + assert values.shape == (len(result.y), len(result.x)) + assert coordinate.startswith("z =") + np.testing.assert_allclose(values, 1.0) + + +def test_slice_selectors_and_prepared_history_are_plane_specific() -> None: + result = _uniform_result() + result = replace( + result, + t=np.array([0.0, 0.5, 1.0]), + u=np.repeat(result.u, 3, axis=0), + v=np.repeat(result.v, 3, axis=0), + w=np.repeat(result.w, 3, axis=0), + pressure=np.repeat(result.pressure, 3, axis=0), + drag_coeff=np.zeros(3), + lift_coeff=np.zeros(3), + side_force_coeff=np.zeros(3), + divergence_l2=np.zeros(3), + max_speed=np.ones(3), + ) + assert slice_index_count(result, "XY") == len(result.z) + assert slice_index_count(result, "XZ") == len(result.y) + assert slice_index_count(result, "YZ") == len(result.x) + assert slice_center_index(result, "XZ") == len(result.y) // 2 + payload = prepare_slice_history( + result, + plane="XZ", + index=slice_center_index(result, "XZ"), + field="Speed", + ) + assert payload.frames.shape == (len(result.t), len(result.z), len(result.x)) + assert payload.coordinate == result.y[len(result.y) // 2] + + +def test_streamline_cache_has_one_entry_per_saved_frame() -> None: + result = _uniform_result() + result = replace( + result, + t=np.array([0.0, 0.5, 1.0]), + u=np.repeat(result.u, 3, axis=0), + v=np.repeat(result.v, 3, axis=0), + w=np.repeat(result.w, 3, axis=0), + ) + cache = build_streamline_cache(result, density=2) + assert cache.seed_density == 2 + assert cache.frame_count == len(result.t) + assert len(cache.frames) == len(result.t) + + +def test_streamline_cache_reuses_one_velocity_interpolator_per_frame() -> None: + result = _uniform_result() + result = replace( + result, + t=np.array([0.0, 0.5, 1.0]), + u=np.repeat(result.u, 3, axis=0), + v=np.repeat(result.v, 3, axis=0), + w=np.repeat(result.w, 3, axis=0), + ) + with patch.object( + result_dialog, + "_velocity_interpolator", + wraps=result_dialog._velocity_interpolator, + ) as interpolator_factory: + cache = build_streamline_cache(result, density=2) + + assert cache.frame_count == len(result.t) + assert interpolator_factory.call_count == len(result.t) + + +def test_log_scaled_arrow_lengths_use_global_minimum_and_maximum() -> None: + magnitude = np.array([0.1, 0.2, 0.5, 1.0]) + lengths = _scaled_arrow_lengths( + magnitude, + 0.1, + 1.0, + minimum=0.04, + maximum=0.16, + ) + + assert np.all(np.diff(lengths) >= 0.0) + assert lengths[0] == 0.04 + assert lengths[-1] == 0.16 + assert np.all(lengths <= 0.16) + + +def test_flow_magnitude_range_spans_all_saved_frames() -> None: + result = _uniform_result() + result = replace( + result, + t=np.array([0.0, 1.0]), + u=np.concatenate((result.u * 0.2, result.u * 5.0), axis=0), + v=np.repeat(result.v, 2, axis=0), + w=np.repeat(result.w, 2, axis=0), + pressure=np.repeat(result.pressure, 2, axis=0), + drag_coeff=np.zeros(2), + lift_coeff=np.zeros(2), + side_force_coeff=np.zeros(2), + divergence_l2=np.zeros(2), + max_speed=np.array([0.2, 5.0]), + ) + + minimum, maximum = _flow_magnitude_range(result, "Velocity vectors", 4) + + assert minimum == 0.2 + assert maximum == 5.0 diff --git a/tests/test_complex_problem_result_dialog_ui.py b/tests/test_complex_problem_result_dialog_ui.py index bbd611a..8dc3f1c 100644 --- a/tests/test_complex_problem_result_dialog_ui.py +++ b/tests/test_complex_problem_result_dialog_ui.py @@ -258,6 +258,7 @@ def test_n_body_close_cleans_every_canvas_attribute() -> None: "module_name", ( "complex_problems.aerodynamics_2d.result_dialog", + "complex_problems.aerodynamics_3d.result_dialog", "complex_problems.antenna_radiation.result_dialog", "complex_problems.coupled_oscillators.result_dialog", "complex_problems.fput_experiment.result_dialog", diff --git a/tests/test_complex_problems_registry.py b/tests/test_complex_problems_registry.py index ebe1815..5ff62be 100644 --- a/tests/test_complex_problems_registry.py +++ b/tests/test_complex_problems_registry.py @@ -17,6 +17,7 @@ def test_registry_contains_expected_plugins() -> None: assert "schrodinger_td" in descriptors assert "antenna_radiation" in descriptors assert "aerodynamics_2d" in descriptors + assert "aerodynamics_3d" in descriptors assert "pipe_flow" in descriptors assert "gravitational_n_body" in descriptors assert "fput_experiment" in descriptors @@ -68,6 +69,13 @@ def test_open_problem_dialog_dispatches_to_aerodynamics_2d() -> None: mock_open.assert_called_once_with(parent) +def test_open_problem_dialog_dispatches_to_aerodynamics_3d() -> None: + parent = object() + with patch("complex_problems.aerodynamics_3d.problem.PROBLEM.open_dialog") as mock_open: + open_problem_dialog("aerodynamics_3d", parent) # type: ignore[arg-type] + mock_open.assert_called_once_with(parent) + + def test_open_problem_dialog_dispatches_to_pipe_flow() -> None: parent = object() with patch("complex_problems.pipe_flow.problem.PROBLEM.open_dialog") as mock_open: