Synchronise Face-Area Fractions Between All Processes

We need a global view of face-area fractions if aquifer connections
happen to be shared between processes.  Add a new helper function,

    BlackoilAquiferModel::computeConnectionAreaFraction()

that performs a collective operation to compute the total face areas
and then defers to the local aquifer objects to compute their face
area fractions.

While here, also split the initialisation of analytic aquifers into
two parts, one for the face area and connection mappings, and one
for the connection depths.  Run the former as part of the object
constructor and the latter as part of 'initQuantities()'.  This
ensures that we can computeConnectionAreaFraction() for all analytic
aquifers before assigning solution quantities from the restart file.

opm/simulators/aquifers/AquiferAnalytical.hpp

@@ -87,17 +87,45 @@ public:
                                           BlackoilIndices::numPhases>;
 
     // Constructor
-    AquiferAnalytical(int aqID,
-                     const std::vector<Aquancon::AquancCell>& connections,
-                     const Simulator& ebosSimulator)
+    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;
     }
 
@@ -226,19 +255,17 @@ protected:
     void initQuantities()
     {
         // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
-        if (!this->solution_set_from_restart_) {
+        if (! this->solution_set_from_restart_) {
             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);
-            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)
+            if (cell_index < 0) {
                 continue;
+            }
 
-            has_active_connection_on_proc = true;
+            auto elemIt = gridView.template begin</*codim=*/ 0>();
+            std::advance(elemIt, cell_index);
+
+            // The global_index is not part of this grid
+            if (elemIt->partitionType() != Dune::InteriorEntity) {
+                continue;
+            }
 
             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;
@@ -325,51 +352,41 @@ protected:
                     break;
                 default:
                     OPM_THROW(std::logic_error,
-                        "Internal error in initialization of aquifer.");
+                              "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,8 +60,7 @@ public:
         , aquifer_data_            (aquifer)
         , connection_flux_         (connections_.size(), Eval{0})
     {
-        const auto connected_face_area = this->initializeConnections();
-        this->area_fraction_ = this->computeFaceAreaFraction(connected_face_area);
+        this->total_face_area_ = this->initializeConnections();
     }
 
     static AquiferConstantFlux serializationTestObject(const Simulator& ebos_simulator)
@@ -72,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;
@@ -150,6 +164,7 @@ private:
     std::vector<int> cellToConnectionIdx_{};
     double flux_rate_{};
     double cumulative_flux_{};
+    double total_face_area_{0.0};
     double area_fraction_{1.0};
 
     double initializeConnections()

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,16 +45,22 @@ 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>
@@ -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,8 +86,9 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::beginTimeStep()
 {
-    for (auto& aquifer : aquifers)
+    for (auto& aquifer : this->aquifers) {
         aquifer->beginTimeStep();
+    }
 }
 
 template <typename TypeTag>
@@ -90,8 +99,9 @@ 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>
@@ -100,8 +110,9 @@ 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>
@@ -113,9 +124,10 @@ template <typename TypeTag>
 void
 BlackoilAquiferModel<TypeTag>::endTimeStep()
 {
-    for (auto& aquifer : aquifers) {
+    using NumAq = AquiferNumerical<TypeTag>;
+
+    for (auto& aquifer : this->aquifers) {
         aquifer->endTimeStep();
-        using NumAq = AquiferNumerical<TypeTag>;
         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