diff --git a/CMakeLists.txt b/CMakeLists.txt index 4ad79d1..8d77084 100644 --- a/CMakeLists.txt +++ b/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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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 $ + 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} $ + 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) diff --git a/bench/setup_spuma.sh b/bench/setup_spuma.sh index fa39d9b..161b32b 100755 --- a/bench/setup_spuma.sh +++ b/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}" diff --git a/src/OpenFOAM/decompose/local_assembly.cuh b/src/OpenFOAM/decompose/local_assembly.cuh index 0981c66..302655f 100644 --- a/src/OpenFOAM/decompose/local_assembly.cuh +++ b/src/OpenFOAM/decompose/local_assembly.cuh @@ -60,6 +60,123 @@ 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; + } + } +} + +// 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& lpatches, + std::vector>& nod, + std::vector>& cv) +{ + 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); + 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(); + + 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) diff --git a/src/OpenFOAM/decompose/scotch_decomposition.cu b/src/OpenFOAM/decompose/scotch_decomposition.cu index a31b0cf..17f8722 100644 --- a/src/OpenFOAM/decompose/scotch_decomposition.cu +++ b/src/OpenFOAM/decompose/scotch_decomposition.cu @@ -3,6 +3,9 @@ #include // scotch.h uses FILE #include #include +#include +#include +#include namespace brae { @@ -67,4 +70,48 @@ std::vector