diff --git a/.cursorrules b/.cursorrules index 0e298f0..f98dc5f 100644 --- a/.cursorrules +++ b/.cursorrules @@ -1,73 +1,139 @@ # Cursor Rules for RBF Incompressible Navier-Stokes MATLAB Codebase ## Project Overview -This is a MATLAB implementation of a mesh-free Radial Basis Function Finite Difference (RBF-FD) method for solving incompressible Navier-Stokes equations around obstacles. The code supports multiple geometries (cylinder, ellipse, rectangle) and uses Polyharmonic Spline (PHS) basis functions with polynomial augmentation. +This is a MATLAB implementation of a mesh-free Radial Basis Function Finite Difference (RBF-FD) method for solving incompressible Navier-Stokes equations around obstacles. The code supports multiple geometries (cylinder, ellipse, rectangle, airfoil, multi-obstacle) and uses Polyharmonic Spline (PHS) basis functions with polynomial augmentation. It includes advanced features like tracer particles, Stokes dynamics, adaptive time stepping, and comprehensive stability diagnostics. ## Core Architecture Principles ### 1. Mathematical Foundation -- **RBF-FD Method**: Uses Polyharmonic Spline basis functions φ(r) = r^m with polynomial augmentation -- **Staggered Grid**: Pressure on P-grid (xy_s), velocity on V-grid (xy) -- **Fractional Step Method**: Adams-Bashforth for advection, Crank-Nicolson for diffusion, projection for incompressibility +- **RBF-FD Method**: Uses Polyharmonic Spline basis functions φ(r) = r^m with polynomial augmentation for enhanced accuracy +- **Staggered Grid**: Pressure on P-grid (xy_s), velocity on V-grid (xy) for proper incompressible flow discretization +- **Fractional Step Method**: Adams-Bashforth for advection, Crank-Nicolson for diffusion, pressure projection for incompressibility - **Mesh-Free Approach**: Uses scattered node distributions with local stencils (k=35 neighbors typically) +- **Adaptive Time Stepping**: CFL-based time step reduction with comprehensive instability detection ### 2. Code Organization ``` -src/ # Core simulation functions -├── geometry/ # Geometry-specific mesh generation -├── RBF_PHS_FD_all.m # Core RBF-FD differentiation matrices -├── NS_2d_fractional_step_PHS.m # Time stepping -├── build_*.m # System assembly functions -└── *.m # Utility functions - -tests/ # Unit tests with golden file validation -lib/distmesh/ # External DistMesh library -config.m # Unified configuration system -simulate.m # Main simulation entry point +src/ # Core simulation functions +├── geometry/ # Geometry-specific mesh generation +│ ├── make_cylinder_geometry.m # Cylinder mesh generation +│ ├── make_ellipse_geometry.m # Ellipse mesh generation +│ ├── make_rectangle_geometry.m # Rectangle mesh generation +│ ├── make_airfoil_geometry.m # NACA airfoil mesh generation +│ ├── make_multi_geometry.m # Multi-obstacle mesh generation +│ ├── distmesh2d_robust.m # Robust DistMesh wrapper +│ ├── dairfoil.m # Airfoil distance function +│ ├── dellipse.m # Ellipse distance function +│ └── drounded_rectangle.m # Rounded rectangle distance function +├── RBF_PHS_FD_all.m # Core RBF-FD differentiation matrices +├── NS_2d_fractional_step_PHS.m # Fractional step time integration +├── build_geometry.m # Unified geometry builder +├── build_stencils.m # RBF stencil construction +├── build_pressure_system.m # Pressure Poisson system +├── build_velocity_operators.m # Velocity diffusion operators +├── build_intergrid_ops.m # P-grid <-> V-grid operators +├── precompute_velocity_solvers.m # LU factorization for velocity +├── init_state.m # Initial condition setup +├── seed_tracers.m # Tracer particle seeding +├── advect_tracers.m # Tracer/Stokes particle advection +├── setup_environment.m # Environment and dependency setup +├── stability_diagnostics.m # Comprehensive stability analysis +├── check_solution_validity.m # Solution validation utilities +├── visualization.m # Plotting and visualization +├── knnsearch.m # K-nearest neighbor search +└── nearest_interp.m # Interpolation utilities + +tests/ # Comprehensive test suite +├── BaseGeometryTest.m # Base class for geometry tests +├── TestGolden*.m # Golden file validation tests +├── TestTracers.m # Tracer particle tests +├── Test*Geometry.m # Geometry-specific tests +├── golden/ # Golden reference files +│ ├── GoldenFileGenerator.m # Utility for generating golden files +│ └── *.mat # Golden reference data files +└── setup_paths.m # Test environment setup + +lib/distmesh/ # External DistMesh library +config.m # Unified configuration system +simulate.m # Main simulation entry point ``` +### 3. Simulation Flow +The main simulation follows this sequence: +1. **Environment Setup**: CI/test detection, path configuration, random seed setting +2. **Configuration Loading**: Unified config system with geometry-specific parameters +3. **Geometry Generation**: Mesh generation using geometry-specific functions +4. **System Assembly**: RBF-FD operators, pressure system, velocity operators +5. **State Initialization**: Initial conditions, tracer seeding, memory allocation +6. **Time Integration**: Fractional step method with adaptive time stepping +7. **Stability Monitoring**: Comprehensive instability detection and diagnostics +8. **Tracer Advection**: Particle dynamics (tracer or Stokes) with obstacle avoidance +9. **Visualization**: Live and final result plotting with velocity heatmaps + ## MATLAB Coding Standards ### 1. Function Documentation - **Always include comprehensive function headers** with: - - Purpose description - - Mathematical background when relevant + - Purpose description and mathematical background - Input/output specifications with dimensions [N x M] - - Algorithm overview + - Algorithm overview and implementation details - References to papers when applicable - Example usage for complex functions ```matlab function [D_all] = RBF_PHS_FD_all(xy1, xy_s, Nearest_Idx, k, m, d) - %RBF_PHS_FD_ALL Generate global RBF-FD differentiation matrices using PHS + polynomials - % - % INPUTS: - % xy1 - Target node coordinates [N1 x 2] where derivatives are evaluated - % xy_s - Source node coordinates [N2 x 2] where function values are known - % Nearest_Idx - Index matrix [N1 x k] of k nearest neighbors for each target node - % k - Stencil size (number of nearest neighbors used in local stencils) - % m - Order of PHS-RBF (r^m, typically m=3 for 2D problems) - % d - Degree of polynomial augmentation (-1=none, 0=constant, 1=linear, 2=quadratic) - % - % OUTPUTS: - % D_all - Cell array containing sparse differentiation matrices: - % {Dx, Dy, L, Dxx, Dyy, Dxy} + %RBF_PHS_FD_ALL Generate global RBF-FD differentiation matrices using PHS + polynomials + % + % This function creates global differentiation matrices using Radial Basis Function + % Finite Difference (RBF-FD) method with Polyharmonic Spline (PHS) basis functions + % augmented with polynomial terms for enhanced accuracy. + % + % INPUTS: + % xy1 - Target node coordinates [N1 x 2] where derivatives are evaluated + % xy_s - Source node coordinates [N2 x 2] where function values are known + % Nearest_Idx - Index matrix [N1 x k] of k nearest neighbors for each target node + % k - Stencil size (number of nearest neighbors used in local stencils) + % m - Order of PHS-RBF (r^m, typically m=3 for 2D problems) + % d - Degree of polynomial augmentation (-1=none, 0=constant, 1=linear, 2=quadratic) + % + % OUTPUTS: + % D_all - Cell array containing sparse differentiation matrices: + % {Dx, Dy, L, Dxx, Dyy, Dxy} where: + % Dx, Dy = first derivatives (d/dx, d/dy) + % L = Laplacian (d^2/dx^2 + d^2/dy^2) + % Dxx,Dyy = second derivatives (d^2/dx^2, d^2/dy^2) + % Dxy = mixed derivative (d^2/dxdy) + % + % REFERENCE: + % [1] T. Chu, O. T. Schmidt, "RBF-FD discretization of the Navier-Stokes + % equations on scattered but staggered nodes", JCP 474, 111756, 2023 ``` ### 2. Variable Naming Conventions -- **Matrices**: Use descriptive names with dimensions implied - - `xy1`, `xy_s` for coordinate arrays - - `D_all` for differentiation matrix collections - - `L_inv` for precomputed inverse operators -- **Indices**: Use `_idx` suffix (e.g., `Nearest_Idx`, `boundary_idx`) -- **Geometry**: Use geometry-specific prefixes (e.g., `rect_width`, `ellipse_a`) -- **Physical quantities**: Use standard notation (e.g., `nu` for viscosity, `dt` for time step) +- **Coordinate arrays**: `xy1` (V-grid), `xy_s` (P-grid), `xy1_s` (complete P-grid) +- **Differentiation matrices**: `D_all` (collection), `Dx`, `Dy`, `L` (Laplacian) +- **Precomputed solvers**: `L_inv_s` (pressure), `L_u_inv`, `L_v_inv` (velocity) +- **Boundary arrays**: `boundary_obs`, `boundary_in`, `boundary_out`, `boundary_y` +- **Geometry parameters**: Use geometry-specific prefixes (`rect_width`, `ellipse_a`, `naca_digits`) +- **Physical quantities**: Standard notation (`nu` for viscosity, `dt` for time step, `Re` for Reynolds) +- **Indices**: Use `_idx` suffix (`Nearest_Idx`, `boundary_obs_ids`) ### 3. Configuration System -- **Use unified config.m**: All parameters centralized in geometry-parameterized config -- **Environment variable support**: Support `GEOMETRY_TYPE` environment variable -- **Geometry-specific sections**: Separate configuration blocks for each geometry type -- **No hardcoded parameters**: All numerical values should be configurable +- **Unified config.m**: All parameters centralized with geometry-specific sections +- **Environment variable support**: `GEOMETRY_TYPE`, `CONFIG_FUNC`, `LOAD_CFG_FILE` +- **Hierarchical structure**: Domain, mesh, simulation, geometry, RBF, tracers, logging sections +- **CI/Test overrides**: Reduced time steps and particle counts for testing +- **Debug controls**: Master debug switch with individual feature toggles + +```matlab +% Configuration loading with environment override +geometry_type = getenv('GEOMETRY_TYPE'); +if isempty(geometry_type) + cfg = config(); % Default (cylinder) +else + cfg = config(geometry_type); % Specific geometry +end +``` ### 4. Code Formatting (MISS_HIT Standards) - **Line length**: Maximum 160 characters @@ -78,140 +144,299 @@ function [D_all] = RBF_PHS_FD_all(xy1, xy_s, Nearest_Idx, k, m, d) ```matlab % Good: Proper line breaking for long function calls -[D_all] = RBF_PHS_FD_all(xy1, xy_s, Nearest_Idx, k, ... - config.rbf.phs_order, config.rbf.poly_degree); - -% Good: Proper operator spacing -result = (a + b) * c / d; +[W(:, j + 1), p0] = NS_2d_fractional_step_PHS(dt, nu, W(:, j - 1), W(:, j), ... + Dy, Dx, L_inv_s, L_u_inv, L_v_inv, ... + L0, L_B, L_B_obs, L_W, L_B_y, ... + length(boundary_s), D0_12_x, D0_12_y, ... + D0_21_x, D0_21_y, Dy_b, Dy_b_1, ... + D0_21_x_obs, D0_21_y_obs, p0, W(:, 1), cfg); ``` ### 5. Error Handling and Validation - **Input validation**: Check dimensions and ranges for critical functions -- **Graceful degradation**: Handle missing files, failed convergence -- **Informative error messages**: Include context and suggested fixes -- **CI/Test compatibility**: Use environment variables to detect test/CI mode +- **Solution validity**: Use `check_solution_validity()` for comprehensive NaN/collapse detection +- **Graceful degradation**: Handle missing files, failed convergence with informative messages +- **CI/Test compatibility**: Use environment variables to detect test/CI mode and adjust behavior ```matlab -% Input validation example -if size(xy1, 2) ~= 2 - error('xy1 must be an [N x 2] coordinate array'); +% Comprehensive solution validation +[is_valid, failure_type, failure_details] = check_solution_validity(W(:, j + 1), p0, 1e-12); +if ~is_valid + fprintf('[CRITICAL] %s detected at step %d\n', failure_type, j); + if cfg.logging.snapshot_on_nan + stability_diagnostics('save_debug_snapshot', cfg, j, xy1, W, p0, dt, Dx, Dy, D0_12_x, D0_12_y, G); + end + error('SIMULATION FAILED: %s', failure_type); end - -% CI detection pattern -isCI = strcmpi(getenv('CI'), 'true'); -isTest = strcmpi(getenv('MATLAB_TEST'), 'true'); ``` ## Geometry System Architecture -### 1. Geometry Independence -- **Unified interface**: All geometries use same function signatures +### 1. Supported Geometries +- **Cylinder**: Circular obstacles with configurable radius +- **Ellipse**: Elliptical obstacles with semi-major/minor axes and rotation +- **Rectangle**: Rectangular obstacles with rounded corners and configurable dimensions +- **Airfoil**: NACA 4-digit airfoils with angle of attack and chord length +- **Multi-obstacle**: Multiple obstacles of different types in the same domain + +### 2. Geometry Independence +- **Unified interface**: All geometries use `make_*_geometry(cfg)` pattern +- **Consistent output**: All return geometry structure `G` with standardized fields - **Parameterized functions**: Geometry-specific parameters passed through config -- **Consistent naming**: Use `obstacle_*` for geometry-agnostic parameters -- **Mesh generation**: Each geometry has dedicated `make_*_geometry.m` function +- **Distance functions**: DistMesh-compatible signed distance functions for each geometry -### 2. Geometry-Specific Implementation -- **Distance functions**: Use DistMesh-compatible distance functions (e.g., `drectangle`, `dcircle`) -- **Boundary conditions**: Geometry-agnostic boundary condition application -- **Mesh refinement**: Geometry-specific refinement parameters in config -- **Robust meshing**: Use `distmesh2d_robust.m` for challenging geometries +### 3. Geometry Structure Output +All geometry functions return a structure `G` containing: +```matlab +G.xy, G.xy_s, G.xt % Velocity nodes, pressure nodes, triangulation +G.boundary_in, G.boundary_out % Inlet/outlet boundary nodes (V-grid) +G.boundary_y, G.boundary_obs % Wall/obstacle boundary nodes (V-grid) +G.boundary_in_s, G.boundary_out_s % Inlet/outlet boundary nodes (P-grid) +G.boundary_y_s, G.boundary_obs_s % Wall/obstacle boundary nodes (P-grid) +G.obs_normals_s % Unit normal vectors at obstacle boundaries +G.fd_obs % Obstacle distance function for tracer collision +% Geometry-specific parameters (radius, a, b, rect_width, naca_digits, etc.) +``` -### 3. Adding New Geometries +### 4. Multi-Obstacle System +- **Obstacle array**: `cfg.geometry.obstacles` contains array of obstacle structures +- **Union distance function**: Combines individual SDFs using recursive `min()` operations +- **Obstacle identification**: `G.boundary_obs_ids` tracks which obstacle each boundary node belongs to +- **Flexible configuration**: Mix different obstacle types (cylinder + ellipse + rectangle + airfoil) + +### 5. Adding New Geometries When adding new geometries, follow this pattern: -1. Add geometry parameters to `config.m` +1. Add geometry parameters to `config.m` in appropriate switch case 2. Create `make_[geometry]_geometry.m` in `src/geometry/` -3. Implement distance function if needed +3. Implement distance function if needed (e.g., `d[geometry].m`) 4. Add test class inheriting from `BaseGeometryTest` -5. Generate golden reference file -6. Update documentation +5. Generate golden reference file using `GoldenFileGenerator` +6. Update documentation and README + +## Tracer Particle System + +### 1. Particle Types +- **Tracer particles**: Massless particles that follow fluid streamlines exactly +- **Stokes particles**: Finite-size particles with one-way drag coupling and optional gravity + +### 2. Tracer Dynamics +```matlab +% Tracer particles (cfg.tracers.dynamics = 'tracer') +x^{n+1} = x^n + 0.5*dt*(u^n + u^{n+1}) % Heun's method + +% Stokes particles (cfg.tracers.dynamics = 'stokes') +dX/dt = V % Position equation +dV/dt = (u(X,t) - V)/tau_p + (1 - rho_f/rho_p)*g % Momentum equation +tau_p = (rho_p * d_p^2) / (18 * mu) % Particle response time +``` + +### 3. Interpolation Methods +- **tri_bary**: Triangular barycentric interpolation (fastest, default) +- **scattered**: MATLAB's scatteredInterpolant (fast, robust) +- **idw**: Inverse distance weighting (original method, slower) + +### 4. Obstacle Avoidance +- **Collision detection**: Uses obstacle distance function `G.fd_obs` +- **Pushback mechanism**: Prevents particles from crossing into obstacles +- **Adaptive step reduction**: Reduces time step when collision detected + +### 5. Boundary Conditions +- **Periodic**: Particles wrap around domain boundaries (default) +- **Absorbing**: Particles are removed when hitting boundaries +- **Reflecting**: Particles bounce off boundaries (future feature) + +## Fractional Step Method Implementation + +### 1. Time Integration Scheme +The `NS_2d_fractional_step_PHS` function implements a three-step fractional step method: + +**Step 1: Advection (Adams-Bashforth)** +```matlab +H_U = -U .* (Dx * U) - V .* (Dy * U); % Nonlinear advection terms +U* = U^n + dt * (3/2 * H^n - 1/2 * H^{n-1}); % Second-order extrapolation +``` + +**Step 2: Diffusion (Crank-Nicolson)** +```matlab +(I - dt*nu/2*∇²) * U** = U* + dt*nu/2*∇²*U^n; % Implicit viscous step +``` + +**Step 3: Pressure Projection** +```matlab +∇²p = ∇·U**/dt; % Pressure Poisson equation +U^{n+1} = U** - dt*∇p; % Velocity correction +``` + +### 2. Boundary Condition Treatment +- **No-slip**: `U = 0`, `V = 0` on obstacle boundaries +- **Slip**: `∂U/∂n = 0`, `V = 0` on wall boundaries +- **Inlet**: Prescribed velocity profile +- **Outlet**: `∂U/∂x = 0`, `∂V/∂x = 0` +- **Pressure**: Neumann conditions derived from momentum equations + +### 3. Adaptive Time Stepping +- **CFL monitoring**: Tracks both convective and viscous CFL conditions +- **Automatic reduction**: Reduces time step when CFL thresholds exceeded +- **Stability recovery**: Recomputes time step with smaller dt +- **Failure detection**: Stops simulation if time step reduction fails to stabilize ## Testing Framework -### 1. Test Structure -- **Inheritance-based**: All geometry tests inherit from `BaseGeometryTest` -- **Golden file validation**: Compare against reference solutions -- **Smoke tests**: Basic execution without errors -- **Configuration tests**: Validate config structure and values +### 1. Test Hierarchy +```matlab +BaseGeometryTest % Abstract base class +├── TestGoldenCylinder % Cylinder golden file tests +├── TestGoldenEllipse % Ellipse golden file tests +├── TestGoldenRectangle % Rectangle golden file tests +├── TestGoldenAirfoil % Airfoil golden file tests +├── TestGoldenMulti % Multi-obstacle golden file tests +├── TestGoldenCylinderTracer % Tracer particle tests +├── TestGoldenCylinderStokes % Stokes particle tests +└── TestTracers % General tracer functionality tests +``` -### 2. Test Implementation Pattern +### 2. Golden File System +- **Naming convention**: `[geometry]_Re100_Nt20_dt[timestep]_seed42.mat` +- **Content structure**: Velocity fields, coordinates, metadata, algorithm info +- **Generation**: Automated using `GoldenFileGenerator.generateGoldenFile()` +- **Validation**: Coordinate and velocity field comparison with configurable tolerances + +### 3. Test Implementation Pattern ```matlab classdef TestGolden[Geometry] < BaseGeometryTest - properties (Constant) - GEOMETRY_TYPE = '[geometry]' - EXPECTED_FIELDS = {'field1', 'field2', ...} + properties (Constant) + GEOMETRY_TYPE = '[geometry]' + EXPECTED_FIELDS = {'field1', 'field2', ...} + end + + methods (Test) + function testMatchesGolden(testCase) + % Load golden reference, run simulation, compare results end + end end ``` -### 3. Golden File Management -- **Consistent naming**: `[geometry]_Re100_Nt20_dt0.01_seed42.mat` -- **Reproducible generation**: Fixed random seed, reduced time steps -- **Metadata inclusion**: Algorithm, geometry, parameters, timestamp -- **Tolerance-based comparison**: Separate tolerances for coordinates and velocities - ### 4. CI Integration -- **Environment detection**: Use `CI` and `MATLAB_TEST` environment variables -- **Headless execution**: Disable plotting in CI/test environments -- **Timeout handling**: Reasonable timeouts for different test types -- **Artifact collection**: Save test results and logs +- **Environment detection**: `CI` and `MATLAB_TEST` environment variables +- **Reduced parameters**: Fewer time steps and particles for faster execution +- **Headless execution**: Automatic plotting disable in CI environments +- **Timeout handling**: Reasonable timeouts for different test complexities + +## RBF-FD Method Implementation + +### 1. Core Algorithm (`RBF_PHS_FD_all`) +- **Local stencils**: k-nearest neighbor stencils for each target node +- **PHS basis**: Polyharmonic spline φ(r) = r^m (typically m=3) +- **Polynomial augmentation**: Adds polynomial terms up to degree d for enhanced accuracy +- **Local scaling**: Coordinate scaling for numerical stability +- **Sparse assembly**: Efficient sparse matrix construction + +### 2. Stencil Configuration +```matlab +% Main interior stencils +cfg.rbf.stencil_size_main = 35; % Interior nodes +cfg.rbf.order_main = 28; % PHS order for main stencils +cfg.rbf.poly_degree_main = 3; % Polynomial degree + +% Boundary stencils (smaller for stability) +cfg.rbf.stencil_size_boundary_obstacle = 20; % Obstacle boundaries +cfg.rbf.stencil_size_boundary_wall = 15; % Wall boundaries +cfg.rbf.order_boundary = 8; % Lower order near boundaries +``` + +### 3. Matrix Assembly +- **Differentiation operators**: Dx, Dy (first derivatives), L (Laplacian), Dxx, Dyy, Dxy (second derivatives) +- **Inter-grid operators**: P-grid ↔ V-grid interpolation and differentiation +- **Boundary operators**: Special treatment near obstacles and domain boundaries +- **Precomputed factorizations**: LU decomposition for repeated solves + +## Stability and Diagnostics + +### 1. Comprehensive Monitoring +- **Solution validity**: NaN/Inf detection, velocity collapse detection +- **CFL conditions**: Both convective (u*dt/h) and viscous (nu*dt/h²) monitoring +- **Instability detection**: Velocity/pressure explosion, gradual instabilities +- **Conservation checks**: Mass conservation (divergence), energy monitoring + +### 2. Diagnostic Functions +```matlab +stability_diagnostics('log_step_stats', ...) % Regular step diagnostics +stability_diagnostics('check_cfl_combined', ...) % CFL condition checking +stability_diagnostics('save_debug_snapshot', ...)% Debug snapshot saving +stability_diagnostics('analyze_instability_cause', ...) % Failure analysis +``` + +### 3. Debug Output Control +- **Master switch**: `cfg.logging.enable` controls all debug features +- **Selective logging**: Individual toggles for different diagnostic types +- **Environment override**: `DEBUG=1` environment variable enables all logging +- **CI compatibility**: Automatic debug disable in test environments ## Performance and Numerical Considerations ### 1. Matrix Operations -- **Sparse matrices**: Use sparse storage for differentiation matrices -- **Precomputed factorizations**: Store LU decompositions for repeated solves -- **Vectorized operations**: Avoid loops where possible, leverage MATLAB's JIT -- **Memory efficiency**: Clear large temporary variables when appropriate +- **Sparse storage**: All differentiation matrices stored as sparse +- **Precomputed factorizations**: LU decompositions cached for velocity and pressure systems +- **Vectorized operations**: Minimize loops, leverage MATLAB's optimized linear algebra +- **Memory management**: Clear large temporary variables, preallocate arrays ### 2. Numerical Stability -- **Condition number monitoring**: Check matrix conditioning in critical operations -- **Scaling**: Use local coordinate scaling in RBF-FD weight computation -- **Convergence criteria**: Appropriate tolerances for DistMesh and time stepping +- **Local coordinate scaling**: Prevents ill-conditioning in RBF weight computation +- **Condition number monitoring**: Track matrix conditioning in critical operations - **Regularization**: Handle near-singular cases in RBF systems +- **Convergence criteria**: Appropriate tolerances for DistMesh and time stepping -### 3. Algorithm-Specific Guidelines -- **RBF-FD stencils**: Use k=35 neighbors, m=3 PHS order, d=3 polynomial degree as defaults -- **Time stepping**: dt=0.01, nu=0.01 (Re=100) for standard validation cases -- **Mesh density**: Balance accuracy vs. computational cost, geometry-dependent -- **Boundary treatment**: Higher-order stencils near obstacles, standard elsewhere +### 3. Algorithm Parameters +```matlab +% Recommended settings for stability and accuracy +cfg.rbf.stencil_size_main = 35; % Main stencil size +cfg.rbf.order_main = 28; % PHS order (r^28) +cfg.rbf.poly_degree_main = 3; % Polynomial degree +cfg.simulation.time_step = 1e-2; % Time step (Re=100) +cfg.simulation.viscosity = 0.01; % Viscosity (Re=100) +cfg.mesh.dist = 0.1; % Base mesh spacing +``` ## Development Workflow ### 1. Code Quality Tools -- **Linting**: Use MATLAB Code Analyzer with project-specific ignore list -- **Formatting**: Use MISS_HIT (mh_style) with 160-character line limit -- **Local formatting**: Run `./format.sh` before commits -- **CI validation**: Automatic linting and formatting checks on PRs +- **Linting**: MATLAB Code Analyzer with project-specific configuration +- **Formatting**: MISS_HIT (mh_style) with 160-character line limit +- **Local scripts**: `./format.sh` for code formatting, `./lint.sh` for linting +- **CI validation**: Automatic checks on pull requests ### 2. Version Control Practices -- **Atomic commits**: Each commit should represent a single logical change -- **Descriptive messages**: Include context and impact of changes -- **Branch naming**: Use descriptive names (e.g., `feature/rectangle-geometry`) +- **Atomic commits**: Each commit represents a single logical change +- **Descriptive messages**: Include context, impact, and reasoning +- **Branch naming**: Use descriptive names (`feature/stokes-particles`, `fix/cfl-detection`) - **Clean history**: Squash related commits before merging -### 3. Documentation Updates -- **README maintenance**: Keep usage examples and feature lists current -- **Function documentation**: Update headers when changing interfaces -- **Configuration documentation**: Document new parameters and their effects -- **Mathematical documentation**: Include relevant equations and references +### 3. Documentation Maintenance +- **Function headers**: Update when changing interfaces or algorithms +- **Configuration docs**: Document new parameters and their effects +- **README updates**: Keep usage examples and feature lists current +- **Mathematical docs**: Include relevant equations and references ## File Organization Rules ### 1. Directory Structure -- **src/**: Core simulation functions only -- **src/geometry/**: Geometry-specific mesh generation and utilities -- **tests/**: All test files, including golden references -- **lib/**: External libraries (DistMesh) -- **old/**: Archive directory for deprecated code (avoid accumulation) - -### 2. File Naming -- **Functions**: Use descriptive names matching primary function (e.g., `build_pressure_system.m`) -- **Classes**: Use PascalCase (e.g., `BaseGeometryTest.m`) -- **Scripts**: Use lowercase with underscores (e.g., `setup_paths.m`) -- **Golden files**: Follow pattern `[geometry]_Re[reynolds]_Nt[steps]_dt[timestep]_seed[seed].mat` - -### 3. Dependencies -- **Minimal external dependencies**: Only DistMesh library required +- **src/**: Core simulation functions and algorithms +- **src/geometry/**: Geometry-specific mesh generation and distance functions +- **tests/**: All test files, base classes, and golden references +- **tests/golden/**: Golden reference files and generation utilities +- **lib/**: External libraries (DistMesh only) +- **debug/**: Debug snapshots and diagnostic output (created at runtime) + +### 2. File Naming Conventions +- **Core functions**: Descriptive names matching primary function (`build_pressure_system.m`) +- **Geometry functions**: `make_[geometry]_geometry.m` pattern +- **Distance functions**: `d[geometry].m` pattern (DistMesh convention) +- **Test classes**: `Test[Feature][Geometry].m` pattern +- **Golden files**: `[geometry]_Re[reynolds]_Nt[steps]_dt[timestep]_seed[seed].mat` + +### 3. Dependencies and Interfaces +- **Minimal dependencies**: Only DistMesh library required - **Self-contained**: All custom functions in src/ directory - **Clear interfaces**: Well-defined function signatures between modules - **Backward compatibility**: Avoid breaking changes to public interfaces @@ -219,55 +444,64 @@ end ## Mathematical and Physical Constraints ### 1. Dimensional Consistency -- **Coordinate systems**: Use consistent units throughout (typically dimensionless) -- **Physical parameters**: Reynolds number, time scales, length scales -- **Boundary conditions**: Ensure physical consistency across geometries -- **Validation cases**: Use established benchmark problems when possible +- **Coordinate systems**: Dimensionless units throughout (domain typically [-8,24] × [-8,8]) +- **Physical parameters**: Reynolds number Re = UL/ν, time scales, length scales +- **Boundary conditions**: Physically consistent across all geometries +- **Validation cases**: Use established benchmark problems (Re=100 cylinder flow) ### 2. Numerical Method Constraints -- **Stencil size**: Minimum k > (d+1)(d+2)/2 + m for stability -- **Time step**: CFL condition and diffusion stability limits -- **Mesh quality**: Avoid highly skewed or degenerate elements -- **Convergence**: Monitor residuals and solution changes +- **Stencil size**: k ≥ (d+1)(d+2)/2 + effective RBF terms for stability +- **Time step limits**: CFL < 0.5 for convection, ν*dt/h² < 0.5 for diffusion +- **Mesh quality**: Avoid highly skewed elements, maintain adequate resolution +- **Convergence monitoring**: Track residuals and solution changes ### 3. Geometry Constraints -- **Obstacle placement**: Ensure adequate domain size relative to obstacle -- **Mesh resolution**: Sufficient points to resolve boundary layers -- **Aspect ratios**: Reasonable geometry proportions for numerical stability -- **Boundary separation**: Adequate spacing between different boundary types +- **Domain sizing**: Ensure 10+ obstacle diameters to boundaries for accurate far-field +- **Mesh resolution**: Sufficient boundary layer resolution (typically 20+ points around obstacle) +- **Aspect ratios**: Reasonable proportions for numerical stability +- **Multi-obstacle spacing**: Adequate separation to avoid mesh generation issues ## Debugging and Troubleshooting Guidelines -### 1. Common Issues -- **Matrix singularity**: Check stencil sizes and polynomial degrees -- **Mesh generation failure**: Adjust DistMesh parameters, use robust version -- **Convergence problems**: Reduce time step, check boundary conditions -- **Memory issues**: Use sparse matrices, clear temporary variables - -### 2. Diagnostic Tools -- **Visualization**: Plot meshes, velocity fields, pressure contours -- **Condition numbers**: Monitor matrix conditioning -- **Residual tracking**: Plot convergence histories +### 1. Common Issues and Solutions +- **Matrix singularity**: Check stencil sizes, reduce polynomial degree, verify mesh quality +- **Mesh generation failure**: Use `distmesh2d_robust`, adjust refinement parameters +- **CFL violations**: Enable adaptive time stepping or reduce initial time step +- **Tracer collision issues**: Increase obstacle margin, check distance function accuracy +- **Memory issues**: Reduce mesh density, use sparse matrices, clear temporary variables + +### 2. Diagnostic Tools and Techniques +- **Mesh visualization**: Plot node distributions, boundary classifications +- **Velocity field plots**: Heatmaps, streamlines, vector fields +- **Stability monitoring**: CFL tracking, condition number monitoring +- **Debug snapshots**: Automatic saving of simulation state on failure - **Profiling**: Use MATLAB profiler for performance bottlenecks -### 3. Validation Strategies -- **Method of manufactured solutions**: Test with known analytical solutions +### 3. Validation and Verification +- **Golden file comparison**: Automated regression testing against reference solutions - **Grid convergence**: Verify solution convergence with mesh refinement -- **Benchmark comparisons**: Compare against established results -- **Conservation checks**: Verify mass conservation and energy balance +- **Benchmark validation**: Compare against established results (cylinder drag, vortex shedding) +- **Conservation checks**: Monitor mass conservation (divergence), energy balance +- **Method of manufactured solutions**: Test with known analytical solutions ## Security and Robustness -### 1. Input Sanitization -- **File path validation**: Check for valid paths and permissions -- **Parameter bounds**: Validate numerical parameters are within reasonable ranges -- **Configuration validation**: Ensure required fields are present and valid -- **Environment variable handling**: Safe parsing of environment inputs - -### 2. Error Recovery -- **Graceful failures**: Provide meaningful error messages and recovery suggestions -- **Resource cleanup**: Ensure proper cleanup of temporary files and variables -- **State preservation**: Maintain simulation state for debugging -- **Logging**: Comprehensive logging for debugging and monitoring - -This comprehensive set of rules should guide development while maintaining the mathematical rigor and computational efficiency required for this RBF-FD Navier-Stokes simulation codebase. +### 1. Input Validation and Sanitization +- **Configuration validation**: Check required fields, parameter ranges, type consistency +- **File path validation**: Verify paths exist and are accessible +- **Environment variable handling**: Safe parsing with fallback defaults +- **Geometry parameter bounds**: Ensure physical reasonableness of obstacle dimensions + +### 2. Error Recovery and Resilience +- **Graceful failure modes**: Informative error messages with suggested fixes +- **Automatic recovery**: Adaptive time stepping, mesh regeneration on failure +- **State preservation**: Save simulation state for post-mortem analysis +- **Resource cleanup**: Proper cleanup of temporary files and variables + +### 3. Testing and Quality Assurance +- **Comprehensive test coverage**: Unit tests, integration tests, golden file validation +- **Continuous integration**: Automated testing on multiple MATLAB versions +- **Performance monitoring**: Track simulation performance and memory usage +- **Documentation testing**: Verify examples and usage instructions work correctly + +This comprehensive set of rules reflects the current state of the RBF-FD Navier-Stokes codebase and should guide development while maintaining mathematical rigor, computational efficiency, and code quality standards. \ No newline at end of file diff --git a/.github/scripts/lint.m b/.github/scripts/lint.m index 964ebc6..ac9fa2d 100644 --- a/.github/scripts/lint.m +++ b/.github/scripts/lint.m @@ -61,6 +61,6 @@ fclose(fid); if hasIssues - error('Code Analyzer found issues. See .github/scripts/lint.log.'); + error('Code Analyzer found issues.'); end end diff --git a/config.m b/config.m index ae52437..edf3043 100644 --- a/config.m +++ b/config.m @@ -150,10 +150,37 @@ config.visualization.plot_tick_x = [-5, 0, 5, 10, 15]; config.visualization.color_axis_range = 1e-0; + %% Tracers (Lagrangian point particles advected by the flow) + config.tracers.enable = true; % Enable tracer particles + config.tracers.num_particles = 100; % Number of tracers (normal runs) + config.tracers.num_particles_ci = 30; % Fewer tracers for CI/tests + + % Particle dynamics (one-way Stokes drag coupling) + config.tracers.dynamics = 'tracer'; % 'tracer' (default, massless tracers) or 'stokes' + config.tracers.rho_f = 1.0; % Fluid density (dimensionless) + config.tracers.rho_p = 1000.0; % Particle density (dimensionless) + config.tracers.diameter = 1e-3; % Particle diameter (dimensionless) + config.tracers.gravity_enable = false; % Enable gravitational/buoyancy force + config.tracers.g = [0, -1]; % Gravitational acceleration vector (dimensionless) + config.tracers.vel_init = 'fluid'; % Initial particle velocity: 'fluid' or 'zero' + + % Particle parameters + config.tracers.k_neighbors = 10; % IDW neighbors for velocity interpolation + config.tracers.idw_power = 2; % IDW power (2 = standard) + config.tracers.obstacle_margin = 3 * config.mesh.boundary_eps; % Seed clearance from obstacles + config.tracers.wall_margin = 3 * config.mesh.boundary_eps; % Seed clearance from walls/domain + config.tracers.max_seed_tries = 20; % Try this many batches before stopping early + config.tracers.integrator = 'heun'; % 'heun' for trapezoidal method + config.tracers.pushback_on_obstacle = true; % Prevent crossing into obstacle during advection + config.tracers.max_pushback_iters = 5; % Limit for step shrinking when crossing obstacle + config.tracers.interp_method = 'tri_bary'; % 'tri_bary' (fastest), 'scattered' (fast), or 'idw' (original) + config.tracers.knn_refresh_interval = 10; % Steps between full KNN rebuild (IDW mode) + config.tracers.periodic = true; % Periodic wrap on domain boundaries + % Live visualization (disabled in CI/tests) config.visualization.live_enable = true; % Show live heatmap during simulation - config.visualization.live_frequency = 10; % Update every N steps - config.visualization.heatmap_nx = 200; % Heatmap grid resolution in x + config.visualization.live_frequency = 50; % Update every N steps + config.visualization.heatmap_nx = 300; % Heatmap grid resolution in x % ny is inferred from domain aspect ratio to keep pixels square-ish %% Logging and Debug Parameters diff --git a/lint.sh b/lint.sh index f806b45..6835cab 100755 --- a/lint.sh +++ b/lint.sh @@ -14,7 +14,6 @@ BLUE='\033[0;34m' NC='\033[0m' # No Color echo -e "${BLUE}[LINT] MATLAB Code Analyzer${NC}" -echo "==========================" # Function to run MATLAB linting run_matlab_lint() { @@ -24,10 +23,6 @@ run_matlab_lint() { return 0 else echo -e "${RED}[FAILED] MATLAB Code Analyzer found issues${NC}" - if [ -f ".github/scripts/lint.log" ]; then - echo -e "${YELLOW}Issues found:${NC}" - cat .github/scripts/lint.log - fi return 1 fi } @@ -42,9 +37,5 @@ else echo "" echo -e "${RED}[FAILED] MATLAB Code Analyzer found issues!${NC}" echo -e "${YELLOW}Please fix the issues above before committing.${NC}" - echo "" - echo -e "${BLUE}Additional checks:${NC}" - echo -e " * For style issues: ${YELLOW}./format.sh${NC} (check/fix formatting)" - echo -e " * For auto-fixes: ${YELLOW}./fix.sh${NC} (attempt automatic fixes)" exit 1 fi diff --git a/lint_fix.sh b/lint_fix.sh deleted file mode 100755 index 6afb075..0000000 --- a/lint_fix.sh +++ /dev/null @@ -1,101 +0,0 @@ -#!/bin/bash -# MATLAB Code Auto-Fixer -# This script automatically fixes MATLAB code issues where possible -# -# Usage: ./lint_fix.sh - -set -e - -# Colors for output -RED='\033[0;31m' -GREEN='\033[0;32m' -YELLOW='\033[1;33m' -BLUE='\033[0;34m' -NC='\033[0m' # No Color - -echo -e "${BLUE}[LINT-FIX] MATLAB Code Auto-Fixer${NC}" -echo "==================================" - -# Function to run MATLAB auto-fixes -run_matlab_fixes() { - echo -e "${YELLOW}Running MATLAB automatic fixes...${NC}" - - # Create a temporary MATLAB script for auto-fixing - cat > temp_fix_script.m << 'EOF' -try - % Analyze code issues in src and tests directories - fprintf('Analyzing code issues...\n'); - issues = codeIssues({'src', 'tests'}, 'IncludeSubfolders', true); - - if height(issues.Issues) == 0 - fprintf('No fixable issues found.\n'); - else - fprintf('Found %d total issues, attempting to fix...\n', height(issues.Issues)); - - % Try to fix issues automatically - % Note: fix() may not work in all MATLAB versions or with all issue types - try - fixedIssues = fix(issues); - catch fixError - fprintf('Auto-fix not available or failed: %s\n', fixError.message); - fprintf('Manual fixes may be required for some issues.\n'); - fixedIssues = []; - end - - if ~isempty(fixedIssues) && height(fixedIssues.Issues) > 0 - fprintf('Successfully fixed %d issues:\n', height(fixedIssues.Issues)); - for i = 1:height(fixedIssues.Issues) - issue = fixedIssues.Issues(i,:); - fprintf(' - %s:%d - %s [%s]\n', issue.FullFilename{1}, issue.LineNumber, issue.Description{1}, issue.CheckID{1}); - end - else - fprintf('No issues could be automatically fixed.\n'); - end - end - - fprintf('MATLAB auto-fix completed.\n'); -catch ME - fprintf('Error during MATLAB auto-fix: %s\n', ME.message); -end - -% Clean up -if exist('temp_fix_script.m', 'file') - delete('temp_fix_script.m'); -end -EOF - - # Run the MATLAB script - if matlab -batch "run('temp_fix_script.m')" -nodesktop -nosplash; then - echo -e "${GREEN}[SUCCESS] MATLAB auto-fix completed!${NC}" - # Clean up the temp script if MATLAB didn't - [ -f "temp_fix_script.m" ] && rm -f temp_fix_script.m - return 0 - else - echo -e "${RED}[FAILED] MATLAB auto-fix encountered errors${NC}" - # Clean up the temp script if MATLAB didn't - [ -f "temp_fix_script.m" ] && rm -f temp_fix_script.m - return 1 - fi -} - - -# Main execution -echo "Attempting to automatically fix MATLAB code issues..." -echo "" - -# Run MATLAB auto-fixes -if run_matlab_fixes; then - echo "" - echo "==================================" - echo -e "${GREEN}[SUCCESS] MATLAB auto-fixes completed! ✓${NC}" - echo -e "${YELLOW}[INFO] Run ./lint.sh to check if all MATLAB issues are resolved${NC}" - echo -e "${YELLOW}[INFO] Run ./format.sh to check/fix any style issues${NC}" - exit 0 -else - echo "" - echo "==================================" - echo -e "${RED}[FAILED] MATLAB auto-fixes failed! ✗${NC}" - echo -e "${YELLOW}[INFO] Some issues may require manual fixing${NC}" - echo -e "${YELLOW}[INFO] Run ./lint.sh to see remaining MATLAB issues${NC}" - exit 1 -fi diff --git a/simulate.m b/simulate.m index 18287ce..d952d69 100644 --- a/simulate.m +++ b/simulate.m @@ -1,5 +1,8 @@ clc; -clear; +% Clear workspace unless running in test/CI mode +if ~strcmpi(getenv('MATLAB_TEST'), 'true') && ~strcmpi(getenv('CI'), 'true') + clear; +end % Add required paths for src functions and lib dependencies scriptDir = fileparts(mfilename('fullpath')); @@ -239,7 +242,7 @@ D0_21_x_obs = P.D0_21_x_obs; D0_21_y_obs = P.D0_21_y_obs; -%% 5) Build inter-grid operators (P-grid ↔ V-grid) +%% 5) Build inter-grid operators (P-grid <-> V-grid) [D0_21_x, D0_21_y, D0_12_x, D0_12_y] = build_intergrid_ops(G, xy, xy1, xy_s, xy1_s, S, cfg); %% 6) Build velocity operators and boundary conditions @@ -306,6 +309,48 @@ end [W, p0] = init_state(xy1, xy1_s, boundary_obs, Nt_alloc); +% Initialize tracers (seed outside obstacles) +tracers = seed_tracers(cfg, G, isCI, isTest); + +% Initialize particle velocities for Stokes dynamics +if strcmpi(cfg.tracers.dynamics, 'stokes') && ~isempty(tracers) + % Extract initial velocity field components + NV = size(xy1, 1); + U0 = W(1:NV, 1); % u-velocity at t=0 + V0 = W(NV + 1:end, 1); % v-velocity at t=0 + + % Initialize particle velocities based on configuration + if strcmpi(cfg.tracers.vel_init, 'fluid') + % Initialize particle velocities to local fluid velocity + Fu = scatteredInterpolant(xy1(:, 1), xy1(:, 2), U0, 'linear', 'nearest'); + Fv = scatteredInterpolant(xy1(:, 1), xy1(:, 2), V0, 'linear', 'nearest'); + tracer_vel = [Fu(tracers(:, 1), tracers(:, 2)), Fv(tracers(:, 1), tracers(:, 2))]; + fprintf('Initialized %d Stokes particles with fluid velocity\n', size(tracers, 1)); + else + % Initialize particle velocities to zero + tracer_vel = zeros(size(tracers)); + fprintf('Initialized %d Stokes particles with zero velocity\n', size(tracers, 1)); + end + + % Report Stokes parameters + nu = cfg.simulation.viscosity; + rho_f = cfg.tracers.rho_f; + rho_p = cfg.tracers.rho_p; + dp = cfg.tracers.diameter; + mu = nu * rho_f; + tau_p = (rho_p * dp^2) / (18 * mu); + fprintf('Stokes parameters: tau_p = %.3e, rho_p/rho_f = %.1f, d_p = %.3e\n', tau_p, rho_p / rho_f, dp); + + if cfg.tracers.gravity_enable + fprintf('Gravity enabled: g = [%.3f, %.3f]\n', cfg.tracers.g(1), cfg.tracers.g(2)); + end +else + tracer_vel = []; + if ~isempty(tracers) + fprintf('Initialized %d tracer particles\n', size(tracers, 1)); + end +end + % Define useful index lengths for boundary handling L_B = length(boundary_obs) + length(boundary_in); % Total special boundaries L_B_y = length(boundary_y); % Wall boundaries @@ -380,6 +425,21 @@ rethrow(ME); end + % Advect tracers with latest velocity field (Heun: uses W(:,j) and W(:,j+1)) + if cfg.tracers.enable && ~isempty(tracers) + W_prev = W(:, j); + W_now = W(:, j + 1); + domain_struct = struct('x_min', x_min, 'x_max', x_max, 'y_min', y_min, 'y_max', y_max); + + if strcmpi(cfg.tracers.dynamics, 'stokes') + % Stokes particle dynamics with velocity evolution + [tracers, tracer_vel] = advect_tracers(tracers, tracer_vel, dt, xy1, W_prev, W_now, cfg, G.fd_obs, domain_struct); + else + % Tracer particle dynamics (original behavior) + tracers = advect_tracers(tracers, dt, xy1, W_prev, W_now, cfg, G.fd_obs, domain_struct); + end + end + % Comprehensive instability detection if cfg.logging.immediate_nan_check || cfg.logging.comprehensive_instability_check instability_detected = false; @@ -435,7 +495,7 @@ instability_type = failure_type; % Immediate notification of detection - fprintf('\n🚨 IMMEDIATE DETECTION at step %d: %s\n', j, failure_type); + fprintf('\n[ALERT] IMMEDIATE DETECTION at step %d: %s\n', j, failure_type); % Find representative bad index based on failure type if contains(failure_type, 'NaN/Inf in velocity') @@ -522,7 +582,7 @@ else fprintf('RECOMMENDATION: Lower CFL thresholds or reduce initial time step\n'); end - elseif strcmp(instability_type, 'Velocity collapse (all velocities → 0)') + elseif strcmp(instability_type, 'Velocity collapse (all velocities -> 0)') fprintf('LIKELY CAUSE: Numerical damping or boundary condition issues\n'); fprintf('RECOMMENDATION: Check boundary conditions and reduce viscosity\n'); elseif contains(instability_type, 'explosion') @@ -643,7 +703,7 @@ [is_valid, failure_type, failure_details] = check_solution_validity(W(:, j + 1), p_for_check, 1e-12); if ~is_valid - fprintf('\n🚨 FATAL ERROR: %s after time step reduction!\n', failure_type); + fprintf('\n[FATAL] ERROR: %s after time step reduction!\n', failure_type); fprintf('Step %d: Time step was reduced from %.3e to %.3e but solution is still corrupted.\n', ... j, dt / cfg.adaptive_dt.reduction_factor, dt); fprintf('This indicates the solution cannot be recovered by time step reduction alone.\n'); @@ -735,7 +795,11 @@ % Live max-|V| heatmap every N steps (interactive only) if doPlot && ~isCI && ~isTest && isfield(cfg.visualization, 'live_enable') && cfg.visualization.live_enable && ... mod(j, cfg.visualization.live_frequency) == 0 - visualization('live_heatmap', cfg, xy1, W(:, j + 1), j, x_min, x_max, y_min, y_max); + if strcmpi(cfg.tracers.dynamics, 'stokes') && exist('tracer_vel', 'var') + visualization('live_heatmap', cfg, xy1, W(:, j + 1), j, x_min, x_max, y_min, y_max, tracers, tracer_vel); + else + visualization('live_heatmap', cfg, xy1, W(:, j + 1), j, x_min, x_max, y_min, y_max, tracers); + end end % Advance simulation time (using current dt, which may have been adapted) @@ -799,7 +863,7 @@ [is_valid, failure_type, failure_details] = check_solution_validity(W(:, j), p_for_check, 1e-12); if ~is_valid - fprintf('🚨 CRITICAL: %s in final solution!\n', failure_type); + fprintf('[CRITICAL] %s in final solution!\n', failure_type); fprintf('The simulation appeared to complete but the solution is corrupted.\n'); % Detailed failure analysis @@ -834,7 +898,7 @@ error('SIMULATION FAILED: %s at step %d. Solution is corrupted.', failure_type, j); else - fprintf('✅ Final solution validation: PASSED\n'); + fprintf('[PASS] Final solution validation: PASSED\n'); fprintf(' Velocity nodes: %d (all valid)\n', failure_details.velocity_nodes); fprintf(' Velocity field: max=%.3e, mean=%.3e\n', failure_details.velocity_max, failure_details.velocity_mean); if ~isempty(p_for_check) @@ -846,4 +910,8 @@ fprintf('=====================================\n\n'); %% 11) Visualization of final results -visualization('final', cfg, doPlot, xy1, W0, Nt, x_min, x_max, y_min, y_max, Dx, Dy); +if strcmpi(cfg.tracers.dynamics, 'stokes') && exist('tracer_vel', 'var') + visualization('final', cfg, doPlot, xy1, W0, Nt, x_min, x_max, y_min, y_max, Dx, Dy, tracers, tracer_vel); +else + visualization('final', cfg, doPlot, xy1, W0, Nt, x_min, x_max, y_min, y_max, Dx, Dy, tracers); +end diff --git a/src/advect_tracers.m b/src/advect_tracers.m new file mode 100644 index 0000000..2094c73 --- /dev/null +++ b/src/advect_tracers.m @@ -0,0 +1,393 @@ +function [tracers, tracer_vel] = advect_tracers(tracers, varargin) + %ADVECT_TRACERS Advect tracers with Heun's method using fluid at t^n and t^{n+1} + % + % This function advances tracer particles using Heun's method (explicit trapezoidal) + % which is second-order accurate and uses the fluid velocity at both the beginning + % and end of the time step. Supports both tracer particles (follow fluid exactly) + % and Stokes particles (one-way coupling with drag and optional gravity). + % + % USAGE (tracer particles - original behavior): + % tracers = advect_tracers(tracers, dt, xy1, W_prev, W_now, cfg, fd_obs, domain) + % + % USAGE (Stokes particles - new physics): + % [tracers, tracer_vel] = advect_tracers(tracers, tracer_vel, dt, xy1, W_prev, W_now, cfg, fd_obs, domain) + % + % INPUTS: + % tracers - [N x 2] tracer positions at t^n + % tracer_vel - [N x 2] particle velocities at t^n (Stokes mode only) + % dt - time step size + % xy1 - [NV x 2] V-grid coordinates (interior + boundary nodes) + % W_prev - [2*NV x 1] velocity field at t^n (stacked [U; V]) + % W_now - [2*NV x 1] velocity field at t^{n+1} (stacked [U; V]) + % cfg - configuration structure (uses cfg.tracers.*) + % fd_obs - obstacle signed distance function (negative inside), or [] + % domain - struct with x_min, x_max, y_min, y_max + % + % OUTPUTS: + % tracers - [N x 2] tracer positions at t^{n+1} + % tracer_vel - [N x 2] particle velocities at t^{n+1} (Stokes mode only) + % + % ALGORITHM: + % Tracer: x^{n+1} = x^n + 0.5*dt*(u^n + u^{n+1}) + % where u^n = u(x^n, t^n) and u^{n+1} = u(x^n + dt*u^n, t^{n+1}) + % Stokes: dX/dt = V, dV/dt = (u(X,t) - V)/tau_p + (1 - rho_f/rho_p)*g + % Integrated with Heun's method using same velocity interpolation + + if isempty(tracers) + tracer_vel = []; + return + end + + % Parse input arguments to determine mode (tracer vs Stokes) + if nargin == 8 + % Tracer mode: advect_tracers(tracers, dt, xy1, W_prev, W_now, cfg, fd_obs, domain) + [dt, xy1, W_prev, W_now, cfg, fd_obs, domain] = deal(varargin{:}); + is_stokes = false; + tracer_vel = []; + elseif nargin == 9 + % Stokes mode: advect_tracers(tracers, tracer_vel, dt, xy1, W_prev, W_now, cfg, fd_obs, domain) + [tracer_vel, dt, xy1, W_prev, W_now, cfg, fd_obs, domain] = deal(varargin{:}); + is_stokes = true; + else + error('Invalid number of arguments. Use 8 for tracer mode or 9 for Stokes mode.'); + end + + % Validate Stokes mode configuration + if is_stokes && ~strcmpi(cfg.tracers.dynamics, 'stokes') + warning('Stokes mode called but cfg.tracers.dynamics is not ''stokes''. Proceeding with Stokes dynamics.'); + end + + % Extract velocity components + NV = size(xy1, 1); + U_prev = W_prev(1:NV); % u-velocity at t^n + V_prev = W_prev(NV + 1:end); % v-velocity at t^n + U_now = W_now(1:NV); % u-velocity at t^{n+1} + V_now = W_now(NV + 1:end); % v-velocity at t^{n+1} + + % IDW interpolation parameters + k = cfg.tracers.k_neighbors; + p = cfg.tracers.idw_power; + eps_w = 1e-12; % Small regularization to avoid division by zero + + % Persistent caches for fast interpolation across calls + persistent FuPrev FvPrev FuNow FvNow cached_NV step_count cached_idx DT tri_nv + if isempty(step_count) + step_count = 0; + end + step_count = step_count + 1; + + % Nested function for velocity interpolation using optimized methods + function vel = interp_vel(pts, U, V) + %INTERP_VEL Fast velocity interpolation at arbitrary points + % + % INPUTS: + % pts - [N x 2] query points + % U,V - [NV x 1] velocity components on V-grid + % + % OUTPUT: + % vel - [N x 2] interpolated velocities [u, v] + + if isempty(pts) + vel = []; + return + end + + % Determine interpolation method + method = 'scattered'; + if isfield(cfg, 'tracers') && isfield(cfg.tracers, 'interp_method') + method = cfg.tracers.interp_method; + end + + switch lower(method) + case 'scattered' + % Micro-optimized: separate interpolants per time level to reduce value updates + % Determine which time level we're interpolating (prev or now) + is_prev_time = isequal(U, U_prev) && isequal(V, V_prev); + + if is_prev_time + % Handle previous time level (t^n) + if isempty(FuPrev) || isempty(FvPrev) || isempty(cached_NV) || cached_NV ~= NV + FuPrev = scatteredInterpolant(xy1(:, 1), xy1(:, 2), U, 'linear', 'none'); + FvPrev = scatteredInterpolant(xy1(:, 1), xy1(:, 2), V, 'linear', 'none'); + cached_NV = NV; + else + FuPrev.Values = U; + FvPrev.Values = V; + end + u = FuPrev(pts(:, 1), pts(:, 2)); + v = FvPrev(pts(:, 1), pts(:, 2)); + else + % Handle current time level (t^{n+1}) + if isempty(FuNow) || isempty(FvNow) || isempty(cached_NV) || cached_NV ~= NV + FuNow = scatteredInterpolant(xy1(:, 1), xy1(:, 2), U, 'linear', 'none'); + FvNow = scatteredInterpolant(xy1(:, 1), xy1(:, 2), V, 'linear', 'none'); + cached_NV = NV; + else + FuNow.Values = U; + FvNow.Values = V; + end + u = FuNow(pts(:, 1), pts(:, 2)); + v = FvNow(pts(:, 1), pts(:, 2)); + end + vel = [u, v]; + + case 'tri_bary' + % Use persistent triangulation (declared at function top level) + if isempty(DT) || isempty(tri_nv) || tri_nv ~= NV + DT = delaunayTriangulation(xy1(:, 1), xy1(:, 2)); + tri_nv = NV; + end + + % Locate triangles and compute barycentric weights + T = pointLocation(DT, pts); % [Nq x 1], NaN if outside hull + inside = ~isnan(T); + Nq = size(pts, 1); + u = zeros(Nq, 1); + v = zeros(Nq, 1); + + if any(inside) + tri = DT.ConnectivityList; % [Nt x 3] + verts = tri(T(inside), :); % [Ni x 3] vertex indices + bc = cartesianToBarycentric(DT, T(inside), pts(inside, :)); % [Ni x 3] + Uv = [U(verts(:, 1)) U(verts(:, 2)) U(verts(:, 3))]; + Vv = [V(verts(:, 1)) V(verts(:, 2)) V(verts(:, 3))]; + u(inside) = sum(bc .* Uv, 2); + v(inside) = sum(bc .* Vv, 2); + end + + if any(~inside) + q = pts(~inside, :); + idx1 = knnsearch(q, xy1, 1); % nearest-neighbor fallback + u(~inside) = U(idx1); + v(~inside) = V(idx1); + end + + vel = [u, v]; + + otherwise % 'idw' (original), but with cached neighbor sets + refresh_int = 10; + if isfield(cfg.tracers, 'knn_refresh_interval') + refresh_int = max(1, cfg.tracers.knn_refresh_interval); + end + + % Rebuild KNN neighbor indices only every refresh_int steps or on size change + need_refresh = isempty(cached_idx) || (size(cached_idx, 1) ~= size(pts, 1)) || (mod(step_count - 1, refresh_int) == 0); + if need_refresh + cached_idx = knnsearch(pts, xy1, k); % [N x k] indices into xy1 + end + + % Compute distances only to cached neighbors (O(N*k)), then IDW weights + Nq = size(pts, 1); + u = zeros(Nq, 1); + v = zeros(Nq, 1); + for i = 1:Nq + nei = cached_idx(i, :); + dx = xy1(nei, 1) - pts(i, 1); + dy = xy1(nei, 2) - pts(i, 2); + d = sqrt(dx .* dx + dy .* dy); + w = 1 ./ max(d, eps_w).^p; + ws = sum(w); + u(i) = (w(:)' * U(nei(:))) / ws; + v(i) = (w(:)' * V(nei(:))) / ws; + end + vel = [u, v]; + end + end + + % Helper function for periodic boundary wrapping + function Pw = wrap_periodic(P, domain) + %WRAP_PERIODIC Apply periodic boundary conditions to particle positions + % + % INPUTS: + % P - [N x 2] particle positions + % domain - struct with x_min, x_max, y_min, y_max + % + % OUTPUT: + % Pw - [N x 2] wrapped positions + + Pw = P; + Lx = domain.x_max - domain.x_min; + Ly = domain.y_max - domain.y_min; + + if Lx > 0 + Pw(:, 1) = domain.x_min + mod(P(:, 1) - domain.x_min, Lx); + end + if Ly > 0 + Pw(:, 2) = domain.y_min + mod(P(:, 2) - domain.y_min, Ly); + end + end + + % Branch on dynamics mode + if is_stokes + % Stokes particle dynamics: dX/dt = V, dV/dt = (u(X,t) - V)/tau_p + (1 - rho_f/rho_p)*g + + % Compute Stokes parameters + nu = cfg.simulation.viscosity; % Kinematic viscosity + rho_f = cfg.tracers.rho_f; % Fluid density + rho_p = cfg.tracers.rho_p; % Particle density + dp = cfg.tracers.diameter; % Particle diameter + mu = nu * rho_f; % Dynamic viscosity + tau_p = (rho_p * dp^2) / (18 * mu); % Particle response time + + % Gravitational acceleration (dimensionless) + g = [0, 0]; + if isfield(cfg.tracers, 'gravity_enable') && cfg.tracers.gravity_enable && ... + isfield(cfg.tracers, 'g') && ~isempty(cfg.tracers.g) + g = cfg.tracers.g; + end + + % Current state + Pn = tracers; % x^n (current positions) + Vn = tracer_vel; % V^n (current particle velocities) + + % Interpolate fluid velocity at current positions and time t^n + un = interp_vel(Pn, U_prev, V_prev); % u(x^n, t^n) + + % Right-hand side at t^n: dV/dt = (u - V)/tau_p + (1 - rho_f/rho_p)*g + buoyancy_factor = 1 - rho_f / rho_p; + Rn = (un - Vn) ./ tau_p + buoyancy_factor .* g; + + % Predictor step + Ppred = Pn + dt * Vn; % Predicted position + Vpred = Vn + dt * Rn; % Predicted velocity + + % Apply periodic wrapping to predictor positions if enabled + if isfield(cfg.tracers, 'periodic') && cfg.tracers.periodic + Ppred = wrap_periodic(Ppred, domain); + end + + % Interpolate fluid velocity at predicted positions and time t^{n+1} + unp1 = interp_vel(Ppred, U_now, V_now); % u(x_pred, t^{n+1}) + + % Right-hand side at t^{n+1} + Rnp1 = (unp1 - Vpred) ./ tau_p + buoyancy_factor .* g; + + % Heun corrector step + Vnew = Vn + 0.5 * dt * (Rn + Rnp1); % Final particle velocity + Pnew = Pn + 0.5 * dt * (Vn + Vnew); % Final position using average velocity + + else + % Tracer particle dynamics (original behavior): follow fluid exactly + Pn = tracers; % x^n (current positions) + Vn = interp_vel(Pn, U_prev, V_prev); % u(x^n, t^n) + + % Predictor step: estimate position at t^{n+1} using velocity at t^n + Ppred = Pn + dt * Vn; + + % Apply periodic wrapping to predictor positions if enabled + if isfield(cfg.tracers, 'periodic') && cfg.tracers.periodic + Ppred = wrap_periodic(Ppred, domain); + end + + % Corrector step: evaluate velocity at predicted position and t^{n+1} + Vnp1 = interp_vel(Ppred, U_now, V_now); % u(x_pred, t^{n+1}) + + % Final Heun update: average of velocities at t^n and t^{n+1} + Pnew = Pn + 0.5 * dt * (Vn + Vnp1); + end + + % Apply boundary conditions + if isfield(cfg.tracers, 'periodic') && cfg.tracers.periodic + % Periodic wrap instead of clamp + Pnew = wrap_periodic(Pnew, domain); + else + % Enforce domain boundaries (clamp to domain limits) + Pnew(:, 1) = min(max(Pnew(:, 1), domain.x_min), domain.x_max); + Pnew(:, 2) = min(max(Pnew(:, 2), domain.y_min), domain.y_max); + end + + % Optional obstacle avoidance via step shrinking + if cfg.tracers.pushback_on_obstacle && ~isempty(fd_obs) + % Check which particles have entered obstacles + dnew = fd_obs(Pnew); + bad = dnew < 0; % Particles inside obstacles (negative distance) + + if any(bad) + % Apply step shrinking to particles that entered obstacles + if is_stokes + base_pos = Pn(bad, :); % Original positions + base_vel = Vn(bad, :); % Original velocities + step_pos = Pnew(bad, :) - base_pos; % Full position step vectors + step_vel = Vnew(bad, :) - base_vel; % Full velocity step vectors + else + base = Pn(bad, :); % Original positions + step = Pnew(bad, :) - base; % Full step vectors + end + scale = 0.5; % Step reduction factor + + % Iteratively shrink steps until particles are outside obstacles + for it = 1:cfg.tracers.max_pushback_iters + if is_stokes + cand_pos = base_pos + scale * step_pos; % Candidate positions with reduced step + cand_vel = base_vel + scale * step_vel; % Candidate velocities with reduced step + d_cand = fd_obs(cand_pos); % Distance to obstacles + else + cand = base + scale * step; % Candidate positions with reduced step + d_cand = fd_obs(cand); % Distance to obstacles + end + ok = d_cand >= 0; % Particles now outside obstacles + + % Accept positions that are now valid + if is_stokes + base_pos(ok, :) = cand_pos(ok, :); + base_vel(ok, :) = cand_vel(ok, :); + % Further shrink steps for particles still inside + step_pos(~ok, :) = scale * step_pos(~ok, :); + step_vel(~ok, :) = scale * step_vel(~ok, :); + else + base(ok, :) = cand(ok, :); + % Further shrink steps for particles still inside + step(~ok, :) = scale * step(~ok, :); + end + + % Stop if all particles are now outside + if all(ok) + break + end + end + + % Update positions (and velocities for Stokes) for particles that were inside obstacles + if is_stokes + Pnew(bad, :) = base_pos; + Vnew(bad, :) = base_vel; + else + Pnew(bad, :) = base; + end + end + end + + % Handle boundary conditions: periodic wrapping or removal + if isfield(cfg.tracers, 'periodic') && cfg.tracers.periodic + % Periodic boundaries: keep all particles (already wrapped) + tracers = Pnew; + if is_stokes + tracer_vel = Vnew; + end + else + % Remove particles that have left or are very close to leaving the domain + x_coords = Pnew(:, 1); + y_coords = Pnew(:, 2); + + % Add a small margin to remove particles before they get stuck at boundaries + boundary_margin = 0.1; % Remove particles within this distance of boundary + + % Check which particles are still well inside the domain + inside_domain = (x_coords >= domain.x_min + boundary_margin) & ... + (x_coords <= domain.x_max - boundary_margin) & ... + (y_coords >= domain.y_min + boundary_margin) & ... + (y_coords <= domain.y_max - boundary_margin); + + % Keep only particles that are inside the domain + tracers = Pnew(inside_domain, :); + if is_stokes + tracer_vel = Vnew(inside_domain, :); + end + + % Report how many particles were removed (only if some were removed) + num_removed = sum(~inside_domain); + if num_removed > 0 + fprintf('Removed %d tracer particles that approached domain boundary\n', num_removed); + end + end + +end diff --git a/src/geometry/make_airfoil_geometry.m b/src/geometry/make_airfoil_geometry.m index 0c113d5..2c861bd 100644 --- a/src/geometry/make_airfoil_geometry.m +++ b/src/geometry/make_airfoil_geometry.m @@ -254,4 +254,7 @@ G.idx_far_boundaries_V = idx_far_boundaries_V; % Velocity nodes far from boundaries G.idx_far_boundaries_P = idx_far_boundaries_P; % Pressure nodes far from boundaries + % Obstacle signed distance function (negative inside airfoil) + G.fd_obs = @(p) dairfoil(p, airfoil_x_center, airfoil_y_center, naca_digits, chord_length, angle_of_attack); + end diff --git a/src/geometry/make_cylinder_geometry.m b/src/geometry/make_cylinder_geometry.m index 4b905d9..0aa2c13 100644 --- a/src/geometry/make_cylinder_geometry.m +++ b/src/geometry/make_cylinder_geometry.m @@ -161,4 +161,7 @@ G.idx_far_boundaries_V = idx_far_boundaries_V; % Velocity nodes far from boundaries G.idx_far_boundaries_P = idx_far_boundaries_P; % Pressure nodes far from boundaries + % Obstacle signed distance function (negative inside obstacle, ~0 on surface) + G.fd_obs = @(p) dcircle(p, 0, 0, radius); + end diff --git a/src/geometry/make_ellipse_geometry.m b/src/geometry/make_ellipse_geometry.m index 4f4d6df..e47b196 100644 --- a/src/geometry/make_ellipse_geometry.m +++ b/src/geometry/make_ellipse_geometry.m @@ -183,4 +183,7 @@ G.idx_far_boundaries_V = idx_far_boundaries_V; % Velocity nodes far from boundaries G.idx_far_boundaries_P = idx_far_boundaries_P; % Pressure nodes far from boundaries + % Obstacle signed distance function (negative inside ellipse) + G.fd_obs = @(p) dellipse(p, 0, 0, a, b); + end diff --git a/src/geometry/make_multi_geometry.m b/src/geometry/make_multi_geometry.m index 771c349..9c57549 100644 --- a/src/geometry/make_multi_geometry.m +++ b/src/geometry/make_multi_geometry.m @@ -162,6 +162,9 @@ G.idx_far_boundaries_V = idx_far_boundaries_V; % Velocity nodes far from boundaries G.idx_far_boundaries_P = idx_far_boundaries_P; % Pressure nodes far from boundaries + % Obstacle union signed distance (negative inside any obstacle) + G.fd_obs = fd_obs_union; + % Store obstacle configuration for reference G.obstacles = obstacles; G.num_obstacles = num_obstacles; diff --git a/src/geometry/make_rectangle_geometry.m b/src/geometry/make_rectangle_geometry.m index b0cd983..567ac34 100644 --- a/src/geometry/make_rectangle_geometry.m +++ b/src/geometry/make_rectangle_geometry.m @@ -226,4 +226,7 @@ G.idx_far_boundaries_V = idx_far_boundaries_V; % Velocity nodes far from boundaries G.idx_far_boundaries_P = idx_far_boundaries_P; % Pressure nodes far from boundaries + % Obstacle signed distance function (negative inside rounded rectangle) + G.fd_obs = @(p) drounded_rectangle(p, rect_x_center, rect_y_center, rect_width, rect_height, corner_radius); + end diff --git a/src/seed_tracers.m b/src/seed_tracers.m new file mode 100644 index 0000000..227593c --- /dev/null +++ b/src/seed_tracers.m @@ -0,0 +1,112 @@ +function tracers = seed_tracers(cfg, G, isCI, isTest) + %SEED_TRACERS Initialize tracer particle positions outside obstacles and walls + % + % This function seeds tracer particles randomly throughout the domain while + % ensuring they are placed outside obstacles (with margin) and away from + % domain boundaries. Uses rejection sampling with batched generation for + % efficiency. + % + % INPUTS: + % cfg - Configuration structure containing tracer parameters + % G - Geometry structure (must contain fd_obs obstacle distance function) + % isCI - Boolean flag indicating CI environment (reduces particle count) + % isTest - Boolean flag indicating test environment (reduces particle count) + % + % OUTPUTS: + % tracers - [N x 2] array of tracer positions [x, y] + % + % ALGORITHM: + % 1. Determine target number of particles based on environment + % 2. Define seeding domain with wall margins + % 3. Use rejection sampling in batches to find valid positions + % 4. Reject points inside obstacles using G.fd_obs distance function + % 5. Return up to target number of valid positions + + if ~cfg.tracers.enable + tracers = []; + return + end + + % Set random seed for reproducible tracer seeding + if isfield(cfg.simulation, 'random_seed') && ~isempty(cfg.simulation.random_seed) + % Use a different seed for tracers to avoid conflicts with mesh generation + tracer_seed = cfg.simulation.random_seed + 1000; % Offset to avoid mesh conflicts + rng(tracer_seed); + end + + % Determine target number of particles + N_target = cfg.tracers.num_particles; + if isCI || isTest + N_target = cfg.tracers.num_particles_ci; + end + + % Define seeding domain with wall margins + x_min = cfg.domain.x_min + cfg.tracers.wall_margin; + x_max = cfg.domain.x_max - cfg.tracers.wall_margin; + y_min = cfg.domain.y_min + cfg.tracers.wall_margin; + y_max = cfg.domain.y_max - cfg.tracers.wall_margin; + + % Safety check for valid domain + if x_min >= x_max || y_min >= y_max + warning('Invalid seeding domain after applying wall margins. No tracers seeded.'); + tracers = []; + return + end + + margin_obs = cfg.tracers.obstacle_margin; + max_tries = max(1, cfg.tracers.max_seed_tries); + + % Initialize output array + tracers = zeros(N_target, 2); + n = 0; % Number of valid tracers found + tries = 0; % Number of batch attempts + + fprintf('Seeding %d tracer particles (target)...\n', N_target); + + while n < N_target && tries < max_tries + % Generate batch of candidate positions + batch_size = max(1000, ceil(1.5 * (N_target - n))); % Oversample for efficiency + xs = x_min + (x_max - x_min) * rand(batch_size, 1); + ys = y_min + (y_max - y_min) * rand(batch_size, 1); + pts = [xs, ys]; + + % Filter out points inside obstacles (with margin) + if isfield(G, 'fd_obs') && ~isempty(G.fd_obs) + d = G.fd_obs(pts); % Negative inside obstacle, positive outside + ok = d > margin_obs; % Keep points outside obstacle with margin + else + ok = true(size(xs)); % No obstacles, all points valid + end + + % Take up to remaining needed particles + n_valid = sum(ok); + n_take = min(n_valid, N_target - n); + + if n_take > 0 + pts_ok = pts(ok, :); + tracers(n + (1:n_take), :) = pts_ok(1:n_take, :); + n = n + n_take; + end + + tries = tries + 1; + + % Progress reporting + if mod(tries, 5) == 0 && n < N_target + fprintf(' Batch %d: found %d/%d tracers (%.1f%% success rate this batch)\n', ... + tries, n, N_target, 100 * n_take / batch_size); + end + end + + % Trim array if we couldn't find enough valid positions + if n < N_target + tracers = tracers(1:n, :); + if n == 0 + warning('Could not seed any valid tracer particles after %d attempts', tries); + else + fprintf('Warning: Only seeded %d/%d tracers after %d attempts\n', n, N_target, tries); + end + else + fprintf('Successfully seeded %d tracer particles in %d batches\n', n, tries); + end + +end diff --git a/src/visualization.m b/src/visualization.m index f3873b9..47f68ed 100644 --- a/src/visualization.m +++ b/src/visualization.m @@ -37,23 +37,33 @@ end end -function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_min, y_max, Dx, Dy) +function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_min, y_max, Dx, Dy, tracers, tracer_vel) %VISUALIZE_FINAL_RESULTS Create visualization of final simulation results % % This function creates plots of the final vorticity field and velocity components. % % INPUTS: - % cfg - Configuration structure - % doPlot - Boolean indicating if plotting should be enabled - % xy1 - Complete velocity grid coordinates - % W_final - Final velocity solution vector [U; V] (single time step) - % ~ - Number of time steps (unused, kept for compatibility) - % x_min - Minimum x-coordinate of domain - % x_max - Maximum x-coordinate of domain - % y_min - Minimum y-coordinate of domain - % y_max - Maximum y-coordinate of domain - % Dx - RBF-FD differentiation matrix for d/dx operator - % Dy - RBF-FD differentiation matrix for d/dy operator + % cfg - Configuration structure + % doPlot - Boolean indicating if plotting should be enabled + % xy1 - Complete velocity grid coordinates + % W_final - Final velocity solution vector [U; V] (single time step) + % ~ - Number of time steps (unused, kept for compatibility) + % x_min - Minimum x-coordinate of domain + % x_max - Maximum x-coordinate of domain + % y_min - Minimum y-coordinate of domain + % y_max - Maximum y-coordinate of domain + % Dx - RBF-FD differentiation matrix for d/dx operator + % Dy - RBF-FD differentiation matrix for d/dy operator + % tracers - (Optional) [N x 2] tracer particle positions + % tracer_vel - (Optional) [N x 2] tracer particle velocities (Stokes mode) + + % Handle optional arguments + if nargin < 12 + tracers = []; + end + if nargin < 13 + tracer_vel = []; + end % Only plot if plotting is enabled if ~doPlot @@ -111,6 +121,22 @@ function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_m colorbar; % Add colorbar for vorticity grid on; grid minor; + + % Overlay tracer particles if available + if ~isempty(tracers) && cfg.tracers.enable + if ~isempty(tracer_vel) && strcmpi(cfg.tracers.dynamics, 'stokes') + % Color Stokes particles by velocity magnitude + vel_mag = sqrt(tracer_vel(:, 1).^2 + tracer_vel(:, 2).^2); + scatter(tracers(:, 1), tracers(:, 2), 12, vel_mag, 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.3, 'DisplayName', sprintf('Stokes particles (%d)', size(tracers, 1))); + fprintf('Overlaid %d Stokes particles colored by velocity magnitude on vorticity plot\n', size(tracers, 1)); + else + scatter(tracers(:, 1), tracers(:, 2), 8, 'w', 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.5, 'DisplayName', sprintf('Tracers (%d)', size(tracers, 1))); + fprintf('Overlaid %d tracer particles on vorticity plot\n', size(tracers, 1)); + end + end + drawnow; % Update display immediately %% Plot 2: U-Velocity Component @@ -139,6 +165,19 @@ function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_m grid on; grid minor; + % Overlay tracer particles if available + if ~isempty(tracers) && cfg.tracers.enable + if ~isempty(tracer_vel) && strcmpi(cfg.tracers.dynamics, 'stokes') + % Color Stokes particles by velocity magnitude + vel_mag = sqrt(tracer_vel(:, 1).^2 + tracer_vel(:, 2).^2); + scatter(tracers(:, 1), tracers(:, 2), 12, vel_mag, 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.3, 'DisplayName', sprintf('Stokes particles (%d)', size(tracers, 1))); + else + scatter(tracers(:, 1), tracers(:, 2), 8, 'w', 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.5, 'DisplayName', sprintf('Tracers (%d)', size(tracers, 1))); + end + end + fprintf('U-velocity field statistics:\n'); fprintf(' Min u-velocity: %.6f\n', min(U_final)); fprintf(' Max u-velocity: %.6f\n', max(U_final)); @@ -172,6 +211,19 @@ function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_m grid on; grid minor; + % Overlay tracer particles if available + if ~isempty(tracers) && cfg.tracers.enable + if ~isempty(tracer_vel) && strcmpi(cfg.tracers.dynamics, 'stokes') + % Color Stokes particles by velocity magnitude + vel_mag = sqrt(tracer_vel(:, 1).^2 + tracer_vel(:, 2).^2); + scatter(tracers(:, 1), tracers(:, 2), 12, vel_mag, 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.3, 'DisplayName', sprintf('Stokes particles (%d)', size(tracers, 1))); + else + scatter(tracers(:, 1), tracers(:, 2), 8, 'w', 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.5, 'DisplayName', sprintf('Tracers (%d)', size(tracers, 1))); + end + end + fprintf('V-velocity field statistics:\n'); fprintf(' Min v-velocity: %.6f\n', min(V_final)); fprintf(' Max v-velocity: %.6f\n', max(V_final)); @@ -205,6 +257,19 @@ function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_m grid on; grid minor; + % Overlay tracer particles if available + if ~isempty(tracers) && cfg.tracers.enable + if ~isempty(tracer_vel) && strcmpi(cfg.tracers.dynamics, 'stokes') + % Color Stokes particles by velocity magnitude + vel_mag = sqrt(tracer_vel(:, 1).^2 + tracer_vel(:, 2).^2); + scatter(tracers(:, 1), tracers(:, 2), 12, vel_mag, 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.3, 'DisplayName', sprintf('Stokes particles (%d)', size(tracers, 1))); + else + scatter(tracers(:, 1), tracers(:, 2), 8, 'w', 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.5, 'DisplayName', sprintf('Tracers (%d)', size(tracers, 1))); + end + end + fprintf('Velocity magnitude statistics:\n'); fprintf(' Min |V|: %.6f\n', min(velocity_magnitude)); fprintf(' Max |V|: %.6f\n', max(velocity_magnitude)); @@ -213,7 +278,7 @@ function visualize_final_results(cfg, doPlot, xy1, W_final, ~, x_min, x_max, y_m drawnow; end -function visualize_live_heatmap(cfg, xy1, W_step, step_idx, x_min, x_max, y_min, y_max) +function visualize_live_heatmap(cfg, xy1, W_step, step_idx, x_min, x_max, y_min, y_max, tracers, tracer_vel) %VISUALIZE_LIVE_HEATMAP Live interpolated velocity magnitude heatmap % % This function creates a live-updating heatmap showing the velocity magnitude @@ -221,11 +286,21 @@ function visualize_live_heatmap(cfg, xy1, W_step, step_idx, x_min, x_max, y_min, % interpolation to create a smooth field over a regular grid. % % INPUTS: - % cfg - Configuration struct containing visualization parameters - % xy1 - Velocity grid coordinates [Nv x 2] - % W_step - Current velocity vector [2*Nv x 1] at this step (stacked [U; V]) - % step_idx - Current time step index (for title) + % cfg - Configuration struct containing visualization parameters + % xy1 - Velocity grid coordinates [Nv x 2] + % W_step - Current velocity vector [2*Nv x 1] at this step (stacked [U; V]) + % step_idx - Current time step index (for title) % x_min, x_max, y_min, y_max - Domain bounds + % tracers - (Optional) [N x 2] tracer particle positions + % tracer_vel - (Optional) [N x 2] tracer particle velocities (Stokes mode) + + % Handle optional arguments + if nargin < 9 + tracers = []; + end + if nargin < 10 + tracer_vel = []; + end % Extract velocity components Nv = size(xy1, 1); @@ -325,6 +400,21 @@ function visualize_live_heatmap(cfg, xy1, W_step, step_idx, x_min, x_max, y_min, grid(ax, 'on'); grid(ax, 'minor'); + % Overlay tracer particles if available + if ~isempty(tracers) && cfg.tracers.enable + hold(ax, 'on'); + if ~isempty(tracer_vel) && strcmpi(cfg.tracers.dynamics, 'stokes') + % Color Stokes particles by velocity magnitude + vel_mag = sqrt(tracer_vel(:, 1).^2 + tracer_vel(:, 2).^2); + scatter(ax, tracers(:, 1), tracers(:, 2), 15, vel_mag, 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.5, 'DisplayName', sprintf('Stokes particles (%d)', size(tracers, 1))); + else + scatter(ax, tracers(:, 1), tracers(:, 2), 12, 'w', 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.8, 'DisplayName', sprintf('Tracers (%d)', size(tracers, 1))); + end + hold(ax, 'off'); + end + else % Update existing figure ax = fig.CurrentAxes; @@ -343,6 +433,28 @@ function visualize_live_heatmap(cfg, xy1, W_step, step_idx, x_min, x_max, y_min, title(ax, sprintf('Velocity Magnitude |V| — step %d', step_idx)); apply_xy_axes(cfg, ax, x_min, x_max, y_min, y_max); clim(ax, [0, cmax]); + + % Update tracer particles (remove old ones first to avoid accumulation) + tracer_scatter = findobj(ax, 'Type', 'scatter'); + % Remove all existing tracer scatter plots + if ~isempty(tracer_scatter) + delete(tracer_scatter); + end + + % Add current tracer positions if enabled + if ~isempty(tracers) && cfg.tracers.enable + hold(ax, 'on'); + if ~isempty(tracer_vel) && strcmpi(cfg.tracers.dynamics, 'stokes') + % Color Stokes particles by velocity magnitude + vel_mag = sqrt(tracer_vel(:, 1).^2 + tracer_vel(:, 2).^2); + scatter(ax, tracers(:, 1), tracers(:, 2), 15, vel_mag, 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.5, 'DisplayName', sprintf('Stokes particles (%d)', size(tracers, 1))); + else + scatter(ax, tracers(:, 1), tracers(:, 2), 12, 'w', 'filled', 'MarkerEdgeColor', 'k', ... + 'LineWidth', 0.8, 'DisplayName', sprintf('Tracers (%d)', size(tracers, 1))); + end + hold(ax, 'off'); + end end % Update display with rate limiting to avoid excessive redraws diff --git a/test_stokes_debug.m b/test_stokes_debug.m new file mode 100644 index 0000000..0519ecb --- /dev/null +++ b/test_stokes_debug.m @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/tests/TestGoldenCylinderStokes.m b/tests/TestGoldenCylinderStokes.m new file mode 100644 index 0000000..bda84cc --- /dev/null +++ b/tests/TestGoldenCylinderStokes.m @@ -0,0 +1,114 @@ +classdef TestGoldenCylinderStokes < BaseGeometryTest + % TESTGOLDENCYLINDERSTOKES Test Stokes particle simulation against golden reference + % + % This class tests the cylinder geometry with Stokes particle dynamics + % against a pre-computed golden reference file to ensure numerical + % consistency and detect regressions. + + properties (Constant) + GEOMETRY_TYPE = 'cylinder' + PARTICLE_MODE = 'stokes' + EXPECTED_FIELDS = {'obstacle_radius'} % Geometry fields for BaseGeometryTest + GOLDEN_FIELDS = {'xy1', 'U', 'V', 'meta', 'tracers_final', 'tracer_vel_final'} % Golden file fields + end + + methods (Static) + + function goldenFile = getGoldenFilePath() + % Get the path to the golden reference file for cylinder Stokes particles + goldenFile = fullfile(fileparts(mfilename('fullpath')), 'golden', ... + 'cylinder_stokes_Re100_Nt10_dt0.01_seed42.mat'); + end + + end + + methods (Test) + + function testMatchesGolden(testCase) + % Test that cylinder Stokes simulation matches golden reference + + goldenFile = TestGoldenCylinderStokes.getGoldenFilePath(); + + if ~exist(goldenFile, 'file') + testCase.assumeFail(sprintf('Golden file not found: %s', goldenFile)); + end + + % Set up environment + setenv('CI', 'true'); + setenv('MATLAB_TEST', 'true'); + setup_paths(); + + % Load golden reference + golden = load(goldenFile); + + % Verify golden file structure + testCase.verifyTrue(isfield(golden, 'gold'), 'Golden file should contain gold structure'); + gold = golden.gold; + + for i = 1:length(testCase.GOLDEN_FIELDS) + field = testCase.GOLDEN_FIELDS{i}; + testCase.verifyTrue(isfield(gold, field), sprintf('Golden file should contain field: %s', field)); + end + + % Run simulation with same parameters as golden file + cfg = config(testCase.GEOMETRY_TYPE); + cfg.tracers.dynamics = testCase.PARTICLE_MODE; + cfg.tracers.num_particles_ci = 20; % Match golden file generation + cfg.simulation.num_time_steps_ci = 10; + cfg.simulation.time_step = 0.01; + + % Save config and run simulation + save('test_config.mat', 'cfg'); + setenv('LOAD_CFG_FILE', 'test_config.mat'); + + simulate; + + % Clean up config file + if exist('test_config.mat', 'file') + delete('test_config.mat'); + end + + % Compare velocity field + n = size(xy1, 1); + U_final = W(1:n, end); + V_final = W(n + 1:end, end); + + % Sort for comparison + [~, idx] = sortrows(xy1, [1, 2]); + xy1_sorted = xy1(idx, :); + U_sorted = U_final(idx); + V_sorted = V_final(idx); + + % Compare coordinates + testCase.verifySize(xy1_sorted, size(gold.xy1), 'Coordinate arrays should have same size'); + testCase.verifyEqual(xy1_sorted, gold.xy1, 'AbsTol', testCase.XY_TOL, ... + 'Coordinates should match within tolerance'); + + % Compare velocity components + testCase.verifySize(U_sorted, size(gold.U), 'U velocity should have same size'); + testCase.verifySize(V_sorted, size(gold.V), 'V velocity should have same size'); + + testCase.verifyEqual(U_sorted, gold.U, 'RelTol', testCase.REL_TOL, 'AbsTol', testCase.ABS_TOL, ... + 'U velocity should match within tolerance'); + testCase.verifyEqual(V_sorted, gold.V, 'RelTol', testCase.REL_TOL, 'AbsTol', testCase.ABS_TOL, ... + 'V velocity should match within tolerance'); + + % Compare tracer positions + if isfield(gold, 'tracers_final') && exist('tracers', 'var') + testCase.verifySize(tracers, size(gold.tracers_final), 'Tracer positions should have same size'); + testCase.verifyEqual(tracers, gold.tracers_final, 'AbsTol', testCase.XY_TOL, ... + 'Tracer positions should match within tolerance'); + end + + % Compare tracer velocities (Stokes mode only) + if isfield(gold, 'tracer_vel_final') && exist('tracer_vel', 'var') + testCase.verifySize(tracer_vel, size(gold.tracer_vel_final), 'Tracer velocities should have same size'); + testCase.verifyEqual(tracer_vel, gold.tracer_vel_final, 'AbsTol', testCase.ABS_TOL, ... + 'Tracer velocities should match within tolerance'); + end + + fprintf('✅ Cylinder Stokes golden test passed\n'); + end + + end +end diff --git a/tests/TestGoldenCylinderTracer.m b/tests/TestGoldenCylinderTracer.m new file mode 100644 index 0000000..b5c78c6 --- /dev/null +++ b/tests/TestGoldenCylinderTracer.m @@ -0,0 +1,107 @@ +classdef TestGoldenCylinderTracer < BaseGeometryTest + % TESTGOLDENCYLINDERTRACER Test tracer particle simulation against golden reference + % + % This class tests the cylinder geometry with tracer particle dynamics + % against a pre-computed golden reference file to ensure numerical + % consistency and detect regressions. + + properties (Constant) + GEOMETRY_TYPE = 'cylinder' + PARTICLE_MODE = 'tracer' + EXPECTED_FIELDS = {'obstacle_radius'} % Geometry fields for BaseGeometryTest + GOLDEN_FIELDS = {'xy1', 'U', 'V', 'meta', 'tracers_final'} % Golden file fields + end + + methods (Static) + + function goldenFile = getGoldenFilePath() + % Get the path to the golden reference file for cylinder tracer particles + goldenFile = fullfile(fileparts(mfilename('fullpath')), 'golden', ... + 'cylinder_tracer_Re100_Nt20_dt0.01_seed42.mat'); + end + + end + + methods (Test) + + function testMatchesGolden(testCase) + % Test that cylinder tracer simulation matches golden reference + + goldenFile = TestGoldenCylinderTracer.getGoldenFilePath(); + + if ~exist(goldenFile, 'file') + testCase.assumeFail(sprintf('Golden file not found: %s', goldenFile)); + end + + % Set up environment + setenv('CI', 'true'); + setenv('MATLAB_TEST', 'true'); + setup_paths(); + + % Load golden reference + golden = load(goldenFile); + + % Verify golden file structure + testCase.verifyTrue(isfield(golden, 'gold'), 'Golden file should contain gold structure'); + gold = golden.gold; + + for i = 1:length(testCase.GOLDEN_FIELDS) + field = testCase.GOLDEN_FIELDS{i}; + testCase.verifyTrue(isfield(gold, field), sprintf('Golden file should contain field: %s', field)); + end + + % Run simulation with same parameters as golden file + cfg = config(testCase.GEOMETRY_TYPE); + cfg.tracers.dynamics = testCase.PARTICLE_MODE; + cfg.tracers.num_particles_ci = 30; % Match golden file generation + cfg.simulation.num_time_steps_ci = 20; + cfg.simulation.time_step = 0.01; + + % Save config and run simulation + save('test_config.mat', 'cfg'); + setenv('LOAD_CFG_FILE', 'test_config.mat'); + + simulate; + + % Clean up config file + if exist('test_config.mat', 'file') + delete('test_config.mat'); + end + + % Compare velocity field + n = size(xy1, 1); + U_final = W(1:n, end); + V_final = W(n + 1:end, end); + + % Sort for comparison + [~, idx] = sortrows(xy1, [1, 2]); + xy1_sorted = xy1(idx, :); + U_sorted = U_final(idx); + V_sorted = V_final(idx); + + % Compare coordinates + testCase.verifySize(xy1_sorted, size(gold.xy1), 'Coordinate arrays should have same size'); + testCase.verifyEqual(xy1_sorted, gold.xy1, 'AbsTol', testCase.XY_TOL, ... + 'Coordinates should match within tolerance'); + + % Compare velocity components + testCase.verifySize(U_sorted, size(gold.U), 'U velocity should have same size'); + testCase.verifySize(V_sorted, size(gold.V), 'V velocity should have same size'); + + testCase.verifyEqual(U_sorted, gold.U, 'RelTol', testCase.REL_TOL, 'AbsTol', testCase.ABS_TOL, ... + 'U velocity should match within tolerance'); + testCase.verifyEqual(V_sorted, gold.V, 'RelTol', testCase.REL_TOL, 'AbsTol', testCase.ABS_TOL, ... + 'V velocity should match within tolerance'); + + % Compare tracer positions + if isfield(gold, 'tracers_final') && exist('tracers', 'var') + testCase.verifySize(tracers, size(gold.tracers_final), 'Tracer positions should have same size'); + testCase.verifyEqual(tracers, gold.tracers_final, 'AbsTol', testCase.XY_TOL, ... + 'Tracer positions should match within tolerance'); + end + + fprintf('✅ Cylinder Tracer golden test passed\n'); + end + + end +end diff --git a/tests/TestTracers.m b/tests/TestTracers.m new file mode 100644 index 0000000..975de5a --- /dev/null +++ b/tests/TestTracers.m @@ -0,0 +1,224 @@ +classdef TestTracers < matlab.unittest.TestCase + % TESTTRACERS Test tracer particle functionality + % + % This class tests the tracer particle system including: + % - Seeding outside obstacles + % - Different interpolation methods (tri_bary, scattered, idw) + % - Periodic vs non-periodic boundary conditions + % - Golden file validation for consistency + + properties (Constant) + % Test tolerances + XY_TOL = 1e-2 % Tolerance for tracer position matching + REL_TOL = 1e-2 % Relative tolerance for tracer comparison + ABS_TOL = 1e-2 % Absolute tolerance for tracer comparison + end + + methods (Static) + + function goldenFile = getGoldenFilePath() + % Get the path to the golden reference file for tracers + goldenFile = fullfile(fileparts(mfilename('fullpath')), 'golden', ... + 'cylinder_tracer_Re100_Nt20_dt0.01_seed42.mat'); + end + + end + + methods (Test) + + function testTracerSeeding(testCase) + % Test that tracers are seeded correctly outside obstacles + + % Set up environment + setenv('CI', 'true'); + setenv('MATLAB_TEST', 'true'); + setup_paths(); + + % Load config and build geometry + cfg = config('cylinder'); + cfg.tracers.enable = true; + cfg.tracers.num_particles_ci = 50; + + G = build_geometry(cfg); + + % Seed tracers + tracers = seed_tracers(cfg, G, true, true); + + % Verify tracers are outside obstacle + distances = G.fd_obs(tracers); + testCase.verifyTrue(all(distances >= 0), ... + 'All tracers should be seeded outside obstacles'); + + % Verify correct number of tracers + testCase.verifyEqual(size(tracers, 1), cfg.tracers.num_particles_ci, ... + 'Should seed the correct number of tracers'); + end + + function testInterpolationMethods(testCase) + % Test different interpolation methods produce reasonable results + + % Set up environment + setenv('CI', 'true'); + setenv('MATLAB_TEST', 'true'); + setup_paths(); + + methods = {'tri_bary', 'scattered', 'idw'}; + results = cell(length(methods), 1); + + % Store testCase reference to avoid clearing issues + tc = testCase; + + for method_idx = 1:length(methods) + % Load config fresh for each test + cfg = config('cylinder'); + cfg.tracers.enable = true; + cfg.tracers.num_particles_ci = 20; + cfg.simulation.num_time_steps_ci = 3; + cfg.tracers.interp_method = methods{method_idx}; + + % Run short simulation in try-catch to protect variables + try + % Run simulation and capture results immediately + simulate; + + % Store results (tracers should exist in workspace after simulate) + if exist('tracers', 'var') + results{method_idx} = tracers; + else + results{method_idx} = MException('TestTracers:NoTracers', 'Tracers not found after simulation'); + end + + catch ME + % Store error for later verification + results{method_idx} = ME; + end + end + + % Verify all methods produced results (after all simulations complete) + for method_idx = 1:length(methods) + if isa(results{method_idx}, 'MException') + tc.verifyFail(sprintf('Simulation failed with %s method: %s', methods{method_idx}, results{method_idx}.message)); + else + tc.verifyTrue(~isempty(results{method_idx}), ... + sprintf('%s method should produce tracer results', methods{method_idx})); + end + end + end + + function testPeriodicBoundaries(testCase) + % Test periodic vs non-periodic boundary conditions + + % Set up environment + setenv('CI', 'true'); + setenv('MATLAB_TEST', 'true'); + setup_paths(); + + % Store testCase reference to avoid clearing issues + tc = testCase; + + % Test periodic boundaries + cfg = config('cylinder'); + cfg.tracers.enable = true; + cfg.tracers.num_particles_ci = 30; + cfg.tracers.periodic = true; + cfg.simulation.num_time_steps_ci = 5; + + periodic_tracers = []; + nonperiodic_tracers = []; + + try + simulate; + if exist('tracers', 'var') + periodic_tracers = tracers; + end + catch ME + tc.verifyFail(sprintf('Periodic simulation failed: %s', ME.message)); + return + end + + % Test non-periodic boundaries + cfg = config('cylinder'); % Reload config + cfg.tracers.enable = true; + cfg.tracers.num_particles_ci = 30; + cfg.tracers.periodic = false; + cfg.simulation.num_time_steps_ci = 5; + + try + simulate; + if exist('tracers', 'var') + nonperiodic_tracers = tracers; + end + catch ME + tc.verifyFail(sprintf('Non-periodic simulation failed: %s', ME.message)); + return + end + + % Verify both modes produced results + tc.verifyTrue(~isempty(periodic_tracers), ... + 'Periodic mode should produce tracer results'); + tc.verifyTrue(~isempty(nonperiodic_tracers), ... + 'Non-periodic mode should produce tracer results'); + + % Periodic mode should preserve more particles (no removal) + tc.verifyGreaterThanOrEqual(size(periodic_tracers, 1), size(nonperiodic_tracers, 1), ... + 'Periodic mode should preserve at least as many particles as non-periodic mode'); + end + + function testMatchesGolden(testCase) + % Test that tracer system matches golden reference + + goldenFile = TestTracers.getGoldenFilePath(); + + if ~exist(goldenFile, 'file') + testCase.assumeFail(sprintf('Golden file not found: %s', goldenFile)); + end + + % Set up environment + setenv('CI', 'true'); + setenv('MATLAB_TEST', 'true'); + setup_paths(); + + % Store testCase reference to avoid clearing issues + tc = testCase; + + % Load golden reference + golden = load(goldenFile); + + % Verify golden file structure + tc.verifyTrue(isfield(golden, 'gold'), 'Golden file should contain gold structure'); + gold = golden.gold; + + % Run simulation with same parameters as golden file + cfg = config('cylinder'); + cfg.tracers.enable = true; + cfg.tracers.dynamics = 'tracer'; % Explicitly set tracer mode + cfg.tracers.num_particles_ci = 30; % Match golden file (30 tracers) + cfg.tracers.interp_method = 'tri_bary'; + cfg.tracers.periodic = true; + cfg.simulation.num_time_steps_ci = 20; % Match golden file (20 steps) + + simulate; + + % Verify tracers exist after simulation + if ~exist('tracers', 'var') + tc.verifyFail('Tracers not found after simulation'); + return + end + + % Compare final tracer positions + tc.verifySize(tracers, size(gold.tracers_final), ... + 'Tracer array should have same size as golden reference'); + + % Compare positions with tolerance + diff = abs(tracers - gold.tracers_final); + max_diff = max(diff(:)); + + tc.verifyLessThanOrEqual(max_diff, tc.XY_TOL, ... + sprintf('Maximum tracer position difference (%.6f) should be within tolerance (%.6f)', ... + max_diff, tc.XY_TOL)); + + fprintf('✅ Tracer golden test passed: max difference = %.6f\n', max_diff); + end + + end +end diff --git a/tests/golden/GoldenFileGenerator.m b/tests/golden/GoldenFileGenerator.m index a52590c..70186c1 100644 --- a/tests/golden/GoldenFileGenerator.m +++ b/tests/golden/GoldenFileGenerator.m @@ -6,13 +6,18 @@ methods (Static) - function generateGoldenFile(geometryType) - % Generate a golden file for the specified geometry + function generateGoldenFile(geometryType, particleMode) + % Generate a golden file for the specified geometry and particle mode % % Parameters: % geometryType - String: 'cylinder', 'ellipse', 'rectangle', or 'airfoil' + % particleMode - String: 'tracer' (default) or 'stokes' - fprintf('=== Creating %s Golden File ===\n', upper(geometryType)); + if nargin < 2 + particleMode = 'tracer'; + end + + fprintf('=== Creating %s %s Golden File ===\n', upper(geometryType), upper(particleMode)); % Set up environment setenv('CI', 'true'); @@ -21,22 +26,43 @@ function generateGoldenFile(geometryType) % Store parameters before simulation (simulate may clear variables) geomType = geometryType; + partMode = particleMode; % Load config for metadata first % Store geometry type before simulate clears variables geometryType = geomType; + particleMode = partMode; cfg = config(geometryType); + % Configure particle mode + cfg.tracers.dynamics = particleMode; + if strcmpi(particleMode, 'stokes') + cfg.tracers.num_particles_ci = 20; % Fewer Stokes particles for faster testing + cfg.simulation.num_time_steps_ci = 10; % Shorter simulation for Stokes + cfg.simulation.time_step = 0.01; % Appropriate time step + fprintf('Configured for Stokes particles: %d particles, %d steps\n', ... + cfg.tracers.num_particles_ci, cfg.simulation.num_time_steps_ci); + else + cfg.tracers.num_particles_ci = 30; % Standard tracer count + cfg.simulation.num_time_steps_ci = 20; % Standard steps + fprintf('Configured for tracer particles: %d particles, %d steps\n', ... + cfg.tracers.num_particles_ci, cfg.simulation.num_time_steps_ci); + end + + % Save config for simulate to use + save('temp_golden_config.mat', 'cfg'); + setenv('LOAD_CFG_FILE', 'temp_golden_config.mat'); + % Run the simulation - fprintf('Running %s simulation...\n', geometryType); - setenv('CONFIG_FUNC', 'config'); - setenv('GEOMETRY_TYPE', geometryType); + fprintf('Running %s %s simulation...\n', geometryType, particleMode); simulate; % Restore variables that may have been cleared by simulate - geometryType = getenv('GEOMETRY_TYPE'); - cfg = config(geometryType); + geometryType = geomType; + particleMode = partMode; + cfg = load('temp_golden_config.mat', 'cfg'); + cfg = cfg.cfg; % Extract the final state % Ensure variables exist from simulate script @@ -52,6 +78,18 @@ function generateGoldenFile(geometryType) fprintf(' Max |U|: %.6f\n', max(abs(U))); fprintf(' Max |V|: %.6f\n', max(abs(V))); + % Extract tracer data if available + tracers_final = []; + tracer_vel_final = []; + if exist('tracers', 'var') && ~isempty(tracers) + tracers_final = tracers; + fprintf(' Number of tracers: %d\n', size(tracers_final, 1)); + if strcmpi(particleMode, 'stokes') && exist('tracer_vel', 'var') && ~isempty(tracer_vel) + tracer_vel_final = tracer_vel; + fprintf(' Max tracer velocity: %.6f\n', max(sqrt(sum(tracer_vel_final.^2, 2)))); + end + end + % Sort nodes for reproducible ordering [~, idx] = sortrows(xy1, [1, 2]); xy1_sorted = xy1(idx, :); @@ -59,13 +97,28 @@ function generateGoldenFile(geometryType) V_sorted = V(idx); % Create metadata - meta.algorithm = sprintf('RBF-FD Navier-Stokes %s flow', geometryType); + meta.algorithm = sprintf('RBF-FD Navier-Stokes %s flow with %s particles', geometryType, particleMode); meta.geometry_type = cfg.geometry.type; + meta.particle_mode = particleMode; meta.reynolds_number = cfg.simulation.reynolds_number; meta.time_step = cfg.simulation.time_step; meta.num_time_steps = cfg.simulation.num_time_steps_ci; meta.random_seed = cfg.simulation.random_seed; + % Add particle-specific metadata + meta.tracers_enable = cfg.tracers.enable; + meta.tracers_dynamics = cfg.tracers.dynamics; + if strcmpi(particleMode, 'stokes') + meta.rho_f = cfg.tracers.rho_f; + meta.rho_p = cfg.tracers.rho_p; + meta.diameter = cfg.tracers.diameter; + meta.gravity_enable = cfg.tracers.gravity_enable; + if cfg.tracers.gravity_enable + meta.gravity = cfg.tracers.g; + end + meta.vel_init = cfg.tracers.vel_init; + end + % Add geometry-specific metadata switch geometryType case 'cylinder' @@ -109,42 +162,64 @@ function generateGoldenFile(geometryType) gold.V = V_sorted; gold.meta = meta; + % Add tracer data if available + if ~isempty(tracers_final) + gold.tracers_final = tracers_final; + if ~isempty(tracer_vel_final) + gold.tracer_vel_final = tracer_vel_final; + end + end + % Save golden file outDir = fullfile('tests', 'golden'); if ~exist(outDir, 'dir') mkdir(outDir); end - outFile = fullfile(outDir, sprintf('%s_Re%d_Nt%d_dt%g_seed%d.mat', ... - geometryType, meta.reynolds_number, meta.num_time_steps, ... + outFile = fullfile(outDir, sprintf('%s_%s_Re%d_Nt%d_dt%g_seed%d.mat', ... + geometryType, particleMode, meta.reynolds_number, meta.num_time_steps, ... meta.time_step, meta.random_seed)); save(outFile, 'gold', '-v7'); - fprintf('=== %s Golden File Created Successfully ===\n', upper(geometryType)); + % Clean up temporary file + if exist('temp_golden_config.mat', 'file') + delete('temp_golden_config.mat'); + end + + fprintf('=== %s %s Golden File Created Successfully ===\n', upper(geometryType), upper(particleMode)); fprintf(' File: %s\n', outFile); fprintf(' Nodes: %d\n', length(gold.U)); fprintf(' Max |U|: %.6f\n', max(abs(gold.U))); fprintf(' Max |V|: %.6f\n', max(abs(gold.V))); fprintf(' Mean U: %.6f\n', mean(gold.U)); fprintf(' Mean V: %.6f\n', mean(gold.V)); + if ~isempty(tracers_final) + fprintf(' Tracers: %d\n', size(tracers_final, 1)); + if ~isempty(tracer_vel_final) + fprintf(' Max tracer vel: %.6f\n', max(sqrt(sum(tracer_vel_final.^2, 2)))); + end + end end function generateAllGoldenFiles() - % generateAllGoldenFiles Generates golden files for all supported geometries + % generateAllGoldenFiles Generates golden files for all supported geometries and particle modes fprintf('=== GENERATING ALL GOLDEN FILES ===\n'); setup_paths(); % Ensure paths are set for all calls geometries = {'cylinder', 'ellipse', 'rectangle', 'airfoil', 'multi'}; + particleModes = {'tracer', 'stokes'}; for i = 1:length(geometries) - try - GoldenFileGenerator.generateGoldenFile(geometries{i}); - catch ME - fprintf('[ERROR] Failed to generate %s golden file: %s\n', geometries{i}, ME.message); + for j = 1:length(particleModes) + try + GoldenFileGenerator.generateGoldenFile(geometries{i}, particleModes{j}); + catch ME + fprintf('[ERROR] Failed to generate %s %s golden file: %s\n', geometries{i}, particleModes{j}, ME.message); + end + % Clear variables for next iteration + clear xy1 W U V xy1_sorted U_sorted V_sorted meta gold tracers tracer_vel tracers_final tracer_vel_final; end - % Clear variables for next iteration - clear xy1 W U V xy1_sorted U_sorted V_sorted meta gold; end fprintf('=== ALL GOLDEN FILES GENERATION COMPLETE ===\n'); end diff --git a/tests/golden/cylinder_stokes_Re100_Nt10_dt0.01_seed42.mat b/tests/golden/cylinder_stokes_Re100_Nt10_dt0.01_seed42.mat new file mode 100644 index 0000000..0d4d765 Binary files /dev/null and b/tests/golden/cylinder_stokes_Re100_Nt10_dt0.01_seed42.mat differ diff --git a/tests/golden/cylinder_tracer_Re100_Nt20_dt0.01_seed42.mat b/tests/golden/cylinder_tracer_Re100_Nt20_dt0.01_seed42.mat new file mode 100644 index 0000000..60cf44f Binary files /dev/null and b/tests/golden/cylinder_tracer_Re100_Nt20_dt0.01_seed42.mat differ