Feat/multi gpu nvshmem

CMakeLists.txt

@@ -61,6 +61,7 @@ list(REMOVE_DUPLICATES BRAE_HEADER_DIRS)
 add_library(brae_core STATIC
     src/parallel/pstream/cf_pstream.cu
     src/parallel/pstream/device_halo.cu
+    src/parallel/pstream/device_reduce.cu
     src/OpenFOAM/meshes/interface/processor_interface.cu
     src/OpenFOAM/db/foam_token_reader.cu
     src/OpenFOAM/meshes/primitive_mesh.cu
@@ -263,10 +264,155 @@ if(BRAE_WITH_NVSHMEM)
                          ${CMAKE_CURRENT_SOURCE_DIR}/validation/kEpsCorrect)
     endforeach()
 
+    # Phase 4b: the distributed kOmegaSST::correct() (2nd RAS model) reproduces np=1
+    add_executable(test_gpu_parallel_komegasst tests/test_gpu_parallel_komegasst.cu)
+    target_include_directories(test_gpu_parallel_komegasst PRIVATE ${BRAE_HEADER_DIRS}
+                               ${CMAKE_CURRENT_SOURCE_DIR}/tests
+                               ${CMAKE_CURRENT_SOURCE_DIR}/src/applications/solvers/simpleFoam)
+    target_link_libraries(test_gpu_parallel_komegasst PRIVATE brae_core)
+    add_test(NAME gpu_parallel_komegasst_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_komegasst>
+                     32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/sst_ref.txt)
+    set_tests_properties(gpu_parallel_komegasst_ref PROPERTIES FIXTURES_SETUP SSTREF)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_komegasst_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_komegasst>
+                         32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/sst_ref.txt)
+        set_tests_properties(gpu_parallel_komegasst_np${_np} PROPERTIES FIXTURES_REQUIRED SSTREF)
+    endforeach()
+
+    # Phase 4b: the CLOSED distributed turbulent SIMPLE loop (momentum+pressure with nuEff=nu+nut, then
+    # kEpsilon correct(), maintained across iterations) reproduces np=1.
+    add_executable(test_gpu_parallel_turbloop tests/test_gpu_parallel_turbloop.cu)
+    target_include_directories(test_gpu_parallel_turbloop PRIVATE ${BRAE_HEADER_DIRS}
+                               ${CMAKE_CURRENT_SOURCE_DIR}/tests
+                               ${CMAKE_CURRENT_SOURCE_DIR}/src/applications/solvers/simpleFoam)
+    target_link_libraries(test_gpu_parallel_turbloop PRIVATE brae_core)
+    add_test(NAME gpu_parallel_turbloop_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_turbloop>
+                     32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/turbloop_ref.txt)
+    set_tests_properties(gpu_parallel_turbloop_ref PROPERTIES FIXTURES_SETUP TURBLOOP)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_turbloop_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_turbloop>
+                         32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/turbloop_ref.txt)
+        set_tests_properties(gpu_parallel_turbloop_np${_np} PROPERTIES FIXTURES_REQUIRED TURBLOOP)
+    endforeach()
+
+    # SST in the closed loop (second field = omega): same gate, sst mode
+    add_test(NAME gpu_parallel_turbloop_sst_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_turbloop>
+                     32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/turbloop_sst_ref.txt "" sst)
+    set_tests_properties(gpu_parallel_turbloop_sst_ref PROPERTIES FIXTURES_SETUP TURBLOOPSST)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_turbloop_sst_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_turbloop>
+                         32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/turbloop_sst_ref.txt "" sst)
+        set_tests_properties(gpu_parallel_turbloop_sst_np${_np} PROPERTIES FIXTURES_REQUIRED TURBLOOPSST)
+    endforeach()
+
+    # Phase 4b: the full distributed kEpsilon::correct() (production + wall functions + eps/k transport + nut)
+    # reproduces np=1. Duct walls are y/z, cut is x -> near-wall cells sit adjacent to the cut.
+    add_executable(test_gpu_parallel_kepsilon tests/test_gpu_parallel_kepsilon.cu)
+    target_include_directories(test_gpu_parallel_kepsilon PRIVATE ${BRAE_HEADER_DIRS}
+                               ${CMAKE_CURRENT_SOURCE_DIR}/tests
+                               ${CMAKE_CURRENT_SOURCE_DIR}/src/applications/solvers/simpleFoam)
+    target_link_libraries(test_gpu_parallel_kepsilon PRIVATE brae_core)
+    add_test(NAME gpu_parallel_kepsilon_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_kepsilon>
+                     32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/keps_ref.txt)
+    set_tests_properties(gpu_parallel_kepsilon_ref PROPERTIES FIXTURES_SETUP KEPSREF)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_kepsilon_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_kepsilon>
+                         32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/keps_ref.txt)
+        set_tests_properties(gpu_parallel_kepsilon_np${_np} PROPERTIES FIXTURES_REQUIRED KEPSREF)
+    endforeach()
+
+    # Phase 4b foundation: the distributed scalar transport solve (deviceParallelScalarTransport) -- the
+    # reusable core of the turbulence step (k/eps/omega/nuTilda) -- reproduces np=1 on a flux-carrying cut.
+    add_executable(test_gpu_parallel_scalar tests/test_gpu_parallel_scalar.cu)
+    target_include_directories(test_gpu_parallel_scalar PRIVATE ${BRAE_HEADER_DIRS}
+                               ${CMAKE_CURRENT_SOURCE_DIR}/tests
+                               ${CMAKE_CURRENT_SOURCE_DIR}/src/applications/solvers/simpleFoam)
+    target_link_libraries(test_gpu_parallel_scalar PRIVATE brae_core)
+    add_test(NAME gpu_parallel_scalar_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_scalar>
+                     32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/scalar_ref.txt)
+    set_tests_properties(gpu_parallel_scalar_ref PROPERTIES FIXTURES_SETUP SCALARREF)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_scalar_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_scalar>
+                         32 16 16 ${CMAKE_CURRENT_BINARY_DIR}/scalar_ref.txt)
+        set_tests_properties(gpu_parallel_scalar_np${_np} PROPERTIES FIXTURES_REQUIRED SCALARREF)
+    endforeach()
+
+    # Phase 4a: the linearUpwind/LUST interface correction on a cut that CARRIES FLUX (the cavity's cut is a
+    # symmetry plane -> phi~0 -> those terms are multiplied by zero and any bug in them is invisible).
+    add_executable(test_gpu_parallel_duct tests/test_gpu_parallel_duct.cu)
+    target_include_directories(test_gpu_parallel_duct PRIVATE ${BRAE_HEADER_DIRS}
+                               ${CMAKE_CURRENT_SOURCE_DIR}/tests
+                               ${CMAKE_CURRENT_SOURCE_DIR}/src/applications/solvers/simpleFoam)
+    target_link_libraries(test_gpu_parallel_duct PRIVATE brae_core)
+    # np=1 writes the reference (no cuts), then np=2/4 must reproduce it. Sequential: np>1 needs the file.
+    add_test(NAME gpu_parallel_duct_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_duct>
+                     32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_ref.txt)
+    set_tests_properties(gpu_parallel_duct_ref PROPERTIES FIXTURES_SETUP DUCTREF)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_duct_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_duct>
+                         32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_ref.txt)
+        set_tests_properties(gpu_parallel_duct_np${_np} PROPERTIES FIXTURES_REQUIRED DUCTREF)
+    endforeach()
+
+    # SHEARED duct: the same flux-carrying cut, but on a NON-ORTHOGONAL mesh (shear 0.5 -> 26.57 deg, =
+    # atan(0.5)). Every other brae test mesh is orthogonal, where the cut-face geometry conventions
+    # (dNei = Cf - C_remote) coincide with their normal-projected variants and cannot be told apart. This is
+    # the only gate that exercises the interface terms on a non-orthogonal cut.
+    add_test(NAME gpu_parallel_duct_skew_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_duct>
+                     32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_skew_ref.txt nolust "" 1 0.5)
+    set_tests_properties(gpu_parallel_duct_skew_ref PROPERTIES FIXTURES_SETUP DUCTSKEW)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_duct_skew_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_duct>
+                         32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_skew_ref.txt nolust "" 1 0.5)
+        set_tests_properties(gpu_parallel_duct_skew_np${_np} PROPERTIES FIXTURES_REQUIRED DUCTSKEW)
+    endforeach()
+
+    # NON-ORTHOGONAL correction on the sheared cut: nonOrthDeltaCoeffs implicit + the explicit corrVec.grad
+    # correction (momentum) + the pEqn source AND face-flux correction. Invisible on an orthogonal mesh
+    # (corrVec == 0), so this gate is the only thing standing between it and silent wrongness on real meshes.
+    add_test(NAME gpu_parallel_duct_nonorth_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_duct>
+                     32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_nonorth_ref.txt nolust "" 1 0.5 nonorth)
+    set_tests_properties(gpu_parallel_duct_nonorth_ref PROPERTIES FIXTURES_SETUP DUCTNO)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_duct_nonorth_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_duct>
+                         32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_nonorth_ref.txt nolust "" 1 0.5 nonorth)
+        set_tests_properties(gpu_parallel_duct_nonorth_np${_np} PROPERTIES FIXTURES_REQUIRED DUCTNO)
+    endforeach()
+
+    # LUST (0.75 linear + 0.25 linearUpwind) over the same flux-carrying cut: it exercises BOTH interface
+    # deferred corrections (linearUpwind's gradient term AND the linear term), so it is the stronger gate.
+    add_test(NAME gpu_parallel_duct_lust_ref
+             COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np 1 $<TARGET_FILE:test_gpu_parallel_duct>
+                     32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_lust_ref.txt lust)
+    set_tests_properties(gpu_parallel_duct_lust_ref PROPERTIES FIXTURES_SETUP DUCTLUST)
+    foreach(_np 2 4)
+        add_test(NAME gpu_parallel_duct_lust_np${_np}
+                 COMMAND ${MPIEXEC_EXECUTABLE} --oversubscribe -np ${_np} $<TARGET_FILE:test_gpu_parallel_duct>
+                         32 16 16 20 ${CMAKE_CURRENT_BINARY_DIR}/duct_lust_ref.txt lust)
+        set_tests_properties(gpu_parallel_duct_lust_np${_np} PROPERTIES FIXTURES_REQUIRED DUCTLUST)
+    endforeach()
+
     # Phase 4a EXIT GATE: the OOM proof -- a mesh too big for one GPU runs split across two. No ctest entry:
     # the mesh size is hardware-specific (it must exceed ONE GPU's VRAM), so it is driven by hand.
     add_executable(test_gpu_parallel_oom tests/test_gpu_parallel_oom.cu)
     target_include_directories(test_gpu_parallel_oom PRIVATE ${BRAE_HEADER_DIRS}
+                               ${CMAKE_CURRENT_SOURCE_DIR}/tests
                                ${CMAKE_CURRENT_SOURCE_DIR}/src/applications/solvers/simpleFoam)
     target_link_libraries(test_gpu_parallel_oom PRIVATE brae_core)
 

bench/setup_spuma.sh

@@ -25,7 +25,7 @@
 #  Env:  SPUMA_SRC (clone target, default $HOME/spuma)   NVARCH (default: autodetect, e.g. 121)
 #        NVHPC_DIR (HPC SDK compilers/bin, default: autodetect)   JOBS (default 16)
 # ============================================================================
-set -u
+set +u   # OpenFOAM etc/bashrc is not nounset-safe; sourcing it under set -u fails
 SPUMA_SRC="${SPUMA_SRC:-$HOME/spuma}"
 NVARCH="${NVARCH:-$(nvidia-smi --query-gpu=compute_cap --format=csv,noheader 2>/dev/null | head -1 | tr -d .)}"
 JOBS="${JOBS:-16}"

src/OpenFOAM/decompose/local_assembly.cuh

@@ -60,6 +60,123 @@ inline std::vector<std::vector<scalar>> computeProcDeltaCoeffs(
     return pdc;
 }
 
+// Per processor face, the two cell-centre->face-centre offset vectors that a linearUpwind deferred correction
+// needs: dOwn = Cf - C[localCell] and dNei = Cf - C[remoteCell]. On an internal face these are dm.dOwn/dNei;
+// across a cut the remote centre lives on the other rank, so exchange it once at setup (geometry is static)
+// exactly as computeProcDeltaCoeffs does. Returned interleaved xyz per face, one vector per processor patch:
+//   dOwn[j][3*i+{0,1,2}], dNei[j][3*i+{0,1,2}]
+inline void computeProcUpwindD(
+    const LocalMesh& lm,
+    const FvGeometry& lg,
+    const std::vector<FvPatch>& lpatches,
+    std::vector<std::vector<scalar>>& dOwn,
+    std::vector<std::vector<scalar>>& dNei)
+{
+    std::vector<const FvPatch*> pp;
+    for (const FvPatch& p : lpatches)
+        if (p.type == "processor") pp.push_back(&p);
+
+    const std::size_t n = pp.size();
+    std::vector<std::vector<scalar>> sendC(n), recvC(n);   // 3 scalars (cell centre) per face
+    for (std::size_t j = 0; j < n; ++j)
+    {
+        const label sz = pp[j]->size;
+        sendC[j].resize(3 * sz);
+        recvC[j].resize(3 * sz);
+        for (label i = 0; i < sz; ++i)
+        {
+            const vector& C = lg.C()[pp[j]->faceCells[i]];
+            sendC[j][3 * i]     = C.x;
+            sendC[j][3 * i + 1] = C.y;
+            sendC[j][3 * i + 2] = C.z;
+        }
+        Pstream::irecv(recvC[j].data(), static_cast<int>(3 * sz), lm.procNbr[j], 0);
+        Pstream::isend(sendC[j].data(), static_cast<int>(3 * sz), lm.procNbr[j], 0);
+    }
+    Pstream::waitAll();
+
+    dOwn.assign(n, {});
+    dNei.assign(n, {});
+    for (std::size_t j = 0; j < n; ++j)
+    {
+        const label sz = pp[j]->size;
+        dOwn[j].resize(3 * sz);
+        dNei[j].resize(3 * sz);
+        for (label i = 0; i < sz; ++i)
+        {
+            const vector& Cf = lg.Cf()[pp[j]->start + i];
+            const vector& Cl = lg.C()[pp[j]->faceCells[i]];
+            const vector  Cr{recvC[j][3 * i], recvC[j][3 * i + 1], recvC[j][3 * i + 2]};
+            const vector dO = Cf - Cl, dN = Cf - Cr;
+            dOwn[j][3 * i] = dO.x; dOwn[j][3 * i + 1] = dO.y; dOwn[j][3 * i + 2] = dO.z;
+            dNei[j][3 * i] = dN.x; dNei[j][3 * i + 1] = dN.y; dNei[j][3 * i + 2] = dN.z;
+        }
+    }
+}
+
+// Per processor face, the two quantities OF's "corrected" laplacian needs, mirroring FvGeometry's internal-face
+// construction EXACTLY (fv_geometry.cu makeInterpolation) with the local cell as owner and the remote as
+// neighbour (delta = C_remote - C_local, unitArea = Sf/magSf outward from local):
+//     nonOrthDeltaCoeffs = 1 / max(unitArea . delta, 0.05*|delta|)     (the IMPLICIT coeff; replaces 1/|delta|)
+//     corrVec            = unitArea - delta*nonOrthDeltaCoeffs         (the EXPLICIT deferred correction vector)
+// The remote centre lives on the other rank, so exchange it once at setup -- geometry is static -- exactly as
+// computeProcDeltaCoeffs does. On an ORTHOGONAL mesh corrVec == 0 and nonOrthDeltaCoeffs == deltaCoeffs, so this
+// collapses to the existing path (which is why an orthogonal mesh cannot tell the two apart).
+// Returned per processor patch: nod[j][i] (scalar), cv[j][3*i+{0,1,2}] (interleaved xyz).
+inline void computeProcNonOrth(
+    const LocalMesh& lm,
+    const FvGeometry& lg,
+    const std::vector<FvPatch>& lpatches,
+    std::vector<std::vector<scalar>>& nod,
+    std::vector<std::vector<scalar>>& cv)
+{
+    std::vector<const FvPatch*> pp;
+    for (const FvPatch& p : lpatches)
+        if (p.type == "processor") pp.push_back(&p);
+
+    const std::size_t n = pp.size();
+    std::vector<std::vector<scalar>> sendC(n), recvC(n);
+    for (std::size_t j = 0; j < n; ++j)
+    {
+        const label sz = pp[j]->size;
+        sendC[j].resize(3 * sz);
+        recvC[j].resize(3 * sz);
+        for (label i = 0; i < sz; ++i)
+        {
+            const vector& C = lg.C()[pp[j]->faceCells[i]];
+            sendC[j][3 * i]     = C.x;
+            sendC[j][3 * i + 1] = C.y;
+            sendC[j][3 * i + 2] = C.z;
+        }
+        Pstream::irecv(recvC[j].data(), static_cast<int>(3 * sz), lm.procNbr[j], 0);
+        Pstream::isend(sendC[j].data(), static_cast<int>(3 * sz), lm.procNbr[j], 0);
+    }
+    Pstream::waitAll();
+
+    nod.assign(n, {});
+    cv.assign(n, {});
+    for (std::size_t j = 0; j < n; ++j)
+    {
+        const label sz = pp[j]->size;
+        nod[j].resize(sz);
+        cv[j].resize(3 * sz);
+        for (label i = 0; i < sz; ++i)
+        {
+            const label  gf = pp[j]->start + i;
+            const vector delta{ recvC[j][3*i]   - lg.C()[pp[j]->faceCells[i]].x,
+                                recvC[j][3*i+1] - lg.C()[pp[j]->faceCells[i]].y,
+                                recvC[j][3*i+2] - lg.C()[pp[j]->faceCells[i]].z };
+            const vector unitArea = lg.Sf()[gf] / lg.magSf()[gf];
+            const scalar d        = 1.0 / std::fmax(dot(unitArea, delta), 0.05 * mag(delta));
+            const vector c        = unitArea - delta * d;
+            nod[j][i]      = d;
+            cv[j][3*i]     = c.x;
+            cv[j][3*i + 1] = c.y;
+            cv[j][3*i + 2] = c.z;
+        }
+    }
+}
+
 // Build the scalar DistributedMatrix structure (diag/upper/lower/interfaceCoeffs) for the momentum
 // div(phi,U) - laplacian(nuEff,U) from the serially-assembled LOCAL FvMatrix (correct internal +
 // real-boundary coeffs; ZERO at processor patches since ProcessorFvPatchField returns zero coeffs)

src/OpenFOAM/decompose/scotch_decomposition.cu

@@ -3,6 +3,9 @@
 #include <cstdio>   // scotch.h uses FILE
 #include <scotch.h>
 #include <stdexcept>
+#include <fstream>
+#include <filesystem>
+#include <string>
 
 namespace brae {
 
@@ -67,4 +70,48 @@ std::vector<label> scotchDecompose(const PrimitiveMesh& m, label nParts)
     return cellToPart;
 }
 
+// Cached SCOTCH decomposition: the cell->partition map is STATIC per (mesh, nParts), but scotchDecompose costs
+// ~1s on a 1M-cell mesh and ran on every launch (a top-3 startup cost). Persist it to
+// cacheDir/.brae_decomp_np<N> and reload when the file is at least as new as cacheDir/owner (mesh unchanged) and
+// the header (magic, nParts, nCells) matches; otherwise recompute and write. Falls back to a plain decompose on
+// any I/O problem, so it can never be wrong -- at worst it recomputes.
+std::vector<label> scotchDecomposeCached(const PrimitiveMesh& m, label nParts, const std::string& cacheDir)
+{
+    namespace fs = std::filesystem;
+    constexpr unsigned kMagic = 0x44454331u;                 // "DEC1"
+    const label nC = m.nCells();
+    const std::string cachePath = cacheDir + "/.brae_decomp_np" + std::to_string(nParts);
+    const std::string ownerPath = cacheDir + "/owner";
+
+    std::error_code ec;
+    const bool fresh = fs::exists(cachePath, ec) && fs::exists(ownerPath, ec)
+                    && fs::last_write_time(cachePath, ec) >= fs::last_write_time(ownerPath, ec);
+    if (fresh)
+    {
+        std::ifstream is(cachePath, std::ios::binary);
+        unsigned magic = 0; label np = 0, n = 0;
+        is.read(reinterpret_cast<char*>(&magic), sizeof magic);
+        is.read(reinterpret_cast<char*>(&np),    sizeof np);
+        is.read(reinterpret_cast<char*>(&n),     sizeof n);
+        if (is && magic == kMagic && np == nParts && n == nC)
+        {
+            std::vector<label> v(nC);
+            is.read(reinterpret_cast<char*>(v.data()), static_cast<std::streamsize>(nC) * sizeof(label));
+            if (is) return v;                                // cache hit
+        }
+    }
+
+    std::vector<label> v = scotchDecompose(m, nParts);       // miss -> compute + persist
+    std::ofstream os(cachePath, std::ios::binary);
+    if (os)
+    {
+        const unsigned magic = kMagic; const label np = nParts, n = nC;
+        os.write(reinterpret_cast<const char*>(&magic), sizeof magic);
+        os.write(reinterpret_cast<const char*>(&np),    sizeof np);
+        os.write(reinterpret_cast<const char*>(&n),     sizeof n);
+        os.write(reinterpret_cast<const char*>(v.data()), static_cast<std::streamsize>(nC) * sizeof(label));
+    }
+    return v;
+}
+
 } // namespace brae