From 654df6fd596e9467eb5bc5583b1b6d1a9a228a31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 15 Jun 2023 14:26:56 +0200 Subject: [PATCH] Well additions for local solves. Also, remove uneeded function updatePerforationIntensiveQuantities(). --- opm/simulators/wells/BlackoilWellModel.hpp | 25 +- .../wells/BlackoilWellModel_impl.hpp | 260 ++++++++++++++---- tests/test_glift1.cpp | 1 - 3 files changed, 236 insertions(+), 50 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index c34469f2c..f3ae7b23c 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -158,12 +158,18 @@ namespace Opm { {} void linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) override; + void linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res); void postSolve(GlobalEqVector& deltaX) override { recoverWellSolutionAndUpdateWellState(deltaX); } + void postSolveDomain(GlobalEqVector& deltaX, const Domain& domain) + { + recoverWellSolutionAndUpdateWellStateDomain(deltaX, domain); + } + ///////////// // ///////////// @@ -273,6 +279,10 @@ namespace Opm { // 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) const; + const SimulatorReportSingle& lastReport() const; void addWellContributions(SparseMatrixAdapter& jacobian) const; @@ -284,7 +294,6 @@ namespace Opm { // called at the beginning of a report step void beginReportStep(const int time_step); - void updatePerforationIntensiveQuantities(); // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation() // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions // twice at the beginning of the time step @@ -295,6 +304,7 @@ namespace Opm { // 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); @@ -331,6 +341,13 @@ namespace Opm { int numLocalNonshutWells() const; + // 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 logPrimaryVars() const; std::vector getPrimaryVarsDomain(const Domain& domain) const; void setPrimaryVarsDomain(const Domain& domain, const std::vector& vars); @@ -416,6 +433,7 @@ namespace Opm { const double dt, DeferredLogger& local_deferredLogger); + // called at the end of a time step void timeStepSucceeded(const double& simulationTime, const double dt); @@ -426,6 +444,10 @@ namespace Opm { // xw to update Well State void recoverWellSolutionAndUpdateWellState(const BVector& x); + // 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); + // setting the well_solutions_ based on well_state. void updatePrimaryVariables(DeferredLogger& deferred_logger); @@ -450,6 +472,7 @@ namespace Opm { 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); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 0c06b1ad5..04ed73c9d 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -38,6 +38,7 @@ #endif #include +#include #include #include @@ -115,6 +116,7 @@ namespace Opm { } } + template void BlackoilWellModel:: @@ -138,34 +140,48 @@ namespace Opm { } } + template void BlackoilWellModel:: linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res) { - if (!param_.matrix_add_well_contributions_) - { - OPM_BEGIN_PARALLEL_TRY_CATCH(); - { - // if the well contributions are not supposed to be included explicitly in - // the matrix, we only apply the vector part of the Schur complement here. - for (const auto& well: well_container_) { - // r = r - duneC_^T * invDuneD_ * resWell_ - well->apply(res); - } - } - OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ", - ebosSimulator_.gridView().comm()); - return; - } - + OPM_BEGIN_PARALLEL_TRY_CATCH(); for (const auto& well: well_container_) { - well->addWellContributions(jacobian); - - // applying the well residual to reservoir residuals + // Modifiy the Jacobian with explicit Schur complement + // contributions if requested. + if (param_.matrix_add_well_contributions_) { + well->addWellContributions(jacobian); + } + // Apply as Schur complement the well residual to reservoir residuals: // r = r - duneC_^T * invDuneD_ * resWell_ well->apply(res); } + OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ", + ebosSimulator_.gridView().comm()); + } + + + template + void + BlackoilWellModel:: + linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res) + { + // 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. + for (const auto& well: well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + // Modifiy the Jacobian with explicit Schur complement + // contributions if requested. + if (param_.matrix_add_well_contributions_) { + well->addWellContributions(jacobian); + } + // Apply as Schur complement the well residual to reservoir residuals: + // r = r - duneC_^T * invDuneD_ * resWell_ + well->apply(res); + } + } } @@ -240,7 +256,6 @@ namespace Opm { beginTimeStep() { OPM_TIMEBLOCK(beginTimeStep); - updatePerforationIntensiveQuantities(); updateAverageFormationFactor(); DeferredLogger local_deferredLogger; switched_prod_groups_.clear(); @@ -914,8 +929,6 @@ namespace Opm { } } - updatePerforationIntensiveQuantities(); - if (iterationIdx == 0 && wellsActive()) { // try-catch is needed here as updateWellControls // contains global communication and has either to @@ -1026,6 +1039,45 @@ namespace Opm { + template + void + BlackoilWellModel:: + assembleDomain(const int iterationIdx, + const double dt, + const Domain& domain) + { + last_report_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + perfTimer.start(); + + { + const int episodeIdx = ebosSimulator_.episodeIndex(); + const auto& network = schedule()[episodeIdx].network(); + if ( !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 (terminal_output_) { + local_deferredLogger.logMessages(); + } + + last_report_.converged = true; + last_report_.assemble_time_well += perfTimer.stop(); + } + template bool @@ -1238,6 +1290,20 @@ namespace Opm { } } + + template + void + BlackoilWellModel:: + assembleWellEqDomain(const double dt, const Domain& domain, DeferredLogger& deferred_logger) + { + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger); + } + } + } + + template void BlackoilWellModel:: @@ -1453,6 +1519,7 @@ namespace Opm { } } + template int BlackoilWellModel:: @@ -1461,28 +1528,49 @@ namespace Opm { return well_container_.size(); } + template void BlackoilWellModel:: recoverWellSolutionAndUpdateWellState(const BVector& x) { - DeferredLogger local_deferredLogger; OPM_BEGIN_PARALLEL_TRY_CATCH(); { + const auto& summary_state = ebosSimulator_.vanguard().summaryState(); for (auto& well : well_container_) { - const auto& summary_state = ebosSimulator_.vanguard().summaryState(); well->recoverWellSolutionAndUpdateWellState(summary_state, x, this->wellState(), local_deferredLogger); } - } OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "recoverWellSolutionAndUpdateWellState() failed: ", terminal_output_, ebosSimulator_.vanguard().grid().comm()); - } + template + void + BlackoilWellModel:: + recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const Domain& domain) + { + // 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; + const auto& summary_state = this->ebosSimulator_.vanguard().summaryState(); + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->recoverWellSolutionAndUpdateWellState(summary_state, x, + this->wellState(), + local_deferredLogger); + } + } + // TODO: avoid losing the logging information that could + // be stored in the local_deferredlogger in a parallel case. + if (terminal_output_) { + local_deferredLogger.logMessages(); + } + } template @@ -1496,6 +1584,84 @@ namespace Opm { } + template + void + BlackoilWellModel:: + initPrimaryVariablesEvaluationDomain(const Domain& domain) const + { + for (auto& well : well_container_) { + if (well_domain_.at(well->name()) == domain.index) { + well->initPrimaryVariablesEvaluation(); + } + } + } + + + + + + + template + ConvergenceReport + BlackoilWellModel:: + getDomainWellConvergence(const Domain& domain, + const std::vector& B_avg) const + { + const auto& summary_state = ebosSimulator_.vanguard().summaryState(); + const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations(); + const bool relax_tolerance = iterationIdx > param_.strict_outer_iter_wells_; + + Opm::DeferredLogger local_deferredLogger; + + ConvergenceReport local_report; + for (const auto& well : well_container_) { + if ((well_domain_.at(well->name()) == domain.index)) { + if (well->isOperableAndSolvable() || well->wellIsStopped()) { + local_report += well->getWellConvergence(summary_state, + this->wellState(), + B_avg, + local_deferredLogger, + iterationIdx > param_.strict_outer_iter_wells_); + } else { + ConvergenceReport report; + using CR = ConvergenceReport; + report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()}); + local_report += report; + } + } + } + + // This function will be called once for each domain on a process, + // and the number of domains on a process will vary, so there is + // no way to communicate here. There is also no need, as a domain + // is local to a single process in our current approach. + // Therefore there is no call to gatherDeferredLogger() or to + // gatherConvergenceReport() below. However, as of now, log messages + // on non-output ranks will be lost here. + // TODO: create the DeferredLogger outside this scope to enable it + // to be gathered and printed on the output rank. + Opm::DeferredLogger global_deferredLogger = local_deferredLogger; + ConvergenceReport report = local_report; + if (terminal_output_) { + global_deferredLogger.logMessages(); + } + + // Log debug messages for NaN or too large residuals. + // TODO: This will as of now only be logged on the output rank. + // In the similar code in getWellConvergence(), all ranks will be + // at the same spot, that does not hold for this per-domain function. + if (terminal_output_) { + for (const auto& f : report.wellFailures()) { + if (f.severity() == ConvergenceReport::Severity::NotANumber) { + OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName()); + } else if (f.severity() == ConvergenceReport::Severity::TooLarge) { + OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName()); + } + } + } + return report; + } + @@ -1641,6 +1807,26 @@ namespace Opm { return { changed_well_group, more_network_update }; } + template + void + BlackoilWellModel:: + updateWellControlsDomain(DeferredLogger& deferred_logger, const Domain& domain) + { + if ( !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 (well_domain_.at(well->name()) == domain.index) { + const auto mode = WellInterface::IndividualOrGroup::Individual; + well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger); + } + } + } + + template void @@ -1944,28 +2130,6 @@ namespace Opm { } } - template - void - BlackoilWellModel:: - updatePerforationIntensiveQuantities() - { - ElementContext elemCtx(ebosSimulator_); - const auto& gridView = ebosSimulator_.gridView(); - - OPM_BEGIN_PARALLEL_TRY_CATCH(); - for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { - elemCtx.updatePrimaryStencil(elem); - int elemIdx = elemCtx.globalSpaceIndex(0, 0); - - if (!is_cell_perforated_[elemIdx]) { - continue; - } - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - } - OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updatePerforationIntensiveQuantities() failed: ", ebosSimulator_.vanguard().grid().comm()); - } - - template typename BlackoilWellModel::WellInterfacePtr BlackoilWellModel:: diff --git a/tests/test_glift1.cpp b/tests/test_glift1.cpp index 4518a37da..2ee603cc0 100644 --- a/tests/test_glift1.cpp +++ b/tests/test_glift1.cpp @@ -143,7 +143,6 @@ BOOST_AUTO_TEST_CASE(G1) int report_step_idx = 0; well_model.beginReportStep(report_step_idx); well_model.beginTimeStep(); - well_model.updatePerforationIntensiveQuantities(); Opm::DeferredLogger deferred_logger; well_model.calculateExplicitQuantities(deferred_logger); well_model.prepareTimeStep(deferred_logger);