Merge pull request #4 from sbryngelson/tracers

Author: sbryngelson

.cursorrules

@@ -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

.github/scripts/lint.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

config.m

@@ -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

lint.sh

@@ -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