From 5586ce63a0f183fc366006dfae1f4fe3ecd6e290 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 9 Mar 2023 18:25:36 +0100 Subject: [PATCH 1/2] Calculate Local Face Area Fraction for Aquifer's Connections This is in preparation of reinitialising the total produced volume from the constant flux aquifer from the restart file in the case that the aquifer's connections are split across multiple MPI ranks. --- .../aquifers/AquiferConstantFlux.hpp | 27 ++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/opm/simulators/aquifers/AquiferConstantFlux.hpp b/opm/simulators/aquifers/AquiferConstantFlux.hpp index ed768a6d0..e3a266393 100644 --- a/opm/simulators/aquifers/AquiferConstantFlux.hpp +++ b/opm/simulators/aquifers/AquiferConstantFlux.hpp @@ -58,7 +58,8 @@ public: , aquifer_data_ (aquifer) , connection_flux_ (connections_.size(), Eval{0}) { - this->initializeConnections(); + const auto connected_face_area = this->initializeConnections(); + this->area_fraction_ = this->computeFaceAreaFraction(connected_face_area); } static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator) @@ -149,9 +150,15 @@ private: std::vector cellToConnectionIdx_{}; double flux_rate_{}; double cumulative_flux_{}; + double area_fraction_{1.0}; + + double initializeConnections() + { + auto connected_face_area = 0.0; + + this->cellToConnectionIdx_ + .resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1); - void initializeConnections() { - this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1); for (std::size_t idx = 0; idx < this->connections_.size(); ++idx) { const auto global_index = this->connections_[idx].global_index; const int cell_index = this->ebos_simulator_.vanguard() @@ -162,11 +169,25 @@ private: } this->cellToConnectionIdx_[cell_index] = idx; + + connected_face_area += this->connections_[idx].effective_facearea; } // TODO: At the moment, we are using the effective_facearea from the // parser. Should we update the facearea here if the grid changed // during the preprocessing? + + return connected_face_area; + } + + double computeFaceAreaFraction(const double connected_face_area) const + { + const auto tot_face_area = this->ebos_simulator_.vanguard() + .grid().comm().sum(connected_face_area); + + return (tot_face_area > 0.0) + ? connected_face_area / tot_face_area + : 0.0; } // TODO: this is a function from AquiferAnalytical From 14a63a4636c3e420654bce815cbfc8dc9b04ff1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 31 Mar 2023 11:47:58 +0200 Subject: [PATCH 2/2] Synchronise Face-Area Fractions Between All Processes We need a global view of face-area fractions if aquifer connections happen to be shared between processes. Add a new helper function, BlackoilAquiferModel::computeConnectionAreaFraction() that performs a collective operation to compute the total face areas and then defers to the local aquifer objects to compute their face area fractions. While here, also split the initialisation of analytic aquifers into two parts, one for the face area and connection mappings, and one for the connection depths. Run the former as part of the object constructor and the latter as part of 'initQuantities()'. This ensures that we can computeConnectionAreaFraction() for all analytic aquifers before assigning solution quantities from the restart file. --- opm/simulators/aquifers/AquiferAnalytical.hpp | 163 ++++++++++-------- .../aquifers/AquiferConstantFlux.hpp | 19 +- opm/simulators/aquifers/AquiferInterface.hpp | 3 + opm/simulators/aquifers/AquiferNumerical.hpp | 8 + .../aquifers/BlackoilAquiferModel.hpp | 2 + .../aquifers/BlackoilAquiferModel_impl.hpp | 53 +++++- 6 files changed, 166 insertions(+), 82 deletions(-) diff --git a/opm/simulators/aquifers/AquiferAnalytical.hpp b/opm/simulators/aquifers/AquiferAnalytical.hpp index 386578ed5..9c53c863d 100644 --- a/opm/simulators/aquifers/AquiferAnalytical.hpp +++ b/opm/simulators/aquifers/AquiferAnalytical.hpp @@ -87,17 +87,45 @@ public: BlackoilIndices::numPhases>; // Constructor - AquiferAnalytical(int aqID, - const std::vector& connections, - const Simulator& ebosSimulator) + AquiferAnalytical(const int aqID, + const std::vector& connections, + const Simulator& ebosSimulator) : AquiferInterface(aqID, ebosSimulator) , connections_(connections) { + this->initializeConnectionMappings(); } // Destructor virtual ~AquiferAnalytical() + {} + + void computeFaceAreaFraction(const std::vector& total_face_area) override { + assert (total_face_area.size() >= static_cast::size_type>(this->aquiferID())); + + const auto tfa = total_face_area[this->aquiferID() - 1]; + const auto eps_sqrt = std::sqrt(std::numeric_limits::epsilon()); + + if (tfa < eps_sqrt) { + this->alphai_.assign(this->size(), Scalar{0}); + } + else { + std::transform(this->faceArea_connected_.begin(), + this->faceArea_connected_.end(), + this->alphai_.begin(), + [tfa](const Scalar area) + { + return area / tfa; + }); + } + + this->area_fraction_ = this->totalFaceArea() / tfa; + } + + double totalFaceArea() const override + { + return this->total_face_area_; } void initFromRestart(const data::Aquifers& aquiferSoln) override @@ -108,8 +136,9 @@ public: this->assignRestartData(xaqPos->second); - this->W_flux_ = xaqPos->second.volume; + this->W_flux_ = xaqPos->second.volume * this->area_fraction_; this->pa0_ = xaqPos->second.initPressure; + this->solution_set_from_restart_ = true; } @@ -226,19 +255,17 @@ protected: void initQuantities() { // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 - if (!this->solution_set_from_restart_) { + if (! this->solution_set_from_restart_) { W_flux_ = Scalar{0}; } - // We next get our connections to the aquifer and initialize these quantities using the initialize_connections - // function - initializeConnections(); - calculateAquiferCondition(); - calculateAquiferConstants(); + this->initializeConnectionDepths(); + this->calculateAquiferCondition(); + this->calculateAquiferConstants(); - pressure_previous_.resize(this->connections_.size(), Scalar{0}); - pressure_current_.resize(this->connections_.size(), Scalar{0}); - Qai_.resize(this->connections_.size(), Scalar{0}); + this->pressure_previous_.resize(this->size(), Scalar{0}); + this->pressure_current_.resize(this->size(), Scalar{0}); + this->Qai_.resize(this->size(), Scalar{0}); } void updateCellPressure(std::vector& pressure_water, @@ -257,54 +284,54 @@ protected: pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value(); } - void initializeConnections() + void initializeConnectionMappings() { - this->cell_depth_.resize(this->size(), this->aquiferDepth()); this->alphai_.resize(this->size(), 1.0); this->faceArea_connected_.resize(this->size(), Scalar{0}); - // Translate the C face tag into the enum used by opm-parser's TransMult class - FaceDir::DirEnum faceDirection; - - bool has_active_connection_on_proc = false; - - // denom_face_areas is the sum of the areas connected to an aquifer - Scalar denom_face_areas{0}; + // total_face_area_ is the sum of the areas connected to an aquifer + this->total_face_area_ = Scalar{0}; this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1); const auto& gridView = this->ebos_simulator_.vanguard().gridView(); for (std::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(); - 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) + if (cell_index < 0) { continue; + } - has_active_connection_on_proc = true; + auto elemIt = gridView.template begin(); + std::advance(elemIt, cell_index); + + // The global_index is not part of this grid + if (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()); - for (const auto& elem : elements(gridView)) { - unsigned cell_index = elemMapper.index(elem); - int idx = this->cellToConnectionIdx_[cell_index]; - // only deal with connections given by the aquifer - if( idx < 0) + // Translate the C face tag into the enum used by opm-parser's TransMult class + FaceDir::DirEnum faceDirection; + + // Get areas for all connections + const auto& elemMapper = this->ebos_simulator_.model().dofMapper(); + for (const auto& elem : elements(gridView)) { + const unsigned cell_index = elemMapper.index(elem); + const int idx = this->cellToConnectionIdx_[cell_index]; + + // Only deal with connections given by the aquifer + if (idx < 0) { continue; + } for (const auto& intersection : intersections(gridView, elem)) { - // only deal with grid boundaries - if (!intersection.boundary()) + // Only deal with grid boundaries + if (! intersection.boundary()) { continue; + } - int insideFaceIdx = intersection.indexInInside(); - switch (insideFaceIdx) { + switch (intersection.indexInInside()) { case 0: faceDirection = FaceDir::XMinus; break; @@ -325,51 +352,41 @@ protected: break; default: OPM_THROW(std::logic_error, - "Internal error in initialization of aquifer."); + "Internal error in initialization of aquifer."); } - if (faceDirection == this->connections_[idx].face_dir) { this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff; break; } } - denom_face_areas += 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::epsilon()); - for (std::size_t idx = 0; idx < this->size(); ++idx) { - // Protect against division by zero NaNs. - this->alphai_.at(idx) = (denom_face_areas < eps_sqrt) - ? Scalar{0} - : this->faceArea_connected_.at(idx) / denom_face_areas; - } - - if (this->solution_set_from_restart_) { - this->rescaleProducedVolume(has_active_connection_on_proc); + this->total_face_area_ += this->faceArea_connected_.at(idx); } } - void rescaleProducedVolume(const bool has_active_connection_on_proc) + void initializeConnectionDepths() { - // Needed in parallel restart to approximate influence of aquifer - // being "owned" by a subset of the parallel processes. If the - // aquifer is fully owned by a single process--i.e., if all cells - // connecting to the aquifer are on a single process--then this_area - // is tot_area on that process and zero elsewhere. + this->cell_depth_.resize(this->size(), this->aquiferDepth()); - const auto this_area = has_active_connection_on_proc - ? std::accumulate(this->alphai_.begin(), - this->alphai_.end(), - Scalar{0}) - : Scalar{0}; + const auto& gridView = this->ebos_simulator_.vanguard().gridView(); + for (std::size_t idx = 0; idx < this->size(); ++idx) { + const int cell_index = this->ebos_simulator_.vanguard() + .compressedIndex(this->connections_[idx].global_index); + if (cell_index < 0) { + continue; + } - const auto tot_area = this->ebos_simulator_.vanguard() - .grid().comm().sum(this_area); + auto elemIt = gridView.template begin(); + std::advance(elemIt, cell_index); - this->W_flux_ *= this_area / tot_area; + // The global_index is not part of this grid + if (elemIt->partitionType() != Dune::InteriorEntity) { + continue; + } + + this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index); + } } // This function is for calculating the aquifer properties from equilibrium state with the reservoir @@ -433,9 +450,13 @@ protected: std::optional Ta0_{}; // initial aquifer temperature Scalar rhow_{}; + Scalar total_face_area_{}; + Scalar area_fraction_{Scalar{1}}; + Eval W_flux_; bool solution_set_from_restart_ {false}; + bool has_active_connection_on_proc_{false}; }; } // namespace Opm diff --git a/opm/simulators/aquifers/AquiferConstantFlux.hpp b/opm/simulators/aquifers/AquiferConstantFlux.hpp index e3a266393..e9c9aa020 100644 --- a/opm/simulators/aquifers/AquiferConstantFlux.hpp +++ b/opm/simulators/aquifers/AquiferConstantFlux.hpp @@ -32,8 +32,10 @@ #include #include +#include #include #include +#include namespace Opm { @@ -58,8 +60,7 @@ public: , aquifer_data_ (aquifer) , connection_flux_ (connections_.size(), Eval{0}) { - const auto connected_face_area = this->initializeConnections(); - this->area_fraction_ = this->computeFaceAreaFraction(connected_face_area); + this->total_face_area_ = this->initializeConnections(); } static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator) @@ -72,6 +73,19 @@ public: virtual ~AquiferConstantFlux() = default; + void computeFaceAreaFraction(const std::vector& total_face_area) override + { + assert (total_face_area.size() >= static_cast::size_type>(this->aquiferID())); + + this->area_fraction_ = this->totalFaceArea() + / total_face_area[this->aquiferID() - 1]; + } + + double totalFaceArea() const override + { + return this->total_face_area_; + } + void updateAquifer(const SingleAquiferFlux& aquifer) { aquifer_data_ = aquifer; @@ -150,6 +164,7 @@ private: std::vector cellToConnectionIdx_{}; double flux_rate_{}; double cumulative_flux_{}; + double total_face_area_{0.0}; double area_fraction_{1.0}; double initializeConnections() diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 69e000b01..3bef5e80a 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -58,6 +58,9 @@ public: virtual data::AquiferData aquiferData() const = 0; + virtual void computeFaceAreaFraction(const std::vector& total_face_area) = 0; + virtual double totalFaceArea() const = 0; + template void addToSource(RateVector& rates, const Context& context, diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp index 8f5fb7600..5b2065ca3 100644 --- a/opm/simulators/aquifers/AquiferNumerical.hpp +++ b/opm/simulators/aquifers/AquiferNumerical.hpp @@ -154,6 +154,14 @@ public: this->cumulative_flux_ = 0.; } + void computeFaceAreaFraction(const std::vector& /*total_face_area*/) override + {} + + double totalFaceArea() const override + { + return 1.0; + } + template void serializeOp(Serializer& serializer) { diff --git a/opm/simulators/aquifers/BlackoilAquiferModel.hpp b/opm/simulators/aquifers/BlackoilAquiferModel.hpp index 879af9a04..0388045ce 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel.hpp @@ -140,6 +140,8 @@ private: createAnalyticAquiferPointer(const AquiferData& aqData, const int aquiferID, std::string_view aqType) const; + + void computeConnectionAreaFraction() const; }; diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp index 1014e5095..41f9badfa 100644 --- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp +++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp @@ -45,16 +45,22 @@ template void BlackoilAquiferModel::initialSolutionApplied() { - for (auto& aquifer : aquifers) + this->computeConnectionAreaFraction(); + + for (auto& aquifer : this->aquifers) { aquifer->initialSolutionApplied(); + } } template void BlackoilAquiferModel::initFromRestart(const data::Aquifers& aquiferSoln) { - for (auto& aquifer : this->aquifers) + this->computeConnectionAreaFraction(); + + for (auto& aquifer : this->aquifers) { aquifer->initFromRestart(aquiferSoln); + } } template @@ -67,6 +73,8 @@ BlackoilAquiferModel::beginEpisode() // SCHEDULE setup in this section it is the beginning of a report step this->createDynamicAquifers(this->simulator_.episodeIndex()); + + this->computeConnectionAreaFraction(); } template @@ -78,8 +86,9 @@ template void BlackoilAquiferModel::beginTimeStep() { - for (auto& aquifer : aquifers) + for (auto& aquifer : this->aquifers) { aquifer->beginTimeStep(); + } } template @@ -90,8 +99,9 @@ BlackoilAquiferModel::addToSource(RateVector& rates, unsigned spaceIdx, unsigned timeIdx) const { - for (auto& aquifer : aquifers) + for (auto& aquifer : this->aquifers) { aquifer->addToSource(rates, context, spaceIdx, timeIdx); + } } template @@ -100,8 +110,9 @@ BlackoilAquiferModel::addToSource(RateVector& rates, unsigned globalSpaceIdx, unsigned timeIdx) const { - for (auto& aquifer : aquifers) + for (auto& aquifer : this->aquifers) { aquifer->addToSource(rates, globalSpaceIdx, timeIdx); + } } template @@ -113,9 +124,10 @@ template void BlackoilAquiferModel::endTimeStep() { - for (auto& aquifer : aquifers) { + using NumAq = AquiferNumerical; + + for (auto& aquifer : this->aquifers) { aquifer->endTimeStep(); - using NumAq = AquiferNumerical; NumAq* num = dynamic_cast(aquifer.get()); if (num) this->simulator_.vanguard().grid().comm().barrier(); @@ -159,8 +171,9 @@ template data::Aquifers BlackoilAquiferModel::aquiferData() const { data::Aquifers data; - for (const auto& aqu : this->aquifers) + for (const auto& aqu : this->aquifers) { data.insert_or_assign(aqu->aquiferID(), aqu->aquiferData()); + } return data; } @@ -170,7 +183,7 @@ template void BlackoilAquiferModel:: serializeOp(Serializer& serializer) { - for (auto& aiPtr : aquifers) { + for (auto& aiPtr : this->aquifers) { auto* ct = dynamic_cast*>(aiPtr.get()); auto* fetp = dynamic_cast*>(aiPtr.get()); auto* num = dynamic_cast*>(aiPtr.get()); @@ -304,4 +317,26 @@ void BlackoilAquiferModel::createDynamicAquifers(const int episode_inde } } +template +void BlackoilAquiferModel::computeConnectionAreaFraction() const +{ + auto maxAquID = + std::accumulate(this->aquifers.begin(), this->aquifers.end(), 0, + [](const int aquID, const auto& aquifer) + { return std::max(aquID, aquifer->aquiferID()); }); + + maxAquID = this->simulator_.vanguard().grid().comm().max(maxAquID); + + auto totalConnArea = std::vector(maxAquID, 0.0); + for (const auto& aquifer : this->aquifers) { + totalConnArea[aquifer->aquiferID() - 1] += aquifer->totalFaceArea(); + } + + this->simulator_.vanguard().grid().comm().sum(totalConnArea.data(), maxAquID); + + for (auto& aquifer : this->aquifers) { + aquifer->computeFaceAreaFraction(totalConnArea); + } +} + } // namespace Opm