following the AQUFLUX design change from opm-common

opm/simulators/aquifers/AquiferConstantFlux.hpp

@@ -37,18 +37,16 @@ public:
     static constexpr int numEq = BlackoilIndices::numEq;
     using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
 
-
+    AquiferConstantFlux(const AquiferFlux& aquifer,
-    // TODO: we need to pass in the previous flux volume
-    AquiferConstantFlux(const std::shared_ptr<AquiferFlux>& aquifer,
                         const std::vector<Aquancon::AquancCell>& connections,
                         const Simulator& ebos_simulator,
                         const double init_cumulative_flux = 0.)
-        : AquiferInterface<TypeTag>(aquifer->id, ebos_simulator)
+        : AquiferInterface<TypeTag>(aquifer.id, ebos_simulator)
          , connections_(connections)
          , aquifer_data_(aquifer)
          , cumulative_flux_(init_cumulative_flux)
     {
-        // flux_volume is the flux volume from previoius running
+        // init_cumulative_flux is the flux volume from previoius running
         this->initializeConnections();
         connection_flux_.resize(this->connections_.size(), {0});
     }
@@ -63,9 +61,7 @@ public:
     }
 
     void initialSolutionApplied() {
-        // Note: we can not do this here
+        // this->initializeConnections();
-        // with the current way, we put the AQUFLUX in the first report step of Schedule
-        // Maybe it might bring some undesiable consequence to remove it from the solution
     }
 
     void beginTimeStep() {
@@ -86,7 +82,7 @@ public:
     data::AquiferData aquiferData() const
     {
         data::AquiferData data;
-        data.aquiferID = this->aquifer_data_->id;
+        data.aquiferID = this->aquifer_data_.id;
         // pressure for constant flux aquifer is 0
         data.pressure = 0.;
         data.fluxRate = 0.;
@@ -113,7 +109,7 @@ public:
             throw std::logic_error("Invalid intensive quantities cache detected in AquiferAnalytical::addToSource()");
         }
 
-        const double fw = this->aquifer_data_->flux;
+        const double fw = this->aquifer_data_.flux;
         // const double m = this->connections_[idx].influx_coeff;
         this->connection_flux_[idx] = fw * this->connections_[idx].effective_facearea;
         rates[BlackoilIndices::conti0EqIdx + compIdx_()]
@@ -128,18 +124,14 @@ public:
 
 private:
     const std::vector<Aquancon::AquancCell> connections_;
-    std::shared_ptr<AquiferFlux> aquifer_data_;
+    AquiferFlux aquifer_data_;
-    // TODO: for simple case, faceArea_connected_ is not needed here, since it is calculated when parsing
-    // But if the grid change, not sure how to handle
-    // std::vector<double> faceArea_connected_;
     std::vector<int> cellToConnectionIdx_;
     std::vector<Eval> connection_flux_;
-    double flux_rate_;
+    double flux_rate_ {};
     double cumulative_flux_ = 0.;
 
 
     void initializeConnections() {
-        // this->faceArea_connected_.resize(this->size(), {0});
 
         this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
         const auto& gridView = this->ebos_simulator_.vanguard().gridView();
@@ -156,9 +148,11 @@ private:
 
             this->cellToConnectionIdx_[cell_index] = idx;
         }
-        // TODO: at the moment, we are using the effective_facearea from the parser. Should we update the facearea here?
+        // 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?
     }
-    // TODO: function from AquiferAnalytical
+
+    // TODO: this is a function from AquiferAnalytical
     int compIdx_() const
     {
         if (this->co2store_())

opm/simulators/aquifers/AquiferInterface.hpp

@@ -56,8 +56,6 @@ public:
     virtual void beginTimeStep() = 0;
     virtual void endTimeStep() = 0;
 
-
-
     virtual data::AquiferData aquiferData() const = 0;
 
     template <class Context>

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -17,8 +17,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-// TODO: remove this include
-#include <iostream>
 
 #include <opm/simulators/aquifers/AquiferConstantFlux.hpp>
 
@@ -59,26 +57,25 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginEpisode()
 {
-    // TODO: not totally sure this is the function should be used.
+    // probably function name beginReportStep() is more appropriate.
     // basically, we want to update the aquifer related information from SCHEDULE setup in this section
     // it is the beginning of a report step
-
     const auto& connections = this->simulator_.vanguard().eclState().aquifer().connections();
     const int report_step  = this->simulator_.episodeIndex();
-    const auto& aqufluxs = this->simulator_.vanguard().schedule()[report_step].aqufluxs;// .aquflu// simulator.vanguard().schedule()[reportStepIdx].events()
+    const auto& aqufluxs = this->simulator_.vanguard().schedule()[report_step].aqufluxs;
     for (const auto& elem : aqufluxs) {
         const int id = elem.first;
         auto find = std::find_if(begin(this->aquifers), end(this->aquifers), [id](auto& v){ return v->aquiferID() == id; });
         if (find == this->aquifers.end()) {
-            // the aquifer id does not exist in aquifers yet
+            // the aquifer id does not exist in BlackoilAquiferModel yet
-            const auto& aquinfo = elem.second;
+            const auto& aquinfo = *(elem.second);
-            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo->id), this->simulator_);
+            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo.id), this->simulator_);
             this->aquifers.push_back(std::move(aqf));
         } else {
             const double prev_cumulative_flux = (*find)->cumulativeFlux();
-            const auto& aquinfo = elem.second;
+            const auto& aquinfo = *(elem.second);
-            // TODO: it should be improved to be something like a update instead of creating a new one
+            // TODO: it should be improved to be something like to update the related information instead of creating a new one
-            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo->id), this->simulator_, prev_cumulative_flux);
+            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo.id), this->simulator_, prev_cumulative_flux);
             *find = std::move(aqf);
         }
     }
@@ -196,6 +193,20 @@ BlackoilAquiferModel<TypeTag>::init()
         aquifers.push_back(std::move(aqf));
     }
 
+    for (const auto& [id, aq] : aquifer.aquflux()) {
+        if ( !aq.active ) continue;
+
+        if (!connections.hasAquiferConnections(id)) {
+            auto msg = fmt::format("No valid connections for constant flux aquifer {}, aquifer {} will be ignored.",
+                                   id, id);
+            OpmLog::warning(msg);
+            continue;
+        }
+
+        auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aq, connections.getConnections(id), this->simulator_);
+        this->aquifers.push_back(std::move(aqf));
+    }
+
     if (aquifer.hasNumericalAquifer()) {
         const auto& num_aquifers = aquifer.numericalAquifers().aquifers();
         for ([[maybe_unused]]const auto& [id, aqu] : num_aquifers) {
@@ -203,9 +214,6 @@ BlackoilAquiferModel<TypeTag>::init()
             aquifers.push_back(std::move(aqf));
         }
     }
-
-    // first time handle constant flux aquifers, which is stored in the schedule. Other aquifer types might also be refactored later
-    // to be able to be updated through SCHEDULE.
 }
 
 template<typename TypeTag>

opm/simulators/utils/readDeck.cpp

@@ -268,6 +268,8 @@ namespace {
                                            errorGuard);
         }
 
+        eclipseState->appendAqufluxSchedule(schedule->getAquiferFluxListEnd());
+
         if (Opm::OpmLog::hasBackend("STDOUT_LOGGER")) {
             // loggers might not be set up!
             setupMessageLimiter((*schedule)[0].message_limits(), "STDOUT_LOGGER");