WIP of supporting AQUFLUX from simulator side

CMakeLists_files.cmake

@@ -308,6 +308,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/aquifers/AquiferAnalytical.hpp
   opm/simulators/aquifers/AquiferCarterTracy.hpp
   opm/simulators/aquifers/AquiferFetkovich.hpp
+  opm/simulators/aquifers/AquiferConstantFlux.hpp
   opm/simulators/aquifers/AquiferInterface.hpp
   opm/simulators/aquifers/AquiferNumerical.hpp
   opm/simulators/aquifers/BlackoilAquiferModel.hpp

opm/simulators/aquifers/AquiferConstantFlux.hpp

@@ -0,0 +1,163 @@
+/*
+  Copyright (C) 2023 Equinor
+
+  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 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/>.
+ */
+
+#ifndef OPM_AQUIFERCONSTANTFLUX_HPP
+#define OPM_AQUIFERCONSTANTFLUX_HPP
+
+#include <opm/simulators/aquifers/AquiferInterface.hpp>
+
+#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
+#include <opm/input/eclipse/EclipseState/Aquifer/AquiferFlux.hpp>
+
+namespace Opm {
+template<typename TypeTag>
+class AquiferConstantFlux : public AquiferInterface<TypeTag> {
+public:
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
+    static constexpr int numEq = BlackoilIndices::numEq;
+    using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
+
+
+    // 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)
+        : AquiferInterface<TypeTag>(aquifer->id, ebos_simulator)
+         , connections_(connections)
+         , aquifer_data_(aquifer)
+    {
+        // flux_volume is the flux volume from previoius running
+        this->initializeConnections();
+        flux_rate_.resize(this->connections_.size(), {0});
+    }
+
+    virtual ~AquiferConstantFlux() = default;
+
+    /* void updateAquifer(const std::shared_ptr<AquiferFlux>& aquifer) {
+        aquifer_data_ = aquifer;
+    } */
+
+    void initFromRestart(const data::Aquifers& /* aquiferSoln */) {
+    }
+
+    void initialSolutionApplied() {
+        // Note: we can not do this here
+        // 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() {
+    }
+
+    void endTimeStep() {
+        for (const auto& q : this->flux_rate_) {
+            this->cumulative_flux_ += Opm::getValue(q) * this->ebos_simulator_.timeStepSize();
+        }
+    }
+
+    data::AquiferData aquiferData() const
+    {
+        data::AquiferData data;
+        data.aquiferID = this->aquifer_data_->id;
+        // pressure for constant flux aquifer is 0
+        data.pressure = 0.;
+        data.fluxRate = 0.;
+        for (const auto& q : this->flux_rate_) {
+            data.fluxRate += q.value();
+        }
+        data.volume = this->cumulative_flux_;
+        // not totally sure whether initPressure matters
+        data.initPressure = 0.;
+        return data;
+    }
+
+    void addToSource(RateVector& rates,
+                     const unsigned cellIdx,
+                     const unsigned timeIdx) {
+        const auto& model = this->ebos_simulator_.model();
+
+        const int idx = this->cellToConnectionIdx_[cellIdx];
+        if (idx < 0)
+            return;
+
+        const auto* intQuantsPtr = model.cachedIntensiveQuantities(cellIdx, timeIdx);
+        if (intQuantsPtr == nullptr) {
+            throw std::logic_error("Invalid intensive quantities cache detected in AquiferAnalytical::addToSource()");
+        }
+
+        const double fw = this->aquifer_data_->flux;
+        // const double m = this->connections_[idx].influx_coeff;
+        this->flux_rate_[idx] = fw * this->connections_[idx].effective_facearea;
+        rates[BlackoilIndices::conti0EqIdx + compIdx_()]
+                += this->flux_rate_[idx] / model.dofTotalVolume(cellIdx);
+    }
+
+    // TODO: repeated function from AquiferAnalytical
+    std::size_t size() const
+    {
+        return this->connections_.size();
+    }
+
+private:
+    const std::vector<Aquancon::AquancCell> connections_;
+    std::shared_ptr<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> 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();
+        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</*codim=*/ 0>();
+            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)
+                continue;
+
+            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: function from AquiferAnalytical
+    int compIdx_() const
+    {
+        if (this->co2store_())
+            return FluidSystem::oilCompIdx;
+
+        return FluidSystem::waterCompIdx;
+    }
+};
+}
+
+#endif //OPM_AQUIFERCONSTANTFLUX_HPP

opm/simulators/aquifers/AquiferInterface.hpp

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

opm/simulators/aquifers/BlackoilAquiferModel.hpp

@@ -122,6 +122,7 @@ protected:
 
     Simulator& simulator_;
 
+    // TODO: possibly better to use unorder_map here for aquifers
     std::vector<std::unique_ptr<AquiferInterface<TypeTag>>> aquifers;
 
     // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -17,6 +17,10 @@
   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>
 
 #include <opm/common/ErrorMacros.hpp>
 #include <stdexcept>
@@ -54,7 +58,30 @@ BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginEpisode()
-{}
+{
+    // TODO: not totally sure this is the function should be used.
+    // 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()
+    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
+            const auto& aquinfo = elem.second;
+            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo->id), this->simulator_);
+            this->aquifers.push_back(std::move(aqf));
+        } else {
+            const auto& aquinfo = elem.second;
+            auto aqf = std::make_unique<AquiferConstantFlux<TypeTag>>(aquinfo, connections.getConnections(aquinfo->id), this->simulator_);
+            *find = std::move(aqf);
+            // create
+        }
+    }
+}
 
 template <typename TypeTag>
 void
@@ -175,6 +202,9 @@ 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/UnsupportedFlowKeywords.cpp

@@ -59,7 +59,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"AQUCHWAT", {true, std::nullopt}},
         {"AQUCWFAC", {true, std::nullopt}},
         {"AQUFET", {true, std::string{"Use the AQUFETP keyword instead"}}},
-        {"AQUFLUX", {true, std::nullopt}},
         {"AQUNNC", {true, std::nullopt}},
         {"AUTOCOAR", {true, std::nullopt}},
         {"AUTOREF", {true, std::nullopt}},