Merge pull request #4576 from bska/face-area-fraction

Calculate Local Face Area Fraction for Aquifer's Connections

opm/simulators/aquifers/AquiferAnalytical.hpp

@@ -87,17 +87,45 @@ public:
                                           BlackoilIndices::numPhases>;
 
     // Constructor
-    AquiferAnalytical(int aqID,
+    AquiferAnalytical(const int aqID,
                       const std::vector<Aquancon::AquancCell>& connections,
                       const Simulator& ebosSimulator)
         : AquiferInterface<TypeTag>(aqID, ebosSimulator)
         , connections_(connections)
     {
+        this->initializeConnectionMappings();
     }
 
     // Destructor
     virtual ~AquiferAnalytical()
+    {}
+
+    void computeFaceAreaFraction(const std::vector<double>& total_face_area) override
     {
+        assert (total_face_area.size() >= static_cast<std::vector<double>::size_type>(this->aquiferID()));
+
+        const auto tfa = total_face_area[this->aquiferID() - 1];
+        const auto eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
+
+        if (tfa < eps_sqrt) {
+            this->alphai_.assign(this->size(), Scalar{0});
+        }
+        else {
+            std::transform(this->faceArea_connected_.begin(),
+                           this->faceArea_connected_.end(),
+                           this->alphai_.begin(),
+                           [tfa](const Scalar area)
+                           {
+                               return area / tfa;
+                           });
+        }
+
+        this->area_fraction_ = this->totalFaceArea() / tfa;
+    }
+
+    double totalFaceArea() const override
+    {
+        return this->total_face_area_;
     }
 
     void initFromRestart(const data::Aquifers& aquiferSoln) override
@@ -108,8 +136,9 @@ public:
 
         this->assignRestartData(xaqPos->second);
 
-        this->W_flux_ = xaqPos->second.volume;
+        this->W_flux_ = xaqPos->second.volume * this->area_fraction_;
         this->pa0_ = xaqPos->second.initPressure;
+
         this->solution_set_from_restart_ = true;
     }
 
@@ -230,15 +259,13 @@ protected:
             W_flux_ = Scalar{0};
         }
 
-        // We next get our connections to the aquifer and initialize these quantities using the initialize_connections
-        // function
-        initializeConnections();
-        calculateAquiferCondition();
-        calculateAquiferConstants();
+        this->initializeConnectionDepths();
+        this->calculateAquiferCondition();
+        this->calculateAquiferConstants();
 
-        pressure_previous_.resize(this->connections_.size(), Scalar{0});
-        pressure_current_.resize(this->connections_.size(), Scalar{0});
-        Qai_.resize(this->connections_.size(), Scalar{0});
+        this->pressure_previous_.resize(this->size(), Scalar{0});
+        this->pressure_current_.resize(this->size(), Scalar{0});
+        this->Qai_.resize(this->size(), Scalar{0});
     }
 
     void updateCellPressure(std::vector<Eval>& pressure_water,
@@ -257,54 +284,54 @@ protected:
         pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
     }
 
-    void initializeConnections()
+    void initializeConnectionMappings()
     {
-        this->cell_depth_.resize(this->size(), this->aquiferDepth());
         this->alphai_.resize(this->size(), 1.0);
         this->faceArea_connected_.resize(this->size(), Scalar{0});
 
-        // Translate the C face tag into the enum used by opm-parser's TransMult class
-        FaceDir::DirEnum faceDirection;
-
-        bool has_active_connection_on_proc = false;
-
-        // denom_face_areas is the sum of the areas connected to an aquifer
-        Scalar denom_face_areas{0};
+        // total_face_area_ is the sum of the areas connected to an aquifer
+        this->total_face_area_ = Scalar{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);
+            if (cell_index < 0) {
+                continue;
+            }
+
             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)
+            // The global_index is not part of this grid
+            if (elemIt->partitionType() != Dune::InteriorEntity) {
                 continue;
-
-            has_active_connection_on_proc = true;
+            }
 
             this->cellToConnectionIdx_[cell_index] = idx;
-            this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
         }
-        // get areas for all connections
-        ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
-        for (const auto& elem : elements(gridView)) {
-            unsigned cell_index = elemMapper.index(elem);
-            int idx = this->cellToConnectionIdx_[cell_index];
 
-            // only deal with connections given by the aquifer
-            if( idx < 0)
+        // Translate the C face tag into the enum used by opm-parser's TransMult class
+        FaceDir::DirEnum faceDirection;
+
+        // Get areas for all connections
+        const auto& elemMapper = this->ebos_simulator_.model().dofMapper();
+        for (const auto& elem : elements(gridView)) {
+            const unsigned cell_index = elemMapper.index(elem);
+            const int idx = this->cellToConnectionIdx_[cell_index];
+
+            // Only deal with connections given by the aquifer
+            if (idx < 0) {
                 continue;
+            }
 
             for (const auto& intersection : intersections(gridView, elem)) {
-                // only deal with grid boundaries
-                if (!intersection.boundary())
+                // Only deal with grid boundaries
+                if (! intersection.boundary()) {
                     continue;
+                }
 
-                int insideFaceIdx  = intersection.indexInInside();
-                switch (insideFaceIdx) {
+                switch (intersection.indexInInside()) {
                 case 0:
                     faceDirection = FaceDir::XMinus;
                     break;
@@ -328,48 +355,38 @@ protected:
                               "Internal error in initialization of aquifer.");
                 }
 
-
                 if (faceDirection == this->connections_[idx].face_dir) {
                     this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
                     break;
                 }
             }
-            denom_face_areas += this->faceArea_connected_.at(idx);
-        }
 
-        const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
-        comm.sum(&denom_face_areas, 1);
-        const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
-        for (std::size_t idx = 0; idx < this->size(); ++idx) {
-            // Protect against division by zero NaNs.
-            this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
-                ? Scalar{0}
-                : this->faceArea_connected_.at(idx) / denom_face_areas;
-        }
-
-        if (this->solution_set_from_restart_) {
-            this->rescaleProducedVolume(has_active_connection_on_proc);
+            this->total_face_area_ += this->faceArea_connected_.at(idx);
         }
     }
 
-    void rescaleProducedVolume(const bool has_active_connection_on_proc)
+    void initializeConnectionDepths()
     {
-        // Needed in parallel restart to approximate influence of aquifer
-        // being "owned" by a subset of the parallel processes.  If the
-        // aquifer is fully owned by a single process--i.e., if all cells
-        // connecting to the aquifer are on a single process--then this_area
-        // is tot_area on that process and zero elsewhere.
+        this->cell_depth_.resize(this->size(), this->aquiferDepth());
 
-        const auto this_area = has_active_connection_on_proc
-            ? std::accumulate(this->alphai_.begin(),
-                              this->alphai_.end(),
-                              Scalar{0})
-            : Scalar{0};
+        const auto& gridView = this->ebos_simulator_.vanguard().gridView();
+        for (std::size_t idx = 0; idx < this->size(); ++idx) {
+            const int cell_index = this->ebos_simulator_.vanguard()
+                .compressedIndex(this->connections_[idx].global_index);
+            if (cell_index < 0) {
+                continue;
+            }
 
-        const auto tot_area = this->ebos_simulator_.vanguard()
-            .grid().comm().sum(this_area);
+            auto elemIt = gridView.template begin</*codim=*/ 0>();
+            std::advance(elemIt, cell_index);
 
-        this->W_flux_ *= this_area / tot_area;
+            // The global_index is not part of this grid
+            if (elemIt->partitionType() != Dune::InteriorEntity) {
+                continue;
+            }
+
+            this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
+        }
     }
 
     // This function is for calculating the aquifer properties from equilibrium state with the reservoir
@@ -433,9 +450,13 @@ protected:
     std::optional<Scalar> Ta0_{}; // initial aquifer temperature
     Scalar rhow_{};
 
+    Scalar total_face_area_{};
+    Scalar area_fraction_{Scalar{1}};
+
     Eval W_flux_;
 
     bool solution_set_from_restart_ {false};
+    bool has_active_connection_on_proc_{false};
 };
 
 } // namespace Opm

opm/simulators/aquifers/AquiferConstantFlux.hpp

@@ -32,8 +32,10 @@
 #include <opm/material/densead/Evaluation.hpp>
 
 #include <algorithm>
+#include <cassert>
 #include <numeric>
 #include <stdexcept>
+#include <vector>
 
 namespace Opm {
 
@@ -58,7 +60,7 @@ public:
         , aquifer_data_            (aquifer)
         , connection_flux_         (connections_.size(), Eval{0})
     {
-        this->initializeConnections();
+        this->total_face_area_ = this->initializeConnections();
     }
 
     static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator)
@@ -71,6 +73,19 @@ public:
 
     virtual ~AquiferConstantFlux() = default;
 
+    void computeFaceAreaFraction(const std::vector<double>& total_face_area) override
+    {
+        assert (total_face_area.size() >= static_cast<std::vector<double>::size_type>(this->aquiferID()));
+
+        this->area_fraction_ = this->totalFaceArea()
+            / total_face_area[this->aquiferID() - 1];
+    }
+
+    double totalFaceArea() const override
+    {
+        return this->total_face_area_;
+    }
+
     void updateAquifer(const SingleAquiferFlux& aquifer)
     {
         aquifer_data_ = aquifer;
@@ -149,9 +164,16 @@ private:
     std::vector<int> cellToConnectionIdx_{};
     double flux_rate_{};
     double cumulative_flux_{};
+    double total_face_area_{0.0};
+    double area_fraction_{1.0};
+
+    double initializeConnections()
+    {
+        auto connected_face_area = 0.0;
+
+        this->cellToConnectionIdx_
+            .resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
 
-    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()
@@ -162,11 +184,25 @@ private:
             }
 
             this->cellToConnectionIdx_[cell_index] = idx;
+
+            connected_face_area += this->connections_[idx].effective_facearea;
         }
 
         // 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?
+
+        return connected_face_area;
+    }
+
+    double computeFaceAreaFraction(const double connected_face_area) const
+    {
+        const auto tot_face_area = this->ebos_simulator_.vanguard()
+            .grid().comm().sum(connected_face_area);
+
+        return (tot_face_area > 0.0)
+            ? connected_face_area / tot_face_area
+            : 0.0;
     }
 
     // TODO: this is a function from AquiferAnalytical

opm/simulators/aquifers/AquiferInterface.hpp

@@ -58,6 +58,9 @@ public:
 
     virtual data::AquiferData aquiferData() const = 0;
 
+    virtual void computeFaceAreaFraction(const std::vector<double>& total_face_area) = 0;
+    virtual double totalFaceArea() const = 0;
+
     template <class Context>
     void addToSource(RateVector& rates,
                      const Context& context,

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -154,6 +154,14 @@ public:
         this->cumulative_flux_ = 0.;
     }
 
+    void computeFaceAreaFraction(const std::vector<double>& /*total_face_area*/) override
+    {}
+
+    double totalFaceArea() const override
+    {
+        return 1.0;
+    }
+
     template<class Serializer>
     void serializeOp(Serializer& serializer)
     {

opm/simulators/aquifers/BlackoilAquiferModel.hpp

@@ -140,6 +140,8 @@ private:
     createAnalyticAquiferPointer(const AquiferData& aqData,
                                  const int          aquiferID,
                                  std::string_view   aqType) const;
+
+    void computeConnectionAreaFraction() const;
 };
 
 

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -45,17 +45,23 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::initialSolutionApplied()
 {
-    for (auto& aquifer : aquifers)
+    this->computeConnectionAreaFraction();
+
+    for (auto& aquifer : this->aquifers) {
         aquifer->initialSolutionApplied();
     }
+}
 
 template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::initFromRestart(const data::Aquifers& aquiferSoln)
 {
-    for (auto& aquifer : this->aquifers)
+    this->computeConnectionAreaFraction();
+
+    for (auto& aquifer : this->aquifers) {
         aquifer->initFromRestart(aquiferSoln);
     }
+}
 
 template <typename TypeTag>
 void
@@ -67,6 +73,8 @@ BlackoilAquiferModel<TypeTag>::beginEpisode()
     // SCHEDULE setup in this section it is the beginning of a report step
 
     this->createDynamicAquifers(this->simulator_.episodeIndex());
+
+    this->computeConnectionAreaFraction();
 }
 
 template <typename TypeTag>
@@ -78,9 +86,10 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginTimeStep()
 {
-    for (auto& aquifer : aquifers)
+    for (auto& aquifer : this->aquifers) {
         aquifer->beginTimeStep();
     }
+}
 
 template <typename TypeTag>
 template <class Context>
@@ -90,9 +99,10 @@ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
                                            unsigned spaceIdx,
                                            unsigned timeIdx) const
 {
-    for (auto& aquifer : aquifers)
+    for (auto& aquifer : this->aquifers) {
         aquifer->addToSource(rates, context, spaceIdx, timeIdx);
     }
+}
 
 template <typename TypeTag>
 void
@@ -100,9 +110,10 @@ BlackoilAquiferModel<TypeTag>::addToSource(RateVector& rates,
                                            unsigned globalSpaceIdx,
                                            unsigned timeIdx) const
 {
-    for (auto& aquifer : aquifers)
+    for (auto& aquifer : this->aquifers) {
         aquifer->addToSource(rates, globalSpaceIdx, timeIdx);
     }
+}
 
 template <typename TypeTag>
 void
@@ -113,9 +124,10 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::endTimeStep()
 {
-    for (auto& aquifer : aquifers) {
-        aquifer->endTimeStep();
     using NumAq = AquiferNumerical<TypeTag>;
+
+    for (auto& aquifer : this->aquifers) {
+        aquifer->endTimeStep();
         NumAq* num = dynamic_cast<NumAq*>(aquifer.get());
         if (num)
             this->simulator_.vanguard().grid().comm().barrier();
@@ -159,8 +171,9 @@ template<typename TypeTag>
 data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const
 {
     data::Aquifers data;
-    for (const auto& aqu : this->aquifers)
+    for (const auto& aqu : this->aquifers) {
         data.insert_or_assign(aqu->aquiferID(), aqu->aquiferData());
+    }
 
     return data;
 }
@@ -170,7 +183,7 @@ template<class Serializer>
 void BlackoilAquiferModel<TypeTag>::
 serializeOp(Serializer& serializer)
 {
-    for (auto& aiPtr : aquifers) {
+    for (auto& aiPtr : this->aquifers) {
         auto* ct = dynamic_cast<AquiferCarterTracy<TypeTag>*>(aiPtr.get());
         auto* fetp = dynamic_cast<AquiferFetkovich<TypeTag>*>(aiPtr.get());
         auto* num = dynamic_cast<AquiferNumerical<TypeTag>*>(aiPtr.get());
@@ -304,4 +317,26 @@ void BlackoilAquiferModel<TypeTag>::createDynamicAquifers(const int episode_inde
     }
 }
 
+template <typename TypeTag>
+void BlackoilAquiferModel<TypeTag>::computeConnectionAreaFraction() const
+{
+    auto maxAquID =
+        std::accumulate(this->aquifers.begin(), this->aquifers.end(), 0,
+                        [](const int aquID, const auto& aquifer)
+                        { return std::max(aquID, aquifer->aquiferID()); });
+
+    maxAquID = this->simulator_.vanguard().grid().comm().max(maxAquID);
+
+    auto totalConnArea = std::vector<double>(maxAquID, 0.0);
+    for (const auto& aquifer : this->aquifers) {
+        totalConnArea[aquifer->aquiferID() - 1] += aquifer->totalFaceArea();
+    }
+
+    this->simulator_.vanguard().grid().comm().sum(totalConnArea.data(), maxAquID);
+
+    for (auto& aquifer : this->aquifers) {
+        aquifer->computeFaceAreaFraction(totalConnArea);
+    }
+}
+
 } // namespace Opm