diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index a91cea4ea..b83443363 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -2379,8 +2379,13 @@ private: this->solventSaturation_[elemIdx] = ssol; } - this->lastRs_[elemIdx] = elemFluidState.Rs(); - this->lastRv_[elemIdx] = elemFluidState.Rv(); + if (! this->lastRs_.empty()) { + this->lastRs_[elemIdx] = elemFluidState.Rs(); + } + + if (! this->lastRv_.empty()) { + this->lastRv_[elemIdx] = elemFluidState.Rv(); + } if constexpr (enablePolymer) this->polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx); diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index 2f53c0eb1..1af0db128 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/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 " - "for Carter-Tracey analytic aquifers"}; + this->fluxValue_ = xaq.volume; + this->rhow_ = this->aquct_data_.waterDensity(); } std::pair @@ -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(); diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index ac53e8ab6..68bb7e7be 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/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) diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 247e1570e..9f0a0104f 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -24,8 +24,6 @@ #include #include -#include -#include #include @@ -35,6 +33,10 @@ #include #include +#include +#include +#include +#include #include #include @@ -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_current_.resize(this->connections_.size(), 0.); - Qai_.resize(this->connections_.size(), 0.0); + pressure_previous_.resize(this->connections_.size(), Scalar{0}); + pressure_current_.resize(this->connections_.size(), Scalar{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(); @@ -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::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 - 0. + ? Scalar{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[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.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(), Scalar{0}); comm.sum(vals, 2); diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp index 06af4689e..087793f1a 100644 --- a/opm/simulators/aquifers/AquiferNumerical.hpp +++ b/opm/simulators/aquifers/AquiferNumerical.hpp @@ -22,9 +22,15 @@ #define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED #include + #include +#include +#include +#include +#include #include +#include namespace Opm { @@ -54,29 +60,49 @@ public: const std::unordered_map& cartesian_to_compressed, const Simulator& ebos_simulator, const int* global_cell) - : id_(aquifer.id()) - , ebos_simulator_(ebos_simulator) - , flux_rate_(0.) - , cumulative_flux_(0.) - , global_cell_(global_cell) - , init_pressure_(aquifer.numCells(), 0.0) + : id_ (aquifer.id()) + , ebos_simulator_ (ebos_simulator) + , flux_rate_ (0.0) + , cumulative_flux_(0.0) + , global_cell_ (global_cell) + , 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(); + 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 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 cell_to_aquifer_cell_idx_; + void checkConnectsToReservoir() + { + ElementContext elem_ctx(this->ebos_simulator_); + auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin(), + this->ebos_simulator_.gridView().template end(), + [&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()) + && "Internal error locating numerical aquifer's connecting cell"); + + this->connects_to_reservoir_ = + elemIt->partitionType() == Dune::InteriorEntity; + } + double calculateAquiferPressure() const { auto capture = std::vector(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(); const auto& elemEndIt = gridView.template end(); 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) { - const auto &face = stencil.interiorFace(face_idx); + for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) { + const auto& face = stencil.interiorFace(face_idx); // dof index - const size_t i = face.interiorIndex(); - const size_t j = face.exteriorIndex(); + const std::size_t i = face.interiorIndex(); + 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 auto &intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0); + const std::size_t up_id = water_flux >= 0.0 ? i : j; + 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; diff --git a/opm/simulators/utils/readDeck.cpp b/opm/simulators/utils/readDeck.cpp index 148087c40..2206457bd 100644 --- a/opm/simulators/utils/readDeck.cpp +++ b/opm/simulators/utils/readDeck.cpp @@ -38,6 +38,7 @@ #include #include #include +#include #include #include @@ -230,36 +231,70 @@ void readDeck(int rank, std::string& deckFilename, std::unique_ptr& d eclipseState = std::make_unique(*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(rst_filename); auto rst_view = std::make_shared(std::move(rst_file), report_step); - const auto rst_state = Opm::RestartIO::RstState::load(std::move(rst_view)); - if (!schedule) - schedule = std::make_unique(*deck, *eclipseState, *parseContext, *errorGuard, python, outputInterval, &rst_state); - udqState = std::make_unique( schedule->operator[](0).udq().params().undefinedValue() ); + + // Note: RstState::load() will just *read* from the grid + // structure, and only do so if the case actually includes + // 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 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(*deck, *eclipseState, + *parseContext, *errorGuard, + python, outputInterval, init_state); + } + + udqState = std::make_unique((*schedule)[0].udq().params().undefinedValue()); udqState->load_rst(rst_state); } else { - if (!schedule) - schedule = std::make_unique(*deck, *eclipseState, *parseContext, *errorGuard, python); - udqState = std::make_unique( schedule->operator[](0).udq().params().undefinedValue() ); + if (!schedule) { + schedule = std::make_unique(*deck, *eclipseState, + *parseContext, *errorGuard, + python); + } + + udqState = std::make_unique((*schedule)[0].udq().params().undefinedValue()); } - if (Opm::OpmLog::hasBackend("STDOUT_LOGGER")) // loggers might not be set up! - { - setupMessageLimiter(schedule->operator[](0).message_limits(), "STDOUT_LOGGER"); + + if (Opm::OpmLog::hasBackend("STDOUT_LOGGER")) { + // loggers might not be set up! + setupMessageLimiter((*schedule)[0].message_limits(), "STDOUT_LOGGER"); } - if (!summaryConfig) - summaryConfig = std::make_unique(*deck, *schedule, eclipseState->fieldProps(), + + if (!summaryConfig) { + summaryConfig = std::make_unique(*deck, *schedule, eclipseState->fieldProps(), eclipseState->aquifer(), *parseContext, *errorGuard); + } Opm::checkConsistentArrayDimensions(*eclipseState, *schedule, *parseContext, *errorGuard); }