fix ci

Author: sbryngelson

config.m

@@ -168,6 +168,15 @@
   config.tracers.knn_refresh_interval = 10;     % Steps between full KNN rebuild (IDW mode)
   config.tracers.periodic = true;               % Periodic wrap on domain boundaries
 
+  % Particle dynamics (one-way Stokes drag coupling)
+  config.tracers.dynamics = 'massless';         % 'massless' (default, current behavior) 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, 0];                    % Gravitational acceleration vector (dimensionless)
+  config.tracers.vel_init = 'fluid';            % Initial particle velocity: 'fluid' or 'zero'
+
   % Live visualization (disabled in CI/tests)
   config.visualization.live_enable = true;               % Show live heatmap during simulation
   config.visualization.live_frequency = 50;              % Update every N steps

example_stokes_particles.m

@@ -0,0 +1,81 @@
+%% Example: One-way Stokes Drag Particle Simulation
+% This script demonstrates the new Stokes particle functionality
+% which implements one-way coupling with drag and optional gravity.
+
+clear; clc;
+
+% Add required paths
+addpath('src'); 
+addpath('src/geometry');
+
+%% Configure Stokes Particles
+cfg = config('cylinder');
+
+% Enable Stokes particle dynamics (default is 'massless')
+cfg.tracers.dynamics = 'stokes';
+
+% Particle properties (dimensionless)
+cfg.tracers.rho_f = 1.0;        % Fluid density
+cfg.tracers.rho_p = 1000.0;     % Particle density (heavy particles)
+cfg.tracers.diameter = 1e-3;    % Particle diameter
+
+% Initial particle velocity
+cfg.tracers.vel_init = 'fluid'; % 'fluid' or 'zero'
+
+% Optional gravity/buoyancy (disabled by default)
+cfg.tracers.gravity_enable = false;
+cfg.tracers.g = [0, -0.1];      % Gravitational acceleration (if enabled)
+
+% Simulation parameters
+cfg.tracers.num_particles = 100;
+cfg.simulation.num_time_steps = 50;
+cfg.simulation.time_step = 0.01;
+
+% Visualization
+cfg.visualization.live_enable = true;
+cfg.visualization.live_frequency = 10;
+
+%% Display 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 Particle Parameters ===\n');
+fprintf('Particle dynamics: %s\n', cfg.tracers.dynamics);
+fprintf('Reynolds number: %.0f\n', 1/nu);
+fprintf('Particle density ratio ρ_p/ρ_f: %.1f\n', rho_p/rho_f);
+fprintf('Particle diameter: %.3e\n', dp);
+fprintf('Response time τ_p: %.3e\n', tau_p);
+fprintf('Time step dt: %.3e\n', cfg.simulation.time_step);
+fprintf('Stokes number St = τ_p/dt: %.2f\n', tau_p/cfg.simulation.time_step);
+
+if cfg.tracers.gravity_enable
+    fprintf('Gravity: [%.3f, %.3f]\n', cfg.tracers.g(1), cfg.tracers.g(2));
+    fprintf('Buoyancy factor (1-ρ_f/ρ_p): %.3f\n', 1 - rho_f/rho_p);
+else
+    fprintf('Gravity: disabled\n');
+end
+fprintf('===================================\n\n');
+
+%% Run Simulation
+fprintf('Running Stokes particle simulation...\n');
+fprintf('- Particles will have their own velocity that relaxes toward fluid velocity\n');
+fprintf('- Drag force: F = 3πμd(u-V) where u=fluid velocity, V=particle velocity\n');
+fprintf('- Particle equation: dV/dt = (u-V)/τ_p + (1-ρ_f/ρ_p)g\n\n');
+
+% Save configuration and run
+save('stokes_config.mat', 'cfg');
+setenv('LOAD_CFG_FILE', 'stokes_config.mat');
+
+% Run the simulation
+simulate;
+
+% Clean up
+delete('stokes_config.mat');
+
+fprintf('\n=== Simulation Complete ===\n');
+fprintf('The visualization shows Stokes particles colored by their velocity magnitude.\n');
+fprintf('Compare this to massless tracers by setting cfg.tracers.dynamics = ''massless''.\n');

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
@@ -309,6 +312,45 @@
 % 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 massless tracers\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
@@ -388,7 +430,14 @@
     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);
-    tracers = advect_tracers(tracers, dt, xy1, W_prev, W_now, cfg, G.fd_obs, domain_struct);
+
+    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
+      % Massless tracer dynamics (original behavior)
+      tracers = advect_tracers(tracers, dt, xy1, W_prev, W_now, cfg, G.fd_obs, domain_struct);
+    end
   end
 
   % Comprehensive instability detection
@@ -446,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')
@@ -533,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')
@@ -654,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');
@@ -746,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, tracers);
+    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)
@@ -810,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
@@ -845,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)
@@ -857,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, tracers);
+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