diff --git a/.github/workflows/verification.yml b/.github/workflows/verification.yml index 86aeed6..c2514eb 100644 --- a/.github/workflows/verification.yml +++ b/.github/workflows/verification.yml @@ -45,3 +45,12 @@ jobs: - name: Test run: ctest --test-dir build/verification --output-on-failure + + canonical-regression: + runs-on: ubuntu-latest + steps: + - name: Check out repository + uses: actions/checkout@v4 + + - name: Run canonical converged case + run: bash scripts/run_canonical_regression.sh diff --git a/CMakeLists.txt b/CMakeLists.txt index 146a29a..358ef6a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.16) -project(LidCavityCPPVerification LANGUAGES CXX) +project(LidCavityCPP LANGUAGES CXX) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) @@ -7,35 +7,49 @@ set(CMAKE_CXX_EXTENSIONS OFF) option(LIDCAVITY_ENABLE_SANITIZERS "Enable AddressSanitizer and UndefinedBehaviorSanitizer" OFF) +function(lidcavity_apply_warnings target) + target_compile_options(${target} PRIVATE + -Wall -Wextra -Wpedantic -Wconversion -Wshadow + ) + if(LIDCAVITY_ENABLE_SANITIZERS AND CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") + target_compile_options(${target} PRIVATE -fsanitize=address,undefined -fno-omit-frame-pointer) + target_link_options(${target} PRIVATE -fsanitize=address,undefined) + endif() +endfunction() + +add_executable(lid_cavity src/lid_cavity.cpp) +lidcavity_apply_warnings(lid_cavity) + add_library(lidcavity_verification src/verification/operators.cpp src/verification/poisson.cpp src/verification/convergence.cpp ) - target_include_directories(lidcavity_verification PUBLIC include) -target_compile_options(lidcavity_verification PRIVATE - -Wall -Wextra -Wpedantic -Wconversion -Wshadow -) - -if(LIDCAVITY_ENABLE_SANITIZERS AND CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") - target_compile_options(lidcavity_verification PRIVATE -fsanitize=address,undefined -fno-omit-frame-pointer) - target_link_options(lidcavity_verification PRIVATE -fsanitize=address,undefined) -endif() +lidcavity_apply_warnings(lidcavity_verification) enable_testing() function(add_lidcavity_test target source) add_executable(${target} ${source}) target_link_libraries(${target} PRIVATE lidcavity_verification) - target_compile_options(${target} PRIVATE -Wall -Wextra -Wpedantic -Wconversion -Wshadow) - if(LIDCAVITY_ENABLE_SANITIZERS AND CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") - target_compile_options(${target} PRIVATE -fsanitize=address,undefined -fno-omit-frame-pointer) - target_link_options(${target} PRIVATE -fsanitize=address,undefined) - endif() + lidcavity_apply_warnings(${target}) add_test(NAME ${target} COMMAND ${target}) endfunction() add_lidcavity_test(test_operators tests/test_operators.cpp) add_lidcavity_test(test_poisson tests/test_poisson.cpp) add_lidcavity_test(test_convergence tests/test_convergence.cpp) + +add_test( + NAME canonical_cavity_regression + COMMAND lid_cavity + --single + --N 32 + --Re 100 + --scheme upwind + --pressure RBSOR + --strict + --no-fields +) +set_tests_properties(canonical_cavity_regression PROPERTIES TIMEOUT 120) diff --git a/README.md b/README.md index 2cf9650..fbf5856 100644 --- a/README.md +++ b/README.md @@ -1,216 +1,223 @@ # Lid-Driven Cavity Flow Solver in C++

- Completed + Phase 2 verified C++17 - GCC 8 and newer - Python post-processing + Staggered MAC grid - C++ build and smoke test + Build and smoke test - MIT License - - Portfolio case study + + Numerical verification + MIT License

-A completed C++17 implementation and parameter study of the two-dimensional lid-driven cavity benchmark. - -This repository is the serial C++ component of a larger CFD comparison project. It solves the same benchmark used by the MATLAB and multi-language implementations so that numerical behavior, result quality, code structure, and runtime can be compared consistently. - -The repository is an educational solver and benchmark study, not a production CFD package. - -## What the code does +A serial C++17 solver for the two-dimensional incompressible lid-driven cavity benchmark. -The solver runs the incompressible lid-driven cavity problem on a structured Cartesian grid. The top wall moves, the other walls are fixed, and the flow develops the characteristic cavity recirculation. +The production solver uses a staggered Marker-and-Cell arrangement, a pseudo-transient projection method, compatible pressure-gradient and divergence operators, and convergence-aware termination. The repository also contains manufactured Poisson verification, discrete-operator tests, sanitizer builds, Ghia centerline comparisons, and configurable parameter studies. -Implemented features: +## What is verified -- serial C++17 solver -- collocated Cartesian grid -- pseudo-transient pressure-correction method -- upwind and central convection schemes -- red-black Gauss-Seidel and red-black SOR pressure solvers -- CSV output for fields, residuals, and study summaries -- Python scripts for plotting fields, residuals, validation, and runtime summaries -- GitHub Actions build, smoke execution, and output verification -- portable filesystem linking for modern compilers and GCC 8 HPC nodes +The current Phase 2 regression set has been exercised with: -The completed study contains 36 configured cases: +- the canonical `N=32`, `Re=100`, upwind, RBSOR case; +- all six `N=32` combinations at `Re=100`, `400`, and `1000` with upwind and central convection using RBSOR; +- the `N=16`, `32`, and `64`, `Re=100`, central, RBSOR grid sequence; +- RBGS/RBSOR manufactured Poisson tests; +- the compatibility test between the selected divergence, gradient, and Laplacian operators; +- GCC and Clang builds in Debug and Release configurations; +- AddressSanitizer and UndefinedBehaviorSanitizer checks. -```text -3 meshes × 3 Reynolds numbers × 2 schemes × 2 pressure solvers -``` +The canonical strict case converges automatically rather than stopping at a configured iteration limit. -## Representative result +## Numerical method -The case shown below uses: +The solver advances the nondimensional incompressible Navier–Stokes equations using: -```text -N = 128 -Re = 1000 -scheme = central -pressure solver = RBSOR -implementation = serial_cpp -``` +1. a staggered MAC grid; +2. an explicit pseudo-time momentum predictor; +3. upwind or central convection differencing; +4. a pressure-correction Poisson equation; +5. RBGS or RBSOR pressure iteration; +6. velocity correction using pressure gradients located consistently with the staggered velocity components; +7. convergence checks based on velocity updates, local divergence, global mass balance, and pressure convergence. -| Streamlines | Velocity magnitude | -|---|---| -| ![Streamlines](assets/figures/re1000_streamlines.svg) | ![Velocity magnitude](assets/figures/re1000_speed.svg) | +Velocity components are stored on cell faces and pressure is stored at cell centers. Cell-centered velocity values are reconstructed for CSV output and post-processing. -## Validation +## Convergence contract -The numerical profiles are compared with the classical Ghia et al. lid-driven cavity data: +A production case is reported as `converged` only when all of the following remain satisfied for the configured number of consecutive iterations: -- `u(y)` on the vertical centerline `x = 0.5` -- `v(x)` on the horizontal centerline `y = 0.5` +- dimensionless velocity-update `Linf` residual; +- dimensionless divergence `Linf` residual; +- dimensionless divergence `L2` residual; +- global boundary mass imbalance; +- successful pressure-Poisson convergence. -| Ghia u comparison | Ghia v comparison | -|---|---| -| ![Ghia u validation](assets/figures/re1000_ghia_u.svg) | ![Ghia v validation](assets/figures/re1000_ghia_v.svg) | +Possible terminal states include: -For each case, the code reports `L2` and `Linf` errors against the benchmark points. +- `converged` +- `max_iterations` +- `pressure_not_converged` +- `stagnated` +- `diverged` +- `non_finite` -The refined-grid central + RBSOR cases are: +Use `--strict` when a script or CI job must fail if any requested case does not converge. -| Re | Case | N | Scheme | Pressure solver | Ghia `u` L2 | Ghia `v` L2 | Runtime [s] | -|---:|---:|---:|---|---|---:|---:|---:| -| 100 | 28 | 128 | central | RBSOR | 0.0031 | 0.0041 | 441.7 | -| 400 | 32 | 128 | central | RBSOR | 0.0539 | 0.0652 | 527.6 | -| 1000 | 36 | 128 | central | RBSOR | 0.1102 | 0.1109 | 647.6 | +## Quick start -Study observations: +On Linux, WSL, or a Linux HPC node: -- all 36 configured cases executed -- 22 cases met the selected Ghia error thresholds -- all 12 cases with `N = 128` met those thresholds -- central differencing produced the best refined-grid agreement -- RBSOR produced similar validation errors to RBGS with lower pressure-solver cost +```bash +bash scripts/run_single.sh +``` -![Ghia error summary](assets/figures/study_ghia_error.svg) +This builds the solver and runs the canonical strict case: -![Pressure solver comparison](assets/figures/study_pressure_solver_iterations.svg) +```text +N = 32 +Re = 100 +scheme = upwind +pressure solver = RBSOR +``` -The selected Ghia limits are comparison thresholds, not a formal verification or grid-independence study. +A successful run ends with `status=converged` and writes CSV files to `results/data/`. -## Convergence interpretation +## Available runs -All full-study cases reached the configured maximum outer-iteration limit. Therefore: +```bash +bash scripts/run_smoke_test.sh # compilation and output check only +bash scripts/run_single.sh # canonical converged regression +bash scripts/run_quick.sh # four fast Re=100 cases +bash scripts/run_medium.sh # six N=32 cases at Re=100/400/1000 +bash scripts/run_grid.sh # N=16/32/64 grid sequence at Re=100 +bash scripts/run_re1000.sh # converged N=32, Re=1000 central case +bash scripts/run_full.sh # complete 36-case configuration +``` -- an executed case is not automatically a converged case -- the reported runtime is the cost of the configured run -- the Ghia error thresholds describe profile agreement, not residual convergence -- the high-Reynolds-number cases require additional convergence tuning +The full mode includes RBGS and grids up to `N=128`; it is intended for a workstation or HPC node. -This distinction is important when comparing the results with other solver implementations. +## Direct command-line use -## Numerical workflow +```bash +bash scripts/build.sh + +bin/lid_cavity \ + --single \ + --N 32 \ + --Re 100 \ + --scheme upwind \ + --pressure RBSOR \ + --strict +``` -The solver advances the nondimensional incompressible Navier-Stokes equations in pseudo-time. Each outer iteration predicts velocity, solves the pressure-correction equation, corrects velocity and pressure, reapplies wall boundary conditions, and records residual information. +Important numerical options include: -More details are available in [`docs/METHODOLOGY.md`](docs/METHODOLOGY.md). +```text +--maxIter +--poisson-maxIter +--tol-velocity +--tol-divergence +--tol-divergence-l2 +--poisson-tol +--alpha-u +--alpha-p +--cfl +--dt-max +--min-iterations +--consecutive-passes +``` -## Running the project +Run `bin/lid_cavity --help` for the complete interface. -On Linux, WSL, or a cluster: +## CMake and tests ```bash -bash scripts/run_smoke_test.sh # compile and run a tiny check -bash scripts/run_single.sh # N=128, Re=1000 example -bash scripts/run_quick.sh # reduced study -bash scripts/run_medium.sh # medium study -bash scripts/run_full.sh # complete 36-case configuration +cmake -S . -B build -DCMAKE_BUILD_TYPE=Release +cmake --build build --parallel +ctest --test-dir build --output-on-failure ``` -The build script compiles the source to an object file, links normally on current compilers, and retries with `-lstdc++fs` only when an older GCC toolchain requires it. The same command therefore works on modern Linux systems and on GCC 8-based cluster nodes. +The CTest suite includes: -Generate the plots with: +- discrete operator verification; +- manufactured Poisson verification and grid refinement; +- convergence-state logic; +- the canonical production-solver regression. + +For sanitizer checks: ```bash -bash scripts/plot_results.sh +cmake -S . -B build/sanitized \ + -DCMAKE_BUILD_TYPE=Debug \ + -DLIDCAVITY_ENABLE_SANITIZERS=ON +cmake --build build/sanitized --parallel +ctest --test-dir build/sanitized --output-on-failure ``` -Outputs are written to: +## Output files + +Each case writes: ```text -results/data/ CSV output files -results/figures/ generated plots +results/data/study_summary_.csv +results/data/case__..._history.csv +results/data/case__..._fields.csv ``` -## Stromboli smoke test — 20 July 2026 +The summary separates: -The repository was compiled and executed successfully on the Stromboli HPC cluster with GCC 8.5.0 after adding the portable filesystem-link fallback. +- execution status; +- iterative convergence; +- pressure convergence; +- divergence metrics; +- runtime; +- Ghia benchmark agreement. -The smoke configuration was intentionally tiny: +The history file stores iteration-by-iteration convergence information. The field file contains cell-center coordinates, velocity, pressure, speed, and vorticity. -| Setting | Value | -|---|---:| -| Grid | `N = 16` | -| Reynolds number | `100` | -| Convection scheme | upwind | -| Pressure solver | RBGS | -| Outer iterations | `20` | -| Runtime | approximately `0.01 s` | +## Ghia benchmark comparison -The smoke run reached the configured `maxIter` limit and did **not** meet the Ghia validation thresholds. That is expected for this deliberately short case. Its purpose is only to verify compilation, execution, argument handling, and CSV output—not numerical convergence. +For `Re=100`, `400`, and `1000`, the solver compares: -The archived log and generated smoke-test data are stored in [`results/stromboli_2026-07-20`](results/stromboli_2026-07-20). +- horizontal velocity `u(y)` on the vertical centerline; +- vertical velocity `v(x)` on the horizontal centerline. -## Continuous integration - -The GitHub Actions workflow runs the existing `scripts/run_smoke_test.sh` path, then verifies: - -- the C++ executable was created -- the smoke-study summary contains exactly one `N = 16`, `Re = 100` case -- at least one convergence-history CSV was generated - -This is a fast build-and-execution check. It does not claim that the full 36-case study is rerun or numerically validated on every commit. +The code reports `L2` and `Linf` errors. These are reference-benchmark comparisons, not experimental validation. ## Repository structure ```text -src/ C++ solver -scripts/ build, run, plot, and clean scripts -postprocess/ Python plotting scripts -assets/ selected README figures -docs/ methodology, running notes, validation, and results -results/data/ full-study CSV output -results/figures/ full-study generated plots -results/stromboli_2026-07-20/ archived HPC smoke test -.github/ build-and-smoke GitHub Actions workflow +src/lid_cavity.cpp production staggered-grid solver +include/lidcavity/ reusable verification interfaces +src/verification/ operator, Poisson, and convergence components +tests/ CTest verification programs +scripts/ build and run helpers +postprocess/ Python plotting scripts +docs/ method and verification notes +results/data/ generated CSV output +.github/workflows/ smoke and numerical-verification CI ``` -## Requirements - -Solver: - -```text -g++ with C++17 support -``` - -Post-processing: - -```bash -python3 -m pip install -r requirements.txt -``` - -WSL is recommended on Windows because the scripts use a Linux-style terminal workflow. - ## Scope and limitations -This completed project records the implemented solver and study as they were configured. Known numerical limitations include: +This is an educational CFD and numerical-verification project, not a production CFD package. + +Current limitations include: -- collocated grid without Rhie-Chow interpolation -- no multigrid pressure solver -- configured maximum-iteration termination in the full study -- high-Reynolds-number cases that need stronger convergence control -- no formal grid-convergence or uncertainty study +- serial CPU execution only; +- explicit pseudo-time momentum advancement; +- no multigrid pressure solver; +- no formal experimental validation; +- the exhaustive 36-case study can be computationally expensive, especially with RBGS and `N=128`. -The natural next research step is to build a stricter verification and convergence protocol around these documented limitations. +Parallel C++, MPI, OpenMP, CUDA, MATLAB, and Python comparisons belong to the separate solver-comparison project and are intentionally not developed in this repository. ## 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/docs/METHODOLOGY.md b/docs/METHODOLOGY.md index 024f05b..421b87f 100644 --- a/docs/METHODOLOGY.md +++ b/docs/METHODOLOGY.md @@ -1,59 +1,115 @@ -# Methodology - -This document explains the numerical workflow used by the completed C++ lid-driven cavity solver. +# Numerical Methodology ## Problem definition -The solver models the classical two-dimensional square lid-driven cavity. The domain is a unit square. The top wall moves with a nondimensional horizontal velocity of `U_lid = 1`, while the left, right, and bottom walls are stationary. No-slip velocity boundary conditions are applied on all walls. - -The Reynolds number is controlled through the kinematic viscosity: +The solver models the classical two-dimensional incompressible lid-driven cavity in a unit square. The top wall moves with nondimensional velocity `U = 1`; the other walls are stationary. The Reynolds number controls the kinematic viscosity: ```text -nu = U_lid * L / Re +nu = U * L / Re ``` -with `L = 1` and `U_lid = 1`. +with `U = 1` and `L = 1`. + +## Staggered-grid arrangement + +The production Phase 2 solver uses a Marker-and-Cell staggered grid: + +- pressure is stored at cell centers; +- horizontal velocity is stored on vertical cell faces; +- vertical velocity is stored on horizontal cell faces. + +This arrangement avoids the pressure checkerboarding problem of the earlier collocated prototype and makes the pressure gradient, velocity correction, and discrete divergence naturally compatible. + +## Projection workflow + +Each outer pseudo-time iteration performs: + +1. calculate a stable pseudo-time step from convection and diffusion limits; +2. predict face velocities from convection, diffusion, and the current pressure field; +3. calculate cell-centered divergence of the predicted velocity; +4. solve the pressure-correction Poisson equation; +5. correct face velocities with pressure-correction gradients; +6. update and normalize pressure; +7. calculate velocity-update, divergence, mass-balance, and pressure metrics; +8. update the convergence state. + +The discrete projection is arranged so that the divergence and pressure-gradient operators compose into the same Laplacian used in the Poisson solve. + +## Momentum discretization + +The code supports: + +- first-order upwind convection; +- second-order central convection; +- second-order central diffusion. -## Numerical model +Tangential no-slip wall conditions are imposed through ghost values. Normal wall velocities are fixed directly on the staggered boundary faces. -The code solves the incompressible Navier-Stokes equations in nondimensional form using a pseudo-transient pressure-correction workflow: +## Pressure Poisson equation -1. initialize velocity and pressure -2. apply lid and wall boundary conditions -3. predict the velocity field -4. solve the pressure-correction Poisson equation -5. correct velocity and pressure -6. record residuals and validation metrics -7. repeat until the iteration limit or stopping criteria are reached +The pressure-correction equation is solved with: -## Spatial discretization +- red-black Gauss-Seidel (`RBGS`); +- red-black successive over-relaxation (`RBSOR`). -The domain is discretized on a structured Cartesian grid. The C++ version uses a collocated storage layout and manual indexing through a flat `std::vector` container. +Homogeneous normal pressure-gradient conditions are represented by the boundary stencil. The right-hand side is projected to zero mean, and the pressure field is normalized to remove the constant null space. -The study supports two convection schemes: +Pressure convergence is measured with a true equation residual. An outer case cannot report `converged` while the pressure solve is failing. -- `upwind`: more dissipative but more stable -- `central`: less dissipative and generally more accurate for the benchmark profiles +## Convergence definition -Diffusion terms are approximated with standard second-order finite differences. +The solver records separate dimensionless quantities: -## Pressure correction +- velocity-update `Linf` residual; +- divergence `Linf` residual; +- divergence `L2` residual; +- global boundary mass imbalance; +- pressure-Poisson relative residual. -The pressure-correction equation is solved iteratively using: +A case reports `converged` only after all configured criteria pass for a required number of consecutive outer iterations and after a minimum iteration count. -- `RBGS`: red-black Gauss-Seidel -- `RBSOR`: red-black Successive Over-Relaxation +Terminal states are: -The recorded study shows that RBSOR substantially reduces the average number of Poisson iterations compared with RBGS. +```text +converged +max_iterations +pressure_not_converged +stagnated +diverged +non_finite +``` + +## Continuation + +Parameter-study modes reuse a converged solution at a lower Reynolds number as the initial state for the next Reynolds number when the grid, convection scheme, and pressure solver are unchanged. This improves the stability and efficiency of the `Re=400` and `Re=1000` cases. + +## Verification + +The repository includes three levels of numerical checking: + +### Operator verification -## Time-step logic +Analytical fields test the gradient, divergence, and Laplacian operators and their discrete compatibility. -The solver uses a pseudo-time step based on convective and diffusive restrictions. This keeps the update conservative enough for the tested Reynolds numbers while allowing the same setup to run across several meshes. +### Poisson verification + +A manufactured Poisson problem is solved on successively refined grids. The tests check error reduction and agreement between RBGS and RBSOR. + +### Production regression + +CTest runs the canonical cavity case: + +```text +N = 32 +Re = 100 +scheme = upwind +pressure solver = RBSOR +``` -## Validation +The regression fails unless the executable reports convergence. -After each case, the code compares the computed centerline velocity profiles with the benchmark data from Ghia et al. The reported values are practical `L2` and `Linf` errors for comparing cases. They do not replace a formal verification or uncertainty study. +## Ghia comparison -## Implementation notes +After each supported Reynolds-number case, the cell-centered solution is interpolated onto the vertical and horizontal centerlines and compared with Ghia et al. The output includes `L2` and `Linf` errors for `u(y)` and `v(x)`. -The code is a completed serial C++17 baseline focused on clarity and reproducibility rather than maximum performance. Accelerated implementations are developed separately in the broader work-in-progress solver-comparison repository so that numerical behavior, runtime, and scalability can be compared transparently. +This is a reference benchmark comparison. It is not experimental validation. diff --git a/docs/PHASE2_LOCAL_TESTS.md b/docs/PHASE2_LOCAL_TESTS.md new file mode 100644 index 0000000..6b9149e --- /dev/null +++ b/docs/PHASE2_LOCAL_TESTS.md @@ -0,0 +1,43 @@ +# Phase 2 Local Test Record + +The following regression sets were run during the Phase 2 integration work. + +## Canonical strict case + +```text +N = 32 +Re = 100 +scheme = upwind +pressure = RBSOR +status = converged +iterations = 2086 +velocity-update Linf = 8.95e-09 +divergence Linf = 5.16e-13 +Ghia u L2 = 1.13e-02 +Ghia v L2 = 8.92e-03 +``` + +## Medium study + +All six `N=32` RBSOR cases converged: + +| Re | Scheme | Iterations | Ghia u L2 | Ghia v L2 | +|---:|---|---:|---:|---:| +| 100 | upwind | 1686 | 0.0113 | 0.0089 | +| 400 | upwind | 2430 | 0.0784 | 0.1020 | +| 1000 | upwind | 4627 | 0.1427 | 0.1956 | +| 100 | central | 1974 | 0.0037 | 0.0029 | +| 400 | central | 3924 | 0.0355 | 0.0439 | +| 1000 | central | 9454 | 0.0842 | 0.0976 | + +## Grid sequence + +The `Re=100`, central, RBSOR sequence converged: + +| N | Iterations | Ghia u L2 | Ghia v L2 | +|---:|---:|---:|---:| +| 16 | 1967 | 0.0205 | 0.0160 | +| 32 | 1974 | 0.0037 | 0.0029 | +| 64 | 3172 | 0.0013 | 0.0037 | + +These timings and iteration counts are configuration- and hardware-dependent. They are recorded as development evidence, while GitHub Actions provides the repeatable acceptance checks. diff --git a/docs/PHASE2_VERIFICATION.md b/docs/PHASE2_VERIFICATION.md index 97fb2ea..b7b9466 100644 --- a/docs/PHASE2_VERIFICATION.md +++ b/docs/PHASE2_VERIFICATION.md @@ -1,78 +1,97 @@ -# Phase 2: Verification and Convergence Contract +# Phase 2: Verification and Convergence -This phase adds a verification layer around the standalone C++ lid-driven-cavity project before the existing production solver is numerically changed. +Phase 2 is now integrated into the standalone C++ lid-driven-cavity solver. -## Why this is needed +## Main production change -The original 36-case study completed its configured executions, but those runs reached their maximum outer-iteration limits. Execution completion, benchmark-profile agreement, and iterative convergence are different results and must remain separate. +The earlier collocated pressure-correction implementation has been replaced in the production path by a staggered Marker-and-Cell solver. Pressure is stored at cell centers, horizontal velocity on vertical faces, and vertical velocity on horizontal faces. + +This gives the solver a compatible pressure-gradient, velocity-correction, and divergence arrangement and removes the pressure-velocity compatibility problem that prevented the original study from reaching strict iterative convergence. ## Convergence contract -A future cavity run may report `converged` only when all of the following hold: +A cavity run reports `converged` only when all of the following hold: -- the velocity-update Linf residual is below its dimensionless tolerance; -- the local divergence Linf residual is below its dimensionless tolerance; -- the local divergence L2 residual is below its dimensionless tolerance; -- the integrated global mass imbalance is below its tolerance; +- dimensionless velocity-update `Linf` residual is below tolerance; +- dimensionless divergence `Linf` residual is below tolerance; +- dimensionless divergence `L2` residual is below tolerance; +- global boundary mass imbalance is below tolerance; - the pressure equation converged for the current outer iteration; -- all conditions remain satisfied for a configured number of consecutive iterations; -- all fields and metrics are finite. +- all conditions remain satisfied for the configured number of consecutive iterations; +- all fields and metrics remain finite. The explicit solver statuses are: -- `running` -- `converged` -- `max_iterations` -- `pressure_not_converged` -- `stagnated` -- `diverged` -- `non_finite` +```text +converged +max_iterations +pressure_not_converged +stagnated +diverged +non_finite +``` -The `max_iterations` status must never be interpreted as convergence. +`max_iterations` is never treated as convergence. -## Verification tests added +## Verification tests ### Operator compatibility -The verification library uses a forward pressure gradient and backward divergence pair. In the interior, their composition is checked against the standard five-point Laplacian: +The verification library checks analytical gradient, divergence, and Laplacian fields, including the discrete compatibility relation: ```text D(G(phi)) = L(phi) ``` -The test also checks constant-field gradients, zero-field divergence, and non-finite-value detection. - ### Manufactured Poisson solution -The Poisson verification problem uses +The independent Poisson verification uses: ```text phi(x,y) = sin(pi x) sin(pi y) +laplacian(phi) = -2 pi^2 sin(pi x) sin(pi y) ``` -with homogeneous Dirichlet boundaries and - -```text -laplacian(phi) = -2 pi^2 sin(pi x) sin(pi y). -``` +The tests cover: -The tests verify: - -- convergence on 17x17, 33x33, and 65x65 grids; -- approximately second-order spatial convergence; -- agreement between RBGS and RBSOR solutions; +- `17x17`, `33x33`, and `65x65` grids; +- approximately second-order error reduction; +- RBGS/RBSOR agreement; - true equation-residual reduction. ### Convergence-state logic -The convergence tracker is tested independently for: +The convergence tracker is tested for: - minimum-iteration protection; - consecutive-pass requirements; -- repeated pressure-solver failures; -- non-finite residuals; +- pressure-solver failure handling; +- non-finite residual handling; - maximum-iteration termination. +### Production regression + +CTest now runs the actual production executable for: + +```text +N = 32 +Re = 100 +scheme = upwind +pressure solver = RBSOR +``` + +The test uses `--strict` and fails unless the case converges. + +## Verified run sets + +The following sets were exercised during Phase 2 development: + +- canonical `N=32`, `Re=100`, upwind, RBSOR; +- six `N=32` cases using `Re=100`, `400`, and `1000`, with upwind and central schemes; +- `N=16`, `32`, and `64` at `Re=100`, central, RBSOR. + +All of these cases reached the configured iterative convergence criteria and passed their selected Ghia centerline thresholds. + ## Running the checks ```bash @@ -87,6 +106,18 @@ cmake --build build/phase2-verification --parallel ctest --test-dir build/phase2-verification --output-on-failure ``` -## Next integration step +## Running the solver + +```bash +bash scripts/run_single.sh +bash scripts/run_medium.sh +bash scripts/run_grid.sh +``` + +The complete 36-case configuration remains available through: + +```bash +bash scripts/run_full.sh +``` -The verification library intentionally does not silently change the existing full-study results. The next step is to integrate the convergence tracker and dimensionless residual definitions into `src/lid_cavity.cpp`, then tune one canonical case (`N=32`, `Re=100`, upwind, RBSOR) before regenerating the 36-case study. +That run contains slower RBGS and `N=128` cases and is intended for a workstation or HPC node. diff --git a/scripts/run_canonical_regression.sh b/scripts/run_canonical_regression.sh new file mode 100644 index 0000000..a44cfe2 --- /dev/null +++ b/scripts/run_canonical_regression.sh @@ -0,0 +1,42 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" +rm -f "$ROOT_DIR/results/data/study_summary_single.csv" + +"$ROOT_DIR/bin/lid_cavity" \ + --single \ + --N 32 \ + --Re 100 \ + --scheme upwind \ + --pressure RBSOR \ + --strict \ + --no-fields + +python3 - "$ROOT_DIR/results/data/study_summary_single.csv" <<'PY' +import csv +import math +import sys + +path = sys.argv[1] +rows = list(csv.DictReader(open(path, encoding="utf-8"))) +if len(rows) != 1: + raise SystemExit(f"Expected one canonical result, found {len(rows)}") +row = rows[0] +if row["Status"] != "converged": + raise SystemExit(f"Canonical case did not converge: {row['Status']}") +if row["ValidationPass"] != "1": + raise SystemExit("Canonical case failed the Ghia benchmark threshold") +if float(row["FinalVelocityLinf"]) > 1.0e-7: + raise SystemExit("Velocity convergence regression") +if float(row["FinalRcDiv"]) > 1.0e-9: + raise SystemExit("Divergence convergence regression") +if int(row["FailedPressureSolves"]) != 0: + raise SystemExit("Pressure solver regression") +for key in ("Runtime_s", "Ghia_u_L2", "Ghia_v_L2"): + if not math.isfinite(float(row[key])): + raise SystemExit(f"Non-finite canonical metric: {key}") +print("Canonical staggered-grid regression passed.") +PY diff --git a/scripts/run_full.sh b/scripts/run_full.sh index c5a5113..2e9b08e 100644 --- a/scripts/run_full.sh +++ b/scripts/run_full.sh @@ -2,8 +2,8 @@ set -euo pipefail ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" -"$ROOT_DIR/scripts/build.sh" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" -# Full C++ study using the same MATLAB meshes/Re/schemes/pressure solvers, with one serial_cpp implementation. -# This can take a long time depending on your CPU. -"$ROOT_DIR/bin/lid_cavity" --mode full +# Complete 36-case convergence study. RBGS and the largest grids are slower. +"$ROOT_DIR/bin/lid_cavity" --mode full --strict diff --git a/scripts/run_grid.sh b/scripts/run_grid.sh new file mode 100644 index 0000000..0ef4bcd --- /dev/null +++ b/scripts/run_grid.sh @@ -0,0 +1,9 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" + +# Re=100 central-difference grid sequence: N=16, 32, 64. +"$ROOT_DIR/bin/lid_cavity" --mode grid --strict diff --git a/scripts/run_medium.sh b/scripts/run_medium.sh index 61e94d9..a216d5e 100644 --- a/scripts/run_medium.sh +++ b/scripts/run_medium.sh @@ -2,8 +2,8 @@ set -euo pipefail ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" -"$ROOT_DIR/scripts/build.sh" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" -# Medium C++ study matching MATLAB meshes/Re choices: -# meshes [32,64], Re [100,400,1000]. -"$ROOT_DIR/bin/lid_cavity" --mode medium +# Six converged N=32 cases at Re=100, 400, and 1000 using upwind and central schemes. +"$ROOT_DIR/bin/lid_cavity" --mode medium --strict diff --git a/scripts/run_quick.sh b/scripts/run_quick.sh index bdbe1e1..481ce0d 100644 --- a/scripts/run_quick.sh +++ b/scripts/run_quick.sh @@ -2,9 +2,8 @@ set -euo pipefail ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" -"$ROOT_DIR/scripts/build.sh" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" -# Same reduced study as MATLAB main_quick.m: -# meshes [32,64], Re [100,400], schemes [upwind,central], -# pressure solvers [RBGS,RBSOR], one C++ implementation [serial_cpp]. -"$ROOT_DIR/bin/lid_cavity" --mode quick +# Four fast converged cases: N=24/32, Re=100, upwind/central, RBSOR. +"$ROOT_DIR/bin/lid_cavity" --mode quick --strict diff --git a/scripts/run_re1000.sh b/scripts/run_re1000.sh new file mode 100644 index 0000000..effab35 --- /dev/null +++ b/scripts/run_re1000.sh @@ -0,0 +1,16 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" + +# Converged higher-Reynolds-number representative case. +"$ROOT_DIR/bin/lid_cavity" \ + --single \ + --N 32 \ + --Re 1000 \ + --scheme central \ + --pressure RBSOR \ + --tol-velocity 1e-7 \ + --strict diff --git a/scripts/run_single.sh b/scripts/run_single.sh index 1ed7d2a..28b41b3 100644 --- a/scripts/run_single.sh +++ b/scripts/run_single.sh @@ -2,8 +2,14 @@ set -euo pipefail ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" -"$ROOT_DIR/scripts/build.sh" +cd "$ROOT_DIR" +bash "$ROOT_DIR/scripts/build.sh" -# Representative README case: -# N=128, Re=1000, central differencing, RBSOR pressure solver, serial C++ implementation. -"$ROOT_DIR/bin/lid_cavity" --single --N 128 --Re 1000 --scheme central --pressure RBSOR --implementation serial_cpp +# Canonical Phase 2 regression case. It should finish with status=converged. +"$ROOT_DIR/bin/lid_cavity" \ + --single \ + --N 32 \ + --Re 100 \ + --scheme upwind \ + --pressure RBSOR \ + --strict diff --git a/src/lid_cavity.cpp b/src/lid_cavity.cpp index 2cea181..ac83e94 100644 --- a/src/lid_cavity.cpp +++ b/src/lid_cavity.cpp @@ -1,926 +1,1177 @@ #include #include -#include #include #include #include #include #include #include +#include #include #include #include #include #include +#include #include namespace fs = std::filesystem; - constexpr double PI = 3.141592653589793238462643383279502884; -struct Matrix { - int n{}; - std::vector a; +struct Field { + int rows{}; + int cols{}; + std::vector values; - Matrix() = default; - explicit Matrix(int n_, double value = 0.0) : n(n_), a(static_cast(n_) * n_, value) {} + Field() = default; + Field(int r, int c, double value = 0.0) + : rows(r), cols(c), values(static_cast(r) * static_cast(c), value) {} - double& operator()(int i, int j) { return a[static_cast(i) * n + j]; } - double operator()(int i, int j) const { return a[static_cast(i) * n + j]; } + double& operator()(int i, int j) { + return values[static_cast(i) * static_cast(cols) + static_cast(j)]; + } + double operator()(int i, int j) const { + return values[static_cast(i) * static_cast(cols) + static_cast(j)]; + } }; struct Config { - double U_lid = 1.0; - double L = 1.0; + double lid_velocity = 1.0; + double length = 1.0; - int maxIter = 4000; - int maxIter_N128_bonus = 3000; - int maxIter_Re1000_bonus = 3000; - int maxIter_central_bonus = 1500; + int max_iterations = 30000; + int minimum_iterations = 200; + int consecutive_passes = 20; + int maximum_pressure_failures = 3; - double tol_mass = 1e-7; - double tol_divergence = 2e-3; - double tol_velocity = 5e-7; + double velocity_tolerance = 1e-8; + double divergence_linf_tolerance = 1e-9; + double divergence_l2_tolerance = 2e-10; double diverged_limit = 1e6; - double cfl = 0.25; - double dt_max = 0.0025; - double dt_min = 1e-6; - - double alpha_u = 0.55; - double alpha_p = 0.20; - - int poisson_maxIter = 2500; - double poisson_tol_abs = 1e-8; - double poisson_tol_rel = 1e-4; - int poisson_check_every = 25; - - std::string sor_omega = "auto"; - double sor_omega_min = 1.15; - double sor_omega_max = 1.90; + double cfl = 0.60; + double dt_max = 0.01; + double dt_min = 1e-8; + double momentum_relaxation = 0.90; + double pressure_relaxation = 1.0; - std::vector meshes = {32, 64, 128}; - std::vector re_list = {100, 400, 1000}; - std::vector schemes = {"upwind", "central"}; - std::vector pressure_solvers = {"RBGS", "RBSOR"}; - std::vector implementations = {"serial_cpp"}; + int poisson_max_iterations = 5000; + int poisson_check_every = 20; + double poisson_absolute_tolerance = 1e-10; + double poisson_relative_tolerance = 1e-9; + double sor_omega = 0.0; - double validation_u_L2_limit_Re100 = 0.030; - double validation_v_L2_limit_Re100 = 0.030; - double validation_u_L2_limit_Re400 = 0.090; - double validation_v_L2_limit_Re400 = 0.120; - double validation_u_L2_limit_Re1000 = 0.160; - double validation_v_L2_limit_Re1000 = 0.180; + int stagnation_window = 1500; + double stagnation_minimum_reduction = 0.005; bool save_fields = true; - std::string results_dir = "results"; - std::string data_dir = "results/data"; + bool strict_exit = false; + std::string data_directory = "results/data"; }; struct PoissonInfo { - int iter = 0; + int iterations = 0; bool converged = false; - double final_true_residual = std::numeric_limits::infinity(); - double final_relative_residual = std::numeric_limits::infinity(); - double final_change = std::numeric_limits::infinity(); + double absolute_residual = std::numeric_limits::infinity(); + double relative_residual = std::numeric_limits::infinity(); double omega = 1.0; - std::vector residual_history; - std::vector change_history; +}; + +struct Residuals { + double velocity_update_linf = 0.0; + double divergence_linf = 0.0; + double divergence_l2 = 0.0; + double global_mass_imbalance = 0.0; +}; + +struct GhiaMetrics { + bool available = false; + bool passed = false; + double u_l2 = std::numeric_limits::quiet_NaN(); + double v_l2 = std::numeric_limits::quiet_NaN(); + double u_linf = std::numeric_limits::quiet_NaN(); + double v_linf = std::numeric_limits::quiet_NaN(); + double u_limit = std::numeric_limits::quiet_NaN(); + double v_limit = std::numeric_limits::quiet_NaN(); }; struct Result { - int N = 0; - int Re = 0; + int case_id = 0; + int cells = 0; + int reynolds = 0; std::string scheme; std::string pressure_solver; - std::string implementation; - - std::vector x, y; - Matrix u, v, p, speed, vorticity; - - std::vector Ru, Rv, Rc_mass, Rc_div, dt, poisson_relative_residual; - std::vector poisson_iters; - std::vector poisson_converged; + std::string status = "max_iterations"; + std::string quality = "needs_improvement"; int iterations = 0; - int localMaxIter = 0; - double runtime = 0.0; - std::string status = "maxIter"; - double final_Ru = 0.0; - double final_Rv = 0.0; - double final_Rc_mass = 0.0; - double final_Rc_div = 0.0; - double avg_poisson_iters = 0.0; - double avg_poisson_relative_residual = 0.0; + int local_max_iterations = 0; + int failed_pressure_solves = 0; + int consecutive_pass_count = 0; + double runtime_seconds = 0.0; + + Field u_face; + Field v_face; + Field pressure; + Field u_center; + Field v_center; + Field speed; + Field vorticity; + std::vector x; + std::vector y; + + std::vector velocity_history; + std::vector divergence_linf_history; + std::vector divergence_l2_history; + std::vector mass_history; + std::vector dt_history; + std::vector poisson_relative_history; + std::vector poisson_iteration_history; + std::vector poisson_converged_history; + + Residuals final_residuals; + double average_poisson_iterations = 0.0; + double average_poisson_relative_residual = 0.0; double pressure_saturation_ratio = 0.0; - int stagnation_counter = 0; + GhiaMetrics ghia; }; -struct GhiaData { - int Re = 0; - std::vector y_u, u, x_v, v; +struct InitialState { + bool available = false; + Field u_face; + Field v_face; + Field pressure; }; -struct Metrics { - bool available = false; - bool pass = false; - double u_L2 = std::numeric_limits::quiet_NaN(); - double v_L2 = std::numeric_limits::quiet_NaN(); - double u_Linf = std::numeric_limits::quiet_NaN(); - double v_Linf = std::numeric_limits::quiet_NaN(); - double u_limit = std::numeric_limits::quiet_NaN(); - double v_limit = std::numeric_limits::quiet_NaN(); +struct GhiaData { + std::vector y_u; + std::vector u; + std::vector x_v; + std::vector v; }; -static std::string lower(std::string s) { - std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return static_cast(std::tolower(c)); }); - return s; +static std::string lower(std::string value) { + std::transform(value.begin(), value.end(), value.begin(), [](unsigned char c) { + return static_cast(std::tolower(c)); + }); + return value; } -static std::string upper(std::string s) { - std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return static_cast(std::toupper(c)); }); - return s; +static std::string upper(std::string value) { + std::transform(value.begin(), value.end(), value.begin(), [](unsigned char c) { + return static_cast(std::toupper(c)); + }); + return value; } -static std::string normalize_implementation(std::string s) { - s = lower(s); - // MATLAB has two implementations: vectorized and loop. This C++ package - // intentionally contains one honest serial loop kernel. For convenience, - // MATLAB labels are accepted as aliases, but CSV output uses serial_cpp. - if (s == "serial" || s == "serial_cpp" || s == "cpp" || s == "loop" || s == "vectorized") { - return "serial_cpp"; +static double maximum_absolute(const Field& field) { + double result = 0.0; + for (const double value : field.values) { + result = std::max(result, std::abs(value)); } - throw std::runtime_error("Unknown implementation: " + s + " (use serial_cpp)"); -} - -static double max_abs(const Matrix& m) { - double r = 0.0; - for (double v : m.a) r = std::max(r, std::abs(v)); - return r; -} - -static double mean_all(const Matrix& m) { - if (m.a.empty()) return 0.0; - return std::accumulate(m.a.begin(), m.a.end(), 0.0) / static_cast(m.a.size()); + return result; } -static bool all_finite(const Matrix& m) { - for (double x : m.a) { - if (!std::isfinite(x)) return false; - } - return true; +static bool all_finite(const Field& field) { + return std::all_of(field.values.begin(), field.values.end(), [](double value) { + return std::isfinite(value); + }); } -static void apply_lid_bc(Matrix& u, Matrix& v, double U_lid) { - const int N = u.n; - for (int j = 0; j < N; ++j) { - u(0, j) = 0.0; // bottom wall - u(N - 1, j) = U_lid; // moving lid +static void remove_mean(Field& field) { + if (field.values.empty()) { + return; } - for (int i = 0; i < N; ++i) { - u(i, 0) = 0.0; // left wall; also sets top-left corner to zero - u(i, N - 1) = 0.0; // right wall; also sets top-right corner to zero + const double mean = std::accumulate(field.values.begin(), field.values.end(), 0.0) + / static_cast(field.values.size()); + for (double& value : field.values) { + value -= mean; } +} - for (int j = 0; j < N; ++j) { - v(0, j) = 0.0; - v(N - 1, j) = 0.0; +static double u_with_wall_ghost(const Field& u, int i, int j, int cells, double lid_velocity) { + if (i < 0) { + return -u(0, j); } - for (int i = 0; i < N; ++i) { - v(i, 0) = 0.0; - v(i, N - 1) = 0.0; + if (i >= cells) { + return 2.0 * lid_velocity - u(cells - 1, j); } + return u(i, j); } -static void apply_pressure_bc(Matrix& p) { - const int N = p.n; - for (int i = 0; i < N; ++i) { - p(i, 0) = p(i, 1); - p(i, N - 1) = p(i, N - 2); +static double v_with_wall_ghost(const Field& v, int i, int j, int cells) { + if (j < 0) { + return -v(i, 0); } - for (int j = 0; j < N; ++j) { - p(0, j) = p(1, j); - p(N - 1, j) = p(N - 2, j); + if (j >= cells) { + return -v(i, cells - 1); } - p(0, 0) = 0.0; -} - -static double compute_dt(const Matrix& u, const Matrix& v, double dx, double dy, double nu, const Config& cfg) { - const double max_vel = std::max({max_abs(u), max_abs(v), cfg.U_lid, 1e-12}); - const double h = std::min(dx, dy); - const double dt_conv = cfg.cfl * h / max_vel; - const double dt_diff = (nu > 0.0) ? 0.25 * h * h / nu : std::numeric_limits::infinity(); - return std::min({dt_conv, dt_diff, cfg.dt_max}); + return v(i, j); } -static Matrix divergence_field(const Matrix& u, const Matrix& v, double dx, double dy) { - const int N = u.n; - Matrix div(N, 0.0); - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - div(i, j) = (u(i, j + 1) - u(i, j - 1)) / (2.0 * dx) - + (v(i + 1, j) - v(i - 1, j)) / (2.0 * dy); +static Field divergence(const Field& u, const Field& v, double dx, double dy, int cells) { + Field result(cells, cells, 0.0); + for (int i = 0; i < cells; ++i) { + for (int j = 0; j < cells; ++j) { + result(i, j) = (u(i, j + 1) - u(i, j)) / dx + + (v(i + 1, j) - v(i, j)) / dy; } } - return div; + return result; } -static Matrix compute_vorticity(const Matrix& u, const Matrix& v, double dx, double dy) { - const int N = u.n; - Matrix omega(N, 0.0); - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - omega(i, j) = (v(i, j + 1) - v(i, j - 1)) / (2.0 * dx) - - (u(i + 1, j) - u(i - 1, j)) / (2.0 * dy); +static double poisson_residual(const Field& pressure, const Field& rhs, double h) { + const int cells = pressure.rows; + double maximum = 0.0; + for (int i = 0; i < cells; ++i) { + for (int j = 0; j < cells; ++j) { + double laplacian = 0.0; + if (i > 0) { + laplacian += pressure(i - 1, j) - pressure(i, j); + } + if (i + 1 < cells) { + laplacian += pressure(i + 1, j) - pressure(i, j); + } + if (j > 0) { + laplacian += pressure(i, j - 1) - pressure(i, j); + } + if (j + 1 < cells) { + laplacian += pressure(i, j + 1) - pressure(i, j); + } + laplacian /= h * h; + maximum = std::max(maximum, std::abs(laplacian - rhs(i, j))); } } - return omega; + return maximum; } -static std::tuple momentum_predictor( - const Matrix& u, const Matrix& v, const Matrix& p, - int Re, const std::string& scheme_in, const Config& cfg +static std::pair solve_pressure_poisson( + Field rhs, + double h, + const std::string& solver_name, + const Config& config ) { - const int N = u.n; - const double dx = cfg.L / static_cast(N - 1); - const double dy = dx; - const double nu = cfg.U_lid * cfg.L / static_cast(Re); - const double dt = compute_dt(u, v, dx, dy, nu, cfg); - const std::string scheme = lower(scheme_in); + const int cells = rhs.rows; + const std::string method = upper(solver_name); + const bool use_sor = method == "RBSOR"; + if (!use_sor && method != "RBGS") { + throw std::runtime_error("Pressure solver must be RBGS or RBSOR"); + } - Matrix u_star = u; - Matrix v_star = v; + const double rhs_mean = std::accumulate(rhs.values.begin(), rhs.values.end(), 0.0) + / static_cast(rhs.values.size()); + double rhs_norm = 0.0; + for (double& value : rhs.values) { + value -= rhs_mean; + rhs_norm = std::max(rhs_norm, std::abs(value)); + } + rhs_norm = std::max(rhs_norm, 1e-30); - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - const double uC = u(i, j); - const double vC = v(i, j); + Field pressure(cells, cells, 0.0); + double omega = use_sor ? config.sor_omega : 1.0; + if (use_sor && omega <= 0.0) { + omega = 2.0 / (1.0 + std::sin(PI / static_cast(cells))); + omega = std::clamp(omega, 1.0, 1.95); + } - const double lap_u = (u(i, j + 1) - 2.0 * u(i, j) + u(i, j - 1)) / (dx * dx) - + (u(i + 1, j) - 2.0 * u(i, j) + u(i - 1, j)) / (dy * dy); - const double lap_v = (v(i, j + 1) - 2.0 * v(i, j) + v(i, j - 1)) / (dx * dx) - + (v(i + 1, j) - 2.0 * v(i, j) + v(i - 1, j)) / (dy * dy); + PoissonInfo info; + info.omega = omega; - double du_dx, du_dy, dv_dx, dv_dy; + for (int iteration = 1; iteration <= config.poisson_max_iterations; ++iteration) { + for (int color = 0; color < 2; ++color) { + for (int i = 0; i < cells; ++i) { + for (int j = 0; j < cells; ++j) { + if (((i + j) & 1) != color) { + continue; + } - if (scheme == "central") { - du_dx = (u(i, j + 1) - u(i, j - 1)) / (2.0 * dx); - du_dy = (u(i + 1, j) - u(i - 1, j)) / (2.0 * dy); - dv_dx = (v(i, j + 1) - v(i, j - 1)) / (2.0 * dx); - dv_dy = (v(i + 1, j) - v(i - 1, j)) / (2.0 * dy); - } else if (scheme == "upwind") { - if (uC >= 0.0) { - du_dx = (u(i, j) - u(i, j - 1)) / dx; - dv_dx = (v(i, j) - v(i, j - 1)) / dx; - } else { - du_dx = (u(i, j + 1) - u(i, j)) / dx; - dv_dx = (v(i, j + 1) - v(i, j)) / dx; - } + double neighbor_sum = 0.0; + int neighbor_count = 0; + if (i > 0) { + neighbor_sum += pressure(i - 1, j); + ++neighbor_count; + } + if (i + 1 < cells) { + neighbor_sum += pressure(i + 1, j); + ++neighbor_count; + } + if (j > 0) { + neighbor_sum += pressure(i, j - 1); + ++neighbor_count; + } + if (j + 1 < cells) { + neighbor_sum += pressure(i, j + 1); + ++neighbor_count; + } - if (vC >= 0.0) { - du_dy = (u(i, j) - u(i - 1, j)) / dy; - dv_dy = (v(i, j) - v(i - 1, j)) / dy; - } else { - du_dy = (u(i + 1, j) - u(i, j)) / dy; - dv_dy = (v(i + 1, j) - v(i, j)) / dy; + const double candidate = (neighbor_sum - rhs(i, j) * h * h) + / static_cast(neighbor_count); + pressure(i, j) = use_sor + ? (1.0 - omega) * pressure(i, j) + omega * candidate + : candidate; } - } else { - throw std::runtime_error("Unknown convection scheme: " + scheme_in); } + } - const double conv_u = uC * du_dx + vC * du_dy; - const double conv_v = uC * dv_dx + vC * dv_dy; - const double dp_dx = (p(i, j + 1) - p(i, j - 1)) / (2.0 * dx); - const double dp_dy = (p(i + 1, j) - p(i - 1, j)) / (2.0 * dy); - - const double u_pred = u(i, j) + dt * (-conv_u - dp_dx + nu * lap_u); - const double v_pred = v(i, j) + dt * (-conv_v - dp_dy + nu * lap_v); + if (iteration % 50 == 0) { + remove_mean(pressure); + } - u_star(i, j) = (1.0 - cfg.alpha_u) * u(i, j) + cfg.alpha_u * u_pred; - v_star(i, j) = (1.0 - cfg.alpha_u) * v(i, j) + cfg.alpha_u * v_pred; + if (iteration == 1 + || iteration % config.poisson_check_every == 0 + || iteration == config.poisson_max_iterations) { + const double absolute = poisson_residual(pressure, rhs, h); + const double relative = absolute / rhs_norm; + info.iterations = iteration; + info.absolute_residual = absolute; + info.relative_residual = relative; + if (absolute <= config.poisson_absolute_tolerance + || relative <= config.poisson_relative_tolerance) { + info.converged = true; + remove_mean(pressure); + return {pressure, info}; + } } } - apply_lid_bc(u_star, v_star, cfg.U_lid); - return {u_star, v_star, dt}; + remove_mean(pressure); + return {pressure, info}; } -static double poisson_true_residual(const Matrix& phi, const Matrix& rhs, double dx, double dy) { - const int N = phi.n; - double res = 0.0; - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - const double lap = (phi(i, j + 1) - 2.0 * phi(i, j) + phi(i, j - 1)) / (dx * dx) - + (phi(i + 1, j) - 2.0 * phi(i, j) + phi(i - 1, j)) / (dy * dy); - res = std::max(res, std::abs(lap - rhs(i, j))); - } - } - return res; +static double compute_time_step( + const Field& u, + const Field& v, + double h, + double viscosity, + const Config& config +) { + const double maximum_velocity = std::max({ + maximum_absolute(u), + maximum_absolute(v), + config.lid_velocity, + 1e-12 + }); + const double convection_limit = config.cfl * h / maximum_velocity; + const double diffusion_limit = 0.24 * h * h / std::max(viscosity, 1e-30); + return std::clamp( + std::min({convection_limit, diffusion_limit, config.dt_max}), + config.dt_min, + config.dt_max + ); } -static std::tuple pressure_poisson( - const Matrix& rhs, double dx, double dy, const std::string& solver_type_in, const Config& cfg +static std::tuple predict_velocity( + const Field& u, + const Field& v, + const Field& pressure, + int reynolds, + const std::string& scheme_name, + const Config& config ) { - const int N = rhs.n; - Matrix phi(N, 0.0); - Matrix rhs2 = rhs; - const std::string solver_type = upper(solver_type_in); - - double interior_sum = 0.0; - int count = 0; - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - interior_sum += rhs2(i, j); - ++count; + const int cells = pressure.rows; + const double h = config.length / static_cast(cells); + const double viscosity = config.lid_velocity * config.length / static_cast(reynolds); + const double dt = compute_time_step(u, v, h, viscosity, config); + const std::string scheme = lower(scheme_name); + + Field u_star = u; + Field v_star = v; + + for (int i = 0; i < cells; ++i) { + for (int j = 1; j < cells; ++j) { + const double u_center = u(i, j); + const double v_at_u = 0.25 * ( + v(i, j - 1) + v(i + 1, j - 1) + + v(i, j) + v(i + 1, j) + ); + + const double u_west = u(i, j - 1); + const double u_east = u(i, j + 1); + const double u_south = u_with_wall_ghost(u, i - 1, j, cells, config.lid_velocity); + const double u_north = u_with_wall_ghost(u, i + 1, j, cells, config.lid_velocity); + + double du_dx = 0.0; + double du_dy = 0.0; + if (scheme == "central") { + du_dx = (u_east - u_west) / (2.0 * h); + du_dy = (u_north - u_south) / (2.0 * h); + } else if (scheme == "upwind") { + du_dx = u_center >= 0.0 + ? (u_center - u_west) / h + : (u_east - u_center) / h; + du_dy = v_at_u >= 0.0 + ? (u_center - u_south) / h + : (u_north - u_center) / h; + } else { + throw std::runtime_error("Convection scheme must be upwind or central"); + } + + const double laplacian = ( + u_east - 2.0 * u_center + u_west + + u_north - 2.0 * u_center + u_south + ) / (h * h); + const double pressure_gradient = (pressure(i, j) - pressure(i, j - 1)) / h; + + u_star(i, j) = u_center + config.momentum_relaxation * dt * ( + -u_center * du_dx - v_at_u * du_dy + - pressure_gradient + viscosity * laplacian + ); } } - const double interior_mean = (count > 0) ? interior_sum / static_cast(count) : 0.0; - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - rhs2(i, j) -= interior_mean; + + for (int i = 1; i < cells; ++i) { + for (int j = 0; j < cells; ++j) { + const double v_center = v(i, j); + const double u_at_v = 0.25 * ( + u(i - 1, j) + u(i - 1, j + 1) + + u(i, j) + u(i, j + 1) + ); + + const double v_south = v(i - 1, j); + const double v_north = v(i + 1, j); + const double v_west = v_with_wall_ghost(v, i, j - 1, cells); + const double v_east = v_with_wall_ghost(v, i, j + 1, cells); + + double dv_dx = 0.0; + double dv_dy = 0.0; + if (scheme == "central") { + dv_dx = (v_east - v_west) / (2.0 * h); + dv_dy = (v_north - v_south) / (2.0 * h); + } else if (scheme == "upwind") { + dv_dx = u_at_v >= 0.0 + ? (v_center - v_west) / h + : (v_east - v_center) / h; + dv_dy = v_center >= 0.0 + ? (v_center - v_south) / h + : (v_north - v_center) / h; + } + + const double laplacian = ( + v_east - 2.0 * v_center + v_west + + v_north - 2.0 * v_center + v_south + ) / (h * h); + const double pressure_gradient = (pressure(i, j) - pressure(i - 1, j)) / h; + + v_star(i, j) = v_center + config.momentum_relaxation * dt * ( + -u_at_v * dv_dx - v_center * dv_dy + - pressure_gradient + viscosity * laplacian + ); } } - const double den = 2.0 * (dx * dx + dy * dy); - double omega; - if (lower(cfg.sor_omega) == "auto") { - omega = 2.0 / (1.0 + std::sin(PI / static_cast(N - 1))); - omega = std::min(std::max(omega, cfg.sor_omega_min), cfg.sor_omega_max); - } else { - omega = std::stod(cfg.sor_omega); + return {u_star, v_star, dt}; +} + +static Residuals calculate_residuals( + const Field& u, + const Field& v, + const Field& old_u, + const Field& old_v, + double h, + const Config& config +) { + Residuals residuals; + for (std::size_t k = 0; k < u.values.size(); ++k) { + residuals.velocity_update_linf = std::max( + residuals.velocity_update_linf, + std::abs(u.values[k] - old_u.values[k]) / config.lid_velocity + ); + } + for (std::size_t k = 0; k < v.values.size(); ++k) { + residuals.velocity_update_linf = std::max( + residuals.velocity_update_linf, + std::abs(v.values[k] - old_v.values[k]) / config.lid_velocity + ); } - double rhs_norm = 1.0; - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - rhs_norm = std::max(rhs_norm, std::abs(rhs2(i, j))); - } + const Field div = divergence(u, v, h, h, u.rows); + double square_sum = 0.0; + for (const double value : div.values) { + const double dimensionless = value / (config.lid_velocity / config.length); + residuals.divergence_linf = std::max(residuals.divergence_linf, std::abs(dimensionless)); + square_sum += dimensionless * dimensionless; } + residuals.divergence_l2 = std::sqrt(square_sum / static_cast(div.values.size())); - PoissonInfo info; - info.omega = omega; - info.residual_history.assign(static_cast(cfg.poisson_maxIter), 0.0); - info.change_history.assign(static_cast(cfg.poisson_maxIter), 0.0); + double boundary_flux = 0.0; + const int cells = u.rows; + for (int i = 0; i < cells; ++i) { + boundary_flux += (u(i, cells) - u(i, 0)) * h; + } + for (int j = 0; j < cells; ++j) { + boundary_flux += (v(cells, j) - v(0, j)) * h; + } + residuals.global_mass_imbalance = std::abs(boundary_flux) + / (config.lid_velocity * config.length); + return residuals; +} - bool converged = false; - double final_res = std::numeric_limits::infinity(); - double final_change = std::numeric_limits::infinity(); - int it = 0; - - for (it = 1; it <= cfg.poisson_maxIter; ++it) { - Matrix phi_old = phi; - - if (solver_type == "JACOBI") { - Matrix phi_new = phi; - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - phi_new(i, j) = ((phi(i + 1, j) + phi(i - 1, j)) * dy * dy - + (phi(i, j + 1) + phi(i, j - 1)) * dx * dx - - rhs2(i, j) * dx * dx * dy * dy) / den; - } - } - phi = std::move(phi_new); - apply_pressure_bc(phi); - } else if (solver_type == "RBGS" || solver_type == "RBSOR") { - for (int color = 0; color <= 1; ++color) { - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - const bool is_red = ((i + 1) + (j + 1)) % 2 == 0; // MATLAB i+j parity - if ((color == 0 && !is_red) || (color == 1 && is_red)) continue; - - const double candidate = ((phi(i + 1, j) + phi(i - 1, j)) * dy * dy - + (phi(i, j + 1) + phi(i, j - 1)) * dx * dx - - rhs2(i, j) * dx * dx * dy * dy) / den; - if (solver_type == "RBSOR") { - phi(i, j) = (1.0 - omega) * phi(i, j) + omega * candidate; - } else { - phi(i, j) = candidate; - } - } - } - apply_pressure_bc(phi); - } - } else { - throw std::runtime_error("Unknown pressure solver: " + solver_type_in); - } +static void build_cell_center_fields(Result& result, const Config& config) { + const int cells = result.cells; + const double h = config.length / static_cast(cells); + result.u_center = Field(cells, cells, 0.0); + result.v_center = Field(cells, cells, 0.0); + result.speed = Field(cells, cells, 0.0); + result.vorticity = Field(cells, cells, 0.0); + result.x.resize(static_cast(cells)); + result.y.resize(static_cast(cells)); + + for (int k = 0; k < cells; ++k) { + result.x[static_cast(k)] = (static_cast(k) + 0.5) * h; + result.y[static_cast(k)] = (static_cast(k) + 0.5) * h; + } - final_change = 0.0; - for (size_t k = 0; k < phi.a.size(); ++k) { - final_change = std::max(final_change, std::abs(phi.a[k] - phi_old.a[k])); + for (int i = 0; i < cells; ++i) { + for (int j = 0; j < cells; ++j) { + result.u_center(i, j) = 0.5 * (result.u_face(i, j) + result.u_face(i, j + 1)); + result.v_center(i, j) = 0.5 * (result.v_face(i, j) + result.v_face(i + 1, j)); + result.speed(i, j) = std::hypot(result.u_center(i, j), result.v_center(i, j)); } - info.change_history[static_cast(it - 1)] = final_change; - - if ((it % cfg.poisson_check_every) == 0 || it == 1 || it == cfg.poisson_maxIter) { - final_res = poisson_true_residual(phi, rhs2, dx, dy); - const double rel_res = final_res / rhs_norm; - info.residual_history[static_cast(it - 1)] = rel_res; - if (final_res < cfg.poisson_tol_abs || rel_res < cfg.poisson_tol_rel) { - converged = true; - break; - } + } + + for (int i = 1; i + 1 < cells; ++i) { + for (int j = 1; j + 1 < cells; ++j) { + const double dv_dx = (result.v_center(i, j + 1) - result.v_center(i, j - 1)) / (2.0 * h); + const double du_dy = (result.u_center(i + 1, j) - result.u_center(i - 1, j)) / (2.0 * h); + result.vorticity(i, j) = dv_dx - du_dy; } } +} - if (it > cfg.poisson_maxIter) { - it = cfg.poisson_maxIter; +static double interpolate(const std::vector& coordinates, const std::vector& values, double query) { + if (coordinates.empty()) { + return std::numeric_limits::quiet_NaN(); } + if (query <= coordinates.front()) { + return values.front(); + } + if (query >= coordinates.back()) { + return values.back(); + } + const auto upper_it = std::upper_bound(coordinates.begin(), coordinates.end(), query); + const std::size_t upper_index = static_cast(std::distance(coordinates.begin(), upper_it)); + const std::size_t lower_index = upper_index - 1; + const double fraction = (query - coordinates[lower_index]) + / (coordinates[upper_index] - coordinates[lower_index]); + return values[lower_index] + fraction * (values[upper_index] - values[lower_index]); +} - if (!converged) { - final_res = poisson_true_residual(phi, rhs2, dx, dy); +static GhiaData ghia_data(int reynolds) { + GhiaData data; + data.y_u = {1.0000,0.9766,0.9688,0.9609,0.9531,0.8516,0.7344,0.6172,0.5000,0.4531,0.2813,0.1719,0.1016,0.0703,0.0625,0.0547,0.0000}; + data.x_v = {1.0000,0.9688,0.9609,0.9531,0.9453,0.9063,0.8594,0.8047,0.5000,0.2344,0.2266,0.1563,0.0938,0.0781,0.0703,0.0625,0.0000}; + if (reynolds == 100) { + data.u = {1.0000,0.84123,0.78871,0.73722,0.68717,0.23151,0.00332,-0.13641,-0.20581,-0.21090,-0.15662,-0.10150,-0.06434,-0.04775,-0.04192,-0.03717,0.0000}; + data.v = {0.0000,-0.05906,-0.07391,-0.08864,-0.10313,-0.16914,-0.22445,-0.24533,0.05454,0.17527,0.17507,0.16077,0.12317,0.10890,0.10091,0.09233,0.0000}; + } else if (reynolds == 400) { + data.u = {1.0000,0.75837,0.68439,0.61756,0.55892,0.29093,0.16256,0.02135,-0.11477,-0.17119,-0.32726,-0.24299,-0.14612,-0.10338,-0.09266,-0.08186,0.0000}; + data.v = {0.0000,-0.12146,-0.15663,-0.19254,-0.22847,-0.23827,-0.44993,-0.38598,0.05186,0.30174,0.30203,0.28124,0.22965,0.20920,0.19713,0.18360,0.0000}; + } else if (reynolds == 1000) { + data.u = {1.0000,0.65928,0.57492,0.51117,0.46604,0.33304,0.18719,0.05702,-0.06080,-0.10648,-0.27805,-0.38289,-0.29730,-0.22220,-0.20196,-0.18109,0.0000}; + data.v = {0.0000,-0.21388,-0.27669,-0.33714,-0.39188,-0.51550,-0.42665,-0.31966,0.02526,0.32235,0.33075,0.37095,0.32627,0.30353,0.29012,0.27485,0.0000}; } + return data; +} - info.iter = it; - info.converged = converged; - info.final_true_residual = final_res; - info.final_relative_residual = final_res / rhs_norm; - info.final_change = final_change; - info.residual_history.resize(static_cast(it)); - info.change_history.resize(static_cast(it)); - return {phi, info}; +static GhiaMetrics compare_with_ghia(const Result& result) { + GhiaMetrics metrics; + const GhiaData reference = ghia_data(result.reynolds); + if (reference.u.empty()) { + return metrics; + } + metrics.available = true; + + const int cells = result.cells; + const int left = (cells - 1) / 2; + const int right = cells / 2; + std::vector y_coordinates; + std::vector u_profile; + std::vector x_coordinates; + std::vector v_profile; + y_coordinates.reserve(static_cast(cells + 2)); + u_profile.reserve(static_cast(cells + 2)); + x_coordinates.reserve(static_cast(cells + 2)); + v_profile.reserve(static_cast(cells + 2)); + + y_coordinates.push_back(0.0); + u_profile.push_back(0.0); + for (int i = 0; i < cells; ++i) { + y_coordinates.push_back(result.y[static_cast(i)]); + u_profile.push_back(0.5 * (result.u_center(i, left) + result.u_center(i, right))); + } + y_coordinates.push_back(1.0); + u_profile.push_back(1.0); + + x_coordinates.push_back(0.0); + v_profile.push_back(0.0); + for (int j = 0; j < cells; ++j) { + x_coordinates.push_back(result.x[static_cast(j)]); + v_profile.push_back(0.5 * (result.v_center(left, j) + result.v_center(right, j))); + } + x_coordinates.push_back(1.0); + v_profile.push_back(0.0); + + double u_square_sum = 0.0; + double v_square_sum = 0.0; + for (std::size_t k = 0; k < reference.y_u.size(); ++k) { + const double error = interpolate(y_coordinates, u_profile, reference.y_u[k]) - reference.u[k]; + u_square_sum += error * error; + metrics.u_linf = std::isnan(metrics.u_linf) ? std::abs(error) : std::max(metrics.u_linf, std::abs(error)); + } + for (std::size_t k = 0; k < reference.x_v.size(); ++k) { + const double error = interpolate(x_coordinates, v_profile, reference.x_v[k]) - reference.v[k]; + v_square_sum += error * error; + metrics.v_linf = std::isnan(metrics.v_linf) ? std::abs(error) : std::max(metrics.v_linf, std::abs(error)); + } + metrics.u_l2 = std::sqrt(u_square_sum / static_cast(reference.y_u.size())); + metrics.v_l2 = std::sqrt(v_square_sum / static_cast(reference.x_v.size())); + + if (result.reynolds == 100) { + metrics.u_limit = 0.035; + metrics.v_limit = 0.035; + } else if (result.reynolds == 400) { + metrics.u_limit = 0.10; + metrics.v_limit = 0.13; + } else { + metrics.u_limit = 0.18; + metrics.v_limit = 0.20; + } + metrics.passed = metrics.u_l2 <= metrics.u_limit && metrics.v_l2 <= metrics.v_limit; + return metrics; } -static std::tuple velocity_residuals( - const Matrix& u, const Matrix& v, const Matrix& u_old, const Matrix& v_old, - double dx, double dy, double U, double L -) { - double Ru = 0.0; - double Rv = 0.0; - for (size_t k = 0; k < u.a.size(); ++k) { - Ru = std::max(Ru, std::abs(u.a[k] - u_old.a[k])); - Rv = std::max(Rv, std::abs(v.a[k] - v_old.a[k])); - } - - Matrix div = divergence_field(u, v, dx, dy); - const double Rc_div = max_abs(div); - const double scale = std::max(U * L, std::numeric_limits::epsilon()); - const double Rc_mass = Rc_div * dx * dy / scale; - return {Ru, Rv, Rc_mass, Rc_div}; +static std::string quality_label(const Result& result) { + if (result.status == "converged" && result.ghia.available && result.ghia.passed) { + return "converged_benchmark_pass"; + } + if (result.status == "converged" && result.ghia.available) { + return "converged_benchmark_needs_improvement"; + } + if (result.status == "converged") { + return "converged_no_benchmark"; + } + if (result.ghia.available && result.ghia.passed) { + return "benchmark_pass_not_converged"; + } + return "needs_improvement"; } -static Result solve_lid_cavity(int N, int Re, std::string scheme, std::string pressure_solver, std::string implementation, const Config& cfg) { +static Result solve_case( + int case_id, + int cells, + int reynolds, + std::string scheme, + std::string pressure_solver, + const Config& config, + const InitialState& initial = {} +) { scheme = lower(scheme); pressure_solver = upper(pressure_solver); - implementation = normalize_implementation(implementation); - - const double dx = cfg.L / static_cast(N - 1); - const double dy = dx; - - int localMaxIter = cfg.maxIter; - if (N >= 128) localMaxIter += cfg.maxIter_N128_bonus; - if (Re >= 1000) localMaxIter += cfg.maxIter_Re1000_bonus; - if (scheme == "central") localMaxIter += cfg.maxIter_central_bonus; - - Matrix u(N, 0.0), v(N, 0.0), p(N, 0.0); - apply_lid_bc(u, v, cfg.U_lid); + const double h = config.length / static_cast(cells); + + Field u(cells, cells + 1, 0.0); + Field v(cells + 1, cells, 0.0); + Field pressure(cells, cells, 0.0); + if (initial.available + && initial.u_face.rows == cells && initial.u_face.cols == cells + 1 + && initial.v_face.rows == cells + 1 && initial.v_face.cols == cells + && initial.pressure.rows == cells && initial.pressure.cols == cells) { + u = initial.u_face; + v = initial.v_face; + pressure = initial.pressure; + } Result result; - result.N = N; - result.Re = Re; + result.case_id = case_id; + result.cells = cells; + result.reynolds = reynolds; result.scheme = scheme; result.pressure_solver = pressure_solver; - result.implementation = implementation; - result.localMaxIter = localMaxIter; - - result.Ru.reserve(localMaxIter); - result.Rv.reserve(localMaxIter); - result.Rc_mass.reserve(localMaxIter); - result.Rc_div.reserve(localMaxIter); - result.dt.reserve(localMaxIter); - result.poisson_iters.reserve(localMaxIter); - result.poisson_relative_residual.reserve(localMaxIter); - result.poisson_converged.reserve(localMaxIter); - - std::string status = "maxIter"; - int stagnation_counter = 0; - double prev_mass = std::numeric_limits::infinity(); - int iter = 0; - - const auto t0 = std::chrono::steady_clock::now(); - - for (iter = 1; iter <= localMaxIter; ++iter) { - Matrix u_old = u; - Matrix v_old = v; - - auto [u_star, v_star, dt] = momentum_predictor(u, v, p, Re, scheme, cfg); - dt = std::max(dt, cfg.dt_min); - - Matrix div_star = divergence_field(u_star, v_star, dx, dy); - Matrix rhs(N, 0.0); - for (int i = 0; i < N; ++i) { - for (int j = 0; j < N; ++j) { - rhs(i, j) = div_star(i, j) / dt; - } + result.local_max_iterations = config.max_iterations; + + int consecutive_passes = 0; + int consecutive_pressure_failures = 0; + std::vector stagnation_history; + stagnation_history.reserve(static_cast(config.stagnation_window)); + + const auto start = std::chrono::steady_clock::now(); + for (int iteration = 1; iteration <= config.max_iterations; ++iteration) { + const Field old_u = u; + const Field old_v = v; + + auto [u_star, v_star, dt] = predict_velocity( + u, v, pressure, reynolds, scheme, config + ); + + Field rhs = divergence(u_star, v_star, h, h, cells); + for (double& value : rhs.values) { + value /= dt; } - auto [p_prime, pinfo] = pressure_poisson(rhs, dx, dy, pressure_solver, cfg); + auto [pressure_correction, poisson] = solve_pressure_poisson( + rhs, h, pressure_solver, config + ); + + if (poisson.converged) { + consecutive_pressure_failures = 0; + } else { + ++consecutive_pressure_failures; + ++result.failed_pressure_solves; + } u = u_star; v = v_star; - for (int i = 1; i < N - 1; ++i) { - for (int j = 1; j < N - 1; ++j) { - const double dpdx = (p_prime(i, j + 1) - p_prime(i, j - 1)) / (2.0 * dx); - const double dpdy = (p_prime(i + 1, j) - p_prime(i - 1, j)) / (2.0 * dy); - u(i, j) = u_star(i, j) - dt * dpdx; - v(i, j) = v_star(i, j) - dt * dpdy; + for (int i = 0; i < cells; ++i) { + for (int j = 1; j < cells; ++j) { + u(i, j) -= dt * (pressure_correction(i, j) - pressure_correction(i, j - 1)) / h; } } - - for (size_t k = 0; k < p.a.size(); ++k) { - p.a[k] += cfg.alpha_p * p_prime.a[k]; + for (int i = 1; i < cells; ++i) { + for (int j = 0; j < cells; ++j) { + v(i, j) -= dt * (pressure_correction(i, j) - pressure_correction(i - 1, j)) / h; + } } - const double p_mean = mean_all(p); - for (double& val : p.a) val -= p_mean; - - apply_lid_bc(u, v, cfg.U_lid); - - auto [Ru, Rv, Rc_mass, Rc_div] = velocity_residuals(u, v, u_old, v_old, dx, dy, cfg.U_lid, cfg.L); - result.Ru.push_back(Ru); - result.Rv.push_back(Rv); - result.Rc_mass.push_back(Rc_mass); - result.Rc_div.push_back(Rc_div); - result.dt.push_back(dt); - result.poisson_iters.push_back(pinfo.iter); - result.poisson_relative_residual.push_back(pinfo.final_relative_residual); - result.poisson_converged.push_back(pinfo.converged); - - if (!all_finite(u) || !all_finite(v) || !all_finite(p) || std::max({Ru, Rv, Rc_div}) > cfg.diverged_limit) { - status = "diverged"; + for (std::size_t k = 0; k < pressure.values.size(); ++k) { + pressure.values[k] += config.pressure_relaxation * pressure_correction.values[k]; + } + remove_mean(pressure); + + const Residuals residuals = calculate_residuals( + u, v, old_u, old_v, h, config + ); + + result.iterations = iteration; + result.final_residuals = residuals; + result.velocity_history.push_back(residuals.velocity_update_linf); + result.divergence_linf_history.push_back(residuals.divergence_linf); + result.divergence_l2_history.push_back(residuals.divergence_l2); + result.mass_history.push_back(residuals.global_mass_imbalance); + result.dt_history.push_back(dt); + result.poisson_iteration_history.push_back(poisson.iterations); + result.poisson_relative_history.push_back(poisson.relative_residual); + result.poisson_converged_history.push_back(poisson.converged); + + if (!all_finite(u) || !all_finite(v) || !all_finite(pressure) + || !std::isfinite(residuals.velocity_update_linf) + || !std::isfinite(residuals.divergence_linf) + || !std::isfinite(residuals.divergence_l2)) { + result.status = "non_finite"; break; } - - if (Rc_mass > 0.995 * prev_mass) { - ++stagnation_counter; - } else { - stagnation_counter = 0; + if (std::max({ + residuals.velocity_update_linf, + residuals.divergence_linf, + residuals.divergence_l2 + }) > config.diverged_limit) { + result.status = "diverged"; + break; } - prev_mass = Rc_mass; - - if (Rc_mass < cfg.tol_mass && std::max(Ru, Rv) < cfg.tol_velocity) { - status = "converged"; + if (consecutive_pressure_failures >= config.maximum_pressure_failures) { + result.status = "pressure_not_converged"; break; } - } - if (iter > localMaxIter) { - iter = localMaxIter; - } - - const auto t1 = std::chrono::steady_clock::now(); - result.runtime = std::chrono::duration(t1 - t0).count(); + const bool passed = poisson.converged + && residuals.velocity_update_linf <= config.velocity_tolerance + && residuals.divergence_linf <= config.divergence_linf_tolerance + && residuals.divergence_l2 <= config.divergence_l2_tolerance + && residuals.global_mass_imbalance <= 1e-12; - result.iterations = iter; - result.status = status; - result.stagnation_counter = stagnation_counter; + if (iteration >= config.minimum_iterations && passed) { + ++consecutive_passes; + } else { + consecutive_passes = 0; + } + result.consecutive_pass_count = consecutive_passes; + if (consecutive_passes >= config.consecutive_passes) { + result.status = "converged"; + break; + } - result.x.resize(N); - result.y.resize(N); - for (int k = 0; k < N; ++k) { - result.x[k] = static_cast(k) * cfg.L / static_cast(N - 1); - result.y[k] = result.x[k]; - } + stagnation_history.push_back(residuals.velocity_update_linf); + if (static_cast(stagnation_history.size()) > config.stagnation_window) { + stagnation_history.erase(stagnation_history.begin()); + } + if (iteration > config.minimum_iterations + && static_cast(stagnation_history.size()) == config.stagnation_window) { + const double first = stagnation_history.front(); + const double last = stagnation_history.back(); + const double relative_reduction = (first - last) / std::max(first, 1e-30); + if (first > config.velocity_tolerance + && relative_reduction < config.stagnation_minimum_reduction) { + result.status = "stagnated"; + break; + } + } - result.u = std::move(u); - result.v = std::move(v); - result.p = std::move(p); - result.speed = Matrix(N, 0.0); - for (int i = 0; i < N; ++i) { - for (int j = 0; j < N; ++j) { - result.speed(i, j) = std::sqrt(result.u(i, j) * result.u(i, j) + result.v(i, j) * result.v(i, j)); + if (iteration % 1000 == 0) { + std::cout << " iter=" << iteration + << " vel=" << std::scientific << residuals.velocity_update_linf + << " div=" << residuals.divergence_linf + << " p=" << poisson.relative_residual << '\n'; } } - result.vorticity = compute_vorticity(result.u, result.v, dx, dy); - result.final_Ru = result.Ru.empty() ? 0.0 : result.Ru.back(); - result.final_Rv = result.Rv.empty() ? 0.0 : result.Rv.back(); - result.final_Rc_mass = result.Rc_mass.empty() ? 0.0 : result.Rc_mass.back(); - result.final_Rc_div = result.Rc_div.empty() ? 0.0 : result.Rc_div.back(); - - if (!result.poisson_iters.empty()) { - result.avg_poisson_iters = std::accumulate(result.poisson_iters.begin(), result.poisson_iters.end(), 0.0) - / static_cast(result.poisson_iters.size()); - result.avg_poisson_relative_residual = std::accumulate(result.poisson_relative_residual.begin(), result.poisson_relative_residual.end(), 0.0) - / static_cast(result.poisson_relative_residual.size()); + const auto end = std::chrono::steady_clock::now(); + result.runtime_seconds = std::chrono::duration(end - start).count(); + result.u_face = std::move(u); + result.v_face = std::move(v); + result.pressure = std::move(pressure); + build_cell_center_fields(result, config); + result.ghia = compare_with_ghia(result); + result.quality = quality_label(result); + + if (!result.poisson_iteration_history.empty()) { + result.average_poisson_iterations = std::accumulate( + result.poisson_iteration_history.begin(), + result.poisson_iteration_history.end(), + 0.0 + ) / static_cast(result.poisson_iteration_history.size()); + result.average_poisson_relative_residual = std::accumulate( + result.poisson_relative_history.begin(), + result.poisson_relative_history.end(), + 0.0 + ) / static_cast(result.poisson_relative_history.size()); int saturated = 0; - for (int pi : result.poisson_iters) if (pi >= cfg.poisson_maxIter) ++saturated; - result.pressure_saturation_ratio = static_cast(saturated) / static_cast(result.poisson_iters.size()); + for (const int pressure_iterations : result.poisson_iteration_history) { + if (pressure_iterations >= config.poisson_max_iterations) { + ++saturated; + } + } + result.pressure_saturation_ratio = static_cast(saturated) + / static_cast(result.poisson_iteration_history.size()); } - return result; } -static GhiaData ghia_data(int Re) { - GhiaData d; - d.Re = Re; - if (Re == 100) { - d.y_u = {1.0000,0.9766,0.9688,0.9609,0.9531,0.8516,0.7344,0.6172,0.5000,0.4531,0.2813,0.1719,0.1016,0.0703,0.0625,0.0547,0.0000}; - d.u = {1.0000,0.84123,0.78871,0.73722,0.68717,0.23151,0.00332,-0.13641,-0.20581,-0.21090,-0.15662,-0.10150,-0.06434,-0.04775,-0.04192,-0.03717,0.0000}; - d.x_v = {1.0000,0.9688,0.9609,0.9531,0.9453,0.9063,0.8594,0.8047,0.5000,0.2344,0.2266,0.1563,0.0938,0.0781,0.0703,0.0625,0.0000}; - d.v = {0.0000,-0.05906,-0.07391,-0.08864,-0.10313,-0.16914,-0.22445,-0.24533,0.05454,0.17527,0.17507,0.16077,0.12317,0.10890,0.10091,0.09233,0.0000}; - } else if (Re == 400) { - d.y_u = {1.0000,0.9766,0.9688,0.9609,0.9531,0.8516,0.7344,0.6172,0.5000,0.4531,0.2813,0.1719,0.1016,0.0703,0.0625,0.0547,0.0000}; - d.u = {1.0000,0.75837,0.68439,0.61756,0.55892,0.29093,0.16256,0.02135,-0.11477,-0.17119,-0.32726,-0.24299,-0.14612,-0.10338,-0.09266,-0.08186,0.0000}; - d.x_v = {1.0000,0.9688,0.9609,0.9531,0.9453,0.9063,0.8594,0.8047,0.5000,0.2344,0.2266,0.1563,0.0938,0.0781,0.0703,0.0625,0.0000}; - d.v = {0.0000,-0.12146,-0.15663,-0.19254,-0.22847,-0.23827,-0.44993,-0.38598,0.05186,0.30174,0.30203,0.28124,0.22965,0.20920,0.19713,0.18360,0.0000}; - } else if (Re == 1000) { - d.y_u = {1.0000,0.9766,0.9688,0.9609,0.9531,0.8516,0.7344,0.6172,0.5000,0.4531,0.2813,0.1719,0.1016,0.0703,0.0625,0.0547,0.0000}; - d.u = {1.0000,0.65928,0.57492,0.51117,0.46604,0.33304,0.18719,0.05702,-0.06080,-0.10648,-0.27805,-0.38289,-0.29730,-0.22220,-0.20196,-0.18109,0.0000}; - d.x_v = {1.0000,0.9688,0.9609,0.9531,0.9453,0.9063,0.8594,0.8047,0.5000,0.2344,0.2266,0.1563,0.0938,0.0781,0.0703,0.0625,0.0000}; - d.v = {0.0000,-0.21388,-0.27669,-0.33714,-0.39188,-0.51550,-0.42665,-0.31966,0.02526,0.32235,0.33075,0.37095,0.32627,0.30353,0.29012,0.27485,0.0000}; - } - return d; +static std::string case_name(const Result& result) { + std::ostringstream stream; + stream << "case_" << std::setw(3) << std::setfill('0') << result.case_id + << "_N" << result.cells + << "_Re" << result.reynolds + << '_' << result.scheme + << '_' << result.pressure_solver + << "_serial_cpp"; + return stream.str(); } -static double interp_linear(const std::vector& x, const std::vector& y, double q) { - if (x.empty()) return std::numeric_limits::quiet_NaN(); - if (q <= x.front()) return y.front(); - if (q >= x.back()) return y.back(); - const auto it = std::upper_bound(x.begin(), x.end(), q); - const size_t k = static_cast(std::distance(x.begin(), it)); - const double x0 = x[k - 1], x1 = x[k]; - const double y0 = y[k - 1], y1 = y[k]; - const double t = (q - x0) / (x1 - x0); - return y0 + t * (y1 - y0); -} - -static Metrics validate_against_ghia(const Result& result, const Config& cfg) { - Metrics m; - GhiaData d = ghia_data(result.Re); - if (d.y_u.empty()) return m; - m.available = true; - - const int N = result.N; - const int mid = static_cast(std::round((N + 1) / 2.0)) - 1; // MATLAB round((N+1)/2) - std::vector u_center(N), v_center(N); - for (int i = 0; i < N; ++i) u_center[i] = result.u(i, mid); - for (int j = 0; j < N; ++j) v_center[j] = result.v(mid, j); - - double sum_u2 = 0.0, sum_v2 = 0.0; - double linf_u = 0.0, linf_v = 0.0; - for (size_t k = 0; k < d.y_u.size(); ++k) { - const double un = interp_linear(result.y, u_center, d.y_u[k]); - const double e = un - d.u[k]; - sum_u2 += e * e; - linf_u = std::max(linf_u, std::abs(e)); - } - for (size_t k = 0; k < d.x_v.size(); ++k) { - const double vn = interp_linear(result.x, v_center, d.x_v[k]); - const double e = vn - d.v[k]; - sum_v2 += e * e; - linf_v = std::max(linf_v, std::abs(e)); - } - - m.u_L2 = std::sqrt(sum_u2 / static_cast(d.y_u.size())); - m.v_L2 = std::sqrt(sum_v2 / static_cast(d.x_v.size())); - m.u_Linf = linf_u; - m.v_Linf = linf_v; - - if (result.Re == 100) { - m.u_limit = cfg.validation_u_L2_limit_Re100; - m.v_limit = cfg.validation_v_L2_limit_Re100; - } else if (result.Re == 400) { - m.u_limit = cfg.validation_u_L2_limit_Re400; - m.v_limit = cfg.validation_v_L2_limit_Re400; - } else if (result.Re == 1000) { - m.u_limit = cfg.validation_u_L2_limit_Re1000; - m.v_limit = cfg.validation_v_L2_limit_Re1000; - } else { - m.u_limit = std::numeric_limits::infinity(); - m.v_limit = std::numeric_limits::infinity(); +static void write_history(const Result& result, const Config& config) { + const fs::path path = fs::path(config.data_directory) / (case_name(result) + "_history.csv"); + std::ofstream output(path); + output << std::setprecision(12); + output << "iter,velocity_update_linf,divergence_linf,divergence_l2,global_mass_imbalance,dt,poisson_iters,poisson_relative_residual,poisson_converged\n"; + for (std::size_t k = 0; k < result.velocity_history.size(); ++k) { + output << (k + 1) << ',' + << result.velocity_history[k] << ',' + << result.divergence_linf_history[k] << ',' + << result.divergence_l2_history[k] << ',' + << result.mass_history[k] << ',' + << result.dt_history[k] << ',' + << result.poisson_iteration_history[k] << ',' + << result.poisson_relative_history[k] << ',' + << (result.poisson_converged_history[k] ? 1 : 0) << '\n'; } - m.pass = (m.u_L2 <= m.u_limit && m.v_L2 <= m.v_limit); - return m; -} - -static std::string quality_label(const Result& result, const Metrics& metrics) { - if (result.status == "converged" && metrics.available && metrics.pass) return "converged_validated"; - if (result.status == "converged" && metrics.available && !metrics.pass) return "converged_not_validated"; - if (result.status == "converged" && !metrics.available) return "converged_no_benchmark"; - if (result.status != "converged" && metrics.available && metrics.pass) return "validated_but_not_converged"; - return "needs_improvement"; } -static void write_field_csv(const Result& r, const std::string& case_name, const Config& cfg) { - fs::create_directories(cfg.data_dir); - const fs::path path = fs::path(cfg.data_dir) / (case_name + "_fields.csv"); - std::ofstream out(path); - out << std::setprecision(12); - out << "i,j,x,y,u,v,p,speed,vorticity\n"; - for (int i = 0; i < r.N; ++i) { - for (int j = 0; j < r.N; ++j) { - out << i << ',' << j << ',' << r.x[j] << ',' << r.y[i] << ',' - << r.u(i, j) << ',' << r.v(i, j) << ',' << r.p(i, j) << ',' - << r.speed(i, j) << ',' << r.vorticity(i, j) << '\n'; +static void write_fields(const Result& result, const Config& config) { + const fs::path path = fs::path(config.data_directory) / (case_name(result) + "_fields.csv"); + std::ofstream output(path); + output << std::setprecision(12); + output << "i,j,x,y,u,v,p,speed,vorticity\n"; + for (int i = 0; i < result.cells; ++i) { + for (int j = 0; j < result.cells; ++j) { + output << i << ',' << j << ',' + << result.x[static_cast(j)] << ',' + << result.y[static_cast(i)] << ',' + << result.u_center(i, j) << ',' + << result.v_center(i, j) << ',' + << result.pressure(i, j) << ',' + << result.speed(i, j) << ',' + << result.vorticity(i, j) << '\n'; } } } -static void write_history_csv(const Result& r, const std::string& case_name, const Config& cfg) { - fs::create_directories(cfg.data_dir); - const fs::path path = fs::path(cfg.data_dir) / (case_name + "_history.csv"); - std::ofstream out(path); - out << std::setprecision(12); - out << "iter,Ru,Rv,Rc_mass,Rc_div,dt,poisson_iters,poisson_relative_residual,poisson_converged\n"; - for (size_t k = 0; k < r.Ru.size(); ++k) { - out << (k + 1) << ',' << r.Ru[k] << ',' << r.Rv[k] << ',' << r.Rc_mass[k] << ',' << r.Rc_div[k] - << ',' << r.dt[k] << ',' << r.poisson_iters[k] << ',' << r.poisson_relative_residual[k] - << ',' << (r.poisson_converged[k] ? 1 : 0) << '\n'; - } +static void write_summary_header(std::ofstream& output) { + output << "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,Ghia_u_L2_Limit,Ghia_v_L2_Limit," + << "FinalVelocityLinf,FinalDivergenceL2,GlobalMassImbalance,FailedPressureSolves,ConsecutivePasses\n"; } -static void write_summary_header(std::ofstream& out) { - out << "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,Ghia_u_L2_Limit,Ghia_v_L2_Limit\n"; +static void write_summary_row(std::ofstream& output, const Result& result) { + output << std::setprecision(12) + << result.case_id << ",serial_cpp," + << result.cells << ',' << result.reynolds << ',' + << result.scheme << ',' << result.pressure_solver << ',' + << result.status << ',' << result.quality << ',' + << result.iterations << ',' << result.local_max_iterations << ',' + << result.final_residuals.velocity_update_linf << ',' + << result.final_residuals.velocity_update_linf << ',' + << result.final_residuals.global_mass_imbalance << ',' + << result.final_residuals.divergence_linf << ',' + << result.runtime_seconds << ',' + << result.average_poisson_iterations << ',' + << result.average_poisson_relative_residual << ',' + << result.pressure_saturation_ratio << ',' + << (result.ghia.available ? 1 : 0) << ',' + << (result.ghia.passed ? 1 : 0) << ',' + << result.ghia.u_l2 << ',' << result.ghia.v_l2 << ',' + << result.ghia.u_linf << ',' << result.ghia.v_linf << ',' + << result.ghia.u_limit << ',' << result.ghia.v_limit << ',' + << result.final_residuals.velocity_update_linf << ',' + << result.final_residuals.divergence_l2 << ',' + << result.final_residuals.global_mass_imbalance << ',' + << result.failed_pressure_solves << ',' + << result.consecutive_pass_count << '\n'; } -static void write_summary_row(std::ofstream& out, int case_id, const Result& r, const Metrics& m, const std::string& quality) { - out << std::setprecision(12); - out << case_id << ',' << r.implementation << ',' << r.N << ',' << r.Re << ',' << r.scheme << ',' << r.pressure_solver << ',' - << r.status << ',' << quality << ',' << r.iterations << ',' << r.localMaxIter << ',' - << r.final_Ru << ',' << r.final_Rv << ',' << r.final_Rc_mass << ',' << r.final_Rc_div << ',' << r.runtime << ',' - << r.avg_poisson_iters << ',' << r.avg_poisson_relative_residual << ',' << r.pressure_saturation_ratio << ',' - << (m.available ? 1 : 0) << ',' << (m.pass ? 1 : 0) << ',' - << m.u_L2 << ',' << m.v_L2 << ',' << m.u_Linf << ',' << m.v_Linf << ',' << m.u_limit << ',' << m.v_limit << '\n'; -} - -static void configure_mode(Config& cfg, const std::string& mode) { - const std::string m = lower(mode); - if (m == "quick") { - 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; - } else if (m == "medium") { - cfg.meshes = {32, 64}; - cfg.re_list = {100, 400, 1000}; - cfg.maxIter = 3500; - cfg.maxIter_N128_bonus = 0; - cfg.poisson_maxIter = 1800; - } else if (m == "full") { - // defaults, equivalent to main.m - } else if (m == "single") { - cfg.meshes = {64}; - cfg.re_list = {100}; - cfg.schemes = {"central"}; - cfg.pressure_solvers = {"RBGS"}; - cfg.implementations = {"serial_cpp"}; - } else if (m == "smoke") { - cfg.meshes = {16}; - cfg.re_list = {100}; - cfg.schemes = {"upwind"}; - cfg.pressure_solvers = {"RBGS"}; - cfg.implementations = {"serial_cpp"}; - cfg.maxIter = 20; - cfg.maxIter_N128_bonus = 0; - cfg.maxIter_Re1000_bonus = 0; - cfg.maxIter_central_bonus = 0; - cfg.poisson_maxIter = 50; - } else { - throw std::runtime_error("Unknown mode: " + mode + " (use quick, medium, full, single, or smoke)"); +static void configure_mode( + const std::string& mode, + Config& config, + std::vector& meshes, + std::vector& reynolds_numbers, + std::vector& schemes, + std::vector& pressure_solvers +) { + const std::string normalized = lower(mode); + if (normalized == "smoke") { + meshes = {16}; + reynolds_numbers = {100}; + schemes = {"upwind"}; + pressure_solvers = {"RBGS"}; + config.max_iterations = 20; + config.minimum_iterations = 100; + config.poisson_max_iterations = 100; + config.save_fields = false; + } else if (normalized == "quick") { + meshes = {24, 32}; + reynolds_numbers = {100}; + schemes = {"upwind", "central"}; + pressure_solvers = {"RBSOR"}; + config.velocity_tolerance = 1e-7; + } else if (normalized == "medium") { + meshes = {32}; + reynolds_numbers = {100, 400, 1000}; + schemes = {"upwind", "central"}; + pressure_solvers = {"RBSOR"}; + config.velocity_tolerance = 1e-7; + } else if (normalized == "grid") { + meshes = {16, 32, 64}; + reynolds_numbers = {100}; + schemes = {"central"}; + pressure_solvers = {"RBSOR"}; + config.velocity_tolerance = 1e-7; + } else if (normalized == "full") { + meshes = {32, 64, 128}; + reynolds_numbers = {100, 400, 1000}; + schemes = {"upwind", "central"}; + pressure_solvers = {"RBGS", "RBSOR"}; + config.velocity_tolerance = 1e-7; + } else if (normalized != "single") { + throw std::runtime_error("Unknown mode: " + mode + " (use smoke, single, quick, medium, grid, or full)"); } } -static void print_usage(const char* exe) { - std::cout << "Usage:\n" - << " " << exe << " --mode quick|medium|full|single|smoke\n" - << " " << exe << " --single --N 64 --Re 100 --scheme central --pressure RBGS --implementation serial_cpp\n\n" - << "Options:\n" - << " --no-fields Do not write full field CSV files\n" - << " --maxIter VALUE Override base outer iterations\n" - << " --poisson-maxIter VALUE Override pressure solver iterations\n" - << "\nImplementation note: this C++ package has one serial_cpp solver. MATLAB labels loop/vectorized are accepted as aliases only.\n"; +static void print_usage(const char* executable) { + std::cout + << "Usage:\n" + << " " << executable << " --mode smoke|single|quick|medium|grid|full\n" + << " " << executable << " --single --N 32 --Re 100 --scheme upwind --pressure RBSOR\n\n" + << "Options:\n" + << " --no-fields Skip full field CSV output\n" + << " --strict Return a non-zero exit code when a case does not converge\n" + << " --maxIter VALUE Maximum outer iterations\n" + << " --poisson-maxIter VALUE Maximum pressure iterations\n" + << " --tol-velocity VALUE Dimensionless velocity-update tolerance\n" + << " --tol-divergence VALUE Dimensionless Linf divergence tolerance\n" + << " --tol-divergence-l2 VALUE Dimensionless L2 divergence tolerance\n" + << " --poisson-tol VALUE Relative pressure residual tolerance\n" + << " --alpha-u VALUE Momentum relaxation factor\n" + << " --alpha-p VALUE Pressure update relaxation factor\n" + << " --cfl VALUE Convective CFL limit\n" + << " --dt-max VALUE Maximum pseudo-time step\n" + << " --min-iterations VALUE Minimum outer iterations\n" + << " --consecutive-passes VALUE Required consecutive converged iterations\n"; } int main(int argc, char** argv) { - Config cfg; + Config config; std::string mode = "quick"; bool explicit_single = false; - int single_N = 64; - int single_Re = 100; - std::string single_scheme = "central"; - std::string single_pressure = "RBGS"; - std::string single_implementation = "serial_cpp"; + int single_cells = 32; + int single_reynolds = 100; + std::string single_scheme = "upwind"; + std::string single_pressure = "RBSOR"; try { for (int i = 1; i < argc; ++i) { - std::string arg = argv[i]; - auto require_value = [&](const std::string& name) -> std::string { - if (i + 1 >= argc) throw std::runtime_error("Missing value for " + name); + const std::string argument = argv[i]; + auto require_value = [&]() -> std::string { + if (i + 1 >= argc) { + throw std::runtime_error("Missing value for " + argument); + } return argv[++i]; }; - if (arg == "--help" || arg == "-h") { + if (argument == "--help" || argument == "-h") { print_usage(argv[0]); return 0; - } else if (arg == "--mode") { - mode = require_value(arg); - } else if (arg == "--single") { - explicit_single = true; + } else if (argument == "--mode") { + mode = require_value(); + } else if (argument == "--single") { mode = "single"; - } else if (arg == "--N") { - single_N = std::stoi(require_value(arg)); explicit_single = true; + } else if (argument == "--N") { + single_cells = std::stoi(require_value()); mode = "single"; - } else if (arg == "--Re") { - single_Re = std::stoi(require_value(arg)); explicit_single = true; + } else if (argument == "--Re") { + single_reynolds = std::stoi(require_value()); mode = "single"; - } else if (arg == "--scheme") { - single_scheme = require_value(arg); explicit_single = true; + } else if (argument == "--scheme") { + single_scheme = require_value(); mode = "single"; - } else if (arg == "--pressure") { - single_pressure = require_value(arg); explicit_single = true; + } else if (argument == "--pressure") { + single_pressure = require_value(); mode = "single"; - } else if (arg == "--implementation") { - single_implementation = require_value(arg); explicit_single = true; - mode = "single"; - } else if (arg == "--no-fields") { - cfg.save_fields = false; - } else if (arg == "--maxIter") { - cfg.maxIter = std::stoi(require_value(arg)); - } else if (arg == "--poisson-maxIter") { - cfg.poisson_maxIter = std::stoi(require_value(arg)); + } else if (argument == "--implementation") { + (void)require_value(); + } else if (argument == "--no-fields") { + config.save_fields = false; + } else if (argument == "--strict") { + config.strict_exit = true; + } else if (argument == "--maxIter") { + config.max_iterations = std::stoi(require_value()); + } else if (argument == "--poisson-maxIter") { + config.poisson_max_iterations = std::stoi(require_value()); + } else if (argument == "--tol-velocity") { + config.velocity_tolerance = std::stod(require_value()); + } else if (argument == "--tol-divergence") { + config.divergence_linf_tolerance = std::stod(require_value()); + } else if (argument == "--tol-divergence-l2") { + config.divergence_l2_tolerance = std::stod(require_value()); + } else if (argument == "--poisson-tol") { + config.poisson_relative_tolerance = std::stod(require_value()); + } else if (argument == "--alpha-u") { + config.momentum_relaxation = std::stod(require_value()); + } else if (argument == "--alpha-p") { + config.pressure_relaxation = std::stod(require_value()); + } else if (argument == "--cfl") { + config.cfl = std::stod(require_value()); + } else if (argument == "--dt-max") { + config.dt_max = std::stod(require_value()); + } else if (argument == "--min-iterations") { + config.minimum_iterations = std::stoi(require_value()); + } else if (argument == "--consecutive-passes") { + config.consecutive_passes = std::stoi(require_value()); } else { - throw std::runtime_error("Unknown argument: " + arg); + throw std::runtime_error("Unknown argument: " + argument); } } - configure_mode(cfg, mode); - if (explicit_single) { - cfg.meshes = {single_N}; - cfg.re_list = {single_Re}; - cfg.schemes = {lower(single_scheme)}; - cfg.pressure_solvers = {upper(single_pressure)}; - cfg.implementations = {normalize_implementation(single_implementation)}; + std::vector meshes; + std::vector reynolds_numbers; + std::vector schemes; + std::vector pressure_solvers; + configure_mode(mode, config, meshes, reynolds_numbers, schemes, pressure_solvers); + if (explicit_single || lower(mode) == "single") { + meshes = {single_cells}; + reynolds_numbers = {single_reynolds}; + schemes = {lower(single_scheme)}; + pressure_solvers = {upper(single_pressure)}; } - fs::create_directories(cfg.data_dir); - const fs::path summary_path = fs::path(cfg.data_dir) / ("study_summary_" + lower(mode) + ".csv"); + fs::create_directories(config.data_directory); + const fs::path summary_path = fs::path(config.data_directory) / ("study_summary_" + lower(mode) + ".csv"); std::ofstream summary(summary_path); write_summary_header(summary); - const int nCases = static_cast(cfg.meshes.size() * cfg.re_list.size() * cfg.schemes.size() - * cfg.pressure_solvers.size() * cfg.implementations.size()); - std::cout << "\nLID-DRIVEN CAVITY C++ SOLVER\n"; - std::cout << "Mode: " << mode << "\n"; - std::cout << "Total simulations: " << nCases << "\n"; - std::cout << "Summary: " << summary_path.string() << "\n\n"; + const int number_of_cases = static_cast( + meshes.size() * reynolds_numbers.size() * schemes.size() * pressure_solvers.size() + ); + std::cout << "\nLID-DRIVEN CAVITY C++ PHASE 2 SOLVER\n" + << "Mode: " << mode << "\n" + << "Total simulations: " << number_of_cases << "\n" + << "Summary: " << summary_path << "\n\n"; int case_id = 0; - for (int N : cfg.meshes) { - for (int Re : cfg.re_list) { - for (const auto& scheme : cfg.schemes) { - for (const auto& pressure_solver : cfg.pressure_solvers) { - for (const auto& implementation : cfg.implementations) { - ++case_id; - std::ostringstream name; - name << "case_" << std::setw(3) << std::setfill('0') << case_id << std::setfill(' ') - << "_N" << N << "_Re" << Re << "_" << lower(scheme) - << "_" << upper(pressure_solver) << "_" << lower(implementation); - const std::string case_name = name.str(); - - std::cout << "[" << std::setw(3) << std::setfill('0') << case_id << std::setfill(' ') - << "] N=" << N << " Re=" << Re << " Scheme=" << scheme - << " Pressure=" << pressure_solver << " Implementation=" << implementation << "\n"; - - Result r = solve_lid_cavity(N, Re, scheme, pressure_solver, implementation, cfg); - Metrics metrics = validate_against_ghia(r, cfg); - const std::string quality = quality_label(r, metrics); - write_summary_row(summary, case_id, r, metrics, quality); - summary.flush(); - write_history_csv(r, case_name, cfg); - if (cfg.save_fields) write_field_csv(r, case_name, cfg); - - std::cout << " status=" << r.status << " quality=" << quality - << " iter=" << r.iterations << "/" << r.localMaxIter - << " Rc_mass=" << std::scientific << std::setprecision(3) << r.final_Rc_mass - << " Rc_div=" << r.final_Rc_div - << " runtime=" << std::fixed << std::setprecision(2) << r.runtime << "s" - << " avgPiter=" << std::setprecision(1) << r.avg_poisson_iters - << " pSat=" << std::setprecision(2) << r.pressure_saturation_ratio << "\n"; - if (metrics.available) { - std::cout << " Ghia L2: u=" << std::scientific << std::setprecision(3) << metrics.u_L2 - << "(limit " << metrics.u_limit << "), v=" << metrics.v_L2 - << "(limit " << metrics.v_limit << "), pass=" << (metrics.pass ? 1 : 0) << "\n"; - } - std::cout << std::defaultfloat; + int failed_cases = 0; + std::map, InitialState> continuation_states; + + for (const int cells : meshes) { + for (const std::string& scheme : schemes) { + for (const std::string& pressure_solver : pressure_solvers) { + for (const int reynolds : reynolds_numbers) { + ++case_id; + const auto key = std::make_tuple(cells, lower(scheme), upper(pressure_solver)); + InitialState initial; + const auto previous = continuation_states.find(key); + if (previous != continuation_states.end()) { + initial = previous->second; + } + + std::cout << '[' << std::setw(3) << std::setfill('0') << case_id << std::setfill(' ') + << "] N=" << cells + << " Re=" << reynolds + << " Scheme=" << scheme + << " Pressure=" << pressure_solver << '\n'; + + Result result = solve_case( + case_id, + cells, + reynolds, + scheme, + pressure_solver, + config, + initial + ); + write_summary_row(summary, result); + summary.flush(); + write_history(result, config); + if (config.save_fields) { + write_fields(result, config); + } + + if (result.status == "converged") { + continuation_states[key] = InitialState{ + true, + result.u_face, + result.v_face, + result.pressure + }; + } else { + ++failed_cases; } + + std::cout << " status=" << result.status + << " quality=" << result.quality + << " iter=" << result.iterations << '/' << result.local_max_iterations + << " vel=" << std::scientific << result.final_residuals.velocity_update_linf + << " div=" << result.final_residuals.divergence_linf + << " runtime=" << std::fixed << std::setprecision(2) << result.runtime_seconds << "s\n"; + if (result.ghia.available) { + std::cout << " Ghia L2: u=" << std::scientific << result.ghia.u_l2 + << " v=" << result.ghia.v_l2 + << " pass=" << (result.ghia.passed ? 1 : 0) << '\n'; + } + std::cout << std::defaultfloat; } } } } - std::cout << "\nFinished. CSV outputs are in " << cfg.data_dir << "\n"; + std::cout << "\nFinished. CSV outputs are in " << config.data_directory << '\n'; + if (config.strict_exit && failed_cases > 0) { + return 2; + } return 0; - } catch (const std::exception& e) { - std::cerr << "ERROR: " << e.what() << "\n"; + } catch (const std::exception& error) { + std::cerr << "ERROR: " << error.what() << '\n'; print_usage(argv[0]); return 1; }