From 80dee99e751e59d28b2e3863fbe3d7589ea28fd1 Mon Sep 17 00:00:00 2001 From: tito Date: Fri, 17 Jul 2026 08:05:47 +0200 Subject: [PATCH 01/13] distributed linearUpwind --- src/OpenFOAM/decompose/local_assembly.cuh | 54 +++++++ .../simpleFoam/parallel_device_foam.cuh | 37 +++-- .../simpleFoam/parallel_device_simple.cuh | 146 +++++++++++++++++- 3 files changed, 222 insertions(+), 15 deletions(-) diff --git a/src/OpenFOAM/decompose/local_assembly.cuh b/src/OpenFOAM/decompose/local_assembly.cuh index 0981c66..d0dcca2 100644 --- a/src/OpenFOAM/decompose/local_assembly.cuh +++ b/src/OpenFOAM/decompose/local_assembly.cuh @@ -60,6 +60,60 @@ inline std::vector> 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& lpatches, + std::vector>& dOwn, + std::vector>& dNei) +{ + std::vector pp; + for (const FvPatch& p : lpatches) + if (p.type == "processor") pp.push_back(&p); + + const std::size_t n = pp.size(); + std::vector> 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(3 * sz), lm.procNbr[j], 0); + Pstream::isend(sendC[j].data(), static_cast(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; + } + } +} + // 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) diff --git a/src/applications/solvers/simpleFoam/parallel_device_foam.cuh b/src/applications/solvers/simpleFoam/parallel_device_foam.cuh index a867974..55b9a95 100644 --- a/src/applications/solvers/simpleFoam/parallel_device_foam.cuh +++ b/src/applications/solvers/simpleFoam/parallel_device_foam.cuh @@ -111,11 +111,12 @@ inline int runParallelDeviceFoam(int argc, char** argv) "`mpirun -np N brae_simpleFoam -case ` for RAS in parallel, or `brae -case ` for " "RAS on a single GPU."); - // fvSchemes div(phi,U). ParallelDeviceSimple implements Gauss upwind (+ the `bounded` -Sp(div(phi),U) - // term). Anything else -- linearUpwind, limitedLinear, LUST, ... -- is a DIFFERENT discretisation, and - // silently substituting upwind produces a converged, plausible, WRONG answer (on the cavity that is a - // ~6% field difference vs the single-GPU solver). Refuse instead, exactly as RAS is refused above. - bool boundedDiv = false; + // fvSchemes div(phi,U). The distributed path implements Gauss upwind, bounded, and linearUpwind (the + // matrix is upwind either way; linearUpwind adds the deferred gradient correction, cut faces included). + // Anything else -- limitedLinear, LUST, linearUpwindV, plain Gauss linear -- is a DIFFERENT + // discretisation, and silently substituting upwind produces a converged, plausible, WRONG answer (on + // the cavity that was a ~6% field difference). Refuse those, exactly as RAS is refused above. + bool boundedDiv = false, linUpwind = false; { std::string divLine; std::istringstream fsch(readFileExpanded(caseDir + "/system/fvSchemes")); @@ -125,15 +126,22 @@ inline int runParallelDeviceFoam(int argc, char** argv) if (!divLine.empty()) { boundedDiv = divLine.find("bounded") != std::string::npos; - const bool upwind = divLine.find("upwind") != std::string::npos; - const bool linearUpwind = divLine.find("linearUpwind") != std::string::npos; - if (!upwind || linearUpwind) + linUpwind = divLine.find("linearUpwind") != std::string::npos; + // NB "linearUpwind" has a capital U, so it does NOT contain the substring "upwind" -- the two + // must be tested separately or a linearUpwind case reads as "no upwind scheme at all". + const bool upwindFamily = (divLine.find("upwind") != std::string::npos) || linUpwind; + // linearUpwindV adds OF's vector direction limiter on top of linearUpwind -- NOT implemented, + // and it contains the substring "linearUpwind", so it must be excluded explicitly. + const bool unsupported = divLine.find("linearUpwindV") != std::string::npos + || divLine.find("limitedLinear") != std::string::npos + || divLine.find("LUST") != std::string::npos; + if (!upwindFamily || unsupported) throw std::runtime_error( "brae -parallel: div(phi,U) scheme '" + divLine + "' is not implemented on the " - "multi-GPU device path (it supports 'Gauss upwind' and 'bounded Gauss upwind'). " - "Substituting upwind would converge to a DIFFERENT answer than the scheme asks for, " - "so this is refused rather than solved wrongly. Use `brae -case ` (single GPU) " - "for the full scheme set."); + "multi-GPU device path (it supports 'Gauss upwind', 'bounded Gauss upwind' and " + "'[bounded] Gauss linearUpwind grad(U)'). Substituting a scheme would converge to a " + "DIFFERENT answer than fvSchemes asks for, so this is refused rather than solved " + "wrongly. Use `brae -case ` (single GPU) for the full scheme set."); } } @@ -179,14 +187,15 @@ inline int runParallelDeviceFoam(int argc, char** argv) std::printf("brae (device, distributed) | case=%s np=%d | laminar | nu=%.3g | %ld cells\n", caseDir.c_str(), nproc, nu, (long)nC); std::printf(" relax U=%.2g p=%.2g | tol p=%.1g U=%.1g | endTime=%d | residualControl=%s | div(phi,U)=%sGauss upwind\n", - relaxU, relaxP, tolP, tolU, endTime, hasRC ? "on" : "off", boundedDiv ? "bounded " : ""); + relaxU, relaxP, tolP, tolU, endTime, hasRC ? "on" : "off", + (std::string(boundedDiv ? "bounded " : "") + (linUpwind ? "Gauss linearUpwind" : "Gauss upwind")).c_str()); } std::vector Ug; std::vector pg; int nIter = 0; { // scope: the solver's symmetric-heap buffers must be freed BEFORE Pstream::finalize - ParallelDeviceSimple solver(P, U0, p0, nu, relaxU, relaxP, tolU, tolP, 2000, boundedDiv); + ParallelDeviceSimple solver(P, U0, p0, nu, relaxU, relaxP, tolU, tolP, 2000, boundedDiv, linUpwind); auto ok = [](scalar res, scalar ctl) { return ctl < 0 || res < ctl; }; bool converged = false; diff --git a/src/applications/solvers/simpleFoam/parallel_device_simple.cuh b/src/applications/solvers/simpleFoam/parallel_device_simple.cuh index 3402924..16bdeef 100644 --- a/src/applications/solvers/simpleFoam/parallel_device_simple.cuh +++ b/src/applications/solvers/simpleFoam/parallel_device_simple.cuh @@ -18,6 +18,7 @@ #include "device_halo.cuh" #include "parallel_simple.cuh" // Partition, distributeFromCells #include "reconstruct.cuh" +#include "local_assembly.cuh" // computeProcUpwindD #include #include @@ -101,6 +102,51 @@ void laplacianInterfaceKernel( atomicAdd(&diag[faceCells[f]], -coeff[f]); } +// linearUpwind deferred correction at a PROCESSOR face -- the distributed analogue of +// linearUpwindCorrKernel (device_simple.cu) and of deviceCyclicAddLinUpwindCorr (no rotation here: a +// processor cut is a plain interior face, the two sides share an orientation). +// +// The matrix stays plain UPWIND; linearUpwind is this explicit source. Per cut face, the local cell is the +// face's owner and phi is the OUTWARD flux, so it is the owner branch of the internal kernel: +// corr[localCell] += phi * (grad(U_comp)[upwind] . d_upwind) +// with the upwind side chosen by sign(phi): +// phi >= 0 -> LOCAL is upwind: grad at the local cell, d = dOwn = Cf - C[local] +// phi < 0 -> REMOTE is upwind: grad at the remote cell (gxN/gyN/gzN, from the halo), d = dNei = Cf - C[remote] +// dOwn/dNei are static geometry, precomputed once by host computeProcUpwindD (the remote centre needs an +// exchange). The caller does relaxSrc -= corr, exactly as for the internal correction. +__global__ +void linUpwindInterfaceKernel( + const label* __restrict__ faceCells, + const scalar* __restrict__ phi, + const scalar* __restrict__ gx, + const scalar* __restrict__ gy, + const scalar* __restrict__ gz, + const scalar* __restrict__ gxN, + const scalar* __restrict__ gyN, + const scalar* __restrict__ gzN, + const scalar* __restrict__ dOwn, + const scalar* __restrict__ dNei, + scalar* __restrict__ corr, + int n) +{ + const int f = blockIdx.x * blockDim.x + threadIdx.x; + if (f >= n) return; + const scalar pf = phi[f]; + const label c = faceCells[f]; + scalar g[3], d[3]; + if (pf >= 0.0) // local upwind: local gradient, local offset + { + g[0] = gx[c]; g[1] = gy[c]; g[2] = gz[c]; + d[0] = dOwn[3*f]; d[1] = dOwn[3*f+1]; d[2] = dOwn[3*f+2]; + } + else // remote upwind: the halo gradient at the remote cell, remote offset + { + g[0] = gxN[f]; g[1] = gyN[f]; g[2] = gzN[f]; + d[0] = dNei[3*f]; d[1] = dNei[3*f+1]; d[2] = dNei[3*f+2]; + } + atomicAdd(&corr[c], pf * (g[0]*d[0] + g[1]*d[1] + g[2]*d[2])); // a cell may own several cut faces +} + // Per-cell sum of |processor off-diagonal|, for fvMatrix::relax's diagonal-dominance term. Host // parallelRelaxMatrix adds |interfaceCoeffs| into sumOff; on device this feeds deviceRelaxDiag's cycSumOff // hook (the same role the cyclic interface's off-diagonal sum plays). @@ -171,6 +217,65 @@ inline void deviceInterfaceOffDiagSum( } } +// Add the processor faces' linearUpwind deferred correction into `corr` (which already holds the internal +// deviceLinearUpwindCorr result). `gx/gy/gz` are grad(U_comp) as CELL fields -- already processor-consistent +// (the halo injected the coupled face value before the gaussGrad). The three exchanges below fetch that same +// gradient at the REMOTE cell, needed only where the remote side is upwind. +// `phiF[i]` is the processor-face flux, `dOwnD/dNeiD[i]` the precomputed offsets (3 per face). +inline void deviceLinUpwindInterface( + DeviceHalo& halo, + const std::vector>& faceCells, + const std::vector>& phiF, + const DeviceBuffer& gx, + const DeviceBuffer& gy, + const DeviceBuffer& gz, + const std::vector>& dOwnD, + const std::vector>& dNeiD, + DeviceBuffer& corr, + cudaStream_t stream = cudaStreamPerThread) +{ + constexpr int TPB = 128; + // The recv buffer is shared and reused, so each component's remote gradient must be COPIED out before the + // next exchange overwrites it (see the hazard note in device_halo.cuh). + const int nI = halo.nInterfaces(); + std::vector> gN[3]; + const DeviceBuffer* gsrc[3] = { &gx, &gy, &gz }; + for (int k = 0; k < 3; ++k) + { + halo.exchange(gsrc[k]->data(), stream); + gN[k].resize(nI); + for (int i = 0; i < nI; ++i) + { + const int n = static_cast(halo.size(i)); + if (n <= 0) continue; + gN[k][i].resize(n); + cudaCheck( + cudaMemcpyAsync(gN[k][i].data(), halo.recvData(i), n * sizeof(scalar), + cudaMemcpyDeviceToDevice, stream), + "linUpwind halo grad copy"); + } + halo.waitExchange(stream); + } + for (int i = 0; i < nI; ++i) + { + const int n = static_cast(halo.size(i)); + if (n <= 0) continue; + detail::linUpwindInterfaceKernel<<<(n + TPB - 1) / TPB, TPB, 0, stream>>>( + faceCells[i].data(), + phiF[i].data(), + gx.data(), + gy.data(), + gz.data(), + gN[0][i].data(), + gN[1][i].data(), + gN[2][i].data(), + dOwnD[i].data(), + dNeiD[i].data(), + corr.data(), + n); + } +} + // Assemble a laplacian(gamma, .) matrix's processor coupling (the pressure equation): fold -coeff into `diag` // and produce `ifCoeff[i]` per interface. `coeffGeo[i] = gammaF*magSf*procDelta` of interface i. inline void deviceLaplacianInterface( @@ -251,7 +356,8 @@ public: scalar tolU, scalar tolP, int maxIter, - bool bounded = false) + bool bounded = false, + bool linearUpwind = false) : P_(part), nu_(nu), relaxU_(relaxU), @@ -260,6 +366,7 @@ public: tolP_(tolP), maxIter_(maxIter), bounded_(bounded), + linearUpwind_(linearUpwind), lnC_(part.Lm.mesh.nCells()), nIf_(part.Lm.mesh.nInternalFaces()), dm_(buildDeviceMesh(part.Lm.mesh, part.lg, part.lp)), @@ -316,6 +423,18 @@ public: for (std::size_t i = 0; i < P_.Lm.procFaceCells.size(); ++i) faceCellsD_[i].copyFrom(P_.Lm.procFaceCells[i]); for (std::size_t i = 0; i < P_.procW.size(); ++i) weightsD_[i].copyFrom(P_.procW[i]); + // linearUpwind needs Cf-C on BOTH sides of a cut; the remote centre needs an exchange, but the + // geometry is static, so do it once here rather than per iteration. + if (linearUpwind_) + { + std::vector> dO, dN; + computeProcUpwindD(P_.Lm, P_.lg, P_.lp, dO, dN); + dOwnD_.resize(dO.size()); + dNeiD_.resize(dN.size()); + for (std::size_t i = 0; i < dO.size(); ++i) dOwnD_[i].copyFrom(dO[i]); + for (std::size_t i = 0; i < dN.size(); ++i) dNeiD_[i].copyFrom(dN[i]); + } + // initial conservative flux phi = flux(U0), internal + boundary (processor faces carry the coupled flux) const SurfaceScalarField phi0 = fvc::flux(U0, P_.Lm.mesh, P_.lg, P_.lp); phiInt_.copyFrom(std::vector(phi0.internal.begin(), phi0.internal.begin() + nIf_)); @@ -359,6 +478,7 @@ private: scalar nu_, relaxU_, relaxP_, tolU_, tolP_; int maxIter_; bool bounded_ = false; + bool linearUpwind_ = false; label lnC_, nIf_, nBnd_ = 0; bool validC_[3] = { true, true, true }; DeviceMesh dm_; @@ -369,6 +489,7 @@ private: std::vector