diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4c2387b3d..4e6aa2418 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -263,6 +263,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_parallelwellinfo.cpp tests/test_partitionCells.cpp tests/test_preconditionerfactory.cpp + tests/test_privarspacking.cpp tests/test_relpermdiagnostics.cpp tests/test_RestartSerialization.cpp tests/test_stoppedwells.cpp @@ -416,6 +417,7 @@ list (APPEND PUBLIC_HEADER_FILES ebos/hdf5serializer.hh ebos/vtkecltracermodule.hh opm/simulators/flow/countGlobalCells.hpp + opm/simulators/flow/priVarsPacking.hpp opm/simulators/flow/BlackoilModelEbos.hpp opm/simulators/flow/BlackoilModelEbosNldd.hpp opm/simulators/flow/BlackoilModelParametersEbos.hpp diff --git a/ebos/FIBlackOilModel.hpp b/ebos/FIBlackOilModel.hpp index 631a26f00..270d6c661 100644 --- a/ebos/FIBlackOilModel.hpp +++ b/ebos/FIBlackOilModel.hpp @@ -58,6 +58,7 @@ public: : BlackOilModel(simulator) { } + void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const { this->invalidateIntensiveQuantitiesCache(timeIdx); @@ -78,17 +79,44 @@ public: } } + void invalidateAndUpdateIntensiveQuantitiesOverlap(unsigned timeIdx) const + { + // loop over all elements + ThreadedEntityIterator threadedElemIt(this->gridView_); +#ifdef _OPENMP +#pragma omp parallel +#endif + { + ElementContext elemCtx(this->simulator_); + auto elemIt = threadedElemIt.beginParallel(); + for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) { + if (elemIt->partitionType() != Dune::OverlapEntity) { + continue; + } + const Element& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + // Mark cache for this element as invalid. + const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx); + for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) { + const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx); + this->setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false); + } + // Update for this element. + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + } + } + } + template void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridSubDomain& gridSubDomain) const { - // loop over all elements... + // loop over all elements in the subdomain using GridViewType = decltype(gridSubDomain.view); ThreadedEntityIterator threadedElemIt(gridSubDomain.view); #ifdef _OPENMP #pragma omp parallel #endif { - ElementContext elemCtx(this->simulator_); auto elemIt = threadedElemIt.beginParallel(); for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) { diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp index 9dc4b6ef9..7032c8a5a 100644 --- a/opm/simulators/flow/BlackoilModelEbosNldd.hpp +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -31,6 +31,7 @@ #include #include +#include #include #include @@ -119,9 +120,11 @@ public: } // Iterate through grid once, setting the seeds of all partitions. + // Note: owned cells only! std::vector count(num_domains, 0); - const auto beg = grid.template leafbegin<0>(); - const auto end = grid.template leafend<0>(); + const auto& gridView = grid.leafGridView(); + const auto beg = gridView.template begin<0, Dune::Interior_Partition>(); + const auto end = gridView.template end<0, Dune::Interior_Partition>(); int cell = 0; for (auto it = beg; it != end; ++it, ++cell) { const int p = partition_vector[cell]; @@ -166,7 +169,8 @@ public: loc_param.linear_solver_reduction_ = 1e-2; } loc_param.linear_solver_print_json_definition_ = false; - domain_linsolvers_.emplace_back(model_.ebosSimulator(), loc_param); + const bool force_serial = true; + domain_linsolvers_.emplace_back(model_.ebosSimulator(), loc_param, force_serial); } assert(int(domains_.size()) == num_domains); @@ -253,6 +257,35 @@ public: model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); } + // Communicate solutions: + // With multiple processes, this process' overlap (i.e. not + // owned) cells' solution values have been modified by local + // solves in the owning processes, and remain unchanged + // here. We must therefore receive the updated solution on the + // overlap cells and update their intensive quantities before + // we move on. + const auto& comm = model_.ebosSimulator().vanguard().grid().comm(); + if (comm.size() > 1) { + const auto* ccomm = model_.ebosSimulator().model().newtonMethod().linearSolver().comm(); + + // Copy numerical values from primary vars. + ccomm->copyOwnerToAll(solution, solution); + + // Copy flags from primary vars. + const std::size_t num = solution.size(); + Dune::BlockVector allmeanings(num); + for (std::size_t ii = 0; ii < num; ++ii) { + allmeanings[ii] = PVUtil::pack(solution[ii]); + } + ccomm->copyOwnerToAll(allmeanings, allmeanings); + for (std::size_t ii = 0; ii < num; ++ii) { + PVUtil::unPack(solution[ii], allmeanings[ii]); + } + + // Update intensive quantities for our overlap values. + model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantitiesOverlap(/*timeIdx=*/0); + } + // Finish with a Newton step. // Note that the "iteration + 100" is a simple way to avoid entering // "if (iteration == 0)" and similar blocks, and also makes it a little @@ -458,7 +491,6 @@ private: std::vector& maxCoeffCell) { const auto& ebosSimulator = model_.ebosSimulator(); - const auto& grid = model_.ebosSimulator().vanguard().grid(); double pvSumLocal = 0.0; double numAquiferPvSumLocal = 0.0; @@ -472,7 +504,6 @@ private: const auto& elemEndIt = gridView.template end(); IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid()); - OPM_BEGIN_PARALLEL_TRY_CATCH(); for (auto elemIt = gridView.template begin(); elemIt != elemEndIt; ++elemIt) @@ -500,7 +531,6 @@ private: model_.getMaxCoeff(cell_idx, intQuants, fs, ebosResid, pvValue, B_avg, R_sum, maxCoeff, maxCoeffCell); } - OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid.comm()); // compute local average in terms of global number of elements const int bSize = B_avg.size(); diff --git a/opm/simulators/flow/priVarsPacking.hpp b/opm/simulators/flow/priVarsPacking.hpp new file mode 100644 index 000000000..778fe54b6 --- /dev/null +++ b/opm/simulators/flow/priVarsPacking.hpp @@ -0,0 +1,54 @@ +/* + Copyright 2023 Total SE + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_PRIVARSPACKING_HEADER_INCLUDED +#define OPM_PRIVARSPACKING_HEADER_INCLUDED + +#include + +namespace Opm { + + namespace PVUtil { + constexpr int fbits = 4; + + template + std::size_t pack(const PV& privar) { + std::size_t m1 = static_cast(privar.primaryVarsMeaningWater()); + std::size_t m2 = static_cast(privar.primaryVarsMeaningPressure()); + std::size_t m3 = static_cast(privar.primaryVarsMeaningGas()); + std::size_t m4 = static_cast(privar.primaryVarsMeaningBrine()); + return m1 + (m2 << fbits*1) + (m3 << fbits*2) + (m4 << fbits*3); + } + + template + void unPack(PV& privar, const std::size_t meanings) { + const std::size_t filter = ((1 << fbits) - 1); + std::size_t m1 = (meanings >> fbits*0) & filter; + std::size_t m2 = (meanings >> fbits*1) & filter; + std::size_t m3 = (meanings >> fbits*2) & filter; + std::size_t m4 = (meanings >> fbits*3) & filter; + privar.setPrimaryVarsMeaningWater(typename PV::WaterMeaning(m1)); + privar.setPrimaryVarsMeaningPressure(typename PV::PressureMeaning(m2)); + privar.setPrimaryVarsMeaningGas(typename PV::GasMeaning(m3)); + privar.setPrimaryVarsMeaningBrine(typename PV::BrineMeaning(m4)); + } + } // namespace PVMeanings +} // namespace Opm + +#endif // OPM_PRIVARSPACKING_HEADER_INCLUDED diff --git a/opm/simulators/linalg/ISTLSolverEbos.cpp b/opm/simulators/linalg/ISTLSolverEbos.cpp index 25b3d2cdd..241df3fb0 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.cpp +++ b/opm/simulators/linalg/ISTLSolverEbos.cpp @@ -124,11 +124,12 @@ void FlexibleSolverInfo::create(const Matrix& matrix, const PropertyTree& prm, std::size_t pressureIndex, std::function trueFunc, + const bool forceSerial, [[maybe_unused]] Comm& comm) { // Write sizes of linear systems on all ranks to debug log. - { + if (!forceSerial) { #if HAVE_MPI auto basic_comm = comm.communicator(); #else diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 6f26d8c07..be1b6c838 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -103,6 +103,7 @@ struct FlexibleSolverInfo const PropertyTree& prm, std::size_t pressureIndex, std::function trueFunc, + const bool forceSerial, Comm& comm); std::unique_ptr solver_; @@ -177,13 +178,14 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, /// \param[in] simulator The opm-models simulator object /// \param[in] parameters Explicit parameters for solver setup, do not /// read them from command line parameters. - ISTLSolverEbos(const Simulator& simulator, const FlowLinearSolverParameters& parameters) + ISTLSolverEbos(const Simulator& simulator, const FlowLinearSolverParameters& parameters, bool forceSerial = false) : simulator_(simulator), iterations_( 0 ), calls_( 0 ), converged_(false), matrix_(nullptr), - parameters_(parameters) + parameters_(parameters), + forceSerial_(forceSerial) { initialize(); } @@ -261,7 +263,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, OPM_TIMEBLOCK(istlSolverEbosPrepare); const bool firstcall = (matrix_ == nullptr); #if HAVE_MPI - if (firstcall) { + if (firstcall && isParallel()) { const std::size_t size = M.N(); detail::copyParValues(parallelInformation_, size, *comm_); } @@ -348,6 +350,8 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, /// \copydoc NewtonIterationBlackoilInterface::parallelInformation const std::any& parallelInformation() const { return parallelInformation_; } + const CommunicationType* comm() const { return comm_.get(); } + protected: #if HAVE_MPI using Comm = Dune::OwnerOverlapCopyCommunication; @@ -377,7 +381,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, bool isParallel() const { #if HAVE_MPI - return comm_->communicator().size() > 1; + return !forceSerial_ && comm_->communicator().size() > 1; #else return false; #endif @@ -403,6 +407,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, prm_, pressureIndex, trueFunc, + forceSerial_, *comm_); } else @@ -502,6 +507,7 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, bool useWellConn_; FlowLinearSolverParameters parameters_; + bool forceSerial_ = false; PropertyTree prm_; std::shared_ptr< CommunicationType > comm_; diff --git a/tests/test_privarspacking.cpp b/tests/test_privarspacking.cpp new file mode 100644 index 000000000..b4cd37680 --- /dev/null +++ b/tests/test_privarspacking.cpp @@ -0,0 +1,118 @@ +/* + Copyright 2023 SINTEF. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include + +#include + +#define BOOST_TEST_MODULE priVarsPacking +#include + + +// Must define a class for testing, using extracts from BlackoilPrimaryVariables, +// but without the typetags. +class PriVarMeaning +{ +public: + + enum class WaterMeaning { + Sw, // water saturation + Rvw, // vaporized water + Rsw, // dissolved gas in water + Disabled, // The primary variable is not used + }; + + enum class PressureMeaning { + Po, // oil pressure + Pg, // gas pressure + Pw, // water pressure + }; + enum class GasMeaning { + Sg, // gas saturation + Rs, // dissolved gas in oil + Rv, // vapporized oil + Disabled, // The primary variable is not used + }; + + enum class BrineMeaning { + Cs, // salt concentration + Sp, // (precipitated) salt saturation + Disabled, // The primary variable is not used + }; + + WaterMeaning primaryVarsMeaningWater() const + { return primaryVarsMeaningWater_; } + void setPrimaryVarsMeaningWater(WaterMeaning newMeaning) + { primaryVarsMeaningWater_ = newMeaning; } + + PressureMeaning primaryVarsMeaningPressure() const + { return primaryVarsMeaningPressure_; } + void setPrimaryVarsMeaningPressure(PressureMeaning newMeaning) + { primaryVarsMeaningPressure_ = newMeaning; } + + GasMeaning primaryVarsMeaningGas() const + { return primaryVarsMeaningGas_; } + void setPrimaryVarsMeaningGas(GasMeaning newMeaning) + { primaryVarsMeaningGas_ = newMeaning; } + + BrineMeaning primaryVarsMeaningBrine() const + { return primaryVarsMeaningBrine_; } + void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning) + { primaryVarsMeaningBrine_ = newMeaning; } + + bool operator==(const PriVarMeaning& other) const + { + return primaryVarsMeaningWater_ == other.primaryVarsMeaningWater_ + && primaryVarsMeaningPressure_ == other.primaryVarsMeaningPressure_ + && primaryVarsMeaningGas_ == other.primaryVarsMeaningGas_ + && primaryVarsMeaningBrine_ == other.primaryVarsMeaningBrine_; + } +private: + WaterMeaning primaryVarsMeaningWater_ = WaterMeaning::Disabled; + PressureMeaning primaryVarsMeaningPressure_ = PressureMeaning::Pw; + GasMeaning primaryVarsMeaningGas_ = GasMeaning::Disabled; + BrineMeaning primaryVarsMeaningBrine_ = BrineMeaning::Disabled; +}; + + + +BOOST_AUTO_TEST_CASE(meanings) +{ + // Test default meanings. + { + PriVarMeaning pv1, pv2; + std::size_t p = Opm::PVUtil::pack(pv1); + Opm::PVUtil::unPack(pv2, p); + BOOST_CHECK(pv1 == pv2); + } + + // Test explicitly set meanings. + { + PriVarMeaning pv1, pv2, pv3; + pv1.setPrimaryVarsMeaningPressure(PriVarMeaning::PressureMeaning::Pw); + pv1.setPrimaryVarsMeaningWater(PriVarMeaning::WaterMeaning::Rvw); + pv1.setPrimaryVarsMeaningGas(PriVarMeaning::GasMeaning::Disabled); + pv1.setPrimaryVarsMeaningBrine(PriVarMeaning::BrineMeaning::Cs); + std::size_t p = Opm::PVUtil::pack(pv1); + Opm::PVUtil::unPack(pv2, p); + BOOST_CHECK(pv1 == pv2); + BOOST_CHECK(!(pv1 == pv3)); + BOOST_CHECK(!(pv2 == pv3)); + } +}