Enh/k omega sst by greole · Pull Request #554 · exasim-project/NeoN

greole · 2026-06-21T06:57:32Z

Motivation

WHAT is it, WHY it is needed. 😘

github-actions · 2026-06-21T07:04:13Z

Thank you for your PR, here are some useful tips:

If this is your first PR, please add yourself to the AUTHORS.md file
Use labels to change the workflow behavior. Set auto-fix to automatically apply fixes, full-ci for a comprehensive set of tests and sanitizer runs, ready-for-review to show that you're done and need a review.
I deployed a test documentation to https://exasim-project.com/NeoN/Build_PR_554
For full-ci builds benchmark results are available at https://github.com/exasim-project/NeoFOAM-BenchmarkData/554
Also check if the GPU pipelines on https://gitlab-ce.lrz.de/greole/neon/-/jobs are successful.

Copilot

Pull request overview

Note

Copilot couldn't run its full agentic review because no GitHub Actions runner was available. Make sure your repository has a runner available to run Copilot's review, or add a copilot-setup-steps.yml file specifying one with the runs-on attribute. See the docs for more details.

This PR adds support for implicit transform boundary conditions (slip/symmetry) via per-component diagonal corrections, extends boundary-condition infrastructure (including flux-aware inletOutlet), and introduces a bounded divergence operator to improve positivity in convection discretizations.

Changes:

Add diagCmpt storage to LinearSystem and update Gauss-Green Laplacian assembly + Ginkgo solvers (serial + distributed) to apply per-component diagonal corrections for implicit slip/symmetry.
Add new volume BCs (slip, inletOutlet) and extend BoundaryContext to carry surface-scalar fields (e.g., face flux phi).
Introduce boundedDiv operator wrapper and add/extend unit tests for BC behavior and implicit-transform solver paths.

Reviewed changes

Copilot reviewed 22 out of 22 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
test/linearAlgebra/ginkgo.cpp	Adds unit test covering implicit transform diag correction in serial Ginkgo solve.
test/finiteVolume/cellCentred/operator/laplacianOperator.cpp	Adds unit test verifying Laplacian assembly populates `diagCmpt` for implicit slip.
test/finiteVolume/cellCentred/boundary/volume/volSymmetry.cpp	Extends symmetry tests to cover deferred vs implicit modes.
test/finiteVolume/cellCentred/boundary/volume/volSlip.cpp	New tests for slip BC (scalar + vector, deferred + implicit).
test/finiteVolume/cellCentred/boundary/volume/volInletOutlet.cpp	New tests for inletOutlet BC with/without `phi` context (scalar + vector).
test/finiteVolume/cellCentred/boundary/volume/CMakeLists.txt	Registers new unit tests.
test/distributed/operator.cpp	Adds distributed unit test for implicit transform diagonal correction solve.
src/linearAlgebra/ginkgo/ginkgoDistributed.cpp	Implements distributed implicit-transform component solve by in-place diagonal edits.
src/linearAlgebra/ginkgo/ginkgo.cpp	Implements serial implicit-transform component solve by in-place diagonal edits.
src/finiteVolume/cellCentred/operators/gaussGreenLaplacian.cpp	Accumulates implicit transform damping into `diagCmpt` instead of shared diagonal.
src/finiteVolume/cellCentred/operators/boundedDiv.cpp	New bounded divergence operator implementation.
src/finiteVolume/cellCentred/boundary/boundaryContext.cpp	Adds support for inserting/reading surface-scalar fields in `BoundaryContext`.
src/CMakeLists.txt	Adds `boundedDiv.cpp` to build.
include/NeoN/linearAlgebra/linearSystem.hpp	Adds `diagCmpt` storage + deep copy/reset/copyToExecutor support.
include/NeoN/finiteVolume/cellCentred/operators/boundedDiv.hpp	Declares bounded divergence operator wrapper.
include/NeoN/finiteVolume/cellCentred/boundary/volumeBoundaryFactory.hpp	Adds `transformImplicit` attribute flag for BCs.
include/NeoN/finiteVolume/cellCentred/boundary/volume/symmetry.hpp	Refactors symmetry to share slip/symmetry implementation and expose implicit mode flag.
include/NeoN/finiteVolume/cellCentred/boundary/volume/slip.hpp	Introduces slip BC (deferred/implicit normal damping).
include/NeoN/finiteVolume/cellCentred/boundary/volume/inletOutlet.hpp	Introduces flux-dependent inletOutlet BC using `BoundaryContext` surface flux.
include/NeoN/finiteVolume/cellCentred/boundary/volume/detail/slipSymmetry.hpp	Shared slip/symmetry implementation + tolerant parsing of `"implicit"`.
include/NeoN/finiteVolume/cellCentred/boundary/boundaryContext.hpp	Adds surface-scalar field slots to `BoundaryContext`.
include/NeoN/finiteVolume/cellCentred/boundary.hpp	Registers slip and inletOutlet boundary types.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

+    parallelFor(
+        exec,
+        {0, nrows},
+        NEON_LAMBDA(const localIdx cell) {
+            Kokkos::atomic_sub(&values[ma.diagIdx(cell)], diagC[cell][I]);
+        },
+        "applyImplicitTransformDiag"
+    );
+    gkoExec->synchronize();


+    parallelFor(
+        exec,
+        {0, nrows},
+        NEON_LAMBDA(const localIdx cell) {
+            Kokkos::atomic_sub(&values[ma.diagIdx(cell)], diagC[cell][I]);
+        },
+        "applyImplicitTransformDiagDist"
+    );
+    gkoExec->synchronize();


+        NEON_LAMBDA(const localIdx cell) {
+            Kokkos::atomic_sub(&values[ma.diagIdx(cell)], diagC[cell][I]);
+        },


+    [[nodiscard]] const std::shared_ptr<Vector<RHSValueType>>& diagCmpt() const
+    {
+        return diagCmpt_;
+    }


    const VolumeField<scalar>& scalarFieldPtr(const std::string& name) const;
    const VolumeField<Vec3>& vectorFieldPtr(const std::string& name) const;
    const VolumeField<Tensor>& tensorFieldPtr(const std::string& name) const;
+    const SurfaceField<scalar>& surfaceScalarField(const std::string& name) const;


+void BoundedDiv<FieldValueType, AssemblyType>::div(
+    la::LinearSystem<AssemblyType, FieldValueType>& ls,
+    const SurfaceField<scalar>& faceFlux,
+    const VolumeField<FieldValueType>& phi,
+    const dsl::Coeff operatorScaling
+) const
+{
+    inner_->div(ls, faceFlux, phi, operatorScaling);
+    applyBoundedDiagInternal<FieldValueType, AssemblyType>(
+        ls, faceFlux, this->mesh_, operatorScaling
+    );
+    applyBoundedDiagBoundary<FieldValueType, AssemblyType>(
+        ls, faceFlux, this->mesh_, operatorScaling
+    );
+    applyBoundedDiagProcBoundary<FieldValueType, AssemblyType>(
+        ls, faceFlux, this->mesh_, operatorScaling
+    );
+}


+        SolverStats stats;
+        solveImplicitTransformComponent<0>(
+            sys, x, exec_, gkoExec_, gkoMtx, factory_, stats, l1Control, values, ma, diagC, nrows
+        );
+        solveImplicitTransformComponent<1>(
+            sys, x, exec_, gkoExec_, gkoMtx, factory_, stats, l1Control, values, ma, diagC, nrows
+        );
+        solveImplicitTransformComponent<2>(
+            sys, x, exec_, gkoExec_, gkoMtx, factory_, stats, l1Control, values, ma, diagC, nrows
+        );


greole · 2026-06-21T19:21:46Z

+    if (sys.diagCmpt() && sys.diagCmpt()->size() > 0)
+    {
+        auto values = const_cast<Vector<scalar>&>(sys.matrix().values()).view();
+        const auto ma = sys.faceToMatrixAddress()->view(sys.matrix().rowOffs().view());
+        auto diagC = sys.diagCmpt()->view();
+        const localIdx nrows = sys.rhs().size();
+        gkoExec_->synchronize();
+
+        SolverStats stats;
+        solveImplicitTransformComponentDist<0>(
+            sys,
+            x,
+            exec_,
+            gkoExec_,
+            comm,
+            gkoMtx,
+            factory_,
+            stats,
+            l1Control,
+            values,
+            ma,
+            diagC,
+            nrows
+        );
+        solveImplicitTransformComponentDist<1>(
+            sys,
+            x,
+            exec_,
+            gkoExec_,
+            comm,
+            gkoMtx,
+            factory_,
+            stats,
+            l1Control,
+            values,
+            ma,
+            diagC,
+            nrows
+        );
+        solveImplicitTransformComponentDist<2>(
+            sys,
+            x,
+            exec_,
+            gkoExec_,
+            comm,
+            gkoMtx,
+            factory_,
+            stats,
+            l1Control,
+            values,
+            ma,
+            diagC,
+            nrows
+        );
+        return stats;
+    }


avoid component wise solve here.

Cherry-picked from exasim-project/NeoN PR #554 (commit 91c3709). Brings the Ginkgo memory-pool handling: the CUDA executor is now backed by NeoN's Umpire device QuickPool via UmpireCudaAllocator, so Ginkgo's per-solve multigrid build/teardown reuses device memory instead of churning cudaMalloc/cudaFree (which fragmented the device heap and OOM'd the second pressure solve on large cases). Serial/CPU/HIP keep the default executor path. Also carries the upstream parse() fix that type-checks the 'solver' any_cast and remaps solver::Ir's inner preconditioner factory to the 'solver' key Ir::parse() reads. Conflict resolution: the upstream commit also refactored the signatures of solveImplicitTransformComponent / ...Dist (slip-BC implicit-transform ginkgo solvers). Those functions do not exist on this branch — PR #554 introduces them separately and they are unused here — so that hunk was dropped; only the memory-pool handling and the parse() fix are applied. ginkgoDistributed.cpp is therefore unchanged. (cherry picked from commit 91c3709)

…on-linearUpwind runs GeometryScheme eagerly computed and stored two SurfaceField<Vec3> (faceDeltaOwner_, faceDeltaNeighbour_) for every mesh in update(), even when the only consumer (linearUpwind) is not selected. On an ~18M-cell mesh that is ~2*nInternalFaces*Vec3 (~2.6 GB) of always-on device memory, preventing the case from fitting on one GPU. Make them lazy/opt-in, mirroring edf9888 (faceFluxCorrection opt-in): - faceDelta* are now mutable std::optional, allocated+filled only via ensureFaceDeltas(). - LinearUpwind's constructor opts in (while the mesh centres are still alive), so cases using Gauss upwind/linear allocate nothing. - update() no longer computes faceDelta* nor frees the mesh centres; reset() (the centre release) is deferred + made idempotent, triggered on first read of any cached geometry field or right after faceDelta* are built, so a late opt-in can still read the centres. - The fast streaming linearUpwind kernel is unchanged. Runs that select linearUpwind keep identical behaviour and field values.

… linear systems Cache one immutable topology-only bundle per mesh — the CSR system sparsity pattern, its FaceToMatrixAddress, and the boundary COO sparsity — in mesh.stencilDB() (mirroring GeometryScheme::readOrCreate) and share it across every LinearSystem built on that mesh. These arrays depend only on mesh topology and were previously rebuilt and duplicated per equation (U, p, nuTilda, ...); only the per-system value/RHS vectors are now allocated fresh. readOrCreateSparsityBundle is keyed on the sparsity types so the CSR system and COO boundary patterns never collide; the MPI off-diagonal/proc-face sparsity stays per-system (v1). Adds neon_test_sharedSparsity proving pointer sharing and value independence. Pure memory/aliasing refactor — no numerical change.

Unify slip and symmetry onto a shared slipSymmetry helper (they apply the same operator, differing only in registered name and in where they may be applied) and add a normal-damping treatment selectable via the opt-in "implicit" dict key. - Deferred (default): refGrad = -deltaCoeffs*(U.n)*n enters the per-component RHS through the existing fixed-gradient assembly, keeping the shared scalar matrix + multi-RHS solve (zero extra memory). Also feeds grad(U) boundary reconstruction, so no tensor BC is needed there. - Implicit (opt-in): flags BoundaryAttributes::transformImplicit; the Vec3 Laplacian assembly accumulates the per-component diagonal weight g|S|*delta*|n_c| into a lazily-allocated LinearSystem::diagCmpt store. The Ginkgo scalar-matrix/Vec3-rhs solve (serial and distributed) solves the three components segregated, subtracting each column's correction from the shared diagonal in place and restoring it -- no matrix copy. Reuses solve_impl/solve_impl_dist, so the l1ScaledResidual criterion works in both modes. Tests: deferred refGrad + implicit attr (volSlip/volSymmetry), diagCmpt assembly (laplacianOperator), serial segregated solve std+l1 (ginkgo), distributed solveDist std+l1 under MPI_SIZE 3 (distributed/operator).

…Wallmodel

greole · 2026-07-22T09:06:35Z

only comment changes, revert.

…ylov solvers Re-adds the SolverLease / findUpdatable / cacheOrUpdateSolver workspace-reuse layer (cachedWorkspace_ + gko::solver::invalidate_and_extract_workspace) dropped from enh/kOmegaSST during a `pull --rebase` onto c5ba297 (reflog HEAD@{10}); the champion build (pre-rebase f0ff3ddbd9) had it, current HEAD had lost it. PBiCGStab/Cg + Jacobi/ILU configs (coupled momentum fallback, k, omega) are not gko::UpdateMatrixValue, so they regenerate every solve; this reclaims the Krylov scratch Workspace after each solve and feeds it into the next generate(matrix, ws). Multigrid (pressure) keeps Strategy-1b in-place update; updatable-but-uncached configs deliberately do NOT reuse (stale coarse sizes -> DimensionMismatch). NOTE: does NOT recover the occDrivAer 0.79->2.6 s/step regression (coupled momentum + this reuse still 2.60 s/step vs champion 0.77) -- that is a lower-level deps change (Kokkos/Ginkgo rebuilt 07-19). Restoration is still correct and worth keeping. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

… (stream-ordered) Restores the champion's fence-free Vector lifecycle (lost in the same rebase). The frees run on the same Kokkos stream as the kernels touching data_, so the device-wide fence(exec_) was redundant -- same rationale as the dsl/solver.hpp fence audit. Perf-neutral on occDrivAer at this operating point (the 0.79->2.6 regression is elsewhere in the rebase-lost assembly diff), but correct to restore. Keeps the nullptr-free guards. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

…esh stencil DB Restores linearUpwindLimitedGrad() (lost in the enh/kOmegaSST rebase): the linearUpwindV deferred correction's CellLimitedGrad operator is a function of the MESH only (GaussGreenGrad factory + GeometryScheme + CellToFaceStencil), yet HEAD reconstructed all of it on EVERY momentum assemble. On the 16M-cell/rank occDrivAer that per-step rebuild was ~90% of the momentum-assembly host cost (profiler: luw.gradOpCtor) -- the bulk of the 0.79->2.6 s/step regression. Caching it in the per-mesh stencil DB (built once, reused every timestep) recovers 2.61 -> 1.07 s/step (2.4x). Correctness unchanged (continuity 2.79e-6, U residuals identical). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

…r solve path createGkoMtxDist re-wrapped the local CSR block as a fresh gko::matrix::Csr on every solve. Csr::create_const with the default (automatical->load_balance) GPU strategy recomputes the srow load-balancing array by scanning the ~O(nnz) sparsity each time -- ~72 ms on the 16M-row local block, per solve, for p/k/omega (~230 ms/step, entirely hidden from the reported Ginkgo "Solve time"). Restore the champion's distributed-matrix cache (cachedDistMtx_ + cachedLocalValPtr_, dropped in a rebase -- present at f0ff3ddbd9, absent on develop): in steady state the topology is fixed and NeoN re-assembles the local block in place, so the local value-buffer pointer is unchanged and the whole dist_mtx wrapper is reused (only the off-diagonal Coo values are refreshed). The pointer guard rebuilds if the buffer was reallocated. Wired to the scalar solveDist path only; the Vec3 momentum paths pass non-persistent locals (the implicit transform-BC branch mutates the diagonal in place, so its wrapper must not be cached). occDrivAer (65M cells, 4xH200, restart 1050): steady 1.105 -> 0.907 s/step (-18%), turbulence phase 287 -> 141 ms (k/omega each -72 ms), pressure pSolve 508 -> 460 ms. Continuity bit-identical (2.735186693e-6), p iters unchanged (8) -- values-only reuse. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

…3 momentum path The scalar cache (3cbdfadbfb) left the Vec3 momentum solveDist rebuilding its distributed matrix every solve, re-running Csr::create_const's load-balancing srow scan (~72 ms). Wire it to cachedDistMtx_/cachedLocalValPtr_ too: the wrapper is a non-owning VIEW over the rank-local value buffer, so it always reflects live values. Safe for all three sub-paths -- the implicit-transform branch shifts the diagonal in place per component but RESTORES it right after each component solve (atomic_sub / atomic_add pair), the buffer is re-assembled in place each step (pointer stable; the guard rebuilds if it changes), and the fused-slip branch applies its shift through FusedDiagShiftMatrix at the operator level (no buffer mutation). Also fix the createGkoMtxDist signature in test/distributed/operator.cpp (was calling the pre-cache overload). occDrivAer (65M, 4xH200, restart 1050): momentum solve 258 -> 182 ms, steady 0.907 -> 0.830 s/step. Continuity bit-identical (2.79e-6), U/p iters unchanged. Cumulative with the two restored caches: 2.6 -> 0.83 s/step (champion 0.78). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Resolves 13 conflicted files. Most were comment-only, where develop's rewordings were taken to follow its "self-contained comments" policy (54600d6) -- this drops the [fence-audit] tags in vector.cpp and the case-specific debugging notes in inletOutlet.hpp / slipSymmetry.hpp. Substantive resolutions, all taking develop's corrected version: * boundedDiv.cpp explicit path: develop divides by cell volume (sumPhi/V * psi * scaling); the branch had dropped the /V believing the volumes cancel. They cancel only for the implicit Sp diagonal, not for the per-unit-volume explicit result. * ginkgo.cpp / ginkgoDistributed.cpp diagonal transform: develop's plain writes replace the branch's atomics. The parallelFor range is {0, nrows} and each iteration writes its own diagIdx(cell), so there is no race. * gaussGreenDivLaplacian.cpp: bounded_ is now always re-derived, so a second read() with a different scheme cannot leave a stale true. * geometryScheme: develop's eager faceDelta* computation in update() wins over the branch's lazy/opt-in version (c20b8c0). The lazy path had a construction-order dependency -- a consumer built after the first read of another cached geometry field found the mesh centres already freed. This gives up that commit's ~2.6 GB saving on large meshes. * src/CMakeLists.txt: develop's NeoN_SRCS list variable, needed for set_source_files_properties(... LANGUAGE CUDA). The branch's explicit source list was a strict subset, so nothing is lost. The branch's ginkgo work is kept in full: solver caching (Strategy 1b), Workspace reuse (Strategy 3), MergedPgm coarseners, mixed-precision inner solves, sellp local matrix format and the fused slip solve. Develop's l1InConfig_ guard (suppressing a double-applied L1 criterion) was merged into those by hand at all three call sites, since develop applied it to code the branch had rewritten around its solver cache. Ginkgo stays pinned at the branch's ba5e6607, which the MergedPgm code requires; develop's 6a3abf8c does not provide gko::UpdateMatrixValue.

greole force-pushed the enh/kOmegaSST branch from 5a09ccd to fa4ab58 Compare June 21, 2026 06:59

greole force-pushed the enh/kOmegaSST branch 2 times, most recently from cb95d4c to 4778a2c Compare June 21, 2026 13:23

greole requested a review from Copilot June 21, 2026 13:24

Copilot AI reviewed Jun 21, 2026

View reviewed changes

Copilot started reviewing on behalf of greole June 21, 2026 13:32 View session

greole force-pushed the enh/kOmegaSST branch from 4778a2c to e18c908 Compare June 21, 2026 19:17

greole commented Jun 21, 2026

View reviewed changes

greole force-pushed the enh/kOmegaSST branch 4 times, most recently from c9d4a73 to 330e30d Compare June 23, 2026 15:38

greole force-pushed the enh/kOmegaSST branch from 5119c04 to 4a7647a Compare June 24, 2026 04:42

greole force-pushed the enh/kOmegaSST branch 2 times, most recently from cfe27e2 to c5ba297 Compare July 13, 2026 06:02

HendriceH and others added 13 commits July 22, 2026 10:54

inline kernels

a7e8443

fixup

4566d18

add inletOutlet

8ea2145

fixup

e07f450

fixup

0dda573

add steady-state

238db2d

fix Gko version

9243870

fixup

f2a6480

fixup

3073e88

add target format

b038515

Gregor Olenik and others added 17 commits July 22, 2026 10:54

update

b8a5144

fully decouple values (own and neig) for pinned values, like in omega…

f9c540c

…Wallmodel

make l1Stop register as config solver

93cbbb3

fix blowup issue on degenerate cells with slip bc

f75d8ce

fix blowup issue on degenerate cells with slip bc

dac962e

bump version

d7e7597

fix gko version

4eeeed6

zero init

f70f129

fix nvcc do_copy trampoline issue in productoin

8f0c11e

bump ginkgo

d5c2581

clean-up

77d86fe

wip implement optional preconditioner/caching reuse

4e2f7c8

add inner precission support

c3c88f2

update

b25e5af

update

577acdf

update

2ecdb92

update

89b3a92

greole force-pushed the enh/kOmegaSST branch from 0b8db25 to 89b3a92 Compare July 22, 2026 08:59

greole commented Jul 22, 2026

View reviewed changes

Comment thread include/NeoN/dsl/expression.hpp

greole Jul 22, 2026

Copy link
Copy Markdown

Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

only comment changes, revert.

greole and others added 7 commits July 22, 2026 12:19

fixup

6b04d90

update patches

dfc6b42

greole force-pushed the enh/kOmegaSST branch from dfc6b42 to 05265f0 Compare July 24, 2026 18:10

update patches

6b7355c

greole force-pushed the enh/kOmegaSST branch from dc9a22c to 6b7355c Compare July 31, 2026 07:41

Navigation Menu

Uh oh!

Enh/k omega sst - #554

Enh/k omega sst#554
greole wants to merge 39 commits into
developfrom
enh/kOmegaSST

greole commented Jun 21, 2026

Uh oh!

github-actions Bot commented Jun 21, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

greole Jun 21, 2026

Uh oh!

greole Jul 22, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Uh oh!

Conversation

greole commented Jun 21, 2026

Motivation

Uh oh!

github-actions Bot commented Jun 21, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

greole Jun 21, 2026

Choose a reason for hiding this comment

Uh oh!

greole Jul 22, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants