diff --git a/ebos/collecttoiorank.cc b/ebos/collecttoiorank.cc index 5ab9c0263..8a8aefaf0 100644 --- a/ebos/collecttoiorank.cc +++ b/ebos/collecttoiorank.cc @@ -37,6 +37,7 @@ #include #endif +#include #include #include #include @@ -627,12 +628,20 @@ public: this->unpackCommon(data, aq); - if (data.aquFet != nullptr) { - this->unpackAquFet(std::move(data.aquFet), aq); + if (auto const* aquFet = data.typeData.get(); + aquFet != nullptr) + { + this->unpackAquFet(*aquFet, aq); } - - if (data.aquCT != nullptr) { - this->unpackCarterTracy(std::move(data.aquCT), aq); + else if (auto const* aquCT = data.typeData.get(); + aquCT != nullptr) + { + this->unpackCarterTracy(*aquCT, aq); + } + else if (auto const* aquNum = data.typeData.get(); + aquNum != nullptr) + { + this->unpackNumericAquifer(*aquNum, aq); } } } @@ -648,14 +657,15 @@ private: aq.volume += data.volume; } - void unpackAquFet(std::shared_ptr aquFet, data::AquiferData& aq) + void unpackAquFet(const data::FetkovichData& aquFet, data::AquiferData& aq) { - if (aq.aquFet == nullptr) { - aq.aquFet = std::move(aquFet); + if (! aq.typeData.is()) { + auto* tData = aq.typeData.create(); + *tData = aquFet; } else { - const auto& src = *aquFet; - auto& dst = *aq.aquFet; + const auto& src = aquFet; + auto& dst = *aq.typeData.getMutable(); dst.initVolume = std::max(dst.initVolume , src.initVolume); dst.prodIndex = std::max(dst.prodIndex , src.prodIndex); @@ -663,19 +673,46 @@ private: } } - void unpackCarterTracy(std::shared_ptr aquCT, data::AquiferData& aq) + void unpackCarterTracy(const data::CarterTracyData& aquCT, data::AquiferData& aq) { - if (aq.aquCT == nullptr) { - aq.aquCT = std::move(aquCT); + if (! aq.typeData.is()) { + auto* tData = aq.typeData.create(); + *tData = aquCT; } else { - const auto& src = *aquCT; - auto& dst = *aq.aquCT; + const auto& src = aquCT; + auto& dst = *aq.typeData.getMutable(); + + dst.timeConstant = std::max(dst.timeConstant , src.timeConstant); + dst.influxConstant = std::max(dst.influxConstant, src.influxConstant); + dst.waterDensity = std::max(dst.waterDensity , src.waterDensity); + dst.waterViscosity = std::max(dst.waterViscosity, src.waterViscosity); dst.dimensionless_time = std::max(dst.dimensionless_time , src.dimensionless_time); dst.dimensionless_pressure = std::max(dst.dimensionless_pressure, src.dimensionless_pressure); } } + + void unpackNumericAquifer(const data::NumericAquiferData& aquNum, data::AquiferData& aq) + { + if (! aq.typeData.is()) { + auto* tData = aq.typeData.create(); + *tData = aquNum; + } + else { + const auto& src = aquNum; + auto& dst = *aq.typeData.getMutable(); + + std::transform(src.initPressure.begin(), + src.initPressure.end(), + dst.initPressure.begin(), + dst.initPressure.begin(), + [](const double p0_1, const double p0_2) + { + return std::max(p0_1, p0_2); + }); + } + } }; diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index e73b15503..2f53c0eb1 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -81,15 +81,14 @@ public: } data.volume = this->W_flux_.value(); data.initPressure = this->pa0_; - data.type = data::AquiferType::CarterTracy; - data.aquCT = std::make_shared(); - data.aquCT->timeConstant = this->aquct_data_.timeConstant(); - data.aquCT->influxConstant = this->aquct_data_.influxConstant(); - data.aquCT->waterDensity = this->aquct_data_.waterDensity(); - data.aquCT->waterViscosity = this->aquct_data_.waterViscosity(); - data.aquCT->dimensionless_time = this->dimensionless_time_; - data.aquCT->dimensionless_pressure = this->dimensionless_pressure_; + auto* aquCT = data.typeData.template create(); + aquCT->timeConstant = this->aquct_data_.timeConstant(); + aquCT->influxConstant = this->aquct_data_.influxConstant(); + aquCT->waterDensity = this->aquct_data_.waterDensity(); + aquCT->waterViscosity = this->aquct_data_.waterViscosity(); + aquCT->dimensionless_time = this->dimensionless_time_; + aquCT->dimensionless_pressure = this->dimensionless_pressure_; return data; } diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index 631f8e0ae..ac53e8ab6 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -27,6 +27,7 @@ along with OPM. If not, see . #include #include +#include namespace Opm { @@ -82,12 +83,11 @@ public: }); data.volume = this->W_flux_.value(); data.initPressure = this->pa0_; - data.type = data::AquiferType::Fetkovich; - data.aquFet = std::make_shared(); - data.aquFet->initVolume = this->aqufetp_data_.initial_watvolume; - data.aquFet->prodIndex = this->aqufetp_data_.prod_index; - data.aquFet->timeConstant = this->aqufetp_data_.timeConstant(); + auto* aquFet = data.typeData.template create(); + aquFet->initVolume = this->aqufetp_data_.initial_watvolume; + aquFet->prodIndex = this->aqufetp_data_.prod_index; + aquFet->timeConstant = this->aqufetp_data_.timeConstant(); return data; } @@ -99,7 +99,7 @@ protected: void assignRestartData(const data::AquiferData& xaq) override { - if (xaq.type != data::AquiferType::Fetkovich) { + if (! xaq.typeData.is()) { throw std::invalid_argument { "Analytic aquifer data for unexpected aquifer " "type passed to Fetkovich aquifer" diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp index 5f5dbf9fd..06af4689e 100644 --- a/opm/simulators/aquifers/AquiferNumerical.hpp +++ b/opm/simulators/aquifers/AquiferNumerical.hpp @@ -24,6 +24,8 @@ #include #include +#include + namespace Opm { template @@ -52,22 +54,22 @@ public: const std::unordered_map& cartesian_to_compressed, const Simulator& ebos_simulator, const int* global_cell) - : id_(aquifer.id()) - , ebos_simulator_(ebos_simulator) - , flux_rate_(0.) - , cumulative_flux_(0.) - , global_cell_(global_cell) + : id_(aquifer.id()) + , ebos_simulator_(ebos_simulator) + , flux_rate_(0.) + , cumulative_flux_(0.) + , global_cell_(global_cell) + , init_pressure_(aquifer.numCells(), 0.0) { this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1); for (size_t idx = 0; idx < aquifer.numCells(); ++idx) { - const auto& cell = *(aquifer.getCellPrt(idx)); - const int global_idx = cell.global_index; - const auto search = cartesian_to_compressed.find(global_idx); + const auto* cell = aquifer.getCellPrt(idx); + // Due to parallelisation, the cell might not exist in the current process + auto search = cartesian_to_compressed.find(cell->global_index); if (search != cartesian_to_compressed.end()) { - const int cell_idx = cartesian_to_compressed.at(global_idx); - this->cell_to_aquifer_cell_idx_[cell_idx] = idx; + this->cell_to_aquifer_cell_idx_[search->second] = idx; } } } @@ -87,19 +89,20 @@ public: data::AquiferData aquiferData() const { data::AquiferData data; - data.aquiferID = this->id_; - data.initPressure = this->init_pressure_; + data.aquiferID = this->aquiferID(); data.pressure = this->pressure_; data.fluxRate = this->flux_rate_; data.volume = this->cumulative_flux_; - data.type = data::AquiferType::Numerical; + + auto* aquNum = data.typeData.template create(); + aquNum->initPressure = this->init_pressure_; + return data; } void initialSolutionApplied() { - this->init_pressure_ = this->calculateAquiferPressure(); - this->pressure_ = this->init_pressure_; + this->pressure_ = this->calculateAquiferPressure(this->init_pressure_); this->flux_rate_ = 0.; this->cumulative_flux_ = 0.; } @@ -115,13 +118,19 @@ private: double flux_rate_; // aquifer influx rate double cumulative_flux_; // cumulative aquifer influx const int* global_cell_; // mapping to global index - double init_pressure_; + std::vector init_pressure_{}; double pressure_; // aquifer pressure // TODO: maybe unordered_map can also do the work to save memory? std::vector cell_to_aquifer_cell_idx_; double calculateAquiferPressure() const + { + auto capture = std::vector(this->init_pressure_.size(), 0.0); + return this->calculateAquiferPressure(capture); + } + + double calculateAquiferPressure(std::vector& cell_pressure) const { double sum_pressure_watervolume = 0.; double sum_watervolume = 0.; @@ -158,11 +167,17 @@ private: const double water_volume = volume * porosity * water_saturation; sum_pressure_watervolume += water_volume * water_pressure_reservoir; sum_watervolume += water_volume; + + cell_pressure[idx] = water_pressure_reservoir; } const auto& comm = this->ebos_simulator_.vanguard().grid().comm(); comm.sum(&sum_pressure_watervolume, 1); comm.sum(&sum_watervolume, 1); + + // Ensure all processes have same notion of the aquifer cells' pressure values. + comm.sum(cell_pressure.data(), cell_pressure.size()); + return sum_pressure_watervolume / sum_watervolume; } diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index bf09ec5c3..80b49b6e4 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -83,14 +83,12 @@ BlackoilAquiferModel::initFromRestart(const data::Aquifers& aquiferSoln template void BlackoilAquiferModel::beginEpisode() -{ -} +{} template void BlackoilAquiferModel::beginIteration() -{ -} +{} template void @@ -131,8 +129,7 @@ BlackoilAquiferModel::addToSource(RateVector& rates, template void BlackoilAquiferModel::endIteration() -{ -} +{} template void @@ -155,11 +152,11 @@ BlackoilAquiferModel::endTimeStep() } } } + template void BlackoilAquiferModel::endEpisode() -{ -} +{} template template @@ -215,12 +212,14 @@ BlackoilAquiferModel::init() } } } + template bool BlackoilAquiferModel::aquiferCarterTracyActive() const { return !aquifers_CarterTracy.empty(); } + template bool BlackoilAquiferModel::aquiferFetkovichActive() const @@ -236,7 +235,8 @@ BlackoilAquiferModel::aquiferNumericalActive() const } template -data::Aquifers BlackoilAquiferModel::aquiferData() const { +data::Aquifers BlackoilAquiferModel::aquiferData() const +{ data::Aquifers data; if (this->aquiferCarterTracyActive()) { for (const auto& aqu : this->aquifers_CarterTracy) { diff --git a/opm/simulators/utils/ParallelRestart.cpp b/opm/simulators/utils/ParallelRestart.cpp index 4bdef0504..6a5f91348 100644 --- a/opm/simulators/utils/ParallelRestart.cpp +++ b/opm/simulators/utils/ParallelRestart.cpp @@ -226,6 +226,11 @@ HANDLE_AS_POD(data::NodeData) HANDLE_AS_POD(data::Rates) HANDLE_AS_POD(data::Segment) +std::size_t packSize(const data::NumericAquiferData& data, Dune::MPIHelper::MPICommunicator comm) +{ + return packSize(data.initPressure, comm); +} + std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunicator comm) { const auto type = 0ull; @@ -238,11 +243,20 @@ std::size_t packSize(const data::AquiferData& data, Dune::MPIHelper::MPICommunic + packSize(data.datumDepth, comm) + packSize(type, comm); - if (data.aquFet != nullptr) { - return base + packSize(*data.aquFet, comm); + if (auto const* aquFet = data.typeData.get(); + aquFet != nullptr) + { + return base + packSize(*aquFet, comm); } - else if (data.aquCT != nullptr) { - return base + packSize(*data.aquCT, comm); + else if (auto const* aquCT = data.typeData.get(); + aquCT != nullptr) + { + return base + packSize(*aquCT, comm); + } + else if (auto const* aquNum = data.typeData.get(); + aquNum != nullptr) + { + return base + packSize(*aquNum, comm); } return base; @@ -513,12 +527,19 @@ void pack(const std::unordered_map& data, std::vector& buffer } } +void pack(const data::NumericAquiferData& data, std::vector& buffer, int& position, + Dune::MPIHelper::MPICommunicator comm) +{ + pack(data.initPressure, buffer, position, comm); +} + void pack(const data::AquiferData& data, std::vector& buffer, int& position, Dune::MPIHelper::MPICommunicator comm) { const auto type = - (data.aquFet != nullptr)*(1ull << 0) - + (data.aquCT != nullptr)*(1ull << 1); + (data.typeData.is() * (1ull << 0)) + + (data.typeData.is() * (1ull << 1)) + + (data.typeData.is() * (1ull << 2)); pack(data.aquiferID, buffer, position, comm); pack(data.pressure, buffer, position, comm); @@ -528,11 +549,20 @@ void pack(const data::AquiferData& data, std::vector& buffer, int& positio pack(data.datumDepth, buffer, position, comm); pack(type, buffer, position, comm); - if (data.aquFet != nullptr) { - pack(*data.aquFet, buffer, position, comm); + if (auto const* aquFet = data.typeData.get(); + aquFet != nullptr) + { + pack(*aquFet, buffer, position, comm); } - else if (data.aquCT != nullptr) { - pack(*data.aquCT, buffer, position, comm); + else if (auto const* aquCT = data.typeData.get(); + aquCT != nullptr) + { + pack(*aquCT, buffer, position, comm); + } + else if (auto const* aquNum = data.typeData.get(); + aquNum != nullptr) + { + pack(*aquNum, buffer, position, comm); } } @@ -865,6 +895,12 @@ void unpack(data::Well& data, std::vector& buffer, int& position, unpack(data.guide_rates, buffer, position, comm); } +void unpack(data::NumericAquiferData& data, std::vector& buffer, int& position, + Dune::MPIHelper::MPICommunicator comm) +{ + unpack(data.initPressure, buffer, position, comm); +} + void unpack(data::AquiferData& data, std::vector& buffer, int& position, Dune::MPIHelper::MPICommunicator comm) { @@ -878,15 +914,17 @@ void unpack(data::AquiferData& data, std::vector& buffer, int& position, unpack(data.datumDepth, buffer, position, comm); unpack(type, buffer, position, comm); - if (type == 1ull) { - data.type = data::AquiferType::Fetkovich; - data.aquFet = std::make_shared(); - unpack(*data.aquFet, buffer, position, comm); + if (type == (1ull << 0)) { + auto* aquFet = data.typeData.create(); + unpack(*aquFet, buffer, position, comm); } - else if (type == 2ull) { - data.type = data::AquiferType::CarterTracy; - data.aquCT = std::make_shared(); - unpack(*data.aquCT, buffer, position, comm); + else if (type == (1ull << 1)) { + auto* aquCT = data.typeData.create(); + unpack(*aquCT, buffer, position, comm); + } + else if (type == (1ull << 2)) { + auto* aquNum = data.typeData.create(); + unpack(*aquNum, buffer, position, comm); } } diff --git a/opm/simulators/utils/ParallelRestart.hpp b/opm/simulators/utils/ParallelRestart.hpp index f4b875445..93835fbf0 100644 --- a/opm/simulators/utils/ParallelRestart.hpp +++ b/opm/simulators/utils/ParallelRestart.hpp @@ -61,6 +61,7 @@ struct GroupData; struct GroupGuideRates; class GuideRateValue; struct NodeData; +struct NumericAquiferData; class Rates; struct Segment; class Solution; @@ -349,6 +350,7 @@ ADD_PACK_PROTOTYPES(data::GroupGuideRates) ADD_PACK_PROTOTYPES(data::GroupData) ADD_PACK_PROTOTYPES(data::NodeData) ADD_PACK_PROTOTYPES(data::GroupAndNetworkValues) +ADD_PACK_PROTOTYPES(data::NumericAquiferData) ADD_PACK_PROTOTYPES(data::Well) ADD_PACK_PROTOTYPES(data::Wells) ADD_PACK_PROTOTYPES(RestartKey) diff --git a/tests/test_ParallelRestart.cpp b/tests/test_ParallelRestart.cpp index 93c15caee..71e0126bd 100644 --- a/tests/test_ParallelRestart.cpp +++ b/tests/test_ParallelRestart.cpp @@ -265,14 +265,14 @@ Opm::data::NodeData getNodeData() Opm::data::AquiferData getFetkovichAquifer(const int aquiferID = 1) { auto aquifer = Opm::data::AquiferData { - aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5, Opm::data::AquiferType::Fetkovich + aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5 }; - aquifer.aquFet = std::make_shared(); + auto* aquFet = aquifer.typeData.create(); - aquifer.aquFet->initVolume = 1.23; - aquifer.aquFet->prodIndex = 45.67; - aquifer.aquFet->timeConstant = 890.123; + aquFet->initVolume = 1.23; + aquFet->prodIndex = 45.67; + aquFet->timeConstant = 890.123; return aquifer; } @@ -280,17 +280,33 @@ Opm::data::AquiferData getFetkovichAquifer(const int aquiferID = 1) Opm::data::AquiferData getCarterTracyAquifer(const int aquiferID = 5) { auto aquifer = Opm::data::AquiferData { - aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5, Opm::data::AquiferType::CarterTracy + aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5 }; - aquifer.aquCT = std::make_shared(); + auto* aquCT = aquifer.typeData.create(); - aquifer.aquCT->timeConstant = 987.65; - aquifer.aquCT->influxConstant = 43.21; - aquifer.aquCT->waterDensity = 1014.5; - aquifer.aquCT->waterViscosity = 0.00318; - aquifer.aquCT->dimensionless_time = 42.0; - aquifer.aquCT->dimensionless_pressure = 2.34; + aquCT->timeConstant = 987.65; + aquCT->influxConstant = 43.21; + aquCT->waterDensity = 1014.5; + aquCT->waterViscosity = 0.00318; + aquCT->dimensionless_time = 42.0; + aquCT->dimensionless_pressure = 2.34; + + return aquifer; +} + +Opm::data::AquiferData getNumericalAquifer(const int aquiferID = 2) +{ + auto aquifer = Opm::data::AquiferData { + aquiferID, 123.456, 56.78, 9.0e10, 290.0, 2515.5 + }; + + auto* aquNum = aquifer.typeData.create(); + + aquNum->initPressure.push_back(1.234); + aquNum->initPressure.push_back(2.345); + aquNum->initPressure.push_back(3.4); + aquNum->initPressure.push_back(9.876); return aquifer; } @@ -363,12 +379,21 @@ BOOST_AUTO_TEST_CASE(dataCarterTracyData) DO_CHECKS(data::CarterTracyData) } +BOOST_AUTO_TEST_CASE(dataNumericAquiferData) +{ + const auto val1 = getNumericalAquifer(); + const auto val2 = PackUnpack(val1); + + DO_CHECKS(data::NumericAquiferData) +} + BOOST_AUTO_TEST_CASE(dataAquifers) { const auto val1 = Opm::data::Aquifers { { 1, getFetkovichAquifer(1) }, { 4, getFetkovichAquifer(4) }, { 5, getCarterTracyAquifer(5) }, + { 9, getNumericalAquifer(9) }, }; const auto val2 = PackUnpack(val1);