Merge pull request #3474 from bska/restart-load-aquifer

Load Analytic Aquifers From Restart File

ebos/eclproblem.hh

@@ -2379,8 +2379,13 @@ private:
                     this->solventSaturation_[elemIdx] = ssol;
             }
 
-            this->lastRs_[elemIdx] = elemFluidState.Rs();
+            if (! this->lastRs_.empty()) {
-            this->lastRv_[elemIdx] = elemFluidState.Rv();
+                this->lastRs_[elemIdx] = elemFluidState.Rs();
+            }
+
+            if (! this->lastRv_.empty()) {
+                this->lastRv_[elemIdx] = elemFluidState.Rv();
+            }
 
             if constexpr (enablePolymer)
                  this->polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -104,10 +104,10 @@ protected:
     Scalar dimensionless_time_{0};
     Scalar dimensionless_pressure_{0};
 
-    void assignRestartData(const data::AquiferData& /* xaq */) override
+    void assignRestartData(const data::AquiferData& xaq) override
     {
-        throw std::runtime_error {"Restart-based initialization not currently supported "
+        this->fluxValue_ = xaq.volume;
-                                  "for Carter-Tracey analytic aquifers"};
+        this->rhow_ = this->aquct_data_.waterDensity();
     }
 
     std::pair<Scalar, Scalar>
@@ -176,6 +176,10 @@ protected:
 
     inline void calculateAquiferCondition() override
     {
+        if (this->solution_set_from_restart_) {
+            return;
+        }
+
         if (! this->aquct_data_.initial_pressure.has_value()) {
             this->aquct_data_.initial_pressure =
                 this->calculateReservoirEquilibrium();

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -107,6 +107,7 @@ protected:
         }
 
         this->aquifer_pressure_ = xaq.pressure;
+        this->rhow_ = this->aqufetp_data_.waterDensity();
     }
 
     inline Eval dpai(int idx)

opm/simulators/aquifers/AquiferInterface.hpp

@@ -24,8 +24,6 @@
 
 #include <opm/common/utility/numeric/linearInterpolation.hpp>
 #include <opm/parser/eclipse/EclipseState/Aquifer/Aquancon.hpp>
-#include <opm/parser/eclipse/EclipseState/Aquifer/AquiferCT.hpp>
-#include <opm/parser/eclipse/EclipseState/Aquifer/Aquifetp.hpp>
 
 #include <opm/output/data/Aquifer.hpp>
 
@@ -35,6 +33,10 @@
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 
 #include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <numeric>
 #include <unordered_map>
 #include <vector>
 
@@ -168,7 +170,7 @@ protected:
     {
         // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
         if (!this->solution_set_from_restart_) {
-            W_flux_ = 0.;
+            W_flux_ = Scalar{0};
         }
 
         // We next get our connections to the aquifer and initialize these quantities using the initialize_connections
@@ -177,9 +179,9 @@ protected:
         calculateAquiferCondition();
         calculateAquiferConstants();
 
-        pressure_previous_.resize(this->connections_.size(), 0.);
+        pressure_previous_.resize(this->connections_.size(), Scalar{0});
-        pressure_current_.resize(this->connections_.size(), 0.);
+        pressure_current_.resize(this->connections_.size(), Scalar{0});
-        Qai_.resize(this->connections_.size(), 0.0);
+        Qai_.resize(this->connections_.size(), Scalar{0});
     }
 
     inline void
@@ -225,16 +227,18 @@ protected:
     {
         this->cell_depth_.resize(this->size(), this->aquiferDepth());
         this->alphai_.resize(this->size(), 1.0);
-        this->faceArea_connected_.resize(this->size(), 0.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.;
+        Scalar denom_face_areas{0};
         this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
         const auto& gridView = this->ebos_simulator_.vanguard().gridView();
-        for (size_t idx = 0; idx < this->size(); ++idx) {
+        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>();
@@ -245,6 +249,8 @@ protected:
             if ( cell_index < 0 || 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);
         }
@@ -308,12 +314,36 @@ protected:
         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 (size_t idx = 0; idx < this->size(); ++idx) {
+        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)
-                ? // Prevent no connection NaNs due to division by zero
+                ? Scalar{0}
-                0.
                 : this->faceArea_connected_.at(idx) / denom_face_areas;
         }
+
+        if (this->solution_set_from_restart_) {
+            this->rescaleProducedVolume(has_active_connection_on_proc);
+        }
+    }
+
+    void rescaleProducedVolume(const bool has_active_connection_on_proc)
+    {
+        // 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.
+
+        const auto this_area = has_active_connection_on_proc
+            ? std::accumulate(this->alphai_.begin(),
+                              this->alphai_.end(),
+                              Scalar{0})
+            : Scalar{0};
+
+        const auto tot_area = this->ebos_simulator_.vanguard()
+            .grid().comm().sum(this_area);
+
+        this->W_flux_ *= this_area / tot_area;
     }
 
     virtual void assignRestartData(const data::AquiferData& xaq) = 0;
@@ -364,8 +394,8 @@ protected:
         const auto& comm = ebos_simulator_.vanguard().grid().comm();
 
         Scalar vals[2];
-        vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.0);
+        vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
-        vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.0);
+        vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
 
         comm.sum(vals, 2);
 

opm/simulators/aquifers/AquiferNumerical.hpp

@@ -22,9 +22,15 @@
 #define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
 
 #include <opm/output/data/Aquifer.hpp>
+
 #include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
 
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <unordered_map>
 #include <utility>
+#include <vector>
 
 namespace Opm
 {
@@ -54,29 +60,49 @@ public:
                      const std::unordered_map<int, int>& cartesian_to_compressed,
                      const Simulator& ebos_simulator,
                      const int* global_cell)
-        : id_(aquifer.id())
+        : id_             (aquifer.id())
-        , ebos_simulator_(ebos_simulator)
+        , ebos_simulator_ (ebos_simulator)
-        , flux_rate_(0.)
+        , flux_rate_      (0.0)
-        , cumulative_flux_(0.)
+        , cumulative_flux_(0.0)
-        , global_cell_(global_cell)
+        , global_cell_    (global_cell)
-        , init_pressure_(aquifer.numCells(), 0.0)
+        , init_pressure_  (aquifer.numCells(), 0.0)
     {
         this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
 
-        for (size_t idx = 0; idx < aquifer.numCells(); ++idx) {
+        auto aquifer_on_process = false;
+        for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
             const auto* cell = aquifer.getCellPrt(idx);
 
             // Due to parallelisation, the cell might not exist in the current process
             auto search = cartesian_to_compressed.find(cell->global_index);
             if (search != cartesian_to_compressed.end()) {
                 this->cell_to_aquifer_cell_idx_[search->second] = idx;
+                aquifer_on_process = true;
             }
         }
+
+        if (aquifer_on_process) {
+            this->checkConnectsToReservoir();
+        }
     }
 
-    void initFromRestart([[maybe_unused]]const data::Aquifers& aquiferSoln)
+    void initFromRestart(const data::Aquifers& aquiferSoln)
     {
-        // NOT handling Restart for now
+        auto xaqPos = aquiferSoln.find(this->aquiferID());
+        if (xaqPos == aquiferSoln.end())
+            return;
+
+        if (this->connects_to_reservoir_) {
+            this->cumulative_flux_ = xaqPos->second.volume;
+        }
+
+        if (const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
+            aqData != nullptr)
+        {
+            this->init_pressure_ = aqData->initPressure;
+        }
+
+        this->solution_set_from_restart_ = true;
     }
 
     void endTimeStep()
@@ -102,6 +128,10 @@ public:
 
     void initialSolutionApplied()
     {
+        if (this->solution_set_from_restart_) {
+            return;
+        }
+
         this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
         this->flux_rate_ = 0.;
         this->cumulative_flux_ = 0.;
@@ -113,17 +143,41 @@ public:
     }
 
 private:
-    const size_t id_;
+    const std::size_t id_;
     const Simulator& ebos_simulator_;
     double flux_rate_; // aquifer influx rate
     double cumulative_flux_; // cumulative aquifer influx
     const int* global_cell_; // mapping to global index
     std::vector<double> init_pressure_{};
     double pressure_; // aquifer pressure
+    bool solution_set_from_restart_ {false};
+    bool connects_to_reservoir_ {false};
 
     // TODO: maybe unordered_map can also do the work to save memory?
     std::vector<int> cell_to_aquifer_cell_idx_;
 
+    void checkConnectsToReservoir()
+    {
+        ElementContext elem_ctx(this->ebos_simulator_);
+        auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
+                                   this->ebos_simulator_.gridView().template end</*codim=*/0>(),
+            [&elem_ctx, this](const auto& elem) -> bool
+        {
+            elem_ctx.updateStencil(elem);
+
+            const auto cell_index = elem_ctx
+                .globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+
+            return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
+        });
+
+        assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
+                && "Internal error locating numerical aquifer's connecting cell");
+
+        this->connects_to_reservoir_ =
+            elemIt->partitionType() == Dune::InteriorEntity;
+    }
+
     double calculateAquiferPressure() const
     {
         auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
@@ -183,21 +237,25 @@ private:
 
     double calculateAquiferFluxRate() const
     {
-        double aquifer_flux = 0.;
+        double aquifer_flux = 0.0;
 
-        ElementContext  elem_ctx(this->ebos_simulator_);
+        if (! this->connects_to_reservoir_) {
+            return aquifer_flux;
+        }
+
+        ElementContext elem_ctx(this->ebos_simulator_);
         const auto& gridView = this->ebos_simulator_.gridView();
         auto elemIt = gridView.template begin</*codim=*/0>();
         const auto& elemEndIt = gridView.template end</*codim=*/0>();
         for (; elemIt != elemEndIt; ++elemIt) {
-            const auto &elem = *elemIt;
+            const auto& elem = *elemIt;
             if (elem.partitionType() != Dune::InteriorEntity) {
                 continue;
             }
             // elem_ctx.updatePrimaryStencil(elem);
             elem_ctx.updateStencil(elem);
 
-            const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+            const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
             const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
             // we only need the first aquifer cell
             if (idx != 0) {
@@ -206,19 +264,19 @@ private:
             elem_ctx.updateAllIntensiveQuantities();
             elem_ctx.updateAllExtensiveQuantities();
 
-            const size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
+            const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
             // const auto &problem = elem_ctx.problem();
-            const auto &stencil = elem_ctx.stencil(0);
+            const auto& stencil = elem_ctx.stencil(0);
             // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
 
-            for (size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
+            for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
-                const auto &face = stencil.interiorFace(face_idx);
+                const auto& face = stencil.interiorFace(face_idx);
                 // dof index
-                const size_t i = face.interiorIndex();
+                const std::size_t i = face.interiorIndex();
-                const size_t j = face.exteriorIndex();
+                const std::size_t j = face.exteriorIndex();
                 // compressed index
                 // const size_t I = stencil.globalSpaceIndex(i);
-                const size_t J = stencil.globalSpaceIndex(j);
+                const std::size_t J = stencil.globalSpaceIndex(j);
 
                 assert(stencil.globalSpaceIndex(i) == cell_index);
 
@@ -227,11 +285,11 @@ private:
                 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
                     continue;
                 }
-                const auto &exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
+                const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
                 const double water_flux = Toolbox::value(exQuants.volumeFlux(waterPhaseIdx));
 
-                const size_t up_id = water_flux >= 0. ? i : j;
+                const std::size_t up_id = water_flux >= 0.0 ? i : j;
-                const auto &intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
+                const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
                 const double invB = Toolbox::value(intQuantsIn.fluidState().invB(waterPhaseIdx));
                 const double face_area = face.area();
                 aquifer_flux += water_flux * invB * face_area;

opm/simulators/utils/readDeck.cpp

@@ -38,6 +38,7 @@
 #include <opm/output/eclipse/RestartIO.hpp>
 #include <opm/io/eclipse/ERst.hpp>
 #include <opm/io/eclipse/RestartFileView.hpp>
+#include <opm/io/eclipse/rst/aquifer.hpp>
 #include <opm/io/eclipse/rst/state.hpp>
 
 #include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -230,36 +231,70 @@ void readDeck(int rank, std::string& deckFilename, std::unique_ptr<Opm::Deck>& d
                 eclipseState = std::make_unique<Opm::EclipseState>(*deck);
 #endif
             }
-            /*
+
-              For the time being initializing wells and groups from the
-              restart file is not possible, but work is underways and it is
-              included here as a switch.
-            */
             const auto& init_config = eclipseState->getInitConfig();
-            if (init_config.restartRequested() && initFromRestart) {
+            if (init_config.restartRequested()) {
+                // Analytic aquifers must always be loaded from the restart
+                // file in restarted runs and the corresponding keywords
+                // (e.g., AQUANCON and AQUCT) do not exist in the input deck
+                // in this case.  In other words, there's no way to check if
+                // there really are analytic aquifers in the run until we
+                // attempt to read the specifications from the restart file.
+                // If the loader determines that there are no analytic
+                // aquifers, then 'EclipseState::loadRestartAquifers()' does
+                // nothing.
                 const int report_step = init_config.getRestartStep();
                 const auto rst_filename = eclipseState->getIOConfig().getRestartFileName( init_config.getRestartRootName(), report_step, false );
                 auto rst_file = std::make_shared<EclIO::ERst>(rst_filename);
                 auto rst_view = std::make_shared<EclIO::RestartFileView>(std::move(rst_file), report_step);
-                const auto rst_state = Opm::RestartIO::RstState::load(std::move(rst_view));
+
-                if (!schedule)
+                // Note: RstState::load() will just *read* from the grid
-                    schedule = std::make_unique<Opm::Schedule>(*deck, *eclipseState, *parseContext, *errorGuard, python, outputInterval, &rst_state);
+                // structure, and only do so if the case actually includes
-                udqState = std::make_unique<Opm::UDQState>( schedule->operator[](0).udq().params().undefinedValue() );
+                // analytic aquifers.  The pointer to the input grid is just
+                // to allow 'nullptr' to signify "don't load aquifers" in
+                // certain unit tests.  Passing an optional<EclipseGrid> is
+                // too expensive however since doing so will create a copy
+                // of the grid inside the optional<>.
+                const auto rst_state = RestartIO::RstState::
+                    load(std::move(rst_view), &eclipseState->getInputGrid());
+
+                eclipseState->loadRestartAquifers(rst_state.aquifers);
+
+                // For the time being initializing wells and groups from the
+                // restart file is not possible.  Work is underway and the
+                // ability is included here contingent on user-level switch
+                // 'initFromRestart' (i.e., setting "--sched-restart=false"
+                // as a command line invocation parameter).
+                const auto* init_state = initFromRestart ? &rst_state : nullptr;
+                if (!schedule) {
+                    schedule = std::make_unique<Schedule>(*deck, *eclipseState,
+                                                          *parseContext, *errorGuard,
+                                                          python, outputInterval, init_state);
+                }
+
+                udqState = std::make_unique<UDQState>((*schedule)[0].udq().params().undefinedValue());
                 udqState->load_rst(rst_state);
             }
             else {
-                if (!schedule)
+                if (!schedule) {
-                    schedule = std::make_unique<Opm::Schedule>(*deck, *eclipseState, *parseContext, *errorGuard, python);
+                    schedule = std::make_unique<Schedule>(*deck, *eclipseState,
-                udqState = std::make_unique<Opm::UDQState>( schedule->operator[](0).udq().params().undefinedValue() );
+                                                          *parseContext, *errorGuard,
+                                                          python);
+                }
+
+                udqState = std::make_unique<UDQState>((*schedule)[0].udq().params().undefinedValue());
             }
 
-            if (Opm::OpmLog::hasBackend("STDOUT_LOGGER")) // loggers might not be set up!
+
-            {
+            if (Opm::OpmLog::hasBackend("STDOUT_LOGGER")) {
-                setupMessageLimiter(schedule->operator[](0).message_limits(), "STDOUT_LOGGER");
+                // loggers might not be set up!
+                setupMessageLimiter((*schedule)[0].message_limits(), "STDOUT_LOGGER");
             }
-            if (!summaryConfig)
+
-                summaryConfig = std::make_unique<Opm::SummaryConfig>(*deck, *schedule, eclipseState->fieldProps(), 
+            if (!summaryConfig) {
+                summaryConfig = std::make_unique<Opm::SummaryConfig>(*deck, *schedule, eclipseState->fieldProps(),
                                                                      eclipseState->aquifer(), *parseContext, *errorGuard);
+            }
 
             Opm::checkConsistentArrayDimensions(*eclipseState, *schedule, *parseContext, *errorGuard);
         }