diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4d8cd0c39d..1e18ed7bc1 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -167,6 +167,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/BlackoilWellModelGasLift.cpp opm/simulators/wells/BlackoilWellModelGeneric.cpp opm/simulators/wells/BlackoilWellModelGuideRates.cpp + opm/simulators/wells/BlackoilWellModelNldd.cpp opm/simulators/wells/BlackoilWellModelRestart.cpp opm/simulators/wells/BlackoilWellModelWBP.cpp opm/simulators/wells/ConnFiltrateData.cpp @@ -972,6 +973,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/BlackoilWellModelGasLift_impl.hpp opm/simulators/wells/BlackoilWellModelGeneric.hpp opm/simulators/wells/BlackoilWellModelGuideRates.hpp + opm/simulators/wells/BlackoilWellModelNldd.hpp + opm/simulators/wells/BlackoilWellModelNldd_impl.hpp opm/simulators/wells/BlackoilWellModelRestart.hpp opm/simulators/wells/BlackoilWellModelWBP.hpp opm/simulators/wells/ConnectionIndexMap.hpp diff --git a/opm/simulators/flow/BlackoilModelNldd.hpp b/opm/simulators/flow/BlackoilModelNldd.hpp index ea00499cb0..4850e7a143 100644 --- a/opm/simulators/flow/BlackoilModelNldd.hpp +++ b/opm/simulators/flow/BlackoilModelNldd.hpp @@ -50,6 +50,8 @@ #include +#include + #include #include @@ -64,6 +66,7 @@ #include #include #include +#include #include #include #include @@ -76,7 +79,8 @@ template class BlackoilModel; /// A NLDD implementation for three-phase black oil. template -class BlackoilModelNldd { +class BlackoilModelNldd +{ public: using ElementContext = GetPropType; using FluidSystem = GetPropType; @@ -98,11 +102,16 @@ class BlackoilModelNldd { //! \param param param Model parameters //! \param compNames Names of the solution components explicit BlackoilModelNldd(BlackoilModel& model) - : model_(model), rank_(model_.simulator().vanguard().grid().comm().rank()) + : model_(model) + , wellModel_(model.wellModel()) + , rank_(model_.simulator().vanguard().grid().comm().rank()) { // Create partitions. const auto& [partition_vector, num_domains] = this->partitionCells(); + // Set nldd handler in main well model + model.wellModel().setNlddAdapter(&wellModel_); + // Scan through partitioning to get correct size for each. std::vector sizes(num_domains, 0); for (const auto& p : partition_vector) { @@ -179,7 +188,7 @@ class BlackoilModelNldd { void prepareStep() { // Setup domain->well mapping. - model_.wellModel().setupDomains(domains_); + wellModel_.setupDomains(domains_); } //! \brief Do one non-linear NLDD iteration. @@ -352,7 +361,6 @@ class BlackoilModelNldd { } private: - //! \brief Solve the equation system for a single domain. std::pair solveDomain(const Domain& domain, @@ -380,9 +388,9 @@ class BlackoilModelNldd { modelSimulator.model().newtonMethod().setIterationIndex(iter); // TODO: we should have a beginIterationLocal function() // only handling the well model for now - modelSimulator.problem().wellModel().assembleDomain(modelSimulator.model().newtonMethod().numIterations(), - modelSimulator.timeStepSize(), - domain); + wellModel_.assemble(modelSimulator.model().newtonMethod().numIterations(), + modelSimulator.timeStepSize(), + domain); // Assemble reservoir locally. report += this->assembleReservoirDomain(domain); report.assemble_time += detailTimer.stop(); @@ -445,9 +453,9 @@ class BlackoilModelNldd { // TODO: we should have a beginIterationLocal function() // only handling the well model for now // Assemble reservoir locally. - modelSimulator.problem().wellModel().assembleDomain(modelSimulator.model().newtonMethod().numIterations(), - modelSimulator.timeStepSize(), - domain); + wellModel_.assemble(modelSimulator.model().newtonMethod().numIterations(), + modelSimulator.timeStepSize(), + domain); report += this->assembleReservoirDomain(domain); report.assemble_time += detailTimer.stop(); @@ -736,7 +744,7 @@ class BlackoilModelNldd { logger, B_avg, residual_norms); - report += model_.wellModel().getDomainWellConvergence(domain, B_avg, logger); + report += wellModel_.getWellConvergence(domain, B_avg, logger); return report; } @@ -815,7 +823,7 @@ class BlackoilModelNldd { const SimulatorTimerInterface& timer, const Domain& domain) { - auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain.index); + auto initial_local_well_primary_vars = wellModel_.getPrimaryVarsDomain(domain.index); auto initial_local_solution = Details::extractVector(solution, domain.cells); auto res = solveDomain(domain, timer, logger, iteration, false); local_report = res.first; @@ -825,7 +833,7 @@ class BlackoilModelNldd { Details::setGlobal(initial_local_solution, domain.cells, solution); model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain); } else { - model_.wellModel().setPrimaryVarsDomain(domain.index, initial_local_well_primary_vars); + wellModel_.setPrimaryVarsDomain(domain.index, initial_local_well_primary_vars); Details::setGlobal(initial_local_solution, domain.cells, solution); model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain); } @@ -840,7 +848,7 @@ class BlackoilModelNldd { const SimulatorTimerInterface& timer, const Domain& domain) { - auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain.index); + auto initial_local_well_primary_vars = wellModel_.getPrimaryVarsDomain(domain.index); auto initial_local_solution = Details::extractVector(solution, domain.cells); auto res = solveDomain(domain, timer, logger, iteration, true); local_report = res.first; @@ -881,7 +889,7 @@ class BlackoilModelNldd { auto local_solution = Details::extractVector(solution, domain.cells); Details::setGlobal(local_solution, domain.cells, locally_solved); } else { - model_.wellModel().setPrimaryVarsDomain(domain.index, initial_local_well_primary_vars); + wellModel_.setPrimaryVarsDomain(domain.index, initial_local_well_primary_vars); Details::setGlobal(initial_local_solution, domain.cells, solution); model_.simulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain); } @@ -996,6 +1004,7 @@ class BlackoilModelNldd { } BlackoilModel& model_; //!< Reference to model + BlackoilWellModelNldd wellModel_; //!< NLDD well model adapter std::vector domains_; //!< Vector of subdomains std::vector> domain_matrices_; //!< Vector of matrix operator for each subdomain std::vector domain_linsolvers_; //!< Vector of linear solvers for each domain diff --git a/opm/simulators/flow/SubDomain.hpp b/opm/simulators/flow/SubDomain.hpp index 27dc0dd354..64a616cd87 100644 --- a/opm/simulators/flow/SubDomain.hpp +++ b/opm/simulators/flow/SubDomain.hpp @@ -57,8 +57,7 @@ namespace Opm /// Representing a part of a grid, in a way suitable for performing /// local solves. - template - struct SubDomain + struct SubDomainIndices { // The index of a subdomain is arbitrary, but can be used by the // solvers to keep track of well locations etc. @@ -72,10 +71,21 @@ namespace Opm std::vector interior; // Enables subdomain solves and linearization using the generic linearization // approach (i.e. FvBaseLinearizer as opposed to TpfaLinearizer). + SubDomainIndices(const int i, std::vector&& c, std::vector&& in) + : index(i), cells(std::move(c)), interior(std::move(in)) + {} + }; + + /// Representing a part of a grid, in a way suitable for performing + /// local solves. + template + struct SubDomain : public SubDomainIndices + { Dune::SubGridPart view; // Constructor that moves from its argument. SubDomain(const int i, std::vector&& c, std::vector&& in, Dune::SubGridPart&& v) - : index(i), cells(std::move(c)), interior(std::move(in)), view(std::move(v)) + : SubDomainIndices(i, std::move(c), std::move(in)) + , view(std::move(v)) {} }; diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp index 82d981b749..d7a9deac9e 100644 --- a/opm/simulators/linalg/WellOperators.hpp +++ b/opm/simulators/linalg/WellOperators.hpp @@ -188,7 +188,8 @@ class DomainWellModelAsLinearOperator : public WellModelAsLinearOperatorwellMod_) { if (this->wellMod_.well_domain().at(well->name()) == domainIndex_) { - this->applySingleWell(x, y, well, this->wellMod_.well_local_cells()[well_index]); + this->applySingleWell(x, y, well, + this->wellMod_.well_local_cells()[well_index]); } ++well_index; } @@ -199,7 +200,10 @@ class DomainWellModelAsLinearOperator : public WellModelAsLinearOperatorwellMod_.addWellPressureEquationsDomain(jacobian, weights, use_well_weights, domainIndex_); + this->wellMod_.addWellPressureEquationsDomain(jacobian, + weights, + use_well_weights, + domainIndex_); } private: diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 2f957ed91f..6e266bb13f 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -30,8 +30,6 @@ #include -#include - #include #include #include @@ -41,7 +39,6 @@ #include #include -#include #include @@ -81,6 +78,9 @@ namespace Opm { +template class BlackoilWellModelNldd; +template class SparseTable; + #if COMPILE_GPU_BRIDGE template class WellContributions; #endif @@ -132,8 +132,6 @@ template class WellContributions; using AverageRegionalPressureType = RegionAverageCalculator:: AverageRegionalPressure >; - using Domain = SubDomain; - explicit BlackoilWellModel(Simulator& simulator); void init(); @@ -248,11 +246,6 @@ template class WellContributions; // Check if well equations is converged. ConvergenceReport getWellConvergence(const std::vector& B_avg, const bool checkWellGroupControls = false) const; - // Check if well equations are converged locally. - ConvergenceReport getDomainWellConvergence(const Domain& domain, - const std::vector& B_avg, - DeferredLogger& local_deferredLogger) const; - const SimulatorReportSingle& lastReport() const; void addWellContributions(SparseMatrixAdapter& jacobian) const; @@ -274,7 +267,6 @@ template class WellContributions; // at the beginning of each time step (Not report step) void prepareTimeStep(DeferredLogger& deferred_logger); void initPrimaryVariablesEvaluation() const; - void initPrimaryVariablesEvaluationDomain(const Domain& domain) const; std::pair updateWellControls(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false); @@ -292,15 +284,23 @@ template class WellContributions; using PressureMatrix = Dune::BCRSMatrix>; - void addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const; - - void addWellPressureEquationsDomain([[maybe_unused]] PressureMatrix& jacobian, - [[maybe_unused]] const BVector& weights, - [[maybe_unused]] const bool use_well_weights, - [[maybe_unused]] const int domainIndex) const; - - + void addWellPressureEquations(PressureMatrix& jacobian, + const BVector& weights, + const bool use_well_weights) const; void addWellPressureEquationsStruct(PressureMatrix& jacobian) const; + void addWellPressureEquationsDomain(PressureMatrix& jacobian, + const BVector& weights, + const bool use_well_weights, + const int domainIndex) const + { + if (!nldd_) { + OPM_THROW(std::logic_error, "Attempt to access NLDD data without a NLDD solver"); + } + return nldd_->addWellPressureEquations(jacobian, + weights, + use_well_weights, + domainIndex); + } /// \brief Get list of local nonshut wells const std::vector& localNonshutWells() const @@ -308,24 +308,32 @@ template class WellContributions; return well_container_; } - // prototype for assemble function for ASPIN solveLocal() - // will try to merge back to assemble() when done prototyping - void assembleDomain(const int iterationIdx, - const double dt, - const Domain& domain); - void updateWellControlsDomain(DeferredLogger& deferred_logger, const Domain& domain); - - void setupDomains(const std::vector& domains); - const SparseTable& well_local_cells() const { - return well_local_cells_; + if (!nldd_) { + OPM_THROW(std::logic_error, "Attempt to access NLDD data without a NLDD solver"); + } + return nldd_->well_local_cells(); } + + const std::map& well_domain() const + { + if (!nldd_) { + OPM_THROW(std::logic_error, "Attempt to access NLDD data without a NLDD solver"); + } + + return nldd_->well_domain(); + } + auto begin() const { return well_container_.begin(); } auto end() const { return well_container_.end(); } bool empty() const { return well_container_.empty(); } - bool addMatrixContributions() const { return param_.matrix_add_well_contributions_; } + bool addMatrixContributions() const + { return param_.matrix_add_well_contributions_; } + + int numStrictIterations() const + { return param_.strict_outer_iter_wells_; } int compressedIndexForInterior(int cartesian_cell_idx) const override { @@ -339,7 +347,16 @@ template class WellContributions; // using the solution x to recover the solution xw for wells and applying // xw to update Well State void recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, - const Domain& domain); + const int domainIdx); + + const Grid& grid() const + { return simulator_.vanguard().grid(); } + + const Simulator& simulator() const + { return simulator_; } + + void setNlddAdapter(BlackoilWellModelNldd* mod) + { nldd_ = mod; } protected: Simulator& simulator_; @@ -385,12 +402,6 @@ template class WellContributions; std::vector B_avg_{}; - // Store the local index of the wells perforated cells in the domain, if using subdomains - SparseTable well_local_cells_; - - const Grid& grid() const - { return simulator_.vanguard().grid(); } - const EquilGrid& equilGrid() const { return simulator_.vanguard().equilGrid(); } @@ -464,7 +475,6 @@ template class WellContributions; int reportStepIndex() const; void assembleWellEq(const double dt, DeferredLogger& deferred_logger); - void assembleWellEqDomain(const double dt, const Domain& domain, DeferredLogger& deferred_logger); void prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger); @@ -495,6 +505,7 @@ template class WellContributions; BlackoilWellModel(Simulator& simulator, const PhaseUsage& pu); BlackoilWellModelGasLift gaslift_; + BlackoilWellModelNldd* nldd_ = nullptr; //!< NLDD well model adapter (not owned) // These members are used to avoid reallocation in specific functions // instead of using local variables. diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index c19d664f7d..dae024afd1 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -2017,35 +2017,6 @@ logPrimaryVars() const OpmLog::debug(os.str()); } -template -std::vector -BlackoilWellModelGeneric:: -getPrimaryVarsDomain(const int domainIdx) const -{ - std::vector ret; - for (const auto& well : this->well_container_generic_) { - if (this->well_domain_.at(well->name()) == domainIdx) { - const auto& pv = well->getPrimaryVars(); - ret.insert(ret.end(), pv.begin(), pv.end()); - } - } - return ret; -} - -template -void BlackoilWellModelGeneric:: -setPrimaryVarsDomain(const int domainIdx, const std::vector& vars) -{ - std::size_t offset = 0; - for (auto& well : this->well_container_generic_) { - if (this->well_domain_.at(well->name()) == domainIdx) { - int num_pri_vars = well->setPrimaryVars(vars.begin() + offset); - offset += num_pri_vars; - } - } - assert(offset == vars.size()); -} - template void BlackoilWellModelGeneric:: reportGroupSwitching(DeferredLogger& local_deferredLogger) const diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index 1eb7c44371..2c2dd4b455 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -118,7 +118,12 @@ class BlackoilWellModelGeneric // whether there exists any multisegment well open on this process bool anyMSWellOpenLocal() const; - const std::vector& eclWells() const { return wells_ecl_; } + const std::vector& eclWells() const + { return wells_ecl_; } + + bool terminalOutput() const + { return terminal_output_; } + const Well& getWellEcl(const std::string& well_name) const; std::vector getLocalWells(const int timeStepIdx) const; const Schedule& schedule() const { return schedule_; } @@ -127,6 +132,9 @@ class BlackoilWellModelGeneric std::vector*> genericWells() const { return {well_container_generic_.begin(), well_container_generic_.end()}; } + std::vector*> genericWells() + { return well_container_generic_; } + /* Immutable version of the currently active wellstate. */ @@ -217,11 +225,6 @@ class BlackoilWellModelGeneric const std::map& wellOpenTimes() const { return well_open_times_; } const std::map& wellCloseTimes() const { return well_close_times_; } - const std::map& well_domain() const - { - return well_domain_; - } - std::vector getCellsForConnections(const Well& well) const; bool reportStepStarts() const { return report_step_starts_; } @@ -236,8 +239,6 @@ class BlackoilWellModelGeneric bool wasDynamicallyShutThisTimeStep(const std::string& well_name) const; void logPrimaryVars() const; - std::vector getPrimaryVarsDomain(const int domainIdx) const; - void setPrimaryVarsDomain(const int domainIdx, const std::vector& vars); template void serializeOp(Serializer& serializer) @@ -537,9 +538,6 @@ class BlackoilWellModelGeneric // Store map of group name and close offending well for output std::map> closed_offending_wells_; - // Keep track of the domain of each well, if using subdomains. - std::map well_domain_; - private: WellInterfaceGeneric* getGenWell(const std::string& well_name); diff --git a/opm/simulators/wells/BlackoilWellModelNldd.cpp b/opm/simulators/wells/BlackoilWellModelNldd.cpp new file mode 100644 index 0000000000..39c8b93bf7 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNldd.cpp @@ -0,0 +1,170 @@ +/* + Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics. + Copyright 2016 - 2018 Equinor ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 Norce AS + + 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 . +*/ + +#include +#include + +namespace Opm { + +template +std::vector +BlackoilWellModelNlddGeneric:: +getPrimaryVarsDomain(const int domainIdx) const +{ + std::vector ret; + for (const auto& well : genWellModel_.genericWells()) { + if (this->well_domain_.at(well->name()) == domainIdx) { + const auto& pv = well->getPrimaryVars(); + ret.insert(ret.end(), pv.begin(), pv.end()); + } + } + return ret; +} + +template +void +BlackoilWellModelNlddGeneric:: +setPrimaryVarsDomain(const int domainIdx, const std::vector& vars) +{ + std::size_t offset = 0; + for (const auto& well : genWellModel_.genericWells()) { + if (this->well_domain_.at(well->name()) == domainIdx) { + int num_pri_vars = well->setPrimaryVars(vars.begin() + offset); + offset += num_pri_vars; + } + } + assert(offset == vars.size()); +} + +template +void +BlackoilWellModelNlddGeneric:: +findWellDomains(const std::vector& domains) +{ + // TODO: This loop nest may be slow for very large numbers of + // domains and wells, but that has not been observed on tests so + // far. Using the partition vector instead would be faster if we + // need to change. + for (const auto& wellPtr : genWellModel_.genericWells()) { + const int first_well_cell = wellPtr->cells().front(); + for (const auto& domain : domains) { + auto cell_present = [domain](const auto cell) + { + return std::binary_search(domain->cells.begin(), + domain->cells.end(), cell); + }; + + if (cell_present(first_well_cell)) { + // Assuming that if the first well cell is found in a domain, + // then all of that well's cells are in that same domain. + this->well_domain_[wellPtr->name()] = domain->index; + + // Verify that all of that well's cells are in that same domain. + for (int well_cell : wellPtr->cells()) { + if (! cell_present(well_cell)) { + OPM_THROW(std::runtime_error, + fmt::format("Well '{}' found on multiple domains.", + wellPtr->name())); + } + } + } + } + } +} + +template +void +BlackoilWellModelNlddGeneric:: +logDomains() const +{ + // Write well/domain info to the DBG file. + const int rank = genWellModel_.comm().rank(); + DeferredLogger local_log; + if (!this->well_domain_.empty()) { + std::ostringstream os; + os << "Well name Rank Domain\n"; + for (const auto& [wname, domain] : this->well_domain_) { + os << wname << std::setw(19 - wname.size()) << rank << std::setw(12) << domain << '\n'; + } + local_log.debug(os.str()); + } + auto global_log = gatherDeferredLogger(local_log, genWellModel_.comm()); + if (genWellModel_.terminalOutput()) { + global_log.logMessages(); + } +} + +template +void +BlackoilWellModelNlddGeneric:: +calcLocalIndices(const std::vector& domains) +{ + well_local_cells_.clear(); + well_local_cells_.reserve(genWellModel_.genericWells().size(), 10); + std::vector local_cells; + for (const auto& well : genWellModel_.genericWells()) { + const auto& global_cells = well->cells(); + const int domain_index = this->well_domain_.at(well->name()); + const auto& domain_cells = domains[domain_index]->cells; + local_cells.resize(global_cells.size()); + + // find the local cell index for each well cell in the domain + // assume domain_cells is sorted + for (size_t i = 0; i < global_cells.size(); ++i) { + auto it = std::lower_bound(domain_cells.begin(), domain_cells.end(), global_cells[i]); + if (it != domain_cells.end() && *it == global_cells[i]) { + local_cells[i] = std::distance(domain_cells.begin(), it); + } else { + OPM_THROW(std::runtime_error, fmt::format("Cell {} not found in domain {}", + global_cells[i], domain_index)); + } + } + well_local_cells_.appendRow(local_cells.begin(), local_cells.end()); + } +} + +template +void +BlackoilWellModelNlddGeneric:: +calcDomains(const std::vector& domains) +{ + const Opm::Parallel::Communication& comm = genWellModel_.comm(); + + OPM_BEGIN_PARALLEL_TRY_CATCH(); + this->findWellDomains(domains); + OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::setupDomains(): " + "well found on multiple domains.", comm); + + // Write well/domain info to the DBG file. + this->logDomains(); + + // Pre-calculate the local cell indices for each well + this->calcLocalIndices(domains); +} + +template class BlackoilWellModelNlddGeneric; + +#if FLOW_INSTANTIATE_FLOAT +template class BlackoilWellModelNlddGeneric; +#endif + +} // namespace Opm diff --git a/opm/simulators/wells/BlackoilWellModelNldd.hpp b/opm/simulators/wells/BlackoilWellModelNldd.hpp new file mode 100644 index 0000000000..fe0025cf85 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNldd.hpp @@ -0,0 +1,141 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 IRIS AS + + 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 . +*/ +#ifndef OPM_BLACKOILWELLMODEL_NLDD_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_NLDD_HEADER_INCLUDED + +#include + +#include + +#include + +#include + +#include +#include +#include + +namespace Opm { + +template +class BlackoilWellModelNlddGeneric +{ +public: + std::vector getPrimaryVarsDomain(const int domainIdx) const; + void setPrimaryVarsDomain(const int domainIdx, const std::vector& vars); + + const SparseTable& well_local_cells() const + { return well_local_cells_; } + + const std::map& well_domain() const + { return well_domain_; } + +protected: + BlackoilWellModelNlddGeneric(BlackoilWellModelGeneric& model) + : genWellModel_(model) + {} + + void calcDomains(const std::vector& domains); + +private: + void logDomains() const; + + void findWellDomains(const std::vector& domains); + + void calcLocalIndices(const std::vector& domains); + + BlackoilWellModelGeneric& genWellModel_; + + // Keep track of the domain of each well + std::map well_domain_{}; + + // Store the local index of the wells perforated cells in the domain + SparseTable well_local_cells_{}; +}; + +/// Class for handling the blackoil well model in a NLDD solver. +template +class BlackoilWellModelNldd : + public BlackoilWellModelNlddGeneric> +{ +public: + // --------- Types --------- + using Grid = GetPropType; + using Scalar = GetPropType; + using PressureMatrix = typename BlackoilWellModel::PressureMatrix; + using BVector = typename BlackoilWellModel::BVector; + + using Domain = SubDomain; + + BlackoilWellModelNldd(BlackoilWellModel& model) + : BlackoilWellModelNlddGeneric(model) + , wellModel_(model) + {} + + void initPrimaryVariablesEvaluation(const Domain& domain) const; + + void addWellPressureEquations(PressureMatrix& jacobian, + const BVector& weights, + const bool use_well_weights, + const int domainIndex) const; + + // prototype for assemble function for ASPIN solveLocal() + // will try to merge back to assemble() when done prototyping + void assemble(const int iterationIdx, + const double dt, + const Domain& domain); + + void updateWellControls(DeferredLogger& deferred_logger, + const Domain& domain); + + void setupDomains(const std::vector& domains); + + // Check if well equations are converged locally. + ConvergenceReport getWellConvergence(const Domain& domain, + const std::vector& B_avg, + DeferredLogger& local_deferredLogger) const; + + // using the solution x to recover the solution xw for wells and applying + // xw to update Well State + void recoverWellSolutionAndUpdateWellState(const BVector& x, + const int domainIdx); + +private: + BlackoilWellModel& wellModel_; + + void assembleWellEq(const double dt, + const Domain& domain, + DeferredLogger& deferred_logger); + + // These members are used to avoid reallocation in specific functions + // instead of using local variables. + // Their state is not relevant between function calls, so they can + // (and must) be mutable, as the functions using them are const. + mutable BVector x_local_; + +}; + +} // namespace Opm + +#include "BlackoilWellModelNldd_impl.hpp" + +#endif diff --git a/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp b/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp new file mode 100644 index 0000000000..6ae4a5d9bb --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp @@ -0,0 +1,227 @@ +/* + Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics. + Copyright 2016 - 2018 Equinor ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 Norce AS + + 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 . +*/ + +#ifndef OPM_BLACKOILWELLMODEL_NLDD_IMPL_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_NLDD_IMPL_HEADER_INCLUDED + +// Improve IDE experience +#ifndef OPM_BLACKOILWELLMODEL_NLDD_HEADER_INCLUDED +#include +#include +#endif + +#include + +namespace Opm { + +template +void +BlackoilWellModelNldd:: +assemble(const int /*iterationIdx*/, + const double dt, + const Domain& domain) +{ + // We assume that calculateExplicitQuantities() and + // prepareTimeStep() have been called already for the entire + // well model, so we do not need to do it here (when + // iterationIdx is 0). + + DeferredLogger local_deferredLogger; + this->updateWellControls(local_deferredLogger, domain); + this->initPrimaryVariablesEvaluation(domain); + this->assembleWellEq(dt, domain, local_deferredLogger); +} + +template +void +BlackoilWellModelNldd:: +assembleWellEq(const double dt, + const Domain& domain, + DeferredLogger& deferred_logger) +{ + for (const auto& well : wellModel_.localNonshutWells()) { + if (this->well_domain().at(well->name()) == domain.index) { + well->assembleWellEq(wellModel_.simulator(), + dt, + wellModel_.wellState(), + wellModel_.groupState(), + deferred_logger); + } + } +} + +template +void +BlackoilWellModelNldd:: +addWellPressureEquations(PressureMatrix& /*jacobian*/, + const BVector& /*weights*/, + const bool /*use_well_weights*/, + const int /*domainIndex*/) const +{ + throw std::logic_error("CPRW is not yet implemented for NLDD subdomains"); + // To fix this function, rdofs should be the size of the domain, and the nw should be the number of wells in the domain + // int nw = this->numLocalWellsEnd(); // should number of wells in the domain + // int rdofs = local_num_cells_; // should be the size of the domain + // for ( int i = 0; i < nw; i++ ) { + // int wdof = rdofs + i; + // jacobian[wdof][wdof] = 1.0;// better scaling ? + // } + + // for ( const auto& well : well_container_ ) { + // if (well_domain_.at(well->name()) == domainIndex) { + // weights should be the size of the domain + // well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState()); + // } + // } +} + +template +void +BlackoilWellModelNldd:: +recoverWellSolutionAndUpdateWellState(const BVector& x, + const int domainIdx) +{ + // Note: no point in trying to do a parallel gathering + // try/catch here, as this function is not called in + // parallel but for each individual domain of each rank. + DeferredLogger local_deferredLogger; + for (const auto& well : wellModel_.localNonshutWells()) { + if (this->well_domain().at(well->name()) == domainIdx) { + const auto& cells = well->cells(); + x_local_.resize(cells.size()); + + for (size_t i = 0; i < cells.size(); ++i) { + x_local_[i] = x[cells[i]]; + } + well->recoverWellSolutionAndUpdateWellState(wellModel_.simulator(), + x_local_, + wellModel_.wellState(), + local_deferredLogger); + } + } + // TODO: avoid losing the logging information that could + // be stored in the local_deferredlogger in a parallel case. + if (wellModel_.terminalOutput()) { + local_deferredLogger.logMessages(); + } +} + +template +void +BlackoilWellModelNldd:: +initPrimaryVariablesEvaluation(const Domain& domain) const +{ + for (auto& well : wellModel_.localNonshutWells()) { + if (this->well_domain().at(well->name()) == domain.index) { + well->initPrimaryVariablesEvaluation(); + } + } +} + +template +ConvergenceReport +BlackoilWellModelNldd:: +getWellConvergence(const Domain& domain, + const std::vector& B_avg, + DeferredLogger& local_deferredLogger) const +{ + const int iterationIdx = wellModel_.simulator().model().newtonMethod().numIterations(); + const bool relax_tolerance = iterationIdx > wellModel_.numStrictIterations(); + + ConvergenceReport report; + for (const auto& well : wellModel_.localNonshutWells()) { + if ((this->well_domain().at(well->name()) == domain.index)) { + if (well->isOperableAndSolvable() || well->wellIsStopped()) { + report += well->getWellConvergence(wellModel_.simulator(), + wellModel_.wellState(), + B_avg, + local_deferredLogger, + relax_tolerance); + } else { + ConvergenceReport xreport; + using CR = ConvergenceReport; + xreport.setWellFailed({CR::WellFailure::Type::Unsolvable, + CR::Severity::Normal, -1, well->name()}); + report += xreport; + } + } + } + + // Log debug messages for NaN or too large residuals. + if (wellModel_.terminalOutput()) { + for (const auto& f : report.wellFailures()) { + if (f.severity() == ConvergenceReport::Severity::NotANumber) { + local_deferredLogger.debug("NaN residual found with phase " + + std::to_string(f.phase()) + + " for well " + f.wellName()); + } else if (f.severity() == ConvergenceReport::Severity::TooLarge) { + local_deferredLogger.debug("Too large residual found with phase " + + std::to_string(f.phase()) + + " for well " + f.wellName()); + } + } + } + return report; +} + +template +void +BlackoilWellModelNldd:: +updateWellControls(DeferredLogger& deferred_logger, + const Domain& domain) +{ + if (!wellModel_.wellsActive()) { + return; + } + + // TODO: decide on and implement an approach to handling of + // group controls, network and similar for domain solves. + + // Check only individual well constraints and communicate. + for (const auto& well : wellModel_.localNonshutWells()) { + if (this->well_domain().at(well->name()) == domain.index) { + constexpr auto mode = WellInterface::IndividualOrGroup::Individual; + well->updateWellControl(wellModel_.simulator(), + mode, + wellModel_.wellState(), + wellModel_.groupState(), + deferred_logger); + } + } +} + +template +void +BlackoilWellModelNldd:: +setupDomains(const std::vector& domains) +{ + std::vector genDomains; + std::transform(domains.begin(), domains.end(), + std::back_inserter(genDomains), + [](const auto& domain) + { return static_cast(&domain); }); + this->calcDomains(genDomains); +} + +} // namespace Opm + +#endif diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 287169ee0e..6c516c48fa 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -41,6 +41,7 @@ #include #include +#include #include #include #include @@ -1449,46 +1450,6 @@ namespace Opm { return well_group_thp_updated; } - template - void - BlackoilWellModel:: - assembleDomain([[maybe_unused]] const int iterationIdx, - const double dt, - const Domain& domain) - { - last_report_ = SimulatorReportSingle(); - Dune::Timer perfTimer; - perfTimer.start(); - - { - const int episodeIdx = simulator_.episodeIndex(); - const auto& network = this->schedule()[episodeIdx].network(); - if (!this->wellsActive() && !network.active()) { - return; - } - } - - // We assume that calculateExplicitQuantities() and - // prepareTimeStep() have been called already for the entire - // well model, so we do not need to do it here (when - // iterationIdx is 0). - - DeferredLogger local_deferredLogger; - updateWellControlsDomain(local_deferredLogger, domain); - initPrimaryVariablesEvaluationDomain(domain); - assembleWellEqDomain(dt, domain, local_deferredLogger); - - // TODO: errors here must be caught higher up, as this method is not called in parallel. - // We will log errors on rank 0, but not other ranks for now. - if (this->terminal_output_) { - local_deferredLogger.logMessages(); - } - - last_report_.converged = true; - last_report_.assemble_time_well += perfTimer.stop(); - } - - template void BlackoilWellModel:: @@ -1500,19 +1461,6 @@ namespace Opm { } - template - void - BlackoilWellModel:: - assembleWellEqDomain(const double dt, const Domain& domain, DeferredLogger& deferred_logger) - { - for (auto& well : well_container_) { - if (this->well_domain_.at(well->name()) == domain.index) { - well->assembleWellEq(simulator_, dt, this->wellState(), this->groupState(), deferred_logger); - } - } - } - - template void BlackoilWellModel:: @@ -1613,31 +1561,6 @@ namespace Opm { } } - template - void - BlackoilWellModel:: - addWellPressureEquationsDomain([[maybe_unused]] PressureMatrix& jacobian, - [[maybe_unused]] const BVector& weights, - [[maybe_unused]] const bool use_well_weights, - [[maybe_unused]] const int domainIndex) const - { - throw std::logic_error("CPRW is not yet implemented for NLDD subdomains"); - // To fix this function, rdofs should be the size of the domain, and the nw should be the number of wells in the domain - // int nw = this->numLocalWellsEnd(); // should number of wells in the domain - // int rdofs = local_num_cells_; // should be the size of the domain - // for ( int i = 0; i < nw; i++ ) { - // int wdof = rdofs + i; - // jacobian[wdof][wdof] = 1.0;// better scaling ? - // } - - // for ( const auto& well : well_container_ ) { - // if (well_domain_.at(well->name()) == domainIndex) { - // weights should be the size of the domain - // well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState()); - // } - // } - } - template void BlackoilWellModel:: addReservoirSourceTerms(GlobalEqVector& residual, @@ -1718,30 +1641,13 @@ namespace Opm { template void BlackoilWellModel:: - recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const Domain& domain) + recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const int domainIdx) { - // Note: no point in trying to do a parallel gathering - // try/catch here, as this function is not called in - // parallel but for each individual domain of each rank. - DeferredLogger local_deferredLogger; - for (auto& well : well_container_) { - if (this->well_domain_.at(well->name()) == domain.index) { - const auto& cells = well->cells(); - x_local_.resize(cells.size()); - - for (size_t i = 0; i < cells.size(); ++i) { - x_local_[i] = x[cells[i]]; - } - well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_, - this->wellState(), - local_deferredLogger); - } - } - // TODO: avoid losing the logging information that could - // be stored in the local_deferredlogger in a parallel case. - if (this->terminal_output_) { - local_deferredLogger.logMessages(); + if (!nldd_) { + OPM_THROW(std::logic_error, "Attempt to call NLDD method without a NLDD solver"); } + + return nldd_->recoverWellSolutionAndUpdateWellState(x, domainIdx); } @@ -1756,68 +1662,6 @@ namespace Opm { } - template - void - BlackoilWellModel:: - initPrimaryVariablesEvaluationDomain(const Domain& domain) const - { - for (auto& well : well_container_) { - if (this->well_domain_.at(well->name()) == domain.index) { - well->initPrimaryVariablesEvaluation(); - } - } - } - - - - - - - template - ConvergenceReport - BlackoilWellModel:: - getDomainWellConvergence(const Domain& domain, - const std::vector& B_avg, - DeferredLogger& local_deferredLogger) const - { - const int iterationIdx = simulator_.model().newtonMethod().numIterations(); - const bool relax_tolerance = iterationIdx > param_.strict_outer_iter_wells_; - - ConvergenceReport report; - for (const auto& well : well_container_) { - if ((this->well_domain_.at(well->name()) == domain.index)) { - if (well->isOperableAndSolvable() || well->wellIsStopped()) { - report += well->getWellConvergence(simulator_, - this->wellState(), - B_avg, - local_deferredLogger, - relax_tolerance); - } else { - ConvergenceReport xreport; - using CR = ConvergenceReport; - xreport.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()}); - report += xreport; - } - } - } - - // Log debug messages for NaN or too large residuals. - if (this->terminal_output_) { - for (const auto& f : report.wellFailures()) { - if (f.severity() == ConvergenceReport::Severity::NotANumber) { - local_deferredLogger.debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName()); - } else if (f.severity() == ConvergenceReport::Severity::TooLarge) { - local_deferredLogger.debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName()); - } - } - } - return report; - } - - - - - template ConvergenceReport BlackoilWellModel:: @@ -1981,29 +1825,6 @@ namespace Opm { return { changed_well_group, more_network_update }; } - template - void - BlackoilWellModel:: - updateWellControlsDomain(DeferredLogger& deferred_logger, const Domain& domain) - { - if ( !this->wellsActive() ) return ; - - // TODO: decide on and implement an approach to handling of - // group controls, network and similar for domain solves. - - // Check only individual well constraints and communicate. - for (const auto& well : well_container_) { - if (this->well_domain_.at(well->name()) == domain.index) { - const auto mode = WellInterface::IndividualOrGroup::Individual; - well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger); - } - } - } - - - - - template void BlackoilWellModel:: @@ -2480,87 +2301,6 @@ namespace Opm { } } } - - - template - void - BlackoilWellModel:: - setupDomains(const std::vector& domains) - { - OPM_BEGIN_PARALLEL_TRY_CATCH(); - // TODO: This loop nest may be slow for very large numbers of - // domains and wells, but that has not been observed on tests so - // far. Using the partition vector instead would be faster if we - // need to change. - for (const auto& wellPtr : this->well_container_) { - const int first_well_cell = wellPtr->cells().front(); - for (const auto& domain : domains) { - auto cell_present = [&domain](const auto cell) - { - return std::binary_search(domain.cells.begin(), - domain.cells.end(), cell); - }; - - if (cell_present(first_well_cell)) { - // Assuming that if the first well cell is found in a domain, - // then all of that well's cells are in that same domain. - this->well_domain_[wellPtr->name()] = domain.index; - - // Verify that all of that well's cells are in that same domain. - for (int well_cell : wellPtr->cells()) { - if (! cell_present(well_cell)) { - OPM_THROW(std::runtime_error, - fmt::format("Well '{}' found on multiple domains.", - wellPtr->name())); - } - } - } - } - } - OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::setupDomains(): well found on multiple domains.", - simulator_.gridView().comm()); - - // Write well/domain info to the DBG file. - const Opm::Parallel::Communication& comm = grid().comm(); - const int rank = comm.rank(); - DeferredLogger local_log; - if (!this->well_domain_.empty()) { - std::ostringstream os; - os << "Well name Rank Domain\n"; - for (const auto& [wname, domain] : this->well_domain_) { - os << wname << std::setw(19 - wname.size()) << rank << std::setw(12) << domain << '\n'; - } - local_log.debug(os.str()); - } - auto global_log = gatherDeferredLogger(local_log, comm); - if (this->terminal_output_) { - global_log.logMessages(); - } - - // Pre-calculate the local cell indices for each well - well_local_cells_.clear(); - well_local_cells_.reserve(well_container_.size(), 10); - std::vector local_cells; - for (const auto& well : well_container_) { - const auto& global_cells = well->cells(); - const int domain_index = this->well_domain_.at(well->name()); - const auto& domain_cells = domains[domain_index].cells; - local_cells.resize(global_cells.size()); - - // find the local cell index for each well cell in the domain - // assume domain_cells is sorted - for (size_t i = 0; i < global_cells.size(); ++i) { - auto it = std::lower_bound(domain_cells.begin(), domain_cells.end(), global_cells[i]); - if (it != domain_cells.end() && *it == global_cells[i]) { - local_cells[i] = std::distance(domain_cells.begin(), it); - } else { - OPM_THROW(std::runtime_error, fmt::format("Cell {} not found in domain {}", - global_cells[i], domain_index)); - } - } - well_local_cells_.appendRow(local_cells.begin(), local_cells.end()); - } - } } // namespace Opm #endif diff --git a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp index 72ff79ed69..473ce40f52 100644 --- a/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp +++ b/opm/simulators/wells/WellConnectionAuxiliaryModule.hpp @@ -128,7 +128,7 @@ class WellConnectionAuxiliaryModule : public BaseAuxiliaryModule void postSolveDomain(GlobalEqVector& deltaX, const Domain& domain) { - model_.recoverWellSolutionAndUpdateWellStateDomain(deltaX, domain); + model_.recoverWellSolutionAndUpdateWellStateDomain(deltaX, domain.index); } template diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index d6df2b7a68..60ea9ee70c 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -39,6 +39,10 @@ namespace Opm { #include +#include + +#include + #include #include #include