Merge pull request #4704 from atgeirr/wells-setprimaryvariables

Add methods for getting and setting primary variables.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -45,6 +45,7 @@
 
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
 #include <opm/simulators/flow/countGlobalCells.hpp>
+#include <opm/simulators/flow/SubDomain.hpp>
 #include <opm/simulators/wells/BlackoilWellModelGeneric.hpp>
 #include <opm/simulators/wells/BlackoilWellModelGuideRates.hpp>
 #include <opm/simulators/wells/GasLiftSingleWell.hpp>
@@ -137,6 +138,8 @@ namespace Opm {
             using AverageRegionalPressureType = RegionAverageCalculator::
                 AverageRegionalPressure<FluidSystem, std::vector<int> >;
 
+            using Domain = SubDomain<Grid>;
+
             BlackoilWellModel(Simulator& ebosSimulator);
 
             void init();
@@ -328,6 +331,12 @@ namespace Opm {
 
             int numLocalNonshutWells() const;
 
+            void logPrimaryVars() const;
+            std::vector<double> getPrimaryVarsDomain(const Domain& domain) const;
+            void setPrimaryVarsDomain(const Domain& domain, const std::vector<double>& vars);
+
+            void setupDomains(const std::vector<Domain>& domains);
+
         protected:
             Simulator& ebosSimulator_;
 
@@ -377,6 +386,9 @@ namespace Opm {
 
             std::vector<Scalar> B_avg_{};
 
+            // Keep track of the domain of each well, if using subdomains.
+            std::map<std::string, int> well_domain_;
+
             const Grid& grid() const
             { return ebosSimulator_.vanguard().grid(); }
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -2094,6 +2094,58 @@ namespace Opm {
     }
 
 
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    logPrimaryVars() const
+    {
+        std::ostringstream os;
+        for (const auto& w : well_container_) {
+            os << w->name() << ":";
+            auto pv = w->getPrimaryVars();
+            for (const double v : pv) {
+                os << ' ' << v;
+            }
+            os << '\n';
+        }
+        OpmLog::debug(os.str());
+    }
+
+
+
+    template <typename TypeTag>
+    std::vector<double>
+    BlackoilWellModel<TypeTag>::
+    getPrimaryVarsDomain(const Domain& domain) const
+    {
+        std::vector<double> ret;
+        for (const auto& well : well_container_) {
+            if (well_domain_.at(well->name()) == domain.index) {
+                const auto& pv = well->getPrimaryVars();
+                ret.insert(ret.end(), pv.begin(), pv.end());
+            }
+        }
+        return ret;
+    }
+
+
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    setPrimaryVarsDomain(const Domain& domain, const std::vector<double>& vars)
+    {
+        std::size_t offset = 0;
+        for (auto& well : well_container_) {
+            if (well_domain_.at(well->name()) == domain.index) {
+                int num_pri_vars = well->setPrimaryVars(vars.begin() + offset);
+                offset += num_pri_vars;
+            }
+        }
+        assert(offset == vars.size());
+    }
+
+
 
     template <typename TypeTag>
     void
@@ -2118,7 +2170,45 @@ namespace Opm {
     }
 
 
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    setupDomains(const std::vector<Domain>& domains)
+    {
+        OPM_BEGIN_PARALLEL_TRY_CATCH();
+        // TODO: This loop nest may be slow for very large numbers
+        // of domains and wells, but that has not been observed on
+        // tests so far. Using the partition vector instead would
+        // be faster if we need to change.
+        for (const auto& wellPtr : this->well_container_) {
+            const int first_well_cell = wellPtr->cells()[0];
+            for (const auto& domain : domains) {
+                const bool found = std::binary_search(domain.cells.begin(), domain.cells.end(), first_well_cell);
+                if (found) {
+                    // Assuming that if the first well cell is found in a domain,
+                    // then all of that well's cells are in that same domain.
+                    well_domain_[wellPtr->name()] = domain.index;
+                    // Verify that all of that well's cells are in that same domain.
+                    for (int well_cell : wellPtr->cells()) {
+                        if (!std::binary_search(domain.cells.begin(), domain.cells.end(), well_cell)) {
+                            OPM_THROW(std::runtime_error, "Well found on multiple domains.");
+                        }
+                    }
+                }
+            }
+        }
+        OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::setupDomains(): well found on multiple domains.",
+                                   ebosSimulator_.gridView().comm());
 
-
+        // Write well/domain info to the DBG file.
+        if (ebosSimulator_.gridView().comm().rank() == 0) {
+            std::ostringstream os;
+            os << "Well name      Domain\n";
+            for (const auto& [wname, domain] : well_domain_) {
+                os << wname << std::setw(21 - wname.size()) << domain << '\n';
+            }
+            OpmLog::debug(os.str());
+        }
+    }
 
 } // namespace Opm

opm/simulators/wells/MultisegmentWell.hpp

@@ -166,6 +166,10 @@ namespace Opm
                                         const double alq_value,
                                         DeferredLogger& deferred_logger) const override;
 
+        std::vector<double> getPrimaryVars() const override;
+
+        int setPrimaryVars(std::vector<double>::const_iterator it) override;
+
     protected:
 
         // regularize msw equation

opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp

@@ -137,10 +137,18 @@ public:
     //! \brief Get WQTotal.
     EvalWell getWQTotal() const;
 
-    //! \brief Returns a const ref to an evaluation.
-    const std::array<EvalWell,numWellEq>& eval(const int idx) const
+    //! \brief Returns a const ref to an array of evaluations.
+    const std::array<EvalWell, numWellEq>& eval(const int idx) const
     { return evaluation_[idx]; }
 
+    //! \brief Returns a value array.
+    const std::array<Scalar, numWellEq>& value(const int idx) const
+    { return value_[idx]; }
+
+    //! \brief Set a value array. Note that this does not also set the corresponding evaluation.
+    void setValue(const int idx, const std::array<Scalar, numWellEq>& val)
+    { value_[idx] = val; }
+
 private:
     //! \brief Handle non-reasonable fractions due to numerical overshoot.
     void processFractions(const int seg);

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1972,4 +1972,41 @@ namespace Opm
     }
 
 
+
+
+    template <typename TypeTag>
+    std::vector<double>
+    MultisegmentWell<TypeTag>::
+    getPrimaryVars() const
+    {
+        const int num_seg = this->numberOfSegments();
+        constexpr int numWellEq = MSWEval::numWellEq;
+        std::vector<double> retval(num_seg * numWellEq);
+        for (int ii = 0; ii < num_seg; ++ii) {
+            const auto& pv = this->primary_variables_.value(ii);
+            std::copy(pv.begin(), pv.end(), retval.begin() + ii * numWellEq);
+        }
+        return retval;
+    }
+
+
+
+
+    template <typename TypeTag>
+    int
+    MultisegmentWell<TypeTag>::
+    setPrimaryVars(std::vector<double>::const_iterator it)
+    {
+        const int num_seg = this->numberOfSegments();
+        constexpr int numWellEq = MSWEval::numWellEq;
+        std::array<double, numWellEq> tmp;
+        for (int ii = 0; ii < num_seg; ++ii) {
+            const auto start = it + num_seg * numWellEq;
+            std::copy(start, start + numWellEq, tmp.begin());
+            this->primary_variables_.setValue(ii, tmp);
+        }
+        return num_seg * numWellEq;
+    }
+
+
 } // namespace Opm

opm/simulators/wells/StandardWell.hpp

@@ -253,6 +253,9 @@ namespace Opm
                                     double* connII,
                                     DeferredLogger& deferred_logger) const;
 
+        std::vector<double> getPrimaryVars() const override;
+
+        int setPrimaryVars(std::vector<double>::const_iterator it) override;
 
     protected:
         bool regularize_;

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -135,7 +135,7 @@ public:
     //! \brief Returns scaled rate for a component.
     EvalWell getQs(const int compIdx) const;
 
-    //! \brief Returns a const ref to an evaluation.
+    //! \brief Returns a value.
     Scalar value(const int idx) const
     { return value_[idx]; }
 
@@ -143,6 +143,10 @@ public:
     const EvalWell& eval(const int idx) const
     { return evaluation_[idx]; }
 
+    //! \brief Set a value. Note that this does not also set the corresponding evaluation.
+    void setValue(const int idx, const Scalar val)
+    { value_[idx] = val; }
+
 private:
     //! \brief Calculate a relaxation factor for producers.
     //! \details To avoid overshoot of the fractions which might result in negative rates.

opm/simulators/wells/StandardWell_impl.hpp

@@ -2543,4 +2543,38 @@ namespace Opm
         const auto mt     = std::accumulate(mobility.begin(), mobility.end(), zero);
         connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
     }
+
+
+
+
+
+    template <typename TypeTag>
+    std::vector<double>
+    StandardWell<TypeTag>::
+    getPrimaryVars() const
+    {
+        const int num_pri_vars = this->primary_variables_.numWellEq();
+        std::vector<double> retval(num_pri_vars);
+        for (int ii = 0; ii < num_pri_vars; ++ii) {
+            retval[ii] = this->primary_variables_.value(ii);
+        }
+        return retval;
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    int
+    StandardWell<TypeTag>::
+    setPrimaryVars(std::vector<double>::const_iterator it)
+    {
+        const int num_pri_vars = this->primary_variables_.numWellEq();
+        for (int ii = 0; ii < num_pri_vars; ++ii) {
+            this->primary_variables_.setValue(ii, it[ii]);
+        }
+        return num_pri_vars;
+    }
+
 } // namespace Opm

opm/simulators/wells/WellInterface.hpp

@@ -328,6 +328,16 @@ public:
         return connectionRates_;
     }
 
+    virtual std::vector<double> getPrimaryVars() const
+    {
+        return {};
+    }
+
+    virtual int setPrimaryVars(std::vector<double>::const_iterator)
+    {
+        return 0;
+    }
+
 protected:
     // simulation parameters
     const ModelParameters& param_;