Merge pull request #4429 from GitPaean/support_aquflux

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/AquiferAnalytical.hpp

@@ -226,7 +226,6 @@ protected:
         return FluidSystem::waterCompIdx;
     }
 
-
     void initQuantities()
     {
         // We reset the cumulative flux at the start of any simulation, so, W_flux = 0

opm/simulators/aquifers/AquiferConstantFlux.hpp

@@ -0,0 +1,147 @@
+/*
+  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>;
+
+    AquiferConstantFlux(const SingleAquiferFlux& aquifer,
+                        const std::vector<Aquancon::AquancCell>& connections,
+                        const Simulator& ebos_simulator)
+        : AquiferInterface<TypeTag>(aquifer.id, ebos_simulator)
+         , connections_(connections)
+         , aquifer_data_(aquifer)
+    {
+        // init_cumulative_flux is the flux volume from previoius running
+        this->initializeConnections();
+        connection_flux_.resize(this->connections_.size(), {0});
+    }
+
+    virtual ~AquiferConstantFlux() = default;
+
+    void updateAquifer(const SingleAquiferFlux& aquifer) {
+        aquifer_data_ = aquifer;
+    }
+
+    void initFromRestart(const data::Aquifers& /* aquiferSoln */) {
+    }
+
+    void initialSolutionApplied() {
+    }
+
+    void beginTimeStep() {
+    }
+
+    void endTimeStep() {
+        this->flux_rate_ = 0.;
+        for (const auto& q : this->connection_flux_) {
+            this->flux_rate_ += Opm::getValue(q);
+        }
+
+        this->cumulative_flux_ += this->flux_rate_ * 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->connection_flux_) {
+            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;
+        this->connection_flux_[idx] = fw * this->connections_[idx].effective_facearea;
+        rates[BlackoilIndices::conti0EqIdx + compIdx_()]
+                += this->connection_flux_[idx] / model.dofTotalVolume(cellIdx);
+    }
+
+private:
+    const std::vector<Aquancon::AquancCell>& connections_;
+    SingleAquiferFlux aquifer_data_;
+    std::vector<int> cellToConnectionIdx_;
+    std::vector<Eval> connection_flux_;
+    double flux_rate_ {};
+    double cumulative_flux_ = 0.;
+
+    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().compressedIndexForInterior(global_index);
+
+            if (cell_index < 0) {
+                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 if
+        //  the grid changed during the preprocessing?
+    }
+
+    // TODO: this is a function from AquiferAnalytical
+    int compIdx_() const
+    {
+        if (this->co2store_())
+            return FluidSystem::oilCompIdx;
+
+        return FluidSystem::waterCompIdx;
+    }
+};
+}
+
+#endif //OPM_AQUIFERCONSTANTFLUX_HPP

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -171,6 +171,11 @@ public:
                this->pressure_ == rhs.pressure_;
     }
 
+    double cumulativeFlux() const
+    {
+        return this->cumulative_flux_;
+    }
+
 private:
     void checkConnectsToReservoir()
     {

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

@@ -18,7 +18,12 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <opm/simulators/aquifers/AquiferConstantFlux.hpp>
+
 #include <opm/common/ErrorMacros.hpp>
+
+#include <algorithm>
+#include <memory>
 #include <stdexcept>
 
 namespace Opm
@@ -54,7 +59,35 @@ BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginEpisode()
-{}
+{
+    // 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;
+    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 BlackoilAquiferModel 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 aqu = dynamic_cast<AquiferConstantFlux<TypeTag>*> (find->get());
+            if (!aqu) {
+                // if the aquifers can return types easily, we might be able to give a better message with type information
+                const auto msg = fmt::format("Aquifer {} is updated with constant flux aquifer keyword AQUFLUX at report step {},"
+                                                       " while it might be specified to be a different type of aquifer before this. We do not support the conversion between"
+                                                       " different types of aquifer.\n", id, report_step);
+                OPM_THROW(std::runtime_error, msg);
+            }
+            aqu->updateAquifer(aquinfo);
+        }
+    }
+}
 
 template <typename TypeTag>
 void
@@ -168,6 +201,21 @@ BlackoilAquiferModel<TypeTag>::init()
         aquifers.push_back(std::move(aqf));
     }
 
+    for (const auto& [id, aq] : aquifer.aquflux()) {
+        // make sure not dummy constant flux aquifers
+        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) {

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}},

opm/simulators/utils/readDeck.cpp

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

parallelTests.cmake

@@ -133,6 +133,14 @@ add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_2aqu
                                      DIR aquifer-num
                                      TEST_ARGS --tolerance-cnv=0.000003 --time-step-control=pid --linear-solver=cpr_trueimpes)
 
+add_test_compare_parallel_simulation(CASENAME aquflux_01
+                                     FILENAME AQUFLUX-01
+                                     SIMULATOR flow
+                                     ABS_TOL ${abs_tol}
+                                     REL_TOL ${coarse_rel_tol_parallel}
+                                     DIR aquifers
+                                     TEST_ARGS --enable-tuning=true)
+
 add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_1aqu
                                      FILENAME 3D_1AQU_3CELLS
                                      SIMULATOR flow

regressionTests.cmake

@@ -186,6 +186,13 @@ add_test_compareECLFiles(CASENAME numerical_aquifer_3d_1aqu
                          DIR aquifer-num
                          TEST_ARGS --tolerance-cnv=0.00003 --time-step-control=pid --linear-solver=cpr_trueimpes)
 
+add_test_compareECLFiles(CASENAME aquflux_01
+                         FILENAME AQUFLUX-01
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR aquifers)
+
 add_test_compareECLFiles(CASENAME spe3
                          FILENAME SPE3CASE1
                          SIMULATOR flow

tests/update_reference_data.sh

@@ -83,6 +83,7 @@ tests[ctaquifer_2d_oilwater]="flow aquifer-oilwater 2D_OW_CTAQUIFER"
 tests[fetkovich_2d]="flow aquifer-fetkovich 2D_FETKOVICHAQUIFER"
 tests[numerical_aquifer_3d_1aqu]="flow aquifer-num 3D_1AQU_3CELLS"
 tests[numerical_aquifer_3d_2aqu]="flow aquifer-num 3D_2AQU_NUM"
+tests[aquflux_01]="flow aquifers AQUFLUX-01"
 tests[msw_2d_h]="flow msw_2d_h 2D_H__"
 tests[msw_3d_hfa]="flow msw_3d_hfa 3D_MSW"
 tests[polymer_oilwater]="flow polymer_oilwater 2D_OILWATER_POLYMER"