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 ed768a6d0..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,7 +60,7 @@ public: , aquifer_data_ (aquifer) , connection_flux_ (connections_.size(), Eval{0}) { - this->initializeConnections(); + this->total_face_area_ = this->initializeConnections(); } static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator) @@ -71,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; @@ -149,9 +164,16 @@ private: std::vector cellToConnectionIdx_{}; double flux_rate_{}; double cumulative_flux_{}; + double total_face_area_{0.0}; + 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 +184,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 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