OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-23 07:53:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4820 from atgeirr/parallel-nldd-fixes

Browse Source 
Fixes to enable using the NLDD nonlinear solver option in MPI parallel runs
This commit is contained in:
Bård Skaflestad 2023-09-01 09:51:06 +02:00 committed by GitHub
commit c6e30fbf3c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 252 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -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

32
ebos/FIBlackOilModel.hpp
View File

@ -58,6 +58,7 @@ public:
: BlackOilModel<TypeTag>(simulator)
{
}
void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
{
this->invalidateIntensiveQuantitiesCache(timeIdx);
@ -78,17 +79,44 @@ public:
}
}
void invalidateAndUpdateIntensiveQuantitiesOverlap(unsigned timeIdx) const
{
// loop over all elements
ThreadedEntityIterator<GridView, /*codim=*/0> 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 <class GridSubDomain>
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<GridViewType, /*codim=*/0> threadedElemIt(gridSubDomain.view);
#ifdef _OPENMP
#pragma omp parallel
#endif
{
ElementContext elemCtx(this->simulator_);
auto elemIt = threadedElemIt.beginParallel();
for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {

42
opm/simulators/flow/BlackoilModelEbosNldd.hpp
View File

@ -31,6 +31,7 @@
#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
#include <opm/simulators/flow/partitionCells.hpp>
#include <opm/simulators/flow/priVarsPacking.hpp>
#include <opm/simulators/flow/SubDomain.hpp>
#include <opm/simulators/linalg/extractMatrix.hpp>
@ -119,9 +120,11 @@ public:
}
// Iterate through grid once, setting the seeds of all partitions.
// Note: owned cells only!
std::vector<int> 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<std::size_t> 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<int>& 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</*codim=*/0>();
IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (auto elemIt = gridView.template begin</*codim=*/0>();
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();

54
opm/simulators/flow/priVarsPacking.hpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PRIVARSPACKING_HEADER_INCLUDED
#define OPM_PRIVARSPACKING_HEADER_INCLUDED
#include <cstddef>
namespace Opm {
namespace PVUtil {
constexpr int fbits = 4;
template <class PV>
std::size_t pack(const PV& privar) {
std::size_t m1 = static_cast<std::size_t>(privar.primaryVarsMeaningWater());
std::size_t m2 = static_cast<std::size_t>(privar.primaryVarsMeaningPressure());
std::size_t m3 = static_cast<std::size_t>(privar.primaryVarsMeaningGas());
std::size_t m4 = static_cast<std::size_t>(privar.primaryVarsMeaningBrine());
return m1 + (m2 << fbits*1) + (m3 << fbits*2) + (m4 << fbits*3);
}
template <class PV>
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

3
opm/simulators/linalg/ISTLSolverEbos.cpp
View File

@ -124,11 +124,12 @@ void FlexibleSolverInfo<Matrix,Vector,Comm>::create(const Matrix& matrix,
const PropertyTree& prm,
std::size_t pressureIndex,
std::function<Vector()> 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

14
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@ -103,6 +103,7 @@ struct FlexibleSolverInfo
const PropertyTree& prm,
std::size_t pressureIndex,
std::function<Vector()> trueFunc,
const bool forceSerial,
Comm& comm);
std::unique_ptr<AbstractSolverType> solver_;
@ -177,13 +178,14 @@ std::unique_ptr<Matrix> 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<Matrix> 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<Matrix> 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<int, int>;
@ -377,7 +381,7 @@ std::unique_ptr<Matrix> 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<Matrix> blockJacobiAdjacency(const Grid& grid,
prm_,
pressureIndex,
trueFunc,
forceSerial_,
*comm_);
}
else
@ -502,6 +507,7 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
bool useWellConn_;
FlowLinearSolverParameters parameters_;
bool forceSerial_ = false;
PropertyTree prm_;
std::shared_ptr< CommunicationType > comm_;

118
tests/test_privarspacking.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/flow/priVarsPacking.hpp>
#define BOOST_TEST_MODULE priVarsPacking
#include <boost/test/unit_test.hpp>
// 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));
}
}