Merge pull request #3094 from bska/wellpi-shut-wells

Compute Well-Level PI For Shut Wells
2026-09-05 04:40:19 -05:00 · 2021-04-28 09:03:57 +02:00
parent 1de5c9539c a6c374a27c
commit 398b284c03
4 changed files with 208 additions and 112 deletions
@@ -410,26 +410,27 @@ namespace Opm {
        protected:
            Simulator& ebosSimulator_;

-            std::vector< Well > wells_ecl_;
-            std::vector< std::vector<PerforationData> > well_perf_data_;
-            std::vector< WellProdIndexCalculator > prod_index_calc_;
+            std::vector< Well > wells_ecl_{};
+            std::vector< std::vector<PerforationData> > well_perf_data_{};
+            std::vector< WellProdIndexCalculator > prod_index_calc_{};
+            std::vector<int> local_shut_wells_{};

            std::vector< ParallelWellInfo > parallel_well_info_;
            std::vector< ParallelWellInfo* > local_parallel_well_info_;

-            bool wells_active_;
+            bool wells_active_{false};

            // a vector of all the wells.
-            std::vector<WellInterfacePtr > well_container_;
+            std::vector<WellInterfacePtr > well_container_{};

-            // map from logically cartesian cell indices to compressed ones
-            // cells not in the interior are not mapped. This deactivates
-            // these for distributed wells and make the distribution non-overlapping.
-            std::vector<int> cartesian_to_compressed_;
+            // Map from logically cartesian cell indices to compressed ones.
+            // Cells not in the interior are not mapped. This deactivates
+            // these for distributed wells and makes the distribution non-overlapping.
+            std::vector<int> cartesian_to_compressed_{};

-            std::vector<bool> is_cell_perforated_;
+            std::vector<bool> is_cell_perforated_{};

-            std::function<bool(const Well&)> not_on_process_;
+            std::function<bool(const Well&)> not_on_process_{};

            void initializeWellProdIndCalculators();
            void initializeWellPerfData();
@@ -440,6 +441,8 @@ namespace Opm {
            // create the well container
            std::vector<WellInterfacePtr > createWellContainer(const int time_step);

+            void inferLocalShutWells();
+
            WellInterfacePtr
            createWellPointer(const int wellID,
                              const int time_step) const;
@@ -453,36 +456,36 @@ namespace Opm {


            const ModelParameters param_;
-            bool terminal_output_;
-            std::vector<int> pvt_region_idx_;
+            bool terminal_output_{false};
+            std::vector<int> pvt_region_idx_{};
            PhaseUsage phase_usage_;
-            size_t global_num_cells_;
+            size_t global_num_cells_{};
            // the number of the cells in the local grid
-            size_t local_num_cells_;
-            double gravity_;
-            std::vector<double> depth_;
-            bool initial_step_;
-            bool report_step_starts_;
+            size_t local_num_cells_{};
+            double gravity_{};
+            std::vector<double> depth_{};
+            bool initial_step_{};
+            bool report_step_starts_{};
            bool glift_debug = false;
-            bool alternative_well_rate_init_;
+            bool alternative_well_rate_init_{};

            std::optional<int> last_run_wellpi_{};

-            std::unique_ptr<RateConverterType> rateConverter_;
-            std::unique_ptr<VFPProperties> vfp_properties_;
+            std::unique_ptr<RateConverterType> rateConverter_{};
+            std::unique_ptr<VFPProperties> vfp_properties_{};

-            SimulatorReportSingle last_report_;
+            SimulatorReportSingle last_report_{};

-            WellTestState wellTestState_;
-            std::unique_ptr<GuideRate> guideRate_;
+            WellTestState wellTestState_{};
+            std::unique_ptr<GuideRate> guideRate_{};

            std::map<std::string, double> node_pressures_{}; // Storing network pressures for output.
            mutable std::unordered_set<std::string> closed_this_step_{};

            // used to better efficiency of calcuation
-            mutable BVector scaleAddRes_;
+            mutable BVector scaleAddRes_{};

-            std::vector< Scalar > B_avg_;
+            std::vector<Scalar> B_avg_{};

            const Grid& grid() const
            { return ebosSimulator_.vanguard().grid(); }
@@ -543,7 +546,10 @@ namespace Opm {

            void calculateEfficiencyFactors(const int reportStepIdx);

+            void calculateProductivityIndexValuesShutWells(const int reportStepIdx, DeferredLogger& deferred_logger);
            void calculateProductivityIndexValues(DeferredLogger& deferred_logger);
+            void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
+                                                  DeferredLogger& deferred_logger);

            // The number of components in the model.
            int numComponents() const;
@@ -1,23 +1,23 @@
 /*
- Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics.
- Copyright 2016 - 2018 Equinor ASA.
- Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
- Copyright 2016 - 2018 Norce AS
+  Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics.
+  Copyright 2016 - 2018 Equinor ASA.
+  Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
+  Copyright 2016 - 2018 Norce AS

- This file is part of the Open Porous Media project (OPM).
+  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 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.
+  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/>.
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */

 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
@@ -27,51 +27,59 @@

 #include <algorithm>
 #include <utility>
+
 #include <fmt/format.h>

 namespace Opm {
    template<typename TypeTag>
    BlackoilWellModel<TypeTag>::
    BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage)
-        : ebosSimulator_(ebosSimulator),
-          phase_usage_(phase_usage),
-          active_wgstate_(phase_usage),
-          last_valid_wgstate_(phase_usage),
-          nupcol_wgstate_(phase_usage)
+        : ebosSimulator_(ebosSimulator)
+        , terminal_output_((ebosSimulator.gridView().comm().rank() == 0) &&
+                           EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput))
+        , phase_usage_(phase_usage)
+        , active_wgstate_(phase_usage)
+        , last_valid_wgstate_(phase_usage)
+        , nupcol_wgstate_(phase_usage)
    {
-        terminal_output_ = false;
-        if (ebosSimulator.gridView().comm().rank() == 0)
-            terminal_output_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput);
-
        // Create the guide rate container.
-        guideRate_.reset(new GuideRate (ebosSimulator_.vanguard().schedule()));
+        this->guideRate_ =
+            std::make_unique<GuideRate>(ebosSimulator_.vanguard().schedule());

        local_num_cells_ = ebosSimulator_.gridView().size(0);
+
        // Number of cells the global grid view
        global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();

        // Set up cartesian mapping.
-        const auto& grid = ebosSimulator_.vanguard().grid();
-        const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
-        setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid),
-                                    cartDims[0]*cartDims[1]*cartDims[2]);
-        auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
-        parallel_well_info_.assign(parallel_wells.begin(), parallel_wells.end());
-        const auto& pwell_info = parallel_well_info_;
-        std::size_t numProcs = ebosSimulator.gridView().comm().size();
-        not_on_process_ = [&pwell_info, numProcs](const Well& well) {
-                if (well.getStatus() == Well::Status::SHUT)
-                    return true;
-                if (numProcs == 1u)
-                    return false;
-                std::pair<std::string, bool> value{well.name(), true}; // false indicate not active!
-                auto candidate = std::lower_bound(pwell_info.begin(),
-                                                  pwell_info.end(),
-                                                  value);
-                return candidate == pwell_info.end() || *candidate != value;
-            };
+        {
+            const auto& grid = this->ebosSimulator_.vanguard().grid();
+            const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
+            setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid),
+                                        cartDims[0] * cartDims[1] * cartDims[2]);

-        alternative_well_rate_init_ = EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
+            auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
+            this->parallel_well_info_.assign(parallel_wells.begin(),
+                                             parallel_wells.end());
+        }
+
+        const auto numProcs = ebosSimulator.gridView().comm().size();
+        this->not_on_process_ = [this, numProcs](const Well& well) {
+            if (numProcs == decltype(numProcs){1})
+                return false;
+
+            // Recall: false indicates NOT active!
+            const auto value = std::make_pair(well.name(), true);
+            auto candidate = std::lower_bound(this->parallel_well_info_.begin(),
+                                              this->parallel_well_info_.end(),
+                                              value);
+
+            return (candidate == this->parallel_well_info_.end())
+                || (*candidate != value);
+        };
+
+        this->alternative_well_rate_init_ =
+            EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
    }

    template<typename TypeTag>
@@ -911,6 +919,30 @@ namespace Opm {



+    template <typename TypeTag>
+    void BlackoilWellModel<TypeTag>::
+    inferLocalShutWells()
+    {
+        this->local_shut_wells_.clear();
+
+        const auto nw = this->numLocalWells();
+
+        auto used = std::vector<bool>(nw, false);
+        for (const auto& wellPtr : this->well_container_) {
+            used[wellPtr->indexOfWell()] = true;
+        }
+
+        for (auto wellID = 0; wellID < nw; ++wellID) {
+            if (! used[wellID]) {
+                this->local_shut_wells_.push_back(wellID);
+            }
+        }
+    }
+
+
+
+
+
    template <typename TypeTag>
    typename BlackoilWellModel<TypeTag>::WellInterfacePtr
    BlackoilWellModel<TypeTag>::
@@ -1586,21 +1618,13 @@ namespace Opm {



-
-    template<typename TypeTag>
+    template <typename TypeTag>
    void
    BlackoilWellModel<TypeTag>::
    calculateProductivityIndexValues(DeferredLogger& deferred_logger)
    {
-        if (! this->localWellsActive()) {
-            return;
-        }
-
        for (const auto& wellPtr : this->well_container_) {
-            wellPtr->updateProductivityIndex(this->ebosSimulator_,
-                                             this->prod_index_calc_[wellPtr->indexOfWell()],
-                                             this->wellState(),
-                                             deferred_logger);
+            this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
        }
    }

@@ -1608,6 +1632,48 @@ namespace Opm {



+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    calculateProductivityIndexValuesShutWells(const int reportStepIdx,
+                                              DeferredLogger& deferred_logger)
+    {
+        // For the purpose of computing PI/II values, it is sufficient to
+        // construct StandardWell instances only.  We don't need to form
+        // well objects that honour the 'isMultisegment()' flag of the
+        // corresponding "this->wells_ecl_[shutWell]".
+
+        for (const auto& shutWell : this->local_shut_wells_) {
+            auto wellPtr = this->template createTypedWellPointer
+                <StandardWell<TypeTag>>(shutWell, reportStepIdx);
+
+            wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
+                          this->local_num_cells_, this->B_avg_);
+
+            this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
+        }
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
+                                     DeferredLogger& deferred_logger)
+    {
+        wellPtr->updateProductivityIndex(this->ebosSimulator_,
+                                         this->prod_index_calc_[wellPtr->indexOfWell()],
+                                         this->wellState(),
+                                         deferred_logger);
+    }
+
+
+
+
+
    template<typename TypeTag>
    void
    BlackoilWellModel<TypeTag>::
@@ -2827,17 +2893,17 @@ namespace Opm {
            auto saved_previous_wgstate = this->prevWGState();
            this->commitWGState();

-            well_container_ = createWellContainer(timeStepIdx);
-            for (auto& well : well_container_) {
-                well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_);
-            }
+            this->well_container_ = this->createWellContainer(timeStepIdx);
+            this->inferLocalShutWells();

-            std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
-            for (auto& well : well_container_) {
-                well->updatePerforatedCell(is_cell_perforated_);
+            for (auto& wellPtr : this->well_container_) {
+                wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
+                              this->local_num_cells_, this->B_avg_);
            }

            this->calculateProductivityIndexValues(local_deferredLogger);
+            this->calculateProductivityIndexValuesShutWells(timeStepIdx, local_deferredLogger);
+
            this->commitWGState(std::move(saved_previous_wgstate));
        }

@@ -271,7 +271,7 @@ namespace Opm
    WellInterface<TypeTag>::
    phaseUsage() const
    {
-        assert(phase_usage_);
+        assert(phase_usage_ != nullptr);

        return *phase_usage_;
    }
@@ -154,7 +154,7 @@ namespace Opm
                const auto& well_info = this->wellMap().at(wname);
                const int connpos = well_info[1];
                const int num_perf_this_well = well_info[2];
-                int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well);
+                const int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well);

                for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) {
                    if (wells_ecl[w].getStatus() == Well::Status::OPEN) {
@@ -225,7 +225,7 @@ namespace Opm
                    }

                    auto it = prevState->wellMap().find(well.name());
-                    if ( it != end )
+                    if (it != end)
                    {
                        const int newIndex = w;
                        const int oldIndex = it->second[ 0 ];
@@ -283,14 +283,16 @@ namespace Opm
                        const int connpos = new_iter->second[1];
                        const int num_perf_this_well = new_iter->second[2];

-                        int num_perf_changed = (num_perf_old_well != num_perf_this_well) ? 1 : 0;
-                        num_perf_changed = parallel_well_info[w]->communication().sum(num_perf_changed);
-                        bool global_num_perf_same = (num_perf_changed == 0);
+                        const int num_perf_changed = parallel_well_info[w]->communication()
+                            .sum(static_cast<int>(num_perf_old_well != num_perf_this_well));
+                        const bool global_num_perf_same = num_perf_changed == 0;

-                        // copy perforation rates when the number of perforations is equal,
-                        // otherwise initialize perfphaserates to well rates divided by the number of perforations.

-                        if( global_num_perf_same )
+                        // copy perforation rates when the number of
+                        // perforations is equal, otherwise initialize
+                        // perfphaserates to well rates divided by the
+                        // number of perforations.
+                        if (global_num_perf_same)
                        {
                            int old_perf_phase_idx = oldPerf_idx_beg *np;
                            for (int perf_phase_idx = connpos*np;
@@ -299,15 +301,16 @@ namespace Opm
                                perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
                            }
                        } else {
-                            int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well);
+                            const int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well);
                            for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) {
                                for (int p = 0; p < np; ++p) {
                                    perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(global_num_perf_this_well);
                                }
                            }
                        }
+
                        // perfPressures
-                        if( global_num_perf_same )
+                        if (global_num_perf_same)
                        {
                            int oldPerf_idx = oldPerf_idx_beg;
                            for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
@@ -315,9 +318,10 @@ namespace Opm
                                perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ];
                            }
                        }
+
                        // perfSolventRates
                        if (pu.has_solvent) {
-                            if( global_num_perf_same )
+                            if (global_num_perf_same)
                            {
                                int oldPerf_idx = oldPerf_idx_beg;
                                for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
@@ -328,10 +332,12 @@ namespace Opm
                        }

                        // polymer injectivity related
-                        // here we did not consider the case that we close some perforation during the running
-                        // and also, wells can be shut and re-opened
+                        //
+                        // here we did not consider the case that we close
+                        // some perforation during the running and also,
+                        // wells can be shut and re-opened
                        if (pu.has_polymermw) {
-                            if( num_perf_old_well == num_perf_this_well )
+                            if (global_num_perf_same)
                            {
                                int oldPerf_idx = oldPerf_idx_beg;
                                for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
@@ -354,10 +360,13 @@ namespace Opm
                        }
                    }

-                    // If in the new step, there is no THP related target/limit anymore, its thp value should be
-                    // set to zero.
-                    const bool has_thp = well.isInjector() ? well.injectionControls(summary_state).hasControl(Well::InjectorCMode::THP)
+                    // If in the new step, there is no THP related
+                    // target/limit anymore, its thp value should be set to
+                    // zero.
+                    const bool has_thp = well.isInjector()
+                        ? well.injectionControls (summary_state).hasControl(Well::InjectorCMode::THP)
                        : well.productionControls(summary_state).hasControl(Well::ProducerCMode::THP);
+
                    if (!has_thp) {
                        thp()[w] = 0.0;
                    }
@@ -382,6 +391,7 @@ namespace Opm
                seg_pressdrop_friction_.assign(nw, 0.);
                seg_pressdrop_acceleration_.assign(nw, 0.);
            }
+
            updateWellsDefaultALQ(wells_ecl);
            do_glift_optimization_ = true;
        }
@@ -816,7 +826,8 @@ namespace Opm

        /// One rate pr well
        double solventWellRate(const int w) const {
-            return parallel_well_info_[w]->sumPerfValues(&perfRateSolvent_[0] + first_perf_index_[w], &perfRateSolvent_[0] + first_perf_index_[w+1]);
+            return parallel_well_info_[w]->sumPerfValues(&perfRateSolvent_[0] + first_perf_index_[w + 0],
+                                                         &perfRateSolvent_[0] + first_perf_index_[w + 1]);
        }

        /// One rate pr well connection.
@@ -825,7 +836,8 @@ namespace Opm

        /// One rate pr well
        double polymerWellRate(const int w) const {
-            return parallel_well_info_[w]->sumPerfValues(&perfRatePolymer_[0] + first_perf_index_[w], &perfRatePolymer_[0] + first_perf_index_[w+1]);
+            return parallel_well_info_[w]->sumPerfValues(&perfRatePolymer_[0] + first_perf_index_[w + 0],
+                                                         &perfRatePolymer_[0] + first_perf_index_[w + 1]);
        }

        /// One rate pr well connection.
@@ -834,7 +846,8 @@ namespace Opm

        /// One rate pr well
        double brineWellRate(const int w) const {
-            return parallel_well_info_[w]->sumPerfValues(&perfRateBrine_[0] + first_perf_index_[w], &perfRateBrine_[0] + first_perf_index_[w+1]);
+            return parallel_well_info_[w]->sumPerfValues(&perfRateBrine_[0] + first_perf_index_[w + 0],
+                                                         &perfRateBrine_[0] + first_perf_index_[w + 1]);
        }

        std::vector<double>& wellReservoirRates()
@@ -980,8 +993,19 @@ namespace Opm
        virtual void shutWell(int well_index) override {
            WellState::shutWell(well_index);
            const int np = numPhases();
-            for (int p = 0; p < np; ++p)
-                this->well_reservoir_rates_[np * well_index + p] = 0;
+
+            auto* resv = &this->well_reservoir_rates_[np*well_index + 0];
+            auto* wpi  = &this->productivity_index_[np*well_index + 0];
+
+            for (int p = 0; p < np; ++p) {
+                resv[p] = 0.0;
+                wpi[p]  = 0.0;
+            }
+
+            const auto first = this->first_perf_index_[well_index + 0]*np;
+            const auto last  = this->first_perf_index_[well_index + 1]*np;
+            std::fill(this->conn_productivity_index_.begin() + first,
+                      this->conn_productivity_index_.begin() + last, 0.0);
        }

        virtual void stopWell(int well_index) override {