Upgrade MATLAB solver to Phase 2 MAC production workflow

README.md

@@ -1,117 +1,89 @@
 # Lid-Driven Cavity Flow Solver in MATLAB
 
-<p align="center">
-  <img src="https://img.shields.io/badge/Status-Completed-brightgreen.svg" alt="Completed">
-  <img src="https://img.shields.io/badge/MATLAB-base%20MATLAB-orange.svg" alt="MATLAB">
-  <img src="https://img.shields.io/badge/License-MIT-lightgrey.svg" alt="MIT License">
-  <a href="https://kandil2001.github.io/">
-    <img src="https://img.shields.io/badge/Portfolio-kandil2001.github.io-2ea44f.svg" alt="Portfolio">
-  </a>
-</p>
+A Phase 2 MATLAB implementation of the two-dimensional incompressible lid-driven cavity benchmark. The production solver now uses a staggered Marker-and-Cell grid and the same strict numerical acceptance logic used by the companion C++ repository.
 
-A completed MATLAB implementation and parameter study of the two-dimensional lid-driven cavity benchmark.
+## Phase 2 production features
 
-This repository is the MATLAB reference implementation for a larger project that compares the same CFD problem across MATLAB, C++, C, Python, OpenMP, MPI, CUDA, and OpenFOAM-oriented workflows. The physical setup is kept consistent so that numerical behavior, accuracy, runtime, implementation style, and scalability can be compared.
-
-## What the project contains
-
-- structured collocated Cartesian grid
-- pseudo-transient pressure-correction algorithm
-- loop-based and vectorized momentum predictors
-- first-order upwind and central convection schemes
+- staggered MAC arrangement for pressure and face velocities
+- projection method with compatible divergence, gradient, and Poisson operators
+- first-order upwind and second-order central convection
 - red-black Gauss-Seidel and red-black SOR pressure solvers
-- Ghia centerline comparison
-- automated field plots, residual histories, validation figures, and CSV summaries
-
-The completed parameter study contains 72 configured combinations:
+- loop-based and vectorized MATLAB momentum predictors
+- strict convergence based on velocity update, divergence `Linf`, divergence `L2`, global mass balance, and pressure convergence
+- required consecutive converged iterations and minimum iteration count
+- Reynolds-number continuation for ordered parameter studies
+- Ghia centerline benchmark comparison
+- standardized summary, history, field, centerline, and MAT outputs
+- named run modes matching the C++ workflow
+
+## Run modes
+
+| Mode | Cases | Purpose |
+|---|---:|---|
+| `single` | 1 | Canonical `N=32`, `Re=100`, upwind, RBSOR regression |
+| `quick` | 2 | Vectorized versus loop implementation on the canonical case |
+| `medium` | 6 | `N=32`, `Re=100/400/1000`, upwind/central, RBSOR |
+| `grid` | 3 | `N=16/32/64`, `Re=100`, central, RBSOR |
+| `re1000` | 1 | Representative `N=64`, `Re=1000`, central, RBSOR case |
+| `full` | 72 | 3 meshes × 3 Reynolds numbers × 2 schemes × 2 pressure solvers × 2 MATLAB implementations |
+
+Open MATLAB in the repository root and run:
 
-```text
-3 meshes × 3 Reynolds numbers × 2 schemes × 2 pressure solvers × 2 implementations
+```matlab
+run_mode('single')
+run_mode('medium')
+run_mode('grid')
+run_mode('re1000')
+run_mode('full')
 ```
 
-## Representative result
-
-This case uses `N = 64`, `Re = 100`, central differencing, RBGS, and the vectorized momentum predictor.
-
-| Flow field | Centerline comparison |
-|---|---|
-| ![Streamlines](assets/figures/case_029_N64_Re100_central_RBGS_vectorized_streamlines.png) | ![Ghia u validation](assets/figures/case_029_N64_Re100_central_RBGS_vectorized_ghia_u.png) |
-| ![Velocity magnitude](assets/figures/case_029_N64_Re100_central_RBGS_vectorized_speed.png) | ![Ghia v validation](assets/figures/case_029_N64_Re100_central_RBGS_vectorized_ghia_v.png) |
+Equivalent scripts are available for Linux systems:
 
-## Numerical approach
-
-The solver advances the nondimensional incompressible Navier-Stokes equations through pseudo-time. Each outer iteration predicts the velocity field, solves a pressure-correction Poisson equation, corrects velocity and pressure, reapplies wall boundary conditions, and records convergence diagnostics.
-
-A detailed description is available in [`docs/METHODOLOGY.md`](docs/METHODOLOGY.md).
-
-## Study observations
+```bash
+bash scripts/run_single.sh
+bash scripts/run_medium.sh
+bash scripts/run_grid.sh
+bash scripts/run_re1000.sh
+bash scripts/run_full.sh
+```
 
-- `44/72` cases met the selected Ghia centerline-error thresholds
-- all `N = 128` cases met the selected validation thresholds
-- coarse high-Reynolds-number cases were less accurate
-- RBSOR reduced pressure-solver iterations and runtime compared with RBGS
-- vectorizing the momentum predictor had a limited effect on total runtime because the pressure solve remained the main cost
+Run the canonical regression test with:
 
-The selected validation thresholds are practical comparison limits, not a substitute for a formal verification, grid-convergence, or uncertainty study. See [`docs/RESULTS.md`](docs/RESULTS.md) for the detailed discussion.
+```bash
+bash scripts/run_tests.sh
+```
 
-![Pressure solver comparison](assets/figures/study_pressure_solver_iterations.png)
+## Output
 
-## Run the project
+Generated data are written to `results/data/`:
 
-Clone the repository, open MATLAB in the repository root, and run one of:
+- `study_summary_<mode>.csv`
+- `<case>_history.csv`
+- `<case>_fields.csv`
+- `<case>_centerlines.csv`
+- `<case>.mat`
 
-```matlab
-main_quick    % reduced check
-main_medium   % study without the N = 128 mesh
-main          % complete 72-case configuration
-```
+Figures are written to `results/figures/`. Large modes save study-level figures by default; the single case also saves per-case flow and validation figures.
 
-Linux shell wrappers are also available:
+## Numerical method
 
-```bash
-bash scripts/run_quick.sh
-bash scripts/run_medium.sh
-bash scripts/run.sh
-```
+The unit-square cavity has a lid velocity of `1`, stationary remaining walls, and viscosity `nu = 1/Re`. Pressure is stored at cell centers, `u` on vertical faces, and `v` on horizontal faces. Each pseudo-time iteration predicts face velocities, solves a pressure-correction Poisson equation, projects the velocities onto a divergence-free field, and evaluates strict convergence metrics.
 
-Generated files are written to `results/data/` and `results/figures/`. Detailed instructions are available in [`docs/RUNNING.md`](docs/RUNNING.md).
-
-## Repository structure
-
-```text
-config/        default solver and study settings
-startup/       path setup and output-folder creation
-core/          solver routines
-studies/       single-case and parameter-study runners
-validation/    Ghia data and error calculations
-post/          plotting and result export
-scripts/       shell wrappers for MATLAB runs
-assets/        selected figures and published summary data
-docs/          methodology, results, validation, and running notes
-results/       generated output; ignored by Git
-```
+See [`docs/METHODOLOGY.md`](docs/METHODOLOGY.md), [`docs/RUNNING.md`](docs/RUNNING.md), and [`docs/PHASE2_VERIFICATION.md`](docs/PHASE2_VERIFICATION.md).
 
 ## Requirements
 
-The project uses base MATLAB scripts and functions. No external MATLAB toolboxes are required for the main solver workflow.
-
-## Scope and limitations
-
-This is a completed educational solver and study, not a replacement for a production CFD package.
-
-Documented limitations include:
+- MATLAB with base language functionality
+- no external toolboxes required for the solver
+- Linux shell only for the optional wrapper scripts
 
-- collocated grid without Rhie-Chow interpolation
-- iterative pressure solver without multigrid acceleration
-- practical validation thresholds rather than a formal verification study
-- high-Reynolds-number cases that require stronger convergence control
-- no uncertainty quantification
+## Scope
 
-The code is kept as the completed MATLAB reference implementation for the broader work-in-progress multi-language comparison project.
+This is an educational and research comparison solver. It is designed for transparent numerical experiments and cross-language comparison, not as a replacement for a production CFD package.
 
 ## Reference
 
-Ghia, U., Ghia, K. N., & Shin, C. T. (1982). *High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method*. Journal of Computational Physics, 48(3), 387-411.
+Ghia, U., Ghia, K. N., & Shin, C. T. (1982). *High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method*. Journal of Computational Physics, 48(3), 387–411.
 
 ## Author
 

config/default_config.m

@@ -1,64 +1,66 @@
 function cfg = default_config()
-%DEFAULT_CONFIG Return the default solver and parameter-study settings.
+%DEFAULT_CONFIG Phase 2 production settings for the MATLAB cavity solver.
 
 cfg.U_lid = 1.0;
 cfg.L = 1.0;
 
-% Outer pressure-correction controls. The local maximum can be increased
-% automatically for the most demanding mesh, Reynolds number, and scheme.
-cfg.maxIter = 4000;
-cfg.maxIter_N128_bonus = 3000;
-cfg.maxIter_Re1000_bonus = 3000;
-cfg.maxIter_central_bonus = 1500;
-
-% Convergence criteria
-cfg.tol_mass = 1e-7;              % normalized mass imbalance
-cfg.tol_divergence = 2e-3;        % diagnostic only
-cfg.tol_velocity = 5e-7;
+% Strict outer convergence definition. A case is converged only when every
+% criterion passes for the requested number of consecutive iterations.
+cfg.maxIter = 30000;
+cfg.minimum_iterations = 200;
+cfg.consecutive_passes = 20;
+cfg.maximum_pressure_failures = 3;
+cfg.tol_velocity_linf = 1e-8;
+cfg.tol_divergence_linf = 1e-9;
+cfg.tol_divergence_l2 = 2e-10;
+cfg.tol_global_mass = 1e-12;
 cfg.diverged_limit = 1e6;
 
-% Pseudo-time controls
-cfg.cfl = 0.25;
-cfg.dt_max = 0.0025;
-cfg.dt_min = 1e-6;
-
-% Velocity and pressure relaxation
-cfg.alpha_u = 0.55;
-cfg.alpha_p = 0.20;
+% Pseudo-time and relaxation controls.
+cfg.cfl = 0.60;
+cfg.dt_max = 0.01;
+cfg.dt_min = 1e-8;
+cfg.alpha_u = 0.90;
+cfg.alpha_p = 1.00;
 
-% Pressure-Poisson controls
-cfg.poisson_maxIter = 2500;
-cfg.poisson_tol_abs = 1e-8;
-cfg.poisson_tol_rel = 1e-4;
-cfg.poisson_check_every = 25;
-
-% SOR controls. With 'auto', omega is estimated from the mesh and clipped
-% to the limits below.
+% Pressure-correction Poisson solver.
+cfg.poisson_maxIter = 5000;
+cfg.poisson_check_every = 20;
+cfg.poisson_tol_abs = 1e-10;
+cfg.poisson_tol_rel = 1e-9;
 cfg.sor_omega = 'auto';
-cfg.sor_omega_min = 1.15;
-cfg.sor_omega_max = 1.90;
-cfg.allow_pressure_maxIter = true;
+cfg.sor_omega_min = 1.0;
+cfg.sor_omega_max = 1.95;
+
+% Stagnation diagnostics.
+cfg.stagnation_window = 1500;
+cfg.stagnation_minimum_reduction = 0.005;
 
-% Full study: 3 meshes x 3 Reynolds numbers x 2 schemes x 2 pressure
-% solvers x 2 implementations = 72 simulations.
+% Study definitions. The production comparison uses the vectorized MATLAB
+% implementation; the full MATLAB study also compares the loop version.
 cfg.meshes = [32, 64, 128];
 cfg.re_list = [100, 400, 1000];
 cfg.schemes = {'upwind','central'};
 cfg.pressure_solvers = {'RBGS','RBSOR'};
 cfg.implementations = {'vectorized','loop'};
+cfg.use_continuation = true;
+cfg.strict = true;
+cfg.progress_every = 1000;
 
-% Practical validation thresholds against the Ghia centreline data.
+% Practical benchmark limits against the Ghia centerline values.
 cfg.validation_u_L2_limit_Re100 = 0.030;
 cfg.validation_v_L2_limit_Re100 = 0.030;
 cfg.validation_u_L2_limit_Re400 = 0.090;
 cfg.validation_v_L2_limit_Re400 = 0.120;
 cfg.validation_u_L2_limit_Re1000 = 0.160;
 cfg.validation_v_L2_limit_Re1000 = 0.180;
 
-% Output
+% Output controls. Per-case figures are disabled for large studies by
+% default because graphics can dominate batch runtime on COMPASS.
 cfg.make_figures = true;
-cfg.figure_every_case = true;
-cfg.results_dir = "results";
-cfg.data_dir = fullfile("results", "data");
-cfg.fig_dir = fullfile("results", "figures");
+cfg.figure_every_case = false;
+cfg.save_fields = true;
+cfg.results_dir = fullfile('results');
+cfg.data_dir = fullfile('results','data');
+cfg.fig_dir = fullfile('results','figures');
 end