diff --git a/README.md b/README.md index e9a6227..be40392 100644 --- a/README.md +++ b/README.md @@ -1,117 +1,89 @@ # Lid-Driven Cavity Flow Solver in MATLAB -

- Completed - MATLAB - MIT License - - Portfolio - -

+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_.csv` +- `_history.csv` +- `_fields.csv` +- `_centerlines.csv` +- `.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 diff --git a/config/default_config.m b/config/default_config.m index cef4d87..3345250 100644 --- a/config/default_config.m +++ b/config/default_config.m @@ -1,53 +1,53 @@ 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; @@ -55,10 +55,12 @@ 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 diff --git a/core/solve_lid_cavity.m b/core/solve_lid_cavity.m index 60a2701..ce44b0c 100644 --- a/core/solve_lid_cavity.m +++ b/core/solve_lid_cavity.m @@ -1,169 +1,548 @@ -function result = solve_lid_cavity(N,Re,scheme,pressure_solver,implementation,cfg) -%SOLVE_LID_CAVITY Solve one lid-driven cavity case with pressure correction. +function result = solve_lid_cavity(N,Re,scheme,pressure_solver,implementation,cfg,initial_state) +%SOLVE_LID_CAVITY Phase 2 staggered-grid projection solver. % -% The output includes the flow fields, residual histories, pressure-solver -% diagnostics, runtime, and stopping status. +% Pressure is stored at cell centers, horizontal velocity on vertical faces, +% and vertical velocity on horizontal faces. The optional initial_state is +% used for Reynolds-number continuation in parameter studies. + +if nargin < 7 + initial_state = struct(); +end -implementation = lower(string(implementation)); scheme = lower(string(scheme)); pressure_solver = upper(string(pressure_solver)); +implementation = lower(string(implementation)); -L = cfg.L; -dx = L/(N-1); -dy = dx; - -% Adaptive maximum iteration count. -localMaxIter = cfg.maxIter; -if N >= 128 - localMaxIter = localMaxIter + cfg.maxIter_N128_bonus; +if scheme ~= "upwind" && scheme ~= "central" + error('Unknown convection scheme: %s',scheme); end -if Re >= 1000 - localMaxIter = localMaxIter + cfg.maxIter_Re1000_bonus; +if pressure_solver ~= "RBGS" && pressure_solver ~= "RBSOR" + error('Unknown pressure solver: %s',pressure_solver); end -if scheme == "central" - localMaxIter = localMaxIter + cfg.maxIter_central_bonus; +if implementation ~= "vectorized" && implementation ~= "loop" + error('Unknown MATLAB implementation: %s',implementation); end -u = zeros(N,N); -v = zeros(N,N); -p = zeros(N,N); -[u,v] = apply_lid_bc(u,v,cfg.U_lid); +h = cfg.L/N; +nu = cfg.U_lid*cfg.L/Re; +maxIter = cfg.maxIter; -Ru = zeros(localMaxIter,1); -Rv = zeros(localMaxIter,1); -Rc_mass = zeros(localMaxIter,1); -Rc_div = zeros(localMaxIter,1); -dt_hist = zeros(localMaxIter,1); -poisson_iters = zeros(localMaxIter,1); -poisson_rel = zeros(localMaxIter,1); -poisson_conv = false(localMaxIter,1); +[u_face,v_face,p] = initialize_state(N,initial_state); +[u_face,v_face] = apply_normal_velocity_bc(u_face,v_face); -tic_total = tic; +velocity_hist = nan(maxIter,1); +div_linf_hist = nan(maxIter,1); +div_l2_hist = nan(maxIter,1); +mass_hist = nan(maxIter,1); +dt_hist = nan(maxIter,1); +poisson_rel_hist = nan(maxIter,1); +poisson_iter_hist = nan(maxIter,1); +poisson_conv_hist = false(maxIter,1); -status = "maxIter"; -stagnation_counter = 0; -prev_mass = inf; +status = "max_iterations"; +consecutive_pass_count = 0; +failed_pressure_solves = 0; -for iter = 1:localMaxIter +start_time = tic; - u_old = u; - v_old = v; +for iter = 1:maxIter + u_old = u_face; + v_old = v_face; + + dt = compute_time_step(u_face,v_face,h,nu,cfg); if implementation == "vectorized" - [u_star,v_star,dt] = momentum_predictor_vectorized(u,v,p,Re,scheme,cfg); - elseif implementation == "loop" - [u_star,v_star,dt] = momentum_predictor_loop(u,v,p,Re,scheme,cfg); + [u_star,v_star] = predict_velocity_vectorized( ... + u_face,v_face,p,h,nu,dt,scheme,cfg); else - error("Unknown implementation: %s", implementation); + [u_star,v_star] = predict_velocity_loop( ... + u_face,v_face,p,h,nu,dt,scheme,cfg); end - dt = max(dt, cfg.dt_min); - - % Pressure correction RHS: div(u*) / dt - div_star = divergence_field(u_star,v_star,dx,dy); - rhs = div_star / dt; + [u_star,v_star] = apply_normal_velocity_bc(u_star,v_star); + div_star = discrete_divergence(u_star,v_star,h); + rhs = div_star/dt; - [p_prime, pinfo] = pressure_poisson(rhs,dx,dy,pressure_solver,cfg); + [p_prime,pinfo] = pressure_poisson_phase2( ... + rhs,h,pressure_solver,cfg); - % Velocity correction, vectorized for both implementations. - u = u_star; - v = v_star; + u_face = u_star; + v_face = v_star; + u_face(:,2:N) = u_star(:,2:N) ... + - dt*(p_prime(:,2:N)-p_prime(:,1:N-1))/h; + v_face(2:N,:) = v_star(2:N,:) ... + - dt*(p_prime(2:N,:)-p_prime(1:N-1,:))/h; + [u_face,v_face] = apply_normal_velocity_bc(u_face,v_face); - C = 2:N-1; - dpdx = (p_prime(C,C+1) - p_prime(C,C-1))/(2*dx); - dpdy = (p_prime(C+1,C) - p_prime(C-1,C))/(2*dy); + p = p + cfg.alpha_p*p_prime; + p = p - mean(p(:)); - u(C,C) = u_star(C,C) - dt*dpdx; - v(C,C) = v_star(C,C) - dt*dpdy; + residuals = calculate_residuals( ... + u_face,v_face,u_old,v_old,h,cfg); - p = p + cfg.alpha_p*p_prime; - p = p - mean(p(:)); % remove arbitrary pressure offset + velocity_hist(iter) = residuals.velocity_linf; + div_linf_hist(iter) = residuals.divergence_linf; + div_l2_hist(iter) = residuals.divergence_l2; + mass_hist(iter) = residuals.global_mass; + dt_hist(iter) = dt; + poisson_rel_hist(iter) = pinfo.relative_residual; + poisson_iter_hist(iter) = pinfo.iterations; + poisson_conv_hist(iter) = pinfo.converged; - [u,v] = apply_lid_bc(u,v,cfg.U_lid); + if ~pinfo.converged + failed_pressure_solves = failed_pressure_solves + 1; + else + failed_pressure_solves = 0; + end - [Ru(iter),Rv(iter),Rc_mass(iter),Rc_div(iter)] = ... - velocity_residuals(u,v,u_old,v_old,dx,dy,cfg.U_lid,cfg.L); + if cfg.progress_every > 0 && mod(iter,cfg.progress_every) == 0 + fprintf(' iter=%d vel=%.6e div=%.6e p=%.6e\n', ... + iter,residuals.velocity_linf,residuals.divergence_linf, ... + pinfo.relative_residual); + end - dt_hist(iter) = dt; - poisson_iters(iter) = pinfo.iter; - poisson_rel(iter) = pinfo.final_relative_residual; - poisson_conv(iter) = pinfo.converged; + if any(~isfinite(u_face(:))) || any(~isfinite(v_face(:))) ... + || any(~isfinite(p(:))) + status = "non_finite"; + break; + end - if any(~isfinite(u(:))) || any(~isfinite(v(:))) || any(~isfinite(p(:))) || ... - max([Ru(iter),Rv(iter),Rc_div(iter)]) > cfg.diverged_limit + if max([residuals.velocity_linf,residuals.divergence_linf, ... + residuals.divergence_l2,pinfo.relative_residual]) ... + > cfg.diverged_limit status = "diverged"; break; end - % Track slow residual reduction as a diagnostic; do not stop early. - if Rc_mass(iter) > 0.995*prev_mass - stagnation_counter = stagnation_counter + 1; + if failed_pressure_solves > cfg.maximum_pressure_failures + status = "pressure_not_converged"; + break; + end + + all_pass = pinfo.converged ... + && residuals.velocity_linf <= cfg.tol_velocity_linf ... + && residuals.divergence_linf <= cfg.tol_divergence_linf ... + && residuals.divergence_l2 <= cfg.tol_divergence_l2 ... + && residuals.global_mass <= cfg.tol_global_mass; + + if iter >= cfg.minimum_iterations && all_pass + consecutive_pass_count = consecutive_pass_count + 1; else - stagnation_counter = 0; + consecutive_pass_count = 0; end - prev_mass = Rc_mass(iter); - % Use the mass and velocity residuals for the outer stopping check. - % Raw divergence is stored as an additional diagnostic. - if Rc_mass(iter) < cfg.tol_mass && max(Ru(iter),Rv(iter)) < cfg.tol_velocity + if consecutive_pass_count >= cfg.consecutive_passes status = "converged"; break; end + + if iter >= cfg.stagnation_window + old_index = iter-cfg.stagnation_window+1; + old_value = velocity_hist(old_index); + if isfinite(old_value) && old_value > 0 + reduction = (old_value-residuals.velocity_linf)/old_value; + if reduction < cfg.stagnation_minimum_reduction ... + && residuals.velocity_linf > 10*cfg.tol_velocity_linf + status = "stagnated"; + break; + end + end + end end -runtime = toc(tic_total); +runtime = toc(start_time); +last = iter; -Ru = Ru(1:iter); -Rv = Rv(1:iter); -Rc_mass = Rc_mass(1:iter); -Rc_div = Rc_div(1:iter); -dt_hist = dt_hist(1:iter); -poisson_iters = poisson_iters(1:iter); -poisson_rel = poisson_rel(1:iter); -poisson_conv = poisson_conv(1:iter); +velocity_hist = velocity_hist(1:last); +div_linf_hist = div_linf_hist(1:last); +div_l2_hist = div_l2_hist(1:last); +mass_hist = mass_hist(1:last); +dt_hist = dt_hist(1:last); +poisson_rel_hist = poisson_rel_hist(1:last); +poisson_iter_hist = poisson_iter_hist(1:last); +poisson_conv_hist = poisson_conv_hist(1:last); -x = linspace(0,L,N); -y = linspace(0,L,N); -[X,Y] = meshgrid(x,y); +[u_center,v_center,x,y] = face_to_cell_center(u_face,v_face,cfg.L); +speed = hypot(u_center,v_center); +vorticity = centered_vorticity(u_center,v_center,h); result.N = N; result.Re = Re; result.scheme = char(scheme); result.pressure_solver = char(pressure_solver); result.implementation = char(implementation); +result.grid_type = 'staggered_MAC'; result.x = x; result.y = y; -result.X = X; -result.Y = Y; - -result.u = u; -result.v = v; +result.u = u_center; +result.v = v_center; result.p = p; -result.speed = sqrt(u.^2 + v.^2); -result.vorticity = compute_vorticity(u,v,dx,dy); - -result.Ru = Ru; -result.Rv = Rv; -result.Rc = Rc_mass; % backward compatible name -result.Rc_mass = Rc_mass; -result.Rc_div = Rc_div; +result.speed = speed; +result.vorticity = vorticity; +result.u_face = u_face; +result.v_face = v_face; + +% Backward-compatible fields used by the existing plotting functions. +result.Ru = velocity_hist; +result.Rv = velocity_hist; +result.Rc = mass_hist; +result.Rc_mass = mass_hist; +result.Rc_div = div_linf_hist; + +result.velocity_update_linf = velocity_hist; +result.divergence_linf = div_linf_hist; +result.divergence_l2 = div_l2_hist; +result.global_mass_imbalance = mass_hist; result.dt = dt_hist; -result.poisson_iters = poisson_iters; -result.poisson_relative_residual = poisson_rel; -result.poisson_converged = poisson_conv; +result.poisson_relative_residual = poisson_rel_hist; +result.poisson_iters = poisson_iter_hist; +result.poisson_converged = poisson_conv_hist; -result.iterations = iter; -result.localMaxIter = localMaxIter; +result.iterations = last; +result.localMaxIter = maxIter; result.runtime = runtime; result.status = char(status); -result.final_Ru = Ru(end); -result.final_Rv = Rv(end); -result.final_Rc = Rc_mass(end); -result.final_Rc_mass = Rc_mass(end); -result.final_Rc_div = Rc_div(end); -result.avg_poisson_iters = mean(poisson_iters); -result.avg_poisson_relative_residual = mean(poisson_rel); -result.pressure_saturation_ratio = mean(poisson_iters >= cfg.poisson_maxIter); -result.stagnation_counter = stagnation_counter; +result.consecutive_pass_count = consecutive_pass_count; +result.failed_pressure_solves = failed_pressure_solves; + +result.final_Ru = velocity_hist(end); +result.final_Rv = velocity_hist(end); +result.final_Rc = mass_hist(end); +result.final_Rc_mass = mass_hist(end); +result.final_Rc_div = div_linf_hist(end); +result.final_velocity_linf = velocity_hist(end); +result.final_divergence_linf = div_linf_hist(end); +result.final_divergence_l2 = div_l2_hist(end); +result.final_global_mass = mass_hist(end); +result.final_poisson_relative_residual = poisson_rel_hist(end); +result.avg_poisson_iters = mean(poisson_iter_hist); +result.avg_poisson_relative_residual = mean(poisson_rel_hist); +result.pressure_saturation_ratio = mean(poisson_iter_hist >= cfg.poisson_maxIter); + +result.continuation_state.available = strcmp(result.status,'converged'); +result.continuation_state.u_face = u_face; +result.continuation_state.v_face = v_face; +result.continuation_state.p = p; +end + +function [u,v,p] = initialize_state(N,initial_state) +u = zeros(N,N+1); +v = zeros(N+1,N); +p = zeros(N,N); + +if ~isstruct(initial_state) || ~isfield(initial_state,'available') ... + || ~initial_state.available + return; +end + +if isfield(initial_state,'u_face') && isequal(size(initial_state.u_face),size(u)) + u = initial_state.u_face; +end +if isfield(initial_state,'v_face') && isequal(size(initial_state.v_face),size(v)) + v = initial_state.v_face; +end +if isfield(initial_state,'p') && isequal(size(initial_state.p),size(p)) + p = initial_state.p; + p = p-mean(p(:)); +end +end + +function [u,v] = apply_normal_velocity_bc(u,v) +u(:,1) = 0; +u(:,end) = 0; +v(1,:) = 0; +v(end,:) = 0; +end + +function dt = compute_time_step(u,v,h,nu,cfg) +max_velocity = max([max(abs(u(:))),max(abs(v(:))),cfg.U_lid,1e-12]); +convection_limit = cfg.cfl*h/max_velocity; +diffusion_limit = 0.24*h*h/max(nu,1e-30); +dt = min([convection_limit,diffusion_limit,cfg.dt_max]); +dt = max(cfg.dt_min,min(cfg.dt_max,dt)); +end + +function [u_star,v_star] = predict_velocity_vectorized( ... + u,v,p,h,nu,dt,scheme,cfg) +N = size(p,1); +u_star = u; +v_star = v; + +% u momentum on vertical faces j=2:N. +u_ext = zeros(N+2,N+1); +u_ext(2:N+1,:) = u; +u_ext(1,:) = -u(1,:); +u_ext(N+2,:) = 2*cfg.U_lid-u(N,:); + +uC = u(:,2:N); +uW = u(:,1:N-1); +uE = u(:,3:N+1); +uS = u_ext(1:N,2:N); +uN = u_ext(3:N+2,2:N); +vAtU = 0.25*(v(1:N,1:N-1)+v(2:N+1,1:N-1) ... + +v(1:N,2:N)+v(2:N+1,2:N)); + +if scheme == "central" + du_dx = (uE-uW)/(2*h); + du_dy = (uN-uS)/(2*h); +else + du_dx = upwind_derivative(uC,uW,uE,uC,h); + du_dy = upwind_derivative(uC,uS,uN,vAtU,h); +end + +lap_u = (uE-2*uC+uW+uN-2*uC+uS)/(h*h); +dp_dx = (p(:,2:N)-p(:,1:N-1))/h; +u_pred = uC+cfg.alpha_u*dt*( ... + -uC.*du_dx-vAtU.*du_dy-dp_dx+nu*lap_u); +u_star(:,2:N) = u_pred; + +% v momentum on horizontal faces i=2:N. +v_ext = zeros(N+1,N+2); +v_ext(:,2:N+1) = v; +v_ext(:,1) = -v(:,1); +v_ext(:,N+2) = -v(:,N); + +vC = v(2:N,:); +vS = v(1:N-1,:); +vN = v(3:N+1,:); +vW = v_ext(2:N,1:N); +vE = v_ext(2:N,3:N+2); +uAtV = 0.25*(u(1:N-1,1:N)+u(1:N-1,2:N+1) ... + +u(2:N,1:N)+u(2:N,2:N+1)); + +if scheme == "central" + dv_dx = (vE-vW)/(2*h); + dv_dy = (vN-vS)/(2*h); +else + dv_dx = upwind_derivative(vC,vW,vE,uAtV,h); + dv_dy = upwind_derivative(vC,vS,vN,vC,h); +end + +lap_v = (vE-2*vC+vW+vN-2*vC+vS)/(h*h); +dp_dy = (p(2:N,:)-p(1:N-1,:))/h; +v_pred = vC+cfg.alpha_u*dt*( ... + -uAtV.*dv_dx-vC.*dv_dy-dp_dy+nu*lap_v); +v_star(2:N,:) = v_pred; +end + +function derivative = upwind_derivative(center,minus,plus,transport,h) +derivative = zeros(size(center)); +pos = transport >= 0; +derivative(pos) = (center(pos)-minus(pos))/h; +derivative(~pos) = (plus(~pos)-center(~pos))/h; +end + +function [u_star,v_star] = predict_velocity_loop( ... + u,v,p,h,nu,dt,scheme,cfg) +N = size(p,1); +u_star = u; +v_star = v; + +for i = 1:N + for j = 2:N + uc = u(i,j); + vw = v(i,j-1); + ve = v(i,j); + vn_w = v(i+1,j-1); + vn_e = v(i+1,j); + v_at_u = 0.25*(vw+ve+vn_w+vn_e); + + u_w = u(i,j-1); + u_e = u(i,j+1); + if i == 1 + u_s = -u(i,j); + else + u_s = u(i-1,j); + end + if i == N + u_n = 2*cfg.U_lid-u(i,j); + else + u_n = u(i+1,j); + end + + if scheme == "central" + du_dx = (u_e-u_w)/(2*h); + du_dy = (u_n-u_s)/(2*h); + else + if uc >= 0 + du_dx = (uc-u_w)/h; + else + du_dx = (u_e-uc)/h; + end + if v_at_u >= 0 + du_dy = (uc-u_s)/h; + else + du_dy = (u_n-uc)/h; + end + end + + lap = (u_e-2*uc+u_w+u_n-2*uc+u_s)/(h*h); + dpdx = (p(i,j)-p(i,j-1))/h; + u_star(i,j) = uc+cfg.alpha_u*dt*( ... + -uc*du_dx-v_at_u*du_dy-dpdx+nu*lap); + end +end + +for i = 2:N + for j = 1:N + vc = v(i,j); + u_at_v = 0.25*(u(i-1,j)+u(i-1,j+1) ... + +u(i,j)+u(i,j+1)); + + v_s = v(i-1,j); + v_n = v(i+1,j); + if j == 1 + v_w = -v(i,j); + else + v_w = v(i,j-1); + end + if j == N + v_e = -v(i,j); + else + v_e = v(i,j+1); + end + + if scheme == "central" + dv_dx = (v_e-v_w)/(2*h); + dv_dy = (v_n-v_s)/(2*h); + else + if u_at_v >= 0 + dv_dx = (vc-v_w)/h; + else + dv_dx = (v_e-vc)/h; + end + if vc >= 0 + dv_dy = (vc-v_s)/h; + else + dv_dy = (v_n-vc)/h; + end + end + + lap = (v_e-2*vc+v_w+v_n-2*vc+v_s)/(h*h); + dpdy = (p(i,j)-p(i-1,j))/h; + v_star(i,j) = vc+cfg.alpha_u*dt*( ... + -u_at_v*dv_dx-vc*dv_dy-dpdy+nu*lap); + end +end +end + +function div = discrete_divergence(u,v,h) +div = (u(:,2:end)-u(:,1:end-1))/h ... + +(v(2:end,:)-v(1:end-1,:))/h; +end + +function [pressure,info] = pressure_poisson_phase2(rhs,h,method,cfg) +N = size(rhs,1); +rhs = rhs-mean(rhs(:)); +rhs_norm = max(max(abs(rhs(:))),1e-30); +pressure = zeros(N,N); + +if method == "RBSOR" + if isnumeric(cfg.sor_omega) + omega = cfg.sor_omega; + else + omega = 2/(1+sin(pi/N)); + omega = min(cfg.sor_omega_max,max(cfg.sor_omega_min,omega)); + end +else + omega = 1.0; +end + +[row_index,col_index] = ndgrid(1:N,1:N); +red = mod(row_index+col_index,2) == 0; +black = ~red; + +info.iterations = 0; +info.converged = false; +info.absolute_residual = inf; +info.relative_residual = inf; +info.omega = omega; + +for k = 1:cfg.poisson_maxIter + pressure = red_black_sweep(pressure,rhs,h,red,omega); + pressure = red_black_sweep(pressure,rhs,h,black,omega); + + if mod(k,50) == 0 + pressure = pressure-mean(pressure(:)); + end + + if k == 1 || mod(k,cfg.poisson_check_every) == 0 ... + || k == cfg.poisson_maxIter + absolute_residual = poisson_residual_linf(pressure,rhs,h); + relative_residual = absolute_residual/rhs_norm; + info.iterations = k; + info.absolute_residual = absolute_residual; + info.relative_residual = relative_residual; + + if absolute_residual <= cfg.poisson_tol_abs ... + || relative_residual <= cfg.poisson_tol_rel + info.converged = true; + pressure = pressure-mean(pressure(:)); + return; + end + end +end + +pressure = pressure-mean(pressure(:)); +end + +function pressure = red_black_sweep(pressure,rhs,h,mask,omega) +[neighbor_sum,neighbor_count] = poisson_neighbors(pressure); +candidate = (neighbor_sum-rhs*h*h)./neighbor_count; +pressure(mask) = (1-omega)*pressure(mask)+omega*candidate(mask); +end + +function [neighbor_sum,neighbor_count] = poisson_neighbors(pressure) +N = size(pressure,1); +neighbor_sum = zeros(N,N); +neighbor_count = zeros(N,N); + +neighbor_sum(2:N,:) = neighbor_sum(2:N,:)+pressure(1:N-1,:); +neighbor_count(2:N,:) = neighbor_count(2:N,:)+1; +neighbor_sum(1:N-1,:) = neighbor_sum(1:N-1,:)+pressure(2:N,:); +neighbor_count(1:N-1,:) = neighbor_count(1:N-1,:)+1; +neighbor_sum(:,2:N) = neighbor_sum(:,2:N)+pressure(:,1:N-1); +neighbor_count(:,2:N) = neighbor_count(:,2:N)+1; +neighbor_sum(:,1:N-1) = neighbor_sum(:,1:N-1)+pressure(:,2:N); +neighbor_count(:,1:N-1) = neighbor_count(:,1:N-1)+1; +end + +function value = poisson_residual_linf(pressure,rhs,h) +[neighbor_sum,neighbor_count] = poisson_neighbors(pressure); +laplacian = (neighbor_sum-neighbor_count.*pressure)/(h*h); +value = max(abs(laplacian(:)-rhs(:))); +end + +function residuals = calculate_residuals(u,v,u_old,v_old,h,cfg) +velocity_change = max([max(abs(u(:)-u_old(:))), ... + max(abs(v(:)-v_old(:)))])/cfg.U_lid; +div = discrete_divergence(u,v,h)/(cfg.U_lid/cfg.L); + +residuals.velocity_linf = velocity_change; +residuals.divergence_linf = max(abs(div(:))); +residuals.divergence_l2 = sqrt(mean(div(:).^2)); + +boundary_flux = sum(u(:,end)-u(:,1))*h ... + +sum(v(end,:)-v(1,:))*h; +residuals.global_mass = abs(boundary_flux)/(cfg.U_lid*cfg.L); +end + +function [u_center,v_center,x,y] = face_to_cell_center(u,v,L) +N = size(u,1); +h = L/N; +u_center = 0.5*(u(:,1:N)+u(:,2:N+1)); +v_center = 0.5*(v(1:N,:)+v(2:N+1,:)); +x = ((1:N)-0.5)*h; +y = ((1:N)-0.5)*h; +end + +function omega = centered_vorticity(u,v,h) +N = size(u,1); +omega = zeros(N,N); +if N > 2 + omega(2:N-1,2:N-1) = ... + (v(2:N-1,3:N)-v(2:N-1,1:N-2))/(2*h) ... + -(u(3:N,2:N-1)-u(1:N-2,2:N-1))/(2*h); +end end diff --git a/docs/METHODOLOGY.md b/docs/METHODOLOGY.md index 768e19c..a2c7f15 100644 --- a/docs/METHODOLOGY.md +++ b/docs/METHODOLOGY.md @@ -1,83 +1,84 @@ -# Methodology +# Numerical Methodology -## Problem setup +## Problem definition -The domain is a unit square filled with an incompressible fluid. The top wall moves from left to right with a nondimensional speed of `1`; the other walls remain stationary. No-slip and no-penetration conditions are applied at all walls. - -The solver uses the nondimensional equations: +The solver models the nondimensional two-dimensional incompressible lid-driven cavity in a unit square. The top wall moves with velocity `U = 1`; the other walls remain stationary. The kinematic viscosity is ```text -Continuity: ∇ · u = 0 -Momentum equation: ∂u/∂t + (u · ∇)u = −∇p + (1/Re)∇²u +nu = U L / Re ``` -The Reynolds number is based on lid velocity and cavity length. - -## Grid and spatial discretization - -Velocity and pressure are stored on a structured Cartesian collocated grid. Spatial derivatives are calculated with finite differences. - -Two convection schemes are available: +with `L = 1`. -- **Upwind:** more stable on coarse meshes, but more diffusive -- **Central:** less diffusive, but more sensitive to mesh resolution +## Staggered MAC grid -Diffusion terms and pressure gradients use central differences. +The Phase 2 production solver uses a Marker-and-Cell arrangement: -## Pressure-correction loop +- pressure at cell centers +- horizontal velocity on vertical cell faces +- vertical velocity on horizontal cell faces -The steady solution is approached through pseudo-time stepping. Each outer iteration: +This removes the checkerboard pressure mode associated with the earlier collocated prototype and makes the pressure gradient, velocity correction, and discrete divergence compatible. -1. applies velocity boundary conditions -2. predicts the intermediate velocity field -3. calculates its divergence -4. solves the pressure-correction Poisson equation -5. corrects velocity and pressure -6. records residuals -7. continues until the stopping criteria or iteration limit is reached +## Projection workflow -The momentum predictor is implemented twice: +Each pseudo-time iteration: -- `momentum_predictor_loop.m` follows the equations cell by cell -- `momentum_predictor_vectorized.m` performs equivalent operations with MATLAB arrays +1. calculates a stable time step from convection and diffusion limits +2. predicts face velocities from convection, diffusion, and pressure gradients +3. calculates cell-centered divergence of the predicted field +4. removes the mean from the Poisson right-hand side +5. solves the pressure-correction Poisson equation +6. corrects face velocities with pressure-correction gradients +7. updates and normalizes pressure +8. evaluates all convergence metrics -Keeping both versions supports runtime comparison while checking numerical consistency. +## Momentum discretization -## Pressure solvers +The code supports first-order upwind and second-order central convection. Diffusion uses second-order central differences. Normal velocities are set directly to zero at solid boundary faces. Tangential no-slip conditions are imposed through ghost values; the upper-wall horizontal-velocity ghost value enforces the moving lid. -The pressure-correction equation can be solved with: +Two MATLAB implementations are available: -- `RBGS`: red-black Gauss-Seidel -- `RBSOR`: red-black Successive Over-Relaxation +- `vectorized`: array-based momentum prediction +- `loop`: explicit cell-by-cell momentum prediction -For RBSOR, the relaxation factor is estimated from grid size and limited by the bounds in `default_config.m`. The pressure solver stops using the relative residual of the Poisson equation. +They share the same pressure solver, projection, convergence checks, and output path. -## Time step and relaxation +## Pressure correction -The pseudo-time step is limited by convection, diffusion, and the configured maximum value. Velocity and pressure under-relaxation factors are also defined in `default_config.m`. +The pressure-correction equation is solved using red-black Gauss-Seidel (`RBGS`) or red-black successive over-relaxation (`RBSOR`). Missing neighbors at boundaries implement homogeneous normal pressure-gradient conditions through the boundary stencil. The pressure correction and right-hand side are normalized to remove the constant null space. -These settings were selected to keep the complete parameter study stable; they are not necessarily optimal for every individual case. +Pressure convergence uses the true discrete equation residual. An outer case cannot report convergence while the pressure solve is failing. -## Residuals +## Strict convergence definition The solver records: -- changes in horizontal and vertical velocity -- normalized mass imbalance -- raw velocity divergence -- pressure-solver iterations and relative residuals +- velocity-update `Linf` +- divergence `Linf` +- divergence `L2` +- global boundary mass imbalance +- pressure-Poisson relative residual -Velocity changes and mass imbalance are used for the outer stopping check. Raw divergence is retained as an additional diagnostic. +A case reports `converged` only when every configured criterion passes for a required number of consecutive iterations and the minimum iteration count has been reached. -## Validation +Terminal states are: + +```text +converged +max_iterations +pressure_not_converged +stagnated +diverged +non_finite +``` -For `Re = 100`, `400`, and `1000`, the computed centerline velocities are interpolated to the sample locations reported by Ghia et al. The code calculates `L2` and maximum errors for: +## Continuation -- `u(y)` at `x = 0.5` -- `v(x)` at `y = 0.5` +Parameter studies reuse a converged lower-Reynolds-number solution when mesh, convection scheme, pressure solver, and MATLAB implementation remain unchanged. This improves stability and reduces startup cost for `Re = 400` and `Re = 1000`. -The limits used to classify cases are stored in `default_config.m` and discussed in [`VALIDATION.md`](VALIDATION.md). +## Benchmark comparison -## Scope +Cell-centered velocities are interpolated to `x = 0.5` and `y = 0.5`, then interpolated again to the Ghia sample coordinates. The code reports `L2` and `Linf` errors for `u(y)` and `v(x)`. -The repository records a completed educational solver and parameter study. It does not include Rhie-Chow interpolation, multigrid acceleration, turbulence modeling, or adaptive meshing. +This is a numerical benchmark comparison, not experimental validation. diff --git a/docs/PHASE2_VERIFICATION.md b/docs/PHASE2_VERIFICATION.md new file mode 100644 index 0000000..8add546 --- /dev/null +++ b/docs/PHASE2_VERIFICATION.md @@ -0,0 +1,48 @@ +# Phase 2 Verification Plan + +## Canonical regression + +The required first check is: + +```text +N = 32 +Re = 100 +scheme = upwind +pressure solver = RBSOR +MATLAB implementation = vectorized +``` + +Run: + +```bash +bash scripts/run_tests.sh +``` + +The test fails unless the case converges, all strict residual limits pass, and the Ghia benchmark limits pass. + +## Production verification modes + +After the regression succeeds: + +1. `medium`: six `N = 32` cases across `Re = 100, 400, 1000` and both convection schemes +2. `grid`: `N = 16, 32, 64` at `Re = 100` +3. `re1000`: representative `N = 64`, `Re = 1000` case +4. `full`: complete 72-case MATLAB comparison + +## Required evidence + +For each case, inspect: + +- terminal status +- velocity-update `Linf` +- divergence `Linf` and `L2` +- global mass imbalance +- pressure residual and pressure iteration count +- Ghia `L2` and `Linf` errors +- runtime + +The standardized CSV outputs allow direct comparison with the C++ Phase 2 results. + +## Current execution status + +The branch implementation was prepared without access to a MATLAB or Octave runtime in the editing environment. The COMPASS regression run is therefore the required numerical acceptance step before merging the Phase 2 branch into `main`. diff --git a/docs/RUNNING.md b/docs/RUNNING.md index ba2fb3d..1fb4618 100644 --- a/docs/RUNNING.md +++ b/docs/RUNNING.md @@ -1,110 +1,56 @@ -# Running the code +# Running the MATLAB Phase 2 Solver -Run the scripts from the repository root so MATLAB can add the required folders correctly. +## MATLAB command window -Generated data and figures are written to: - -```text -results/data/ -results/figures/ -``` - -The folders remain in the repository, but their generated contents are ignored by Git. - -## Run modes - -### Quick check +Start MATLAB in the repository root and run one mode: ```matlab -main_quick +run_mode('single') +run_mode('quick') +run_mode('medium') +run_mode('grid') +run_mode('re1000') +run_mode('full') ``` -Use this after cloning the repository or changing a solver function. It runs a reduced study while comparing both convection schemes, pressure solvers, and implementations. - -### Intermediate study - -```matlab -main_medium -``` +The equivalent script entry points are `main_single`, `main_quick`, `main_medium`, `main_grid`, `main_re1000`, and `main`. -This runs more cases without the `N = 128` mesh. - -### Complete configured study - -```matlab -main -``` - -This runs all 72 configured combinations. It can take a long time because some cases require thousands of outer iterations and many pressure iterations. - -Linux shell wrappers are included in `scripts/`: +## Linux batch execution ```bash +bash scripts/run_single.sh bash scripts/run_quick.sh bash scripts/run_medium.sh -bash scripts/run.sh +bash scripts/run_grid.sh +bash scripts/run_re1000.sh +bash scripts/run_full.sh ``` -They require the `matlab` command to be available in the shell path. - -## Running one editable case - -Open and run: +The wrapper uses `matlab -batch` and writes a timestamped log in `logs/`. When MATLAB is not on `PATH`, set its executable explicitly: -```matlab -studies/run_single_case.m +```bash +MATLAB_BIN=/path/to/matlab bash scripts/run_single.sh ``` -The mesh, Reynolds number, convection scheme, pressure solver, and implementation are defined near the top of the file. - -`studies/run_representative_case.m` runs the `N = 64`, `Re = 100`, central-differencing case shown in the README. +## Recommended COMPASS sequence -## Changing the setup +Run the modes sequentially: -Most settings are collected in: - -```text -config/default_config.m +```bash +bash scripts/run_tests.sh +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 ``` -This includes: - -- meshes and Reynolds numbers -- convergence tolerances -- relaxation factors -- pressure-solver limits -- validation limits -- figure-output settings - -The quick and medium scripts override some values to reduce runtime. - -## Output files - -The parameter study creates: - -- one `.mat` file per case -- CSV and MAT summary tables -- residual and pressure-residual histories -- velocity, pressure, vorticity, streamline, and vector plots -- Ghia centerline-comparison plots -- study-level runtime and error comparisons - -Selected figures and the published summary are copied to `assets/` so they remain visible on GitHub. - -## Common issues - -### MATLAB cannot find a function - -Start MATLAB in the repository root and run one of the main scripts. The scripts call `startup/setup_project.m`, which adds the required folders to the path. - -### A case reaches `maxIter` - -The outer stopping criteria were not met before the configured limit. Check the residual history, mass imbalance, validation error, and flow plots together. A low validation error alone does not prove numerical convergence. +Do not start multiple modes in the same repository simultaneously because they share `results/data/` and `results/figures/`. -### The complete study takes too long +## Outputs -Start with `main_quick`. The `N = 128` and RBGS cases are the most expensive in the included setup. +Each mode writes `study_summary_.csv`. Each case writes convergence history, flattened fields, centerline profiles, and a MATLAB result file. Strict mode raises a MATLAB error when any requested case does not reach the numerical convergence definition. -### Generated figures do not appear on GitHub +## Runtime note -The contents of `results/` are ignored by Git. Copy only figures intended for publication into `assets/figures/` and reference those files in the README. +The full 72-case study includes the loop implementation and the `N = 128` mesh. It can be substantially slower than the production verification modes. Use `medium` and `grid` first to verify the environment and solver behavior. diff --git a/main.m b/main.m index c9e35e3..b4c4c97 100644 --- a/main.m +++ b/main.m @@ -1,36 +1,3 @@ -% MAIN Run the full 72-case MATLAB parameter study. - +% MAIN Run the complete 72-case MATLAB Phase 2 study. clear; clc; close all; -addpath("startup"); -setup_project(); - -cfg = default_config(); - -fprintf("\nFULL STUDY SELECTED\n"); -fprintf("Meshes: %s\n", mat2str(cfg.meshes)); -fprintf("Reynolds numbers: %s\n", mat2str(cfg.re_list)); -fprintf("Schemes: %s\n", strjoin(string(cfg.schemes), ", ")); -fprintf("Pressure solvers: %s\n", strjoin(string(cfg.pressure_solvers), ", ")); -fprintf("Implementations: %s\n", strjoin(string(cfg.implementations), ", ")); - -nCases = numel(cfg.meshes) * numel(cfg.re_list) * numel(cfg.schemes) * ... - numel(cfg.pressure_solvers) * numel(cfg.implementations); - -fprintf("Total simulations: %d\n", nCases); -fprintf("Base max outer iterations per simulation: %d\n", cfg.maxIter); -fprintf("N=128 bonus iterations: %d\n", cfg.maxIter_N128_bonus); -fprintf("Re=1000 bonus iterations: %d\n", cfg.maxIter_Re1000_bonus); -fprintf("Central scheme bonus iterations: %d\n", cfg.maxIter_central_bonus); -fprintf("Max pressure iterations per outer iteration: %d\n\n", cfg.poisson_maxIter); - -T = run_parametric_study(cfg); - -disp(" "); -disp("Finished. Summary table:"); -disp(T); - -writetable(T, fullfile("results", "data", "study_summary.csv")); -plot_study_summary(T, cfg); - -disp(" "); -disp("Results saved in results/data and results/figures."); +run_mode('full'); diff --git a/main_grid.m b/main_grid.m new file mode 100644 index 0000000..b9934af --- /dev/null +++ b/main_grid.m @@ -0,0 +1,3 @@ +% MAIN_GRID Run the N=16/32/64 grid sequence at Re=100. +clear; clc; close all; +run_mode('grid'); diff --git a/main_medium.m b/main_medium.m index bfb2365..307251d 100644 --- a/main_medium.m +++ b/main_medium.m @@ -1,26 +1,3 @@ -% MAIN_MEDIUM Run the MATLAB study without the N = 128 mesh. -% 48 simulations: [32,64] x [100,400,1000] x 2 schemes x 2 pressure solvers x 2 implementations. - +% MAIN_MEDIUM Run the six production verification cases. clear; clc; close all; -addpath("startup"); -setup_project(); - -cfg = default_config(); -cfg.meshes = [32, 64]; -cfg.re_list = [100, 400, 1000]; -cfg.maxIter = 3500; -cfg.maxIter_N128_bonus = 0; -cfg.poisson_maxIter = 1800; - -fprintf("\nMEDIUM STUDY SELECTED\n"); - -nCases = numel(cfg.meshes) * numel(cfg.re_list) * numel(cfg.schemes) * ... - numel(cfg.pressure_solvers) * numel(cfg.implementations); - -fprintf("Total simulations: %d\n", nCases); -fprintf("Base max outer iterations per simulation: %d\n\n", cfg.maxIter); - -T = run_parametric_study(cfg); -disp(T); -writetable(T, fullfile("results", "data", "study_summary_medium.csv")); -plot_study_summary(T, cfg); +run_mode('medium'); diff --git a/main_quick.m b/main_quick.m index 0c4584a..e638a38 100644 --- a/main_quick.m +++ b/main_quick.m @@ -1,27 +1,3 @@ -% MAIN_QUICK Run a reduced MATLAB study for a fast solver check. - +% MAIN_QUICK Compare vectorized and loop MATLAB implementations. clear; clc; close all; -addpath("startup"); -setup_project(); - -cfg = default_config(); -cfg.meshes = [32, 64]; -cfg.re_list = [100, 400]; -cfg.maxIter = 2000; -cfg.maxIter_N128_bonus = 0; -cfg.maxIter_Re1000_bonus = 0; -cfg.maxIter_central_bonus = 500; -cfg.poisson_maxIter = 1200; - -fprintf("\nQUICK STUDY SELECTED\n"); - -nCases = numel(cfg.meshes) * numel(cfg.re_list) * numel(cfg.schemes) * ... - numel(cfg.pressure_solvers) * numel(cfg.implementations); - -fprintf("Total simulations: %d\n", nCases); -fprintf("Base max outer iterations per simulation: %d\n\n", cfg.maxIter); - -T = run_parametric_study(cfg); -disp(T); -writetable(T, fullfile("results", "data", "study_summary_quick.csv")); -plot_study_summary(T, cfg); +run_mode('quick'); diff --git a/main_re1000.m b/main_re1000.m new file mode 100644 index 0000000..fc21d77 --- /dev/null +++ b/main_re1000.m @@ -0,0 +1,3 @@ +% MAIN_RE1000 Run the representative Re=1000 production case. +clear; clc; close all; +run_mode('re1000'); diff --git a/main_single.m b/main_single.m new file mode 100644 index 0000000..581057b --- /dev/null +++ b/main_single.m @@ -0,0 +1,3 @@ +% MAIN_SINGLE Run the canonical Phase 2 regression case. +clear; clc; close all; +run_mode('single'); diff --git a/post/plot_residuals.m b/post/plot_residuals.m index 4c12967..12a940b 100644 --- a/post/plot_residuals.m +++ b/post/plot_residuals.m @@ -1,36 +1,28 @@ function plot_residuals(result,cfg,case_name) -%PLOT_RESIDUALS Saves residual history. +%PLOT_RESIDUALS Save Phase 2 convergence and pressure histories. figure('Visible','off'); -semilogy(result.Ru,'LineWidth',1.5); hold on; -semilogy(result.Rv,'LineWidth',1.5); -semilogy(result.Rc_mass,'LineWidth',1.5); -if isfield(result,'Rc_div') - semilogy(result.Rc_div,'--','LineWidth',1.0); -end -grid on; +semilogy(result.velocity_update_linf,'LineWidth',1.5); hold on; +semilogy(result.divergence_linf,'LineWidth',1.3); +semilogy(result.divergence_l2,'LineWidth',1.3); +semilogy(max(result.global_mass_imbalance,eps),'--','LineWidth',1.0); +hold off; grid on; xlabel('Outer iteration'); -ylabel('Residual'); -if isfield(result,'Rc_div') - legend('R_u','R_v','R_c mass','R_c raw div','Location','best'); -else - legend('R_u','R_v','R_c','Location','best'); -end -title(sprintf('Residuals: N=%d Re=%d %s %s %s', ... - result.N,result.Re,result.scheme,result.pressure_solver,result.implementation), ... - 'Interpreter','none'); - -save_current_figure(fullfile(cfg.fig_dir, case_name + "_residuals")); +ylabel('Dimensionless residual'); +legend('Velocity update L_\infty','Divergence L_\infty', ... + 'Divergence L_2','Global mass imbalance','Location','best'); +title(sprintf('Convergence: N=%d Re=%d %s %s %s', ... + result.N,result.Re,result.scheme,result.pressure_solver, ... + result.implementation),'Interpreter','none'); +save_current_figure(fullfile(cfg.fig_dir,case_name + "_residuals")); close; -if isfield(result,'poisson_relative_residual') - figure('Visible','off'); - semilogy(result.poisson_relative_residual,'LineWidth',1.5); - grid on; - xlabel('Outer iteration'); - ylabel('Pressure Poisson relative residual'); - title('Pressure correction residual'); - save_current_figure(fullfile(cfg.fig_dir, case_name + "_pressure_poisson_residual")); - close; -end +figure('Visible','off'); +semilogy(result.poisson_relative_residual,'LineWidth',1.5); +grid on; +xlabel('Outer iteration'); +ylabel('Pressure Poisson relative residual'); +title('Pressure-correction convergence'); +save_current_figure(fullfile(cfg.fig_dir,case_name + "_pressure_poisson_residual")); +close; end diff --git a/post/plot_study_summary.m b/post/plot_study_summary.m index 97a0df0..d36a34a 100644 --- a/post/plot_study_summary.m +++ b/post/plot_study_summary.m @@ -1,104 +1,99 @@ function plot_study_summary(T,cfg) -%PLOT_STUDY_SUMMARY Generates study-level plots from the summary table. -% -% Uses base MATLAB only. No Statistics Toolbox functions. +%PLOT_STUDY_SUMMARY Generate compact Phase 2 study-level figures. if isempty(T) return; end -% Runtime comparison: loop vs vectorized +if ~exist(cfg.fig_dir,'dir') + mkdir(cfg.fig_dir); +end + +% Runtime by implementation. figure('Visible','off'); -impl = unique(string(T.Implementation), 'stable'); +impl = unique(string(T.Implementation),'stable'); hold on; for k = 1:numel(impl) mask = string(T.Implementation) == impl(k); values = T.Runtime_s(mask); - x = k * ones(size(values)); - plot(x, values, 'o', 'MarkerSize', 6, 'DisplayName', char(impl(k))); - plot(k, mean(values), 'x', 'MarkerSize', 10, 'LineWidth', 2, 'HandleVisibility','off'); + plot(k*ones(size(values)),values,'o','MarkerSize',6, ... + 'DisplayName',char(impl(k))); + plot(k,mean(values),'x','MarkerSize',10,'LineWidth',2, ... + 'HandleVisibility','off'); end -hold off; -grid on; -xlim([0.5, numel(impl)+0.5]); -set(gca, 'XTick', 1:numel(impl), 'XTickLabel', cellstr(impl)); +hold off; grid on; +xlim([0.5,numel(impl)+0.5]); +set(gca,'XTick',1:numel(impl),'XTickLabel',cellstr(impl)); ylabel('Runtime [s]'); -title('Runtime: loop vs vectorized'); -save_current_figure(fullfile(cfg.fig_dir, "study_runtime_implementation")); +title('MATLAB implementation runtime'); +save_current_figure(fullfile(cfg.fig_dir,'study_runtime_implementation')); close; -% Pressure solver iterations +% Pressure solver effort. figure('Visible','off'); -ps = unique(string(T.PressureSolver), 'stable'); +solvers = unique(string(T.PressureSolver),'stable'); hold on; -for k = 1:numel(ps) - mask = string(T.PressureSolver) == ps(k); +for k = 1:numel(solvers) + mask = string(T.PressureSolver) == solvers(k); values = T.AvgPoissonIterations(mask); - x = k * ones(size(values)); - plot(x, values, 'o', 'MarkerSize', 6, 'DisplayName', char(ps(k))); - plot(k, mean(values), 'x', 'MarkerSize', 10, 'LineWidth', 2, 'HandleVisibility','off'); + plot(k*ones(size(values)),values,'o','MarkerSize',6, ... + 'DisplayName',char(solvers(k))); + plot(k,mean(values),'x','MarkerSize',10,'LineWidth',2, ... + 'HandleVisibility','off'); end -hold off; -grid on; -xlim([0.5, numel(ps)+0.5]); -set(gca, 'XTick', 1:numel(ps), 'XTickLabel', cellstr(ps)); +hold off; grid on; +xlim([0.5,numel(solvers)+0.5]); +set(gca,'XTick',1:numel(solvers),'XTickLabel',cellstr(solvers)); ylabel('Average pressure iterations'); title('Pressure solver comparison'); -save_current_figure(fullfile(cfg.fig_dir, "study_pressure_solver_iterations")); +save_current_figure(fullfile(cfg.fig_dir,'study_pressure_solver_iterations')); close; -% Continuity mass residual by case -if any(strcmp(T.Properties.VariableNames,'FinalRcMass')) - figure('Visible','off'); - semilogy(T.CaseID, T.FinalRcMass, 'o-','LineWidth',1.2); - grid on; - xlabel('Case ID'); - ylabel('Final normalized mass residual'); - title('Final continuity residual by case'); - save_current_figure(fullfile(cfg.fig_dir, "study_final_mass_residual")); - close; -end +% Final divergence and velocity-update diagnostics. +figure('Visible','off'); +semilogy(T.CaseID,T.FinalVelocityLinf,'o-','LineWidth',1.2); hold on; +semilogy(T.CaseID,T.FinalDivergenceLinf,'s-','LineWidth',1.2); +semilogy(T.CaseID,T.FinalDivergenceL2,'^-','LineWidth',1.2); +hold off; grid on; +xlabel('Case ID'); +ylabel('Final dimensionless residual'); +legend('Velocity update L_\infty','Divergence L_\infty', ... + 'Divergence L_2','Location','best'); +title('Final convergence metrics'); +save_current_figure(fullfile(cfg.fig_dir,'study_final_residuals')); +close; -% Validation error vs mesh +% Ghia error vs mesh. valid = ~isnan(T.Ghia_u_L2); if any(valid) - figure('Visible','off'); - hold on; - - schemes = unique(string(T.Scheme(valid)), 'stable'); + figure('Visible','off'); hold on; + schemes = unique(string(T.Scheme(valid)),'stable'); for k = 1:numel(schemes) - mask = valid & string(T.Scheme)==schemes(k); - plot(T.N(mask), T.Ghia_u_L2(mask), 'o-', ... - 'LineWidth', 1.5, ... - 'MarkerSize', 6, ... - 'DisplayName', char(schemes(k))); + mask = valid & string(T.Scheme) == schemes(k); + plot(T.N(mask),T.Ghia_u_L2(mask),'o-','LineWidth',1.5, ... + 'MarkerSize',6,'DisplayName',char(schemes(k))); end - - hold off; - grid on; - xlabel('Mesh size N'); + hold off; grid on; + xlabel('Mesh cells per direction, N'); ylabel('L2 error in u centerline'); legend('Location','best'); - title('Mesh / scheme validation error vs Ghia'); - save_current_figure(fullfile(cfg.fig_dir, "study_ghia_error")); + title('Ghia centerline error'); + save_current_figure(fullfile(cfg.fig_dir,'study_ghia_error')); close; end -% Quality summary -if any(strcmp(T.Properties.VariableNames,'Quality')) - figure('Visible','off'); - qualities = unique(string(T.Quality), 'stable'); - counts = zeros(size(qualities)); - for k = 1:numel(qualities) - counts(k) = sum(string(T.Quality)==qualities(k)); - end - bar(counts); - grid on; - set(gca,'XTick',1:numel(qualities),'XTickLabel',cellstr(qualities)); - xtickangle(30); - ylabel('Number of cases'); - title('Case quality classification'); - save_current_figure(fullfile(cfg.fig_dir, "study_quality_summary")); - close; +% Quality summary. +figure('Visible','off'); +qualities = unique(string(T.Quality),'stable'); +counts = zeros(size(qualities)); +for k = 1:numel(qualities) + counts(k) = sum(string(T.Quality) == qualities(k)); end +bar(counts); grid on; +set(gca,'XTick',1:numel(qualities),'XTickLabel',cellstr(qualities)); +xtickangle(30); +ylabel('Number of cases'); +title('Case quality classification'); +save_current_figure(fullfile(cfg.fig_dir,'study_quality_summary')); +close; end diff --git a/post/plot_validation.m b/post/plot_validation.m index 9253c3d..2f86905 100644 --- a/post/plot_validation.m +++ b/post/plot_validation.m @@ -6,27 +6,37 @@ function plot_validation(result,cfg,case_name) return; end -mid = round((result.N+1)/2); +u_centerline = zeros(numel(result.y),1); +for i = 1:numel(result.y) + u_centerline(i) = interp1(result.x,result.u(i,:),0.5, ... + 'linear','extrap'); +end + +v_centerline = zeros(numel(result.x),1); +for j = 1:numel(result.x) + v_centerline(j) = interp1(result.y,result.v(:,j),0.5, ... + 'linear','extrap'); +end figure('Visible','off'); -plot(result.u(:,mid), result.y, 'LineWidth', 1.5); hold on; -plot(data.u, data.y_u, 'o', 'LineWidth', 1.2); +plot(u_centerline,result.y,'LineWidth',1.5); hold on; +plot(data.u,data.y_u,'o','LineWidth',1.2); grid on; xlabel('u velocity at x=0.5'); ylabel('y'); legend('Solver','Ghia et al.','Location','best'); -title(sprintf('Vertical centerline validation, Re=%d', result.Re)); -save_current_figure(fullfile(cfg.fig_dir, case_name + "_ghia_u")); +title(sprintf('Vertical centerline benchmark, Re=%d',result.Re)); +save_current_figure(fullfile(cfg.fig_dir,case_name + "_ghia_u")); close; figure('Visible','off'); -plot(result.x, result.v(mid,:), 'LineWidth', 1.5); hold on; -plot(data.x_v, data.v, 'o', 'LineWidth', 1.2); +plot(result.x,v_centerline,'LineWidth',1.5); hold on; +plot(data.x_v,data.v,'o','LineWidth',1.2); grid on; xlabel('x'); ylabel('v velocity at y=0.5'); legend('Solver','Ghia et al.','Location','best'); -title(sprintf('Horizontal centerline validation, Re=%d', result.Re)); -save_current_figure(fullfile(cfg.fig_dir, case_name + "_ghia_v")); +title(sprintf('Horizontal centerline benchmark, Re=%d',result.Re)); +save_current_figure(fullfile(cfg.fig_dir,case_name + "_ghia_v")); close; end diff --git a/post/write_case_outputs.m b/post/write_case_outputs.m new file mode 100644 index 0000000..13bba44 --- /dev/null +++ b/post/write_case_outputs.m @@ -0,0 +1,56 @@ +function write_case_outputs(result,metrics,quality,cfg,case_name) +%WRITE_CASE_OUTPUTS Save standardized MAT and CSV files for one case. + +if ~exist(cfg.data_dir,'dir') + mkdir(cfg.data_dir); +end + +save(fullfile(cfg.data_dir,case_name + ".mat"), ... + 'result','metrics','quality','-v7.3'); + +iteration = (1:result.iterations)'; +history = table(iteration, ... + result.velocity_update_linf(:), ... + result.divergence_linf(:), ... + result.divergence_l2(:), ... + result.global_mass_imbalance(:), ... + result.dt(:), ... + result.poisson_relative_residual(:), ... + result.poisson_iters(:), ... + result.poisson_converged(:), ... + 'VariableNames',{'Iteration','VelocityUpdateLinf', ... + 'DivergenceLinf','DivergenceL2','GlobalMassImbalance', ... + 'Dt','PoissonRelativeResidual','PoissonIterations', ... + 'PoissonConverged'}); +writetable(history,fullfile(cfg.data_dir,case_name + "_history.csv")); + +if cfg.save_fields + [X,Y] = meshgrid(result.x,result.y); + fields = table(X(:),Y(:),result.u(:),result.v(:),result.p(:), ... + result.speed(:),result.vorticity(:), ... + 'VariableNames',{'x','y','u','v','p','speed','vorticity'}); + writetable(fields,fullfile(cfg.data_dir,case_name + "_fields.csv")); +end + +u_centerline = interpolate_vertical_centerline(result.u,result.x); +v_centerline = interpolate_horizontal_centerline(result.v,result.y); +profiles = table(result.x(:),result.y(:),u_centerline(:),v_centerline(:), ... + 'VariableNames',{'x','y','u_at_x_0p5','v_at_y_0p5'}); +writetable(profiles,fullfile(cfg.data_dir,case_name + "_centerlines.csv")); +end + +function values = interpolate_vertical_centerline(field,x) +N = size(field,1); +values = zeros(N,1); +for i = 1:N + values(i) = interp1(x,field(i,:),0.5,'linear','extrap'); +end +end + +function values = interpolate_horizontal_centerline(field,y) +N = size(field,2); +values = zeros(N,1); +for j = 1:N + values(j) = interp1(y,field(:,j),0.5,'linear','extrap'); +end +end diff --git a/run_mode.m b/run_mode.m new file mode 100644 index 0000000..905b9e9 --- /dev/null +++ b/run_mode.m @@ -0,0 +1,33 @@ +function T = run_mode(mode) +%RUN_MODE Run one named Phase 2 MATLAB study. + +if nargin < 1 + mode = 'single'; +end +mode = lower(string(mode)); + +project_root = fileparts(mfilename('fullpath')); +original_folder = pwd; +cleanup = onCleanup(@() cd(original_folder)); %#ok +cd(project_root); + +addpath('startup'); +setup_project(); +cfg = default_config(); + +% Keep the canonical run visual, while larger batch studies prioritize data. +if mode == "single" + cfg.make_figures = true; + cfg.figure_every_case = true; +else + cfg.make_figures = true; + cfg.figure_every_case = false; +end + +T = run_parametric_study(cfg,mode); + +disp(' '); +disp('Finished. Summary table:'); +disp(T); +fprintf('Results saved in %s\n',cfg.data_dir); +end diff --git a/run_tests.m b/run_tests.m new file mode 100644 index 0000000..b642b02 --- /dev/null +++ b/run_tests.m @@ -0,0 +1,8 @@ +function run_tests() +%RUN_TESTS Execute the MATLAB Phase 2 regression suite. +project_root = fileparts(mfilename('fullpath')); +addpath(fullfile(project_root,'startup')); +setup_project(); +addpath(fullfile(project_root,'tests')); +test_phase2_regression(); +end diff --git a/scripts/run.sh b/scripts/run.sh index 17e0667..bc355a5 100644 --- a/scripts/run.sh +++ b/scripts/run.sh @@ -1,5 +1,4 @@ -#!/bin/bash -set -e -cd "$(dirname "$0")/.." -mkdir -p results/data results/figures -matlab -nodisplay -nosplash -r "main; exit" +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" full diff --git a/scripts/run_full.sh b/scripts/run_full.sh new file mode 100644 index 0000000..bc355a5 --- /dev/null +++ b/scripts/run_full.sh @@ -0,0 +1,4 @@ +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" full diff --git a/scripts/run_grid.sh b/scripts/run_grid.sh new file mode 100644 index 0000000..e31861b --- /dev/null +++ b/scripts/run_grid.sh @@ -0,0 +1,4 @@ +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" grid diff --git a/scripts/run_medium.sh b/scripts/run_medium.sh index 2842eaf..0167775 100644 --- a/scripts/run_medium.sh +++ b/scripts/run_medium.sh @@ -1,5 +1,4 @@ -#!/bin/bash -set -e -cd "$(dirname "$0")/.." -mkdir -p results/data results/figures -matlab -nodisplay -nosplash -r "main_medium; exit" +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" medium diff --git a/scripts/run_mode.sh b/scripts/run_mode.sh new file mode 100644 index 0000000..b6c82ea --- /dev/null +++ b/scripts/run_mode.sh @@ -0,0 +1,30 @@ +#!/usr/bin/env bash +set -euo pipefail + +MODE="${1:-single}" +ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +MATLAB_BIN="${MATLAB_BIN:-matlab}" + +case "$MODE" in + single|quick|medium|grid|re1000|full) ;; + *) + echo "Unknown mode: $MODE" >&2 + echo "Use: single, quick, medium, grid, re1000, or full" >&2 + exit 2 + ;; +esac + +if ! command -v "$MATLAB_BIN" >/dev/null 2>&1; then + echo "MATLAB executable not found: $MATLAB_BIN" >&2 + echo "Load the MATLAB module or set MATLAB_BIN to its executable." >&2 + exit 127 +fi + +mkdir -p "$ROOT/logs" +LOG="$ROOT/logs/matlab_${MODE}_$(date +%Y%m%d_%H%M%S).log" + +echo "Project: $ROOT" +echo "Mode: $MODE" +echo "Log: $LOG" + +"$MATLAB_BIN" -batch "cd('$ROOT'); run_mode('$MODE');" 2>&1 | tee "$LOG" diff --git a/scripts/run_quick.sh b/scripts/run_quick.sh index d558031..39f2660 100644 --- a/scripts/run_quick.sh +++ b/scripts/run_quick.sh @@ -1,5 +1,4 @@ -#!/bin/bash -set -e -cd "$(dirname "$0")/.." -mkdir -p results/data results/figures -matlab -nodisplay -nosplash -r "main_quick; exit" +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" quick diff --git a/scripts/run_re1000.sh b/scripts/run_re1000.sh new file mode 100644 index 0000000..4d12d00 --- /dev/null +++ b/scripts/run_re1000.sh @@ -0,0 +1,4 @@ +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" re1000 diff --git a/scripts/run_single.sh b/scripts/run_single.sh new file mode 100644 index 0000000..cdb27b4 --- /dev/null +++ b/scripts/run_single.sh @@ -0,0 +1,4 @@ +#!/usr/bin/env bash +set -euo pipefail +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +bash "$SCRIPT_DIR/run_mode.sh" single diff --git a/scripts/run_tests.sh b/scripts/run_tests.sh new file mode 100644 index 0000000..1679018 --- /dev/null +++ b/scripts/run_tests.sh @@ -0,0 +1,5 @@ +#!/usr/bin/env bash +set -euo pipefail +ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +MATLAB_BIN="${MATLAB_BIN:-matlab}" +"$MATLAB_BIN" -batch "cd('$ROOT'); run_tests();" diff --git a/studies/run_parametric_study.m b/studies/run_parametric_study.m index 247ee52..806cfa1 100644 --- a/studies/run_parametric_study.m +++ b/studies/run_parametric_study.m @@ -1,105 +1,120 @@ -function T = run_parametric_study(cfg) -%RUN_PARAMETRIC_STUDY Runs mesh, Re, scheme, solver, implementation study. - -rows = {}; - -case_id = 0; - -for iN = 1:numel(cfg.meshes) - N = cfg.meshes(iN); - - for iR = 1:numel(cfg.re_list) - Re = cfg.re_list(iR); - - for is = 1:numel(cfg.schemes) - scheme = cfg.schemes{is}; - - for ip = 1:numel(cfg.pressure_solvers) - pressure_solver = cfg.pressure_solvers{ip}; - - for ii = 1:numel(cfg.implementations) - implementation = cfg.implementations{ii}; - - case_id = case_id + 1; - case_name = sprintf("case_%03d_N%d_Re%d_%s_%s_%s", ... - case_id,N,Re,scheme,pressure_solver,implementation); - - fprintf("\n[%03d] N=%d Re=%d Scheme=%s Pressure=%s Implementation=%s\n", ... - case_id,N,Re,scheme,pressure_solver,implementation); - - result = solve_lid_cavity(N,Re,scheme,pressure_solver,implementation,cfg); - metrics = validate_against_ghia(result,cfg); - - quality = "not_validated"; - if result.status == "converged" && metrics.available && metrics.pass - quality = "converged_validated"; - elseif result.status == "converged" && metrics.available && ~metrics.pass - quality = "converged_not_validated"; - elseif result.status == "converged" && ~metrics.available - quality = "converged_no_benchmark"; - elseif result.status ~= "converged" && metrics.available && metrics.pass - quality = "validated_but_not_converged"; - else - quality = "needs_improvement"; - end - - fprintf(" status=%s quality=%s iter=%d/%d Rc_mass=%.3e Rc_div=%.3e runtime=%.2fs avgPiter=%.1f pSat=%.2f\n", ... - result.status,quality,result.iterations,result.localMaxIter,result.final_Rc_mass,result.final_Rc_div, ... - result.runtime,result.avg_poisson_iters,result.pressure_saturation_ratio); - - if metrics.available - fprintf(" Ghia L2: u=%.3e(limit %.3e), v=%.3e(limit %.3e), pass=%d\n", ... - metrics.u_L2,metrics.u_limit,metrics.v_L2,metrics.v_limit,metrics.pass); - end - - save(fullfile(cfg.data_dir, case_name + ".mat"), "result", "metrics", "quality"); - - if cfg.make_figures && cfg.figure_every_case - plot_residuals(result,cfg,case_name); - plot_fields(result,cfg,case_name); - plot_validation(result,cfg,case_name); - end - - rows(end+1,:) = { ... - case_id, ... - string(implementation), ... - N, ... - Re, ... - string(scheme), ... - string(pressure_solver), ... - string(result.status), ... - string(quality), ... - result.iterations, ... - result.localMaxIter, ... - result.final_Ru, ... - result.final_Rv, ... - result.final_Rc_mass, ... - result.final_Rc_div, ... - result.runtime, ... - result.avg_poisson_iters, ... - result.avg_poisson_relative_residual, ... - result.pressure_saturation_ratio, ... - metrics.available, ... - metrics.pass, ... - metrics.u_L2, ... - metrics.v_L2, ... - metrics.u_Linf, ... - metrics.v_Linf, ... - metrics.u_limit, ... - metrics.v_limit}; - end - end - end +function T = run_parametric_study(cfg,mode) +%RUN_PARAMETRIC_STUDY Run a Phase 2 mode with continuation and strict checks. + +if nargin < 2 + mode = 'full'; +end +mode = lower(string(mode)); +cases = study_definition(mode,cfg); +rows = cell(height(cases),29); +continuation = containers.Map('KeyType','char','ValueType','any'); + +fprintf('\nLID-DRIVEN CAVITY MATLAB PHASE 2 SOLVER\n'); +fprintf('Mode: %s\n',mode); +fprintf('Total simulations: %d\n',height(cases)); +fprintf('Summary: %s\n\n', ... + fullfile(cfg.data_dir,"study_summary_" + mode + ".csv")); + +for k = 1:height(cases) + N = cases.N(k); + Re = cases.Re(k); + scheme = string(cases.Scheme(k)); + pressure_solver = string(cases.PressureSolver(k)); + implementation = string(cases.Implementation(k)); + + continuation_key = sprintf('N%d_%s_%s_%s',N, ... + char(scheme),char(pressure_solver),char(implementation)); + initial_state = struct('available',false); + if cfg.use_continuation && isKey(continuation,continuation_key) + initial_state = continuation(continuation_key); + end + + case_name = sprintf('case_%03d_N%d_Re%d_%s_%s_%s', ... + k,N,Re,char(scheme),char(pressure_solver),char(implementation)); + + fprintf('[%03d] N=%d Re=%d Scheme=%s Pressure=%s Implementation=%s\n', ... + k,N,Re,scheme,pressure_solver,implementation); + + result = solve_lid_cavity(N,Re,scheme,pressure_solver, ... + implementation,cfg,initial_state); + metrics = validate_against_ghia(result,cfg); + quality = classify_quality(result,metrics); + + if strcmp(result.status,'converged') + continuation(continuation_key) = result.continuation_state; + elseif isKey(continuation,continuation_key) + remove(continuation,continuation_key); + end + + fprintf([' status=%s quality=%s iter=%d/%d vel=%.6e ' ... + 'div=%.6e runtime=%.2fs\n'], ... + result.status,quality,result.iterations,result.localMaxIter, ... + result.final_velocity_linf,result.final_divergence_linf, ... + result.runtime); + if metrics.available + fprintf(' Ghia L2: u=%.3e v=%.3e pass=%d\n', ... + metrics.u_L2,metrics.v_L2,metrics.pass); + end + + write_case_outputs(result,metrics,quality,cfg,case_name); + + if cfg.make_figures && cfg.figure_every_case + plot_residuals(result,cfg,case_name); + plot_fields(result,cfg,case_name); + plot_validation(result,cfg,case_name); end + + rows(k,:) = { ... + k,N,Re,scheme,pressure_solver,implementation, ... + string(result.status),string(quality),result.iterations, ... + result.localMaxIter,result.consecutive_pass_count, ... + result.failed_pressure_solves,result.final_velocity_linf, ... + result.final_divergence_linf,result.final_divergence_l2, ... + result.final_global_mass,result.final_poisson_relative_residual, ... + result.runtime,result.avg_poisson_iters, ... + result.avg_poisson_relative_residual, ... + result.pressure_saturation_ratio,metrics.available,metrics.pass, ... + metrics.u_L2,metrics.v_L2,metrics.u_Linf,metrics.v_Linf, ... + metrics.u_limit,metrics.v_limit}; end -T = cell2table(rows, 'VariableNames', { ... - 'CaseID','Implementation','N','Re','Scheme','PressureSolver','Status','Quality', ... - 'Iterations','LocalMaxIter','FinalRu','FinalRv','FinalRcMass','FinalRcDiv','Runtime_s', ... - 'AvgPoissonIterations','AvgPoissonRelResidual','PressureSaturationRatio', ... - 'HasGhia','ValidationPass','Ghia_u_L2','Ghia_v_L2','Ghia_u_Linf','Ghia_v_Linf', ... +T = cell2table(rows,'VariableNames',{ ... + 'CaseID','N','Re','Scheme','PressureSolver','Implementation', ... + 'Status','Quality','Iterations','LocalMaxIter', ... + 'ConsecutivePassCount','FailedPressureSolves', ... + 'FinalVelocityLinf','FinalDivergenceLinf','FinalDivergenceL2', ... + 'FinalGlobalMass','FinalPoissonRelativeResidual','Runtime_s', ... + 'AvgPoissonIterations','AvgPoissonRelResidual', ... + 'PressureSaturationRatio','HasGhia','ValidationPass', ... + 'Ghia_u_L2','Ghia_v_L2','Ghia_u_Linf','Ghia_v_Linf', ... 'Ghia_u_L2_Limit','Ghia_v_L2_Limit'}); -writetable(T, fullfile(cfg.data_dir, "study_summary.csv")); -save(fullfile(cfg.data_dir, "study_summary.mat"), "T"); +summary_path = fullfile(cfg.data_dir,"study_summary_" + mode + ".csv"); +writetable(T,summary_path); +save(fullfile(cfg.data_dir,"study_summary_" + mode + ".mat"),'T'); + +if cfg.make_figures + plot_study_summary(T,cfg); +end + +if cfg.strict && any(T.Status ~= "converged") + failed = T(T.Status ~= "converged", ... + {'CaseID','N','Re','Scheme','PressureSolver','Implementation','Status'}); + disp(failed); + error('Phase 2 strict mode: %d case(s) did not converge.',height(failed)); +end +end + +function quality = classify_quality(result,metrics) +if strcmp(result.status,'converged') && metrics.available && metrics.pass + quality = "converged_benchmark_pass"; +elseif strcmp(result.status,'converged') && metrics.available + quality = "converged_benchmark_fail"; +elseif strcmp(result.status,'converged') + quality = "converged_no_benchmark"; +elseif metrics.available && metrics.pass + quality = "benchmark_pass_not_converged"; +else + quality = "needs_improvement"; +end end diff --git a/studies/run_single_case.m b/studies/run_single_case.m index b70d700..fa6e704 100644 --- a/studies/run_single_case.m +++ b/studies/run_single_case.m @@ -1,28 +1,35 @@ function result = run_single_case() -%RUN_SINGLE_CASE Fast test case for debugging. +%RUN_SINGLE_CASE Run the canonical Phase 2 case interactively. -addpath(fullfile("..", "startup")); -if exist("setup_project", "file") ~= 2 - addpath("startup"); -end +project_root = fileparts(fileparts(mfilename('fullpath'))); +addpath(fullfile(project_root,'startup')); setup_project(); cfg = default_config(); -cfg.maxIter = 500; cfg.make_figures = true; +cfg.figure_every_case = true; -N = 32; -Re = 100; -scheme = 'upwind'; -pressure_solver = 'RBSOR'; -implementation = 'vectorized'; +result = solve_lid_cavity(32,100,'upwind','RBSOR','vectorized',cfg); +metrics = validate_against_ghia(result,cfg); -result = solve_lid_cavity(N,Re,scheme,pressure_solver,implementation,cfg); +if strcmp(result.status,'converged') && metrics.available && metrics.pass + quality = "converged_benchmark_pass"; +elseif strcmp(result.status,'converged') + quality = "converged_benchmark_fail"; +else + quality = "needs_improvement"; +end -case_name = sprintf("single_N%d_Re%d_%s_%s_%s",N,Re,scheme,pressure_solver,implementation); +case_name = "single_N32_Re100_upwind_RBSOR_vectorized"; +write_case_outputs(result,metrics,quality,cfg,case_name); plot_residuals(result,cfg,case_name); plot_fields(result,cfg,case_name); plot_validation(result,cfg,case_name); disp(result); +disp(metrics); + +if cfg.strict && ~strcmp(result.status,'converged') + error('Canonical Phase 2 case did not converge: %s',result.status); +end end diff --git a/studies/study_definition.m b/studies/study_definition.m new file mode 100644 index 0000000..82a8a24 --- /dev/null +++ b/studies/study_definition.m @@ -0,0 +1,63 @@ +function cases = study_definition(mode,cfg) +%STUDY_DEFINITION Return the ordered Phase 2 case matrix for one run mode. + +mode = lower(string(mode)); +rows = {}; + +switch mode + case "single" + rows = add_cases(rows,32,100,{'upwind'},{'RBSOR'},{'vectorized'}); + + case "quick" + % Compare MATLAB implementations on the canonical case. + rows = add_cases(rows,32,100,{'upwind'},{'RBSOR'}, ... + {'vectorized','loop'}); + + case "medium" + % Same six production cases used by the C++ Phase 2 verification. + rows = add_cases(rows,32,[100,400,1000], ... + {'upwind','central'},{'RBSOR'},{'vectorized'}); + + case "grid" + rows = add_cases(rows,[16,32,64],100, ... + {'central'},{'RBSOR'},{'vectorized'}); + + case "re1000" + rows = add_cases(rows,64,1000, ... + {'central'},{'RBSOR'},{'vectorized'}); + + case "full" + % Complete MATLAB study, including loop/vectorized comparison. + rows = add_cases(rows,cfg.meshes,cfg.re_list,cfg.schemes, ... + cfg.pressure_solvers,cfg.implementations); + + otherwise + error('Unknown study mode: %s',mode); +end + +case_ids = num2cell((1:size(rows,1))'); +rows = [case_ids,rows]; +cases = cell2table(rows,'VariableNames', ... + {'CaseID','N','Re','Scheme','PressureSolver','Implementation'}); +end + +function rows = add_cases(rows,meshes,reynolds,schemes,solvers,implementations) +% Order Reynolds number inside a fixed numerical configuration so that +% continuation can reuse the lower-Re solution. +for iN = 1:numel(meshes) + for iS = 1:numel(schemes) + for iP = 1:numel(solvers) + for iI = 1:numel(implementations) + for iR = 1:numel(reynolds) + rows(end+1,:) = { ... + meshes(iN), ... + reynolds(iR), ... + string(schemes{iS}), ... + string(solvers{iP}), ... + string(implementations{iI})}; %#ok + end + end + end + end +end +end diff --git a/tests/test_phase2_regression.m b/tests/test_phase2_regression.m new file mode 100644 index 0000000..d628c17 --- /dev/null +++ b/tests/test_phase2_regression.m @@ -0,0 +1,33 @@ +function test_phase2_regression() +%TEST_PHASE2_REGRESSION Run the canonical strict-convergence regression. + +project_root = fileparts(fileparts(mfilename('fullpath'))); +addpath(fullfile(project_root,'startup')); +setup_project(); + +cfg = default_config(); +cfg.make_figures = false; +cfg.figure_every_case = false; +cfg.save_fields = false; +cfg.strict = true; + +result = solve_lid_cavity(32,100,'upwind','RBSOR','vectorized',cfg); +metrics = validate_against_ghia(result,cfg); + +assert(strcmp(result.status,'converged'), ... + 'Canonical case did not converge: %s',result.status); +assert(result.final_velocity_linf <= cfg.tol_velocity_linf, ... + 'Velocity-update tolerance failed.'); +assert(result.final_divergence_linf <= cfg.tol_divergence_linf, ... + 'Divergence Linf tolerance failed.'); +assert(result.final_divergence_l2 <= cfg.tol_divergence_l2, ... + 'Divergence L2 tolerance failed.'); +assert(result.final_global_mass <= cfg.tol_global_mass, ... + 'Global mass tolerance failed.'); +assert(metrics.available && metrics.pass, ... + 'Canonical Ghia benchmark comparison failed.'); + +fprintf(['PASS: canonical Phase 2 case converged in %d iterations, ' ... + 'u_L2=%.3e, v_L2=%.3e\n'], ... + result.iterations,metrics.u_L2,metrics.v_L2); +end diff --git a/validation/validate_against_ghia.m b/validation/validate_against_ghia.m index bb4bff7..a6d4467 100644 --- a/validation/validate_against_ghia.m +++ b/validation/validate_against_ghia.m @@ -1,12 +1,11 @@ function metrics = validate_against_ghia(result,cfg) -%VALIDATE_AGAINST_GHIA Compares centerline profiles to Ghia data if available. +%VALIDATE_AGAINST_GHIA Compare interpolated centerlines with Ghia data. if nargin < 2 cfg = default_config(); end data = ghia_data(result.Re); - metrics.available = ~isempty(data); metrics.u_L2 = NaN; metrics.v_L2 = NaN; @@ -20,20 +19,23 @@ return; end -N = result.N; -mid = round((N+1)/2); - -% Numerical profiles. u at x=0.5 as function of y, v at y=0.5 as function of x. -u_center = result.u(:,mid); -v_center = result.v(mid,:); +u_centerline = zeros(numel(result.y),1); +for i = 1:numel(result.y) + u_centerline(i) = interp1(result.x,result.u(i,:),0.5, ... + 'linear','extrap'); +end -% Interpolate numerical values onto Ghia sample points. -u_num = interp1(result.y, u_center, data.y_u, 'linear', 'extrap'); -v_num = interp1(result.x, v_center, data.x_v, 'linear', 'extrap'); +v_centerline = zeros(numel(result.x),1); +for j = 1:numel(result.x) + v_centerline(j) = interp1(result.y,result.v(:,j),0.5, ... + 'linear','extrap'); +end -eu = u_num(:) - data.u(:); -ev = v_num(:) - data.v(:); +u_num = interp1(result.y,u_centerline,data.y_u,'linear','extrap'); +v_num = interp1(result.x,v_centerline,data.x_v,'linear','extrap'); +eu = u_num(:)-data.u(:); +ev = v_num(:)-data.v(:); metrics.u_L2 = sqrt(mean(eu.^2)); metrics.v_L2 = sqrt(mean(ev.^2)); metrics.u_Linf = max(abs(eu)); @@ -54,5 +56,6 @@ metrics.v_limit = inf; end -metrics.pass = metrics.u_L2 <= metrics.u_limit && metrics.v_L2 <= metrics.v_limit; +metrics.pass = metrics.u_L2 <= metrics.u_limit ... + && metrics.v_L2 <= metrics.v_limit; end