diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.cpp b/opm/simulators/linalg/FlowLinearSolverParameters.cpp index 724076038a..1b3641f980 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.cpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.cpp @@ -72,6 +72,7 @@ void FlowLinearSolverParameters::init(bool cprRequestedInDataFile) opencl_platform_id_ = Parameters::Get(); opencl_ilu_parallel_ = Parameters::Get(); linear_solver_accelerator_ = Parameters::linearSolverAcceleratorTypeFromCLI(); + cpr_weights_thread_parallel_ = Parameters::Get(); if (linear_solver_accelerator_ == Parameters::LinearSolverAcceleratorType::GPU) { if (!Parameters::IsSet()) { @@ -176,6 +177,10 @@ void FlowLinearSolverParameters::registerParameters() "Note that the verification is not exhaustive, " "and some configurations will not verify, but will work in practice. " "Usage: --verify-gpu-aware-mpi=[true|false]. "); + Parameters::Register + ("Enable OpenMP thread parallelization of CPR weight calculation. " + "This can improve performance for large models but is disabled by default." + "Usage: --cpr-weights-thread-parallel=[true|false]. "); Parameters::SetDefault(0); } @@ -207,6 +212,7 @@ void FlowLinearSolverParameters::reset() linear_solver_accelerator_ = Parameters::LinearSolverAcceleratorType::CPU; gpu_aware_mpi_ = false; verify_gpu_aware_mpi_ = false; + cpr_weights_thread_parallel_ = false; } } // namespace Opm diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.hpp b/opm/simulators/linalg/FlowLinearSolverParameters.hpp index 0fba75be68..2cbf6e6568 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.hpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.hpp @@ -88,6 +88,7 @@ struct OpenclPlatformId { static constexpr int value = 0; }; struct OpenclIluParallel { static constexpr bool value = true; }; // note: false should only be used in debug struct GpuAwareMpi { static constexpr bool value = false; }; struct VerifyGpuAwareMpi { static constexpr bool value = false; }; +struct CprWeightsThreadParallel { static constexpr bool value = false; }; } // namespace Opm::Parameters namespace Opm { @@ -120,6 +121,7 @@ struct FlowLinearSolverParameters Parameters::LinearSolverAcceleratorType linear_solver_accelerator_; bool gpu_aware_mpi_; bool verify_gpu_aware_mpi_; + bool cpr_weights_thread_parallel_; FlowLinearSolverParameters() { reset(); } diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 7a6574a25e..7a5223f7ba 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -580,19 +580,21 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, if (preconditionerType == "cpr" || preconditionerType == "cprt" || preconditionerType == "cprw" || preconditionerType == "cprwt") { const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt"; + const bool enableThreadParallel = this->parameters_[0].cpr_weights_thread_parallel_; const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s); if (weightsType == "quasiimpes") { // weights will be created as default in the solver // assignment p = pressureIndex prevent compiler warning about // capturing variable with non-automatic storage duration - weightsCalculator = [matrix, transpose, pressIndex]() { + weightsCalculator = [matrix, transpose, pressIndex, enableThreadParallel]() { return Amg::getQuasiImpesWeights(matrix, pressIndex, - transpose); + transpose, + enableThreadParallel); }; } else if ( weightsType == "trueimpes" ) { weightsCalculator = - [this, pressIndex] + [this, pressIndex, enableThreadParallel] { Vector weights(rhs_->size()); ElementContext elemCtx(simulator_); @@ -600,13 +602,14 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, weights, elemCtx, simulator_.model(), - *element_chunks_ + *element_chunks_, + enableThreadParallel ); return weights; }; } else if (weightsType == "trueimpesanalytic" ) { weightsCalculator = - [this, pressIndex] + [this, pressIndex, enableThreadParallel] { Vector weights(rhs_->size()); ElementContext elemCtx(simulator_); @@ -614,7 +617,8 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, weights, elemCtx, simulator_.model(), - *element_chunks_ + *element_chunks_, + enableThreadParallel ); return weights; }; diff --git a/opm/simulators/linalg/getQuasiImpesWeights.hpp b/opm/simulators/linalg/getQuasiImpesWeights.hpp index 8c0731bcd3..3b40b3fd0b 100644 --- a/opm/simulators/linalg/getQuasiImpesWeights.hpp +++ b/opm/simulators/linalg/getQuasiImpesWeights.hpp @@ -58,7 +58,7 @@ namespace Details namespace Amg { template - void getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose, Vector& weights) + void getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose, Vector& weights, bool enable_thread_parallel) { using VectorBlockType = typename Vector::block_type; using MatrixBlockType = typename Matrix::block_type; @@ -72,10 +72,10 @@ namespace Amg VectorBlockType bweights; MatrixBlockType diag_block_transpose; - // Use OpenMP to parallelize over matrix rows - #ifdef _OPENMP - #pragma omp parallel for private(diag_block, bweights, diag_block_transpose) - #endif + // Use OpenMP to parallelize over matrix rows (runtime controlled via if clause) +#ifdef _OPENMP +#pragma omp parallel for private(diag_block, bweights, diag_block_transpose) if(enable_thread_parallel) +#endif for (int row_idx = 0; row_idx < static_cast(A.N()); ++row_idx) { diag_block = MatrixBlockType(0.0); // Find diagonal block for this row @@ -102,10 +102,10 @@ namespace Amg } template - Vector getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose) + Vector getQuasiImpesWeights(const Matrix& matrix, const int pressureVarIndex, const bool transpose, bool enable_thread_parallel) { Vector weights(matrix.N()); - getQuasiImpesWeights(matrix, pressureVarIndex, transpose, weights); + getQuasiImpesWeights(matrix, pressureVarIndex, transpose, weights, enable_thread_parallel); return weights; } @@ -152,7 +152,8 @@ namespace Amg void getTrueImpesWeights(int pressureVarIndex, Vector& weights, const ElementContext& elemCtx, const Model& model, - const ElementChunksType& element_chunks) + const ElementChunksType& element_chunks, + bool enable_thread_parallel) { using VectorBlockType = typename Vector::block_type; using Matrix = typename std::decay_t; @@ -171,9 +172,9 @@ namespace Amg Dune::FieldVector storage; OPM_BEGIN_PARALLEL_TRY_CATCH(); - #ifdef _OPENMP - #pragma omp parallel for private(block, bweights, block_transpose, storage) - #endif +#ifdef _OPENMP +#pragma omp parallel for private(block, bweights, block_transpose, storage) if(enable_thread_parallel) +#endif for (const auto& chunk : element_chunks) { const std::size_t thread_id = ThreadManager::threadId(); ElementContext localElemCtx(elemCtx.simulator()); @@ -217,7 +218,8 @@ namespace Amg Vector& weights, const ElementContext& elemCtx, const Model& model, - const ElementChunksType& element_chunks) + const ElementChunksType& element_chunks, + bool enable_thread_parallel) { // The sequential residual is a linear combination of the // mass balance residuals, with coefficients equal to (for @@ -239,11 +241,11 @@ namespace Amg const auto& solution = model.solution(/*timeIdx*/ 0); VectorBlockType bweights; - // Use OpenMP to parallelize over element chunks + // Use OpenMP to parallelize over element chunks (runtime controlled via if clause) OPM_BEGIN_PARALLEL_TRY_CATCH(); - #ifdef _OPENMP - #pragma omp parallel for private(bweights) - #endif +#ifdef _OPENMP +#pragma omp parallel for private(bweights) if(enable_thread_parallel) +#endif for (const auto& chunk : element_chunks) { // Each thread gets a unique copy of elemCtx diff --git a/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp b/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp index 33c06a4888..68c89c86e3 100644 --- a/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp +++ b/opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp @@ -375,14 +375,16 @@ private: const_cast(&m_cpuWeights[0][0]), m_cpuWeights.dim()); m_weights.emplace(m_cpuWeights); + const bool enableThreadParallel = m_parameters.cpr_weights_thread_parallel_; // CPU implementation wrapped for GPU - weightsCalculator = [this]() -> GPUVector& { + weightsCalculator = [this, enableThreadParallel]() -> GPUVector& { // Use the CPU implementation to calculate the weights ElementContext elemCtx(m_simulator); Amg::getTrueImpesWeights(pressureIndex, m_cpuWeights, elemCtx, m_simulator.model(), - m_element_chunks); + m_element_chunks, + enableThreadParallel); // Copy CPU vector to GPU vector using main stream and asynchronous transfer m_weights->copyFromHostAsync(m_cpuWeights); @@ -394,14 +396,16 @@ private: m_pinnedWeightsMemory = std::make_unique>( const_cast(&m_cpuWeights[0][0]), m_cpuWeights.dim()); m_weights.emplace(m_cpuWeights); + const bool enableThreadParallel = m_parameters.cpr_weights_thread_parallel_; // CPU implementation wrapped for GPU - weightsCalculator = [this]() -> GPUVector& { + weightsCalculator = [this, enableThreadParallel]() -> GPUVector& { // Use the CPU implementation to calculate the weights ElementContext elemCtx(m_simulator); Amg::getTrueImpesWeightsAnalytic(pressureIndex, m_cpuWeights, elemCtx, m_simulator.model(), - m_element_chunks); + m_element_chunks, + enableThreadParallel); // Copy CPU vector to GPU vector using main stream and asynchronous transfer m_weights->copyFromHostAsync(m_cpuWeights); diff --git a/tests/gpuistl/test_GpuPressureTransferPolicy.cpp b/tests/gpuistl/test_GpuPressureTransferPolicy.cpp index e14c5607b9..471b677893 100644 --- a/tests/gpuistl/test_GpuPressureTransferPolicy.cpp +++ b/tests/gpuistl/test_GpuPressureTransferPolicy.cpp @@ -132,7 +132,7 @@ testQuasiImpesWeights(int pressureVarIndex) // Calculate quasiimpes weights for CPU auto cpuWeights = Opm::Amg::getQuasiImpesWeights::MatrixType, typename TestFixture::VectorType>( - fixture.matrix, pressureVarIndex, transpose); + fixture.matrix, pressureVarIndex, transpose, false); // Calculate quasiimpes weights for GPU auto gpuMatrix = fixture.createGpuMatrix(); @@ -165,7 +165,7 @@ testCoarseMatrixCreation(int pressureVarIndex) // Calculate weights for CPU auto cpuWeights = Opm::Amg::getQuasiImpesWeights::MatrixType, typename TestFixture::VectorType>( - fixture.matrix, pressureVarIndex, transpose); + fixture.matrix, pressureVarIndex, transpose, false); // Calculate quasiimpes weights for GPU auto gpuMatrix = fixture.createGpuMatrix(); @@ -234,7 +234,7 @@ testRestriction(int pressureVarIndex) // Calculate weights for CPU auto cpuWeights = Opm::Amg::getQuasiImpesWeights::MatrixType, typename TestFixture::VectorType>( - fixture.matrix, pressureVarIndex, transpose); + fixture.matrix, pressureVarIndex, transpose, false); // Calculate quasiimpes weights for GPU auto gpuMatrix = fixture.createGpuMatrix(); @@ -307,7 +307,7 @@ testProlongation(int pressureVarIndex) // Calculate weights and create policies auto cpuWeights = Opm::Amg::getQuasiImpesWeights::MatrixType, typename TestFixture::VectorType>( - fixture.matrix, pressureVarIndex, transpose); + fixture.matrix, pressureVarIndex, transpose, false); // Calculate quasiimpes weights for GPU auto gpuMatrix = fixture.createGpuMatrix(); diff --git a/tests/test_blackoil_amg.cpp b/tests/test_blackoil_amg.cpp index 86d6473871..2e3ee68109 100644 --- a/tests/test_blackoil_amg.cpp +++ b/tests/test_blackoil_amg.cpp @@ -295,7 +295,7 @@ void runBlackoilAmgLaplace() Opm::PropertyTree prm; prm.put("type", "amg"s); std::function weights = [&mat]() { - return Opm::Amg::getQuasiImpesWeights(mat, 0, false); + return Opm::Amg::getQuasiImpesWeights(mat, 0, false, false); }; auto amg = Opm::PreconditionerFactory::create(fop, prm, weights, comm); diff --git a/tests/test_flexiblesolver.cpp b/tests/test_flexiblesolver.cpp index a648ed05c4..4781405111 100644 --- a/tests/test_flexiblesolver.cpp +++ b/tests/test_flexiblesolver.cpp @@ -63,9 +63,9 @@ testSolver(const Opm::PropertyTree& prm, const std::string& matrix_filename, con if(prm.get("preconditioner.type") == "cprt"){ transpose = true; } - auto wc = [&matrix, transpose]() + std::function wc = [&matrix, transpose]() { - return Opm::Amg::getQuasiImpesWeights(matrix, 1, transpose); + return Opm::Amg::getQuasiImpesWeights(matrix, 1, transpose, false); }; using SeqOperatorType = Dune::MatrixAdapter; diff --git a/tests/test_ghostlastmatrixadapter.cpp b/tests/test_ghostlastmatrixadapter.cpp index 7ce57068c7..63ac689720 100644 --- a/tests/test_ghostlastmatrixadapter.cpp +++ b/tests/test_ghostlastmatrixadapter.cpp @@ -292,7 +292,7 @@ int main(int argc, char** argv) Opm::PropertyTree prm; prm.put("type", "Jac"s); std::function weights = [&mat]() { - return Opm::Amg::getQuasiImpesWeights(mat, 0, false); + return Opm::Amg::getQuasiImpesWeights(mat, 0, false, false); }; auto fullPreconditioner = Opm::PreconditionerFactory::create(op, prm, weights, comm); diff --git a/tests/test_preconditionerfactory.cpp b/tests/test_preconditionerfactory.cpp index 31c0d05b1f..2ec020ef0b 100644 --- a/tests/test_preconditionerfactory.cpp +++ b/tests/test_preconditionerfactory.cpp @@ -95,9 +95,9 @@ testPrec(const Opm::PropertyTree& prm, const std::string& matrix_filename, const if(prm.get("preconditioner.type") == "cprt"){ transpose = true; } - auto wc = [&matrix, transpose]() + std::function wc = [&matrix, transpose]() { - return Opm::Amg::getQuasiImpesWeights(matrix, 1, transpose); + return Opm::Amg::getQuasiImpesWeights(matrix, 1, transpose, false); }; auto prec = PrecFactory::create(op, prm.get_child("preconditioner"), wc, 1);