OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2987 from akva2/parallel_aquifers

Browse Source 
Add supports for aquifers in parallel
This commit is contained in:
Arne Morten Kvarving
2020-12-21 13:53:24 +01:00
committed by GitHub
parent eb4aabc71a b812088519
commit 560a064ea7
7 changed files with 205 additions and 194 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
compareECLFiles.cmake
View File

@@ -945,4 +945,20 @@ if(MPI_FOUND)
ABS_TOL ${abs_tol_parallel}
REL_TOL ${rel_tol_parallel}
TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
add_test_compare_parallel_simulation(CASENAME fetkovich_2d
FILENAME 2D_FETKOVICHAQUIFER
SIMULATOR flow
ABS_TOL ${abs_tol_parallel}
REL_TOL ${rel_tol_parallel}
DIR aquifer-fetkovich
TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
add_test_compare_parallel_simulation(CASENAME ctaquifer_2d_oilwater
FILENAME 2D_OW_CTAQUIFER
SIMULATOR flow
ABS_TOL ${abs_tol_parallel}
REL_TOL ${rel_tol_parallel}
DIR aquifer-oilwater
TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
endif()

74
ebos/collecttoiorank.hh
View File

@@ -708,19 +708,80 @@ public:
}
}
};
class PackUnPackAquiferData : public P2PCommunicatorType::DataHandleInterface
{
const Opm::data::Aquifers& localAquiferData_;
data::Aquifers& globalAquiferData_;
public:
PackUnPackAquiferData(const Opm::data::Aquifers& localAquiferData,
Opm::data::Aquifers& globalAquiferData,
bool isIORank)
: localAquiferData_(localAquiferData)
, globalAquiferData_(globalAquiferData)
{
if (isIORank) {
MessageBufferType buffer;
pack(0, buffer);
// pass a dummyLink to satisfy virtual class
int dummyLink = -1;
unpack(dummyLink, buffer);
}
}
// pack all data associated with link
void pack(int link, MessageBufferType& buffer)
{
// we should only get one link
if (link != 0)
throw std::logic_error("link in method pack is not 0 as expected");
int size = localAquiferData_.size();
buffer.write(size);
for (const auto& [key, data] : localAquiferData_) {
buffer.write(key);
data.write(buffer);
}
}
// unpack all data associated with link
void unpack(int /*link*/, MessageBufferType& buffer)
{
int size;
buffer.read(size);
for (int i = 0; i < size; ++i) {
int key;
buffer.read(key);
data::AquiferData data;
data.read(buffer);
auto& aq = globalAquiferData_[key];
aq.aquiferID = std::max(aq.aquiferID, data.aquiferID);
aq.pressure = std::max(aq.pressure, data.pressure);
aq.initPressure = std::max(aq.initPressure, data.initPressure);
aq.datumDepth = std::max(aq.datumDepth, data.datumDepth);
aq.fluxRate += data.fluxRate;
aq.volume += data.volume;
}
}
};
// gather solution to rank 0 for EclipseWriter
void collect(const Opm::data::Solution& localCellData,
const std::map<std::pair<std::string, int>, double>& localBlockData,
const std::map<std::size_t, double>& localWBPData,
const Opm::data::Wells& localWellData,
const Opm::data::GroupAndNetworkValues& localGroupAndNetworkData)
const Opm::data::GroupAndNetworkValues& localGroupAndNetworkData,
const Opm::data::Aquifers& localAquiferData)
{
globalCellData_ = {};
globalBlockData_.clear();
globalWBPData_.clear();
globalWellData_.clear();
globalGroupAndNetworkData_.clear();
globalAquiferData_.clear();
// index maps only have to be build when reordering is needed
if(!needsReordering && !isParallel())
@@ -765,11 +826,18 @@ public:
this->isIORank()
};
PackUnPackAquiferData packUnpackAquiferData {
localAquiferData,
this->globalAquiferData_,
this->isIORank()
};
toIORankComm_.exchange(packUnpackCellData);
toIORankComm_.exchange(packUnpackWellData);
toIORankComm_.exchange(packUnpackGroupAndNetworkData);
toIORankComm_.exchange(packUnpackBlockData);
toIORankComm_.exchange(packUnpackWBPData);
toIORankComm_.exchange(packUnpackAquiferData);
#ifndef NDEBUG
// mkae sure every process is on the same page
@@ -792,6 +860,9 @@ public:
const Opm::data::GroupAndNetworkValues& globalGroupAndNetworkData() const
{ return globalGroupAndNetworkData_; }
const Opm::data::Aquifers& globalAquiferData() const
{ return globalAquiferData_; }
bool isIORank() const
{ return toIORankComm_.rank() == ioRank; }
@@ -839,6 +910,7 @@ protected:
std::map<std::size_t, double> globalWBPData_;
Opm::data::Wells globalWellData_;
Opm::data::GroupAndNetworkValues globalGroupAndNetworkData_;
Opm::data::Aquifers globalAquiferData_;
std::vector<int> localIdxToGlobalIdx_;
/// \brief sorted list of cartesian indices present-
///

12
ebos/eclwriter.hh
View File

@@ -295,7 +295,9 @@ public:
eclOutputModule_->getBlockData(),
eclOutputModule_->getWBPData(),
localWellData,
localGroupAndNetworkData);
localGroupAndNetworkData,
localAquiferData);
std::map<std::string, double> miscSummaryData;
std::map<std::string, std::vector<double>> regionData;
@@ -334,8 +336,9 @@ public:
? this->collectToIORank_.globalGroupAndNetworkData()
: localGroupAndNetworkData;
// Aquifer can not be parallel running yet
const auto& aquiferData = localAquiferData;
const auto& aquiferData = this->collectToIORank_.isParallel()
? this->collectToIORank_.globalAquiferData()
: localAquiferData;
const auto& blockData
= this->collectToIORank_.isParallel()
@@ -410,7 +413,8 @@ public:
eclOutputModule_->getBlockData(),
eclOutputModule_->getWBPData(),
localWellData,
localGroupAndNetworkData);
localGroupAndNetworkData,
{});
}
if (this->collectToIORank_.isIORank()) {

94
opm/simulators/aquifers/AquiferCarterTracy.hpp
View File

@@ -63,6 +63,9 @@ public:
for (const auto& q : this->Qai_) {
this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
}
this->fluxValue_ = this->W_flux_.value();
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&this->fluxValue_, 1);
}
Opm::data::AquiferData aquiferData() const
@@ -87,94 +90,7 @@ protected:
Scalar beta_; // Influx constant
// TODO: it is possible it should be a AD variable
Scalar mu_w_; // water viscosity
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer
inline void initializeConnections() override
{
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
// Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection;
// denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas = 0.;
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < this->size(); ++idx) {
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
if (cell_index < 0) //the global_index is not part of this grid
continue;
this->cellToConnectionIdx_[cell_index] = idx;
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get default areas for all intersections
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
// only deal with connections given by the aquifer
if( idx < 0)
continue;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// only deal with grid boundaries
if (!intersection.boundary())
continue;
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = Opm::FaceDir::XMinus;
break;
case 1:
faceDirection = Opm::FaceDir::XPlus;
break;
case 2:
faceDirection = Opm::FaceDir::YMinus;
break;
case 3:
faceDirection = Opm::FaceDir::YPlus;
break;
case 4:
faceDirection = Opm::FaceDir::ZMinus;
break;
case 5:
faceDirection = Opm::FaceDir::ZPlus;
break;
default:
OPM_THROW(Opm::NumericalIssue,
"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined");
}
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
}
}
}
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (size_t idx = 0; idx < this->size(); ++idx) {
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
? // Prevent no connection NaNs due to division by zero
0.
: (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx)) / denom_face_areas;
}
}
Scalar fluxValue_; // value of flux
void assignRestartData(const data::AquiferData& /* xaq */) override
{
@@ -205,7 +121,7 @@ protected:
Scalar PItdprime = 0.;
Scalar PItd = 0.;
getInfluenceTableValues(PItd, PItdprime, td_plus_dt);
a = 1.0 / this->Tc_ * ((beta_ * dpai(idx)) - (this->W_flux_.value() * PItdprime)) / (PItd - td * PItdprime);
a = 1.0 / this->Tc_ * ((beta_ * dpai(idx)) - (this->fluxValue_ * PItdprime)) / (PItd - td * PItdprime);
b = beta_ / (this->Tc_ * (PItd - td * PItdprime));
}

95
opm/simulators/aquifers/AquiferFetkovich.hpp
View File

@@ -64,8 +64,8 @@ public:
{
for (const auto& q : this->Qai_) {
this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
aquifer_pressure_ = aquiferPressure();
}
aquifer_pressure_ = aquiferPressure();
}
Opm::data::AquiferData aquiferData() const
@@ -92,97 +92,6 @@ protected:
const Aquifetp::AQUFETP_data aqufetp_data_;
Scalar aquifer_pressure_; // aquifer
inline void initializeConnections() override
{
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
// Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection;
// denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas = 0.;
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
for (size_t idx = 0; idx < this->size(); ++idx) {
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
if (cell_index < 0) //the global_index is not part of this grid
continue;
this->cellToConnectionIdx_[cell_index] = idx;
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get areas for all connections
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
// only deal with connections given by the aquifer
if( idx < 0)
continue;
if (!this->connections_[idx].influx_coeff.first) { // influx_coeff is defaulted
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// only deal with grid boundaries
if (!intersection.boundary())
continue;
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = Opm::FaceDir::XMinus;
break;
case 1:
faceDirection = Opm::FaceDir::XPlus;
break;
case 2:
faceDirection = Opm::FaceDir::YMinus;
break;
case 3:
faceDirection = Opm::FaceDir::YPlus;
break;
case 4:
faceDirection = Opm::FaceDir::ZMinus;
break;
case 5:
faceDirection = Opm::FaceDir::ZPlus;
break;
default:
OPM_THROW(Opm::NumericalIssue,
"Initialization of Aquifer Fetkovich problem. Make sure faceTag is correctly defined"); }
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
break;
}
}
} else {
this->faceArea_connected_.at(idx) = this->connections_[idx].influx_coeff.second;
}
denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
}
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (size_t idx = 0; idx < this->size(); ++idx) {
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
? // Prevent no connection NaNs due to division by zero
0.
: (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx)) / denom_face_areas;
}
}
void assignRestartData(const data::AquiferData& xaq) override
{
if (xaq.type != data::AquiferType::Fetkovich) {
@@ -204,6 +113,8 @@ protected:
inline Scalar aquiferPressure()
{
Scalar Flux = this->W_flux_.value();
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&Flux, 1);
Scalar pa_ = this->pa0_ - Flux / (aqufetp_data_.C_t * aqufetp_data_.V0);
return pa_;
}

104
opm/simulators/aquifers/AquiferInterface.hpp
View File

@@ -233,7 +233,100 @@ protected:
bool solution_set_from_restart_ {false};
virtual void initializeConnections() = 0;
void initializeConnections()
{
this->cell_depth_.resize(this->size(), this->aquiferDepth());
this->alphai_.resize(this->size(), 1.0);
this->faceArea_connected_.resize(this->size(), 0.0);
// Translate the C face tag into the enum used by opm-parser's TransMult class
Opm::FaceDir::DirEnum faceDirection;
// denom_face_areas is the sum of the areas connected to an aquifer
Scalar denom_face_areas = 0.;
this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
const auto& gridView = this->ebos_simulator_.vanguard().gridView();
for (size_t idx = 0; idx < this->size(); ++idx) {
const auto global_index = this->connections_[idx].global_index;
const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
auto elemIt = gridView.template begin</*codim=*/ 0>();
if (cell_index > 0)
std::advance(elemIt, cell_index);
//the global_index is not part of this grid
if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
continue;
this->cellToConnectionIdx_[cell_index] = idx;
this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
}
// get areas for all connections
ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
unsigned cell_index = elemMapper.index(elem);
int idx = this->cellToConnectionIdx_[cell_index];
// only deal with connections given by the aquifer
if( idx < 0)
continue;
auto isIt = gridView.ibegin(elem);
const auto& isEndIt = gridView.iend(elem);
for (; isIt != isEndIt; ++ isIt) {
// store intersection, this might be costly
const auto& intersection = *isIt;
// only deal with grid boundaries
if (!intersection.boundary())
continue;
int insideFaceIdx = intersection.indexInInside();
switch (insideFaceIdx) {
case 0:
faceDirection = Opm::FaceDir::XMinus;
break;
case 1:
faceDirection = Opm::FaceDir::XPlus;
break;
case 2:
faceDirection = Opm::FaceDir::YMinus;
break;
case 3:
faceDirection = Opm::FaceDir::YPlus;
break;
case 4:
faceDirection = Opm::FaceDir::ZMinus;
break;
case 5:
faceDirection = Opm::FaceDir::ZPlus;
break;
default:
OPM_THROW(std::logic_error,
"Internal error in initialization of aquifer.");
}
if (faceDirection == this->connections_[idx].face_dir) {
this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
break;
}
}
denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
}
const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
comm.sum(&denom_face_areas, 1);
const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
for (size_t idx = 0; idx < this->size(); ++idx) {
this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
? // Prevent no connection NaNs due to division by zero
0.
: (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx)) / denom_face_areas;
}
}
virtual void assignRestartData(const data::AquiferData& xaq) = 0;
@@ -278,9 +371,12 @@ protected:
}
// We take the average of the calculated equilibrium pressures.
const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.);
const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha;
return aquifer_pres_avg;
const auto& comm = ebos_simulator_.vanguard().grid().comm();
Scalar vals[2];
vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.);
vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.);
comm.sum(vals, 2);
return vals[1] / vals[0];
}
// This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer

4
opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
View File

@@ -168,10 +168,6 @@ BlackoilAquiferModel<TypeTag>::init()
return;
}
const auto& comm = this->simulator_.vanguard().gridView().comm();
if (comm.size() > 1)
throw std::runtime_error("Aquifers currently do not work in parallel.");
// Get all the carter tracy aquifer properties data and put it in aquifers vector
const auto& connections = aquifer.connections();
for (const auto& aq : aquifer.ct()) {