From 8fc2b8394890912edb0b67b9aed23d7fa9b32f23 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Jul 2023 15:37:50 +0200 Subject: [PATCH 1/7] BlackoilModelEbos: make initialLinearization public to make it possible to call externally --- opm/simulators/flow/BlackoilModelEbos.hpp | 91 ++++++++++++----------- 1 file changed, 46 insertions(+), 45 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 109e5f7da..681c17118 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -394,6 +394,52 @@ namespace Opm { } + void initialLinearization(SimulatorReportSingle& report, + const int iteration, + const int minIter, + const SimulatorTimerInterface& timer) + { + // ----------- Set up reports and timer ----------- + failureReport_ = SimulatorReportSingle(); + Dune::Timer perfTimer; + + perfTimer.start(); + report.total_linearizations = 1; + + // ----------- Assemble ----------- + try { + report += assembleReservoir(timer, iteration); + report.assemble_time += perfTimer.stop(); + } + catch (...) { + report.assemble_time += perfTimer.stop(); + failureReport_ += report; + throw; // continue throwing the stick + } + + // ----------- Check if converged ----------- + std::vector residual_norms; + perfTimer.reset(); + perfTimer.start(); + // the step is not considered converged until at least minIter iterations is done + { + auto convrep = getConvergence(timer, iteration, residual_norms); + report.converged = convrep.converged() && iteration > minIter; + ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); + convergence_reports_.back().report.push_back(std::move(convrep)); + + // Throw if any NaN or too large residual found. + if (severity == ConvergenceReport::Severity::NotANumber) { + OPM_THROW(NumericalProblem, "NaN residual found!"); + } else if (severity == ConvergenceReport::Severity::TooLarge) { + OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); + } + } + report.update_time += perfTimer.stop(); + residual_norms_history_.push_back(residual_norms); + } + + /// Called once per nonlinear iteration. /// This model will perform a Newton-Raphson update, changing reservoir_state /// and well_state. It will also use the nonlinear_solver to do relaxation of @@ -1677,51 +1723,6 @@ namespace Opm { } } - void initialLinearization(SimulatorReportSingle& report, - const int iteration, - const int minIter, - const SimulatorTimerInterface& timer) - { - // ----------- Set up reports and timer ----------- - failureReport_ = SimulatorReportSingle(); - Dune::Timer perfTimer; - - perfTimer.start(); - report.total_linearizations = 1; - - // ----------- Assemble ----------- - try { - report += assembleReservoir(timer, iteration); - report.assemble_time += perfTimer.stop(); - } - catch (...) { - report.assemble_time += perfTimer.stop(); - failureReport_ += report; - throw; // continue throwing the stick - } - - // ----------- Check if converged ----------- - std::vector residual_norms; - perfTimer.reset(); - perfTimer.start(); - // the step is not considered converged until at least minIter iterations is done - { - auto convrep = getConvergence(timer, iteration, residual_norms); - report.converged = convrep.converged() && iteration > minIter; - ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); - convergence_reports_.back().report.push_back(std::move(convrep)); - - // Throw if any NaN or too large residual found. - if (severity == ConvergenceReport::Severity::NotANumber) { - OPM_THROW(NumericalProblem, "NaN residual found!"); - } else if (severity == ConvergenceReport::Severity::TooLarge) { - OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!"); - } - } - report.update_time += perfTimer.stop(); - residual_norms_history_.push_back(residual_norms); - } - double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } double drMaxRel() const { return param_.dr_max_rel_; } From cee8dc4c6e412cadaea1eaa6e5ec5906e6865285 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Jul 2023 15:37:50 +0200 Subject: [PATCH 2/7] BlackoilModelEbos: make getMaxCoeff public to make it possible to call externally --- opm/simulators/flow/BlackoilModelEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 681c17118..8c1d33864 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -1606,7 +1606,6 @@ namespace Opm { ebosSimulator_.problem().endEpisode(); } - private: template void getMaxCoeff(const unsigned cell_idx, const IntensiveQuantities& intQuants, @@ -1723,6 +1722,7 @@ namespace Opm { } } + private: double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } double drMaxRel() const { return param_.dr_max_rel_; } From 6f726049ce3c163cad764c02971d769f829b6f90 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Jul 2023 15:41:57 +0200 Subject: [PATCH 3/7] BlackoilModelEbos: add accessor to linear solve setup time to allow for external usage --- opm/simulators/flow/BlackoilModelEbos.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 8c1d33864..24c17cc71 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -888,6 +888,13 @@ namespace Opm { } + // Obtain reference to linear solver setup time + double& linearSolveSetupTime() + { + return linear_solve_setup_time_; + } + + void solveJacobianSystemDomain(const Domain& domain, BVector& global_x) { Dune::Timer perfTimer; From 0c8446eed3ce99efb557f836996e8828341b8988 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 Jul 2023 10:25:06 +0200 Subject: [PATCH 4/7] BlackoilModelEbos: add acccessor for model parameters --- opm/simulators/flow/BlackoilModelEbos.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 24c17cc71..a947abf3b 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -1729,6 +1729,12 @@ namespace Opm { } } + //! \brief Returns const reference to model parameters. + const ModelParameters& param() const + { + return param_; + } + private: double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } From 184c6128cd9027b9d04fd2daadce6c4bd5620984 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 Jul 2023 10:25:06 +0200 Subject: [PATCH 5/7] BlackoilModelEbos: add acccessor for component names --- opm/simulators/flow/BlackoilModelEbos.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index a947abf3b..def6eee88 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -1735,6 +1735,12 @@ namespace Opm { return param_; } + //! \brief Returns const reference to component names. + const ComponentName& compNames() const + { + return compNames_; + } + private: double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } From 241c1d327931d4d63665702a44acba0f982ae556 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 Jul 2023 09:58:04 +0200 Subject: [PATCH 6/7] BlackoilModelEbos: use value in localAccumulatedReports cheap to copy and aids refactoring --- opm/simulators/flow/BlackoilModelEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index def6eee88..48cbb2639 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -1550,7 +1550,7 @@ namespace Opm { { return failureReport_; } /// return the statistics if the nonlinearIteration() method failed - const SimulatorReportSingle& localAccumulatedReports() const + SimulatorReportSingle localAccumulatedReports() const { return local_reports_accumulated_; } const std::vector& stepReports() const From 8e7de83218aa2edb623a9ec886262f0c252ab295 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 4 Jul 2023 15:15:29 +0200 Subject: [PATCH 7/7] BlackoilModelEbos: put Nldd solver in separate class --- CMakeLists_files.cmake | 1 + opm/simulators/flow/BlackoilModelEbos.hpp | 738 +--------------- opm/simulators/flow/BlackoilModelEbosNldd.hpp | 810 ++++++++++++++++++ 3 files changed, 827 insertions(+), 722 deletions(-) create mode 100644 opm/simulators/flow/BlackoilModelEbosNldd.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 9aeadc38e..39cf7ca25 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -360,6 +360,7 @@ list (APPEND PUBLIC_HEADER_FILES ebos/eclinterregflows.hh opm/simulators/flow/countGlobalCells.hpp opm/simulators/flow/BlackoilModelEbos.hpp + opm/simulators/flow/BlackoilModelEbosNldd.hpp opm/simulators/flow/BlackoilModelParametersEbos.hpp opm/simulators/flow/Banners.hpp opm/simulators/flow/ConvergenceOutputConfiguration.hpp diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 48cbb2639..4f4a1b541 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -34,21 +34,15 @@ #include -#include -#include - #include #include #include #include +#include #include -#include #include #include -#include -#include -#include #include #include #include @@ -67,8 +61,10 @@ #include #include #include +#include #include -#include +#include +#include #include #include #include @@ -151,6 +147,7 @@ struct LinearSolverSplice { } // namespace Opm::Properties namespace Opm { + /// A model implementation for three-phase black oil. /// /// The simulator is capable of handling three-phase problems @@ -209,9 +206,6 @@ namespace Opm { using BVector = Dune::BlockVector; using ComponentName = ::Opm::ComponentName; - using Domain = SubDomain; - - using ISTLSolverType = ISTLSolverEbos; // --------- Public methods --------- @@ -246,7 +240,7 @@ namespace Opm { if (terminal_output) { OpmLog::info("Using Non-Linear Domain Decomposition solver (nldd)."); } - setupSubDomains(); + nlddSolver_ = std::make_unique>(*this); } else if (param_.nonlinear_solver_ == "newton") { if (terminal_output) { OpmLog::info("Using Newton nonlinear solver."); @@ -265,91 +259,6 @@ namespace Opm { { return ebosSimulator_.vanguard().eclState(); } - - void setupSubDomains() - { - // Create partitions. - const auto& [partition_vector, num_domains] = - partitionCells(this->grid_, - this->ebosSimulator_.vanguard().schedule().getWellsatEnd(), - this->param_.local_domain_partition_method_, - this->param_.num_local_domains_, - this->param_.local_domain_partition_imbalance_); - - // Scan through partitioning to get correct size for each. - std::vector sizes(num_domains, 0); - for (const auto& p : partition_vector) { - ++sizes[p]; - } - - // Set up correctly sized vectors of entity seeds and of indices for each partition. - using EntitySeed = typename Grid::template Codim<0>::EntitySeed; - std::vector> seeds(num_domains); - std::vector> partitions(num_domains); - for (int domain = 0; domain < num_domains; ++domain) { - seeds[domain].resize(sizes[domain]); - partitions[domain].resize(sizes[domain]); - } - - // Iterate through grid once, setting the seeds of all partitions. - std::vector count(num_domains, 0); - const auto beg = grid_.template leafbegin<0>(); - const auto end = grid_.template leafend<0>(); - int cell = 0; - for (auto it = beg; it != end; ++it, ++cell) { - const int p = partition_vector[cell]; - seeds[p][count[p]] = it->seed(); - partitions[p][count[p]] = cell; - ++count[p]; - } - assert(count == sizes); - - // Create the domains. - for (int index = 0; index < num_domains; ++index) { - std::vector interior(partition_vector.size(), false); - for (int ix : partitions[index]) { - interior[ix] = true; - } - - Dune::SubGridPart view { - ebosSimulator_.vanguard().grid(), - std::move(seeds[index]) - }; - - this->domains_.emplace_back(index, - std::move(partitions[index]), - std::move(interior), - std::move(view)); - } - - // Set up container for the local system matrices. - domain_matrices_.resize(num_domains); - - // Set up container for the local linear solvers. - for (int index = 0; index < num_domains; ++index) { - // TODO: The ISTLSolverEbos constructor will make - // parallel structures appropriate for the full grid - // only. This must be addressed before going parallel. - FlowLinearSolverParameters param; - param.template init(ebosSimulator_.vanguard().eclState().getSimulationConfig().useCPR()); - // Override solver type with umfpack if small domain. - // Otherwise hardcode to ILU0 - if (domains_[index].cells.size() < 200) { - param.linsolver_ = "umfpack"; - } else { - param.linsolver_ = "ilu0"; - param.linear_solver_reduction_ = 1e-2; - } - param.linear_solver_print_json_definition_ = false; - domain_linsolvers_.emplace_back(ebosSimulator_, param); - } - - assert(int(domains_.size()) == num_domains); - } - - - - /// Called once before each time step. /// \param[in] timer simulation timer SimulatorReportSingle prepareStep(const SimulatorTimerInterface& timer) @@ -383,9 +292,8 @@ namespace Opm { //updateEquationsScaling(); } - if (!domains_.empty()) { - // Setup domain->well mapping. - wellModel().setupDomains(domains_); + if (nlddSolver_) { + nlddSolver_->prepareStep(); } report.pre_post_time += perfTimer.stop(); @@ -466,7 +374,7 @@ namespace Opm { return nonlinearIterationNewton(iteration, timer, nonlinear_solver); } if (param_.nonlinear_solver_ == "nldd") { - return nonlinearIterationNldd(iteration, timer, nonlinear_solver); + return nlddSolver_->nonlinearIterationNldd(iteration, timer, nonlinear_solver); } else { return nonlinearIterationNewton(iteration, timer, nonlinear_solver); } @@ -557,213 +465,6 @@ namespace Opm { } - - - template - SimulatorReportSingle nonlinearIterationNldd(const int iteration, - const SimulatorTimerInterface& timer, - NonlinearSolverType& nonlinear_solver) - { - // ----------- Set up reports and timer ----------- - SimulatorReportSingle report; - Dune::Timer perfTimer; - - this->initialLinearization(report, iteration, nonlinear_solver.minIter(), timer); - - if (report.converged) { - return report; - } - - // ----------- If not converged, do an NLDD iteration ----------- - - auto& solution = ebosSimulator().model().solution(0); - auto initial_solution = solution; - auto locally_solved = initial_solution; - - // ----------- Decide on an ordering for the domains ----------- - const auto domain_order = this->getSubdomainOrder(); - - // ----------- Solve each domain separately ----------- - std::vector domain_reports(domains_.size()); - for (const int domain_index : domain_order) { - const auto& domain = domains_[domain_index]; - SimulatorReportSingle local_report; - switch (param_.local_solve_approach_) { - case DomainSolveApproach::Jacobi: - solveDomainJacobi(solution, locally_solved, local_report, - iteration, timer, domain); - break; - default: - case DomainSolveApproach::GaussSeidel: - solveDomainGaussSeidel(solution, locally_solved, local_report, - iteration, timer, domain); - break; - } - // This should have updated the global matrix to be - // dR_i/du_j evaluated at new local solutions for - // i == j, at old solution for i != j. - if (!local_report.converged) { - // TODO: more proper treatment, including in parallel. - OpmLog::debug("Convergence failure in domain " + std::to_string(domain.index)); - } - domain_reports[domain.index] = local_report; - } - - // Log summary of local solve convergence to DBG file. - { - int num_converged = 0; - SimulatorReportSingle rep; - for (const auto& dr : domain_reports) { - if (dr.converged) { - ++num_converged; - } - rep += dr; - } - std::ostringstream os; - os << fmt::format("Local solves finished. Converged for {}/{} domains.\n", - num_converged, domain_reports.size()); - rep.reportFullyImplicit(os, nullptr); - OpmLog::debug(os.str()); - local_reports_accumulated_ += rep; - } - - if (param_.local_solve_approach_ == DomainSolveApproach::Jacobi) { - solution = locally_solved; - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); - } - - // Finish with a Newton step. - // Note that the "iteration + 100" is a simple way to avoid entering - // "if (iteration == 0)" and similar blocks, and also makes it a little - // easier to spot the iteration residuals in the DBG file. A more sophisticated - // approach can be done later. - auto rep = nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver); - report += rep; - if (rep.converged) { - report.converged = true; - } - return report; - } - - - - - std::pair - solveDomain(const Domain& domain, - const SimulatorTimerInterface& timer, - [[maybe_unused]] const int global_iteration, - const bool initial_assembly_required = false) - { - SimulatorReportSingle report; - Dune::Timer solveTimer; - solveTimer.start(); - Dune::Timer detailTimer; - - ebosSimulator_.model().newtonMethod().setIterationIndex(0); - - // When called, if assembly has already been performed - // with the initial values, we only need to check - // for local convergence. Otherwise, we must do a local - // assembly. - int iter = 0; - if (initial_assembly_required) { - detailTimer.start(); - ebosSimulator_.model().newtonMethod().setIterationIndex(iter); - // TODO: we should have a beginIterationLocal function() - // only handling the well model for now - ebosSimulator_.problem().wellModel().assembleDomain(ebosSimulator_.model().newtonMethod().numIterations(), - ebosSimulator_.timeStepSize(), - domain); - // Assemble reservoir locally. - report += assembleReservoirDomain(domain); - report.assemble_time += detailTimer.stop(); - } - detailTimer.reset(); - detailTimer.start(); - std::vector resnorms; - auto convreport = getDomainConvergence(domain, timer, 0, resnorms); - if (convreport.converged()) { - // TODO: set more info, timing etc. - report.converged = true; - return { report, convreport }; - } - - // We have already assembled for the first iteration, - // but not done the Schur complement for the wells yet. - detailTimer.reset(); - detailTimer.start(); - wellModel().linearizeDomain(domain, - ebosSimulator().model().linearizer().jacobian(), - ebosSimulator().model().linearizer().residual()); - const double tt1 = detailTimer.stop(); - report.assemble_time += tt1; - report.assemble_time_well += tt1; - - // Local Newton loop. - const int max_iter = param_.max_local_solve_iterations_; - do { - // Solve local linear system. - // Note that x has full size, we expect it to be nonzero only for in-domain cells. - const int nc = grid_.size(0); - BVector x(nc); - detailTimer.reset(); - detailTimer.start(); - solveJacobianSystemDomain(domain, x); - wellModel().postSolveDomain(x, domain); - report.linear_solve_time += detailTimer.stop(); - report.linear_solve_setup_time += linear_solve_setup_time_; - report.total_linear_iterations = linearIterationsLastSolve(); - - // Update local solution. // TODO: x is still full size, should we optimize it? - detailTimer.reset(); - detailTimer.start(); - updateDomainSolution(domain, x); - report.update_time += detailTimer.stop(); - - // Assemble well and reservoir. - detailTimer.reset(); - detailTimer.start(); - ++iter; - ebosSimulator_.model().newtonMethod().setIterationIndex(iter); - // TODO: we should have a beginIterationLocal function() - // only handling the well model for now - // Assemble reservoir locally. - ebosSimulator_.problem().wellModel().assembleDomain(ebosSimulator_.model().newtonMethod().numIterations(), - ebosSimulator_.timeStepSize(), - domain); - report += assembleReservoirDomain(domain); - report.assemble_time += detailTimer.stop(); - - // Check for local convergence. - detailTimer.reset(); - detailTimer.start(); - convreport = getDomainConvergence(domain, timer, iter, resnorms); - - // apply the Schur complement of the well model to the - // reservoir linearized equations - detailTimer.reset(); - detailTimer.start(); - wellModel().linearizeDomain(domain, - ebosSimulator().model().linearizer().jacobian(), - ebosSimulator().model().linearizer().residual()); - const double tt2 = detailTimer.stop(); - report.assemble_time += tt2; - report.assemble_time_well += tt2; - } while (!convreport.converged() && iter <= max_iter); - - ebosSimulator_.problem().endIteration(); - - report.converged = convreport.converged(); - report.total_newton_iterations = iter; - report.total_linearizations = iter; - report.total_time = solveTimer.stop(); - // TODO: set more info, timing etc. - return { report, convreport }; - } - - - - /// Called once after each time step. /// In this class, this function does nothing. /// \param[in] timer simulation timer @@ -789,14 +490,6 @@ namespace Opm { return wellModel().lastReport(); } - /// Assemble the residual and Jacobian of the nonlinear system. - SimulatorReportSingle assembleReservoirDomain(const Domain& domain) - { - // -------- Mass balance equations -------- - ebosSimulator_.model().linearizer().linearizeDomain(domain); - return wellModel().lastReport(); - } - // compute the "relative" change of the solution between time steps double relativeChange() const { @@ -895,35 +588,6 @@ namespace Opm { } - void solveJacobianSystemDomain(const Domain& domain, BVector& global_x) - { - Dune::Timer perfTimer; - perfTimer.start(); - - const Mat& main_matrix = ebosSimulator_.model().linearizer().jacobian().istlMatrix(); - if (domain_matrices_[domain.index]) { - Details::copySubMatrix(main_matrix, domain.cells, *domain_matrices_[domain.index]); - } else { - domain_matrices_[domain.index] = std::make_unique(Details::extractMatrix(main_matrix, domain.cells)); - } - auto& jac = *domain_matrices_[domain.index]; - auto res = Details::extractVector(ebosSimulator_.model().linearizer().residual(), domain.cells); - auto x = res; - - // set initial guess - global_x = 0.0; - x = 0.0; - - auto& linsolver = domain_linsolvers_[domain.index]; - - linsolver.prepare(jac, res); - linear_solve_setup_time_ = perfTimer.stop(); - linsolver.setResidual(res); - linsolver.solve(x); - - Details::setGlobal(x, domain.cells, global_x); - } - /// Solve the Jacobian system Jx = r where J is the Jacobian and /// r is the residual. void solveJacobianSystem(BVector& x) @@ -949,28 +613,6 @@ namespace Opm { } - - /// Apply an update to the primary variables. - void updateDomainSolution(const Domain& domain, const BVector& dx) - { - auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod(); - SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0); - - ebosNewtonMethod.update_(/*nextSolution=*/solution, - /*curSolution=*/solution, - /*update=*/dx, - /*resid=*/dx, - domain.cells); // the update routines of the black - // oil model do not care about the - // residual - - // if the solution is updated, the intensive quantities need to be recalculated - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); - } - - - - /// Apply an update to the primary variables. void updateSolution(const BVector& dx) { @@ -1113,95 +755,6 @@ namespace Opm { } - // Get reservoir quantities on this process needed for convergence calculations. - std::pair localDomainConvergenceData(const Domain& domain, - std::vector& R_sum, - std::vector& maxCoeff, - std::vector& B_avg, - std::vector& maxCoeffCell) - { - double pvSumLocal = 0.0; - double numAquiferPvSumLocal = 0.0; - const auto& ebosModel = ebosSimulator_.model(); - const auto& ebosProblem = ebosSimulator_.problem(); - - const auto& ebosResid = ebosSimulator_.model().linearizer().residual(); - - ElementContext elemCtx(ebosSimulator_); - const auto& gridView = domain.view; - const auto& elemEndIt = gridView.template end(); - IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid()); - OPM_BEGIN_PARALLEL_TRY_CATCH(); - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; - ++elemIt) - { - if (elemIt->partitionType() != Dune::InteriorEntity) { - continue; - } - const auto& elem = *elemIt; - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = intQuants.fluidState(); - - const auto pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * - ebosModel.dofTotalVolume(cell_idx); - pvSumLocal += pvValue; - - if (isNumericalAquiferCell(elem)) - { - numAquiferPvSumLocal += pvValue; - } - - this->getMaxCoeff(cell_idx, intQuants, fs, ebosResid, pvValue, - B_avg, R_sum, maxCoeff, maxCoeffCell); - } - - OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid_.comm()); - - // compute local average in terms of global number of elements - const int bSize = B_avg.size(); - for ( int i = 0; i& B_avg, double dt) - { - double errorPV{}; - const auto& ebosModel = ebosSimulator_.model(); - const auto& ebosProblem = ebosSimulator_.problem(); - const auto& ebosResid = ebosSimulator_.model().linearizer().residual(); - - for (const int cell_idx : domain.cells) - { - const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx ); - const auto& cellResidual = ebosResid[cell_idx]; - bool cnvViolated = false; - - for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) - { - using std::fabs; - Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue; - cnvViolated = cnvViolated || (fabs(CNV) > param_.tolerance_cnv_); - } - - if (cnvViolated) - { - errorPV += pvValue; - } - } - return errorPV; - } - - /// \brief Compute the total pore volume of cells violating CNV that are not part /// of a numerical aquifer. double computeCnvErrorPv(const std::vector& B_avg, double dt) @@ -1216,7 +769,6 @@ namespace Opm { IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid()); OPM_BEGIN_PARALLEL_TRY_CATCH(); - for (const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) { // Skip cells of numerical Aquifer @@ -1250,115 +802,6 @@ namespace Opm { } - ConvergenceReport getDomainReservoirConvergence(const double reportTime, - const double dt, - const int iteration, - const Domain& domain, - std::vector& B_avg, - std::vector& residual_norms) - { - typedef std::vector< Scalar > Vector; - - const int numComp = numEq; - Vector R_sum(numComp, 0.0 ); - Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() ); - std::vector maxCoeffCell(numComp, -1); - const auto [ pvSum, numAquiferPvSum] - = localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell); - - auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt); - cnvErrorPvFraction /= (pvSum - numAquiferPvSum); - - const double tol_mb = param_.local_tolerance_scaling_mb_ * param_.tolerance_mb_; - // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0. - // For each iteration, we need to determine whether to use the relaxed CNV tolerance. - // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0. - const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.min_strict_cnv_iter_; - // Tighter bound for local convergence should increase the - // likelyhood of: local convergence => global convergence - const double tol_cnv = param_.local_tolerance_scaling_cnv_ - * (use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_); - - // Finish computation - std::vector CNV(numComp); - std::vector mass_balance_residual(numComp); - for ( int compIdx = 0; compIdx < numComp; ++compIdx ) - { - CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx]; - mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum; - residual_norms.push_back(CNV[compIdx]); - } - - // Create convergence report. - ConvergenceReport report{reportTime}; - using CR = ConvergenceReport; - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] }; - CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance, - CR::ReservoirFailure::Type::Cnv }; - double tol[2] = { tol_mb, tol_cnv }; - for (int ii : {0, 1}) { - if (std::isnan(res[ii])) { - report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); - if ( terminal_output_ ) { - OpmLog::debug("NaN residual for " + compNames_.name(compIdx) + " equation."); - } - } else if (res[ii] > maxResidualAllowed()) { - report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); - if ( terminal_output_ ) { - OpmLog::debug("Too large residual for " + compNames_.name(compIdx) + " equation."); - } - } else if (res[ii] < 0.0) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); - if ( terminal_output_ ) { - OpmLog::debug("Negative residual for " + compNames_.name(compIdx) + " equation."); - } - } else if (res[ii] > tol[ii]) { - report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); - } - } - } - - // Output of residuals. - if ( terminal_output_ ) - { - // Only rank 0 does print to std::cout - if (iteration == 0) { - std::string msg = fmt::format("Domain {}, size {}, containing cell {}\n| Iter", - domain.index, domain.cells.size(), domain.cells[0]); - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - msg += " MB("; - msg += compNames_.name(compIdx)[0]; - msg += ") "; - } - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - msg += " CNV("; - msg += compNames_.name(compIdx)[0]; - msg += ") "; - } - OpmLog::debug(msg); - } - std::ostringstream ss; - ss << "| "; - const std::streamsize oprec = ss.precision(3); - const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); - ss << std::setw(4) << iteration; - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - ss << std::setw(11) << mass_balance_residual[compIdx]; - } - for (int compIdx = 0; compIdx < numComp; ++compIdx) { - ss << std::setw(11) << CNV[compIdx]; - } - ss.precision(oprec); - ss.flags(oflags); - OpmLog::debug(ss.str()); - } - - return report; - } - - - void updateTUNING(const Tuning& tuning) { param_.tolerance_mb_ = tuning.XXXMBE; if ( terminal_output_ ) { @@ -1475,22 +918,6 @@ namespace Opm { return report; } - ConvergenceReport getDomainConvergence(const Domain& domain, - const SimulatorTimerInterface& timer, - const int iteration, - std::vector& residual_norms) - { - std::vector B_avg(numEq, 0.0); - auto report = getDomainReservoirConvergence(timer.simulationTimeElapsed(), - timer.currentStepLength(), - iteration, - domain, - B_avg, - residual_norms); - report += wellModel().getDomainWellConvergence(domain, B_avg); - return report; - } - /// Compute convergence based on total mass balance (tol_mb) and maximum /// residual mass balance (tol_cnv). /// \param[in] timer simulation timer @@ -1551,7 +978,10 @@ namespace Opm { /// return the statistics if the nonlinearIteration() method failed SimulatorReportSingle localAccumulatedReports() const - { return local_reports_accumulated_; } + { + return nlddSolver_ ? nlddSolver_->localAccumulatedReports() + : SimulatorReportSingle{}; + } const std::vector& stepReports() const { @@ -1575,7 +1005,6 @@ namespace Opm { ModelParameters param_; SimulatorReportSingle failureReport_; - SimulatorReportSingle local_reports_accumulated_; // Well Model BlackoilWellModel& well_model_; @@ -1591,9 +1020,8 @@ namespace Opm { std::vector convergence_reports_; ComponentName compNames_{}; - std::vector domains_; - std::vector> domain_matrices_; - std::vector domain_linsolvers_; + + std::unique_ptr> nlddSolver_; //!< Non-linear DD solver public: /// return the StandardWells object @@ -1748,141 +1176,7 @@ namespace Opm { double maxResidualAllowed() const { return param_.max_residual_allowed_; } double linear_solve_setup_time_; - //! \brief Returns subdomain ordered according to method and ordering measure. - std::vector getSubdomainOrder() - { - const auto& solution = ebosSimulator().model().solution(0); - - std::vector domain_order(domains_.size()); - switch (param_.local_solve_approach_) { - case DomainSolveApproach::GaussSeidel: { - switch (param_.local_domain_ordering_) { - case DomainOrderingMeasure::AveragePressure: { - // Use average pressures to order domains. - std::vector> avgpress_per_domain(domains_.size()); - for (const auto& domain : domains_) { - double press_sum = 0.0; - for (const int c : domain.cells) { - press_sum += solution[c][Indices::pressureSwitchIdx]; - } - const double avgpress = press_sum / domain.cells.size(); - avgpress_per_domain[domain.index] = std::make_pair(avgpress, domain.index); - } - // Lexicographical sort by pressure, then index. - std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end()); - // Reverse since we want high-pressure regions solved first. - std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end()); - for (size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = avgpress_per_domain[ii].second; - } - break; - } - case DomainOrderingMeasure::Residual: { - // Use maximum residual to order domains. - const auto& residual = ebosSimulator().model().linearizer().residual(); - const int num_vars = residual[0].size(); - std::vector> maxres_per_domain(domains_.size()); - for (const auto& domain : domains_) { - double maxres = 0.0; - for (const int c : domain.cells) { - for (int ii = 0; ii < num_vars; ++ii) { - maxres = std::max(maxres, std::fabs(residual[c][ii])); - } - } - maxres_per_domain[domain.index] = std::make_pair(maxres, domain.index); - } - // Lexicographical sort by pressure, then index. - std::sort(maxres_per_domain.begin(), maxres_per_domain.end()); - // Reverse since we want high-pressure regions solved first. - std::reverse(maxres_per_domain.begin(), maxres_per_domain.end()); - for (size_t ii = 0; ii < domains_.size(); ++ii) { - domain_order[ii] = maxres_per_domain[ii].second; - } - } - break; - } - break; - } - - case DomainSolveApproach::Jacobi: - default: - std::iota(domain_order.begin(), domain_order.end(), 0); - break; - } - - return domain_order; - } - - template - void solveDomainJacobi(GlobalEqVector& solution, - GlobalEqVector& locally_solved, - SimulatorReportSingle& local_report, - const int iteration, - const SimulatorTimerInterface& timer, - const Domain& domain) - { - auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); - auto initial_local_solution = Details::extractVector(solution, domain.cells); - auto res = solveDomain(domain, timer, iteration); - local_report = res.first; - if (local_report.converged) { - auto local_solution = Details::extractVector(solution, domain.cells); - Details::setGlobal(local_solution, domain.cells, locally_solved); - Details::setGlobal(initial_local_solution, domain.cells, solution); - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); - } else { - wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); - Details::setGlobal(initial_local_solution, domain.cells, solution); - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); - } - } - - template - void solveDomainGaussSeidel(GlobalEqVector& solution, - GlobalEqVector& locally_solved, - SimulatorReportSingle& local_report, - const int iteration, - const SimulatorTimerInterface& timer, - const Domain& domain) - { - auto initial_local_well_primary_vars = wellModel().getPrimaryVarsDomain(domain); - auto initial_local_solution = Details::extractVector(solution, domain.cells); - auto res = solveDomain(domain, timer, iteration); - local_report = res.first; - if (!local_report.converged) { - // We look at the detailed convergence report to evaluate - // if we should accept the unconverged solution. - const auto& convrep = res.second; - // We do not accept a solution if the wells are unconverged. - if (!convrep.wellFailed()) { - // Calculare the sums of the mb and cnv failures. - double mb_sum = 0.0; - double cnv_sum = 0.0; - for (const auto& rc : convrep.reservoirConvergence()) { - if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) { - mb_sum += rc.value(); - } else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) { - cnv_sum += rc.value(); - } - } - // If not too high, we overrule the convergence failure. - const double acceptable_local_mb_sum = 1e-3; - const double acceptable_local_cnv_sum = 1.0; - if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) { - local_report.converged = true; - OpmLog::debug("Accepting solution in unconverged domain " + std::to_string(domain.index)); - } - } - } - if (local_report.converged) { - auto local_solution = Details::extractVector(solution, domain.cells); - Details::setGlobal(local_solution, domain.cells, locally_solved); - } else { - wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); - Details::setGlobal(initial_local_solution, domain.cells, solution); - ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); - } - } + SimulatorReportSingle local_reports_accumulated_; //!< Accumulated convergence report for subdomain solvers public: std::vector wasSwitched_; diff --git a/opm/simulators/flow/BlackoilModelEbosNldd.hpp b/opm/simulators/flow/BlackoilModelEbosNldd.hpp new file mode 100644 index 000000000..dbf24025c --- /dev/null +++ b/opm/simulators/flow/BlackoilModelEbosNldd.hpp @@ -0,0 +1,810 @@ +/* + Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2014, 2015 Statoil ASA. + Copyright 2015 NTNU + Copyright 2015, 2016, 2017 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_BLACKOILMODELEBOS_NLDD_HEADER_INCLUDED +#define OPM_BLACKOILMODELEBOS_NLDD_HEADER_INCLUDED + +#include + +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + +template class BlackoilModelEbos; + +/// A NLDD implementation for three-phase black oil. +template +class BlackoilModelEbosNldd { +public: + using ElementContext = GetPropType; + using FluidSystem = GetPropType; + using Grid = GetPropType; + using Indices = GetPropType; + using ModelParameters = BlackoilModelParametersEbos; + using Scalar = GetPropType; + using SolutionVector = GetPropType; + + using BVector = typename BlackoilModelEbos::BVector; + using Domain = SubDomain; + using ISTLSolverType = ISTLSolverEbos; + using Mat = typename BlackoilModelEbos::Mat; + + static constexpr int numEq = Indices::numEq; + + //! \brief The constructor sets up the subdomains. + //! \param model BlackOil model to solve for + //! \param param param Model parameters + //! \param compNames Names of the solution components + BlackoilModelEbosNldd(BlackoilModelEbos& model) + : model_(model) + { + const auto& grid = model_.ebosSimulator().vanguard().grid(); + const auto& schedule = model_.ebosSimulator().vanguard().schedule(); + + // Create partitions. + const auto& [partition_vector, num_domains] = + partitionCells(grid, + schedule.getWellsatEnd(), + model_.param().local_domain_partition_method_, + model_.param().num_local_domains_, + model_.param().local_domain_partition_imbalance_); + + // Scan through partitioning to get correct size for each. + std::vector sizes(num_domains, 0); + for (const auto& p : partition_vector) { + ++sizes[p]; + } + + // Set up correctly sized vectors of entity seeds and of indices for each partition. + using EntitySeed = typename Grid::template Codim<0>::EntitySeed; + std::vector> seeds(num_domains); + std::vector> partitions(num_domains); + for (int domain = 0; domain < num_domains; ++domain) { + seeds[domain].resize(sizes[domain]); + partitions[domain].resize(sizes[domain]); + } + + // Iterate through grid once, setting the seeds of all partitions. + std::vector count(num_domains, 0); + const auto beg = grid.template leafbegin<0>(); + const auto end = grid.template leafend<0>(); + int cell = 0; + for (auto it = beg; it != end; ++it, ++cell) { + const int p = partition_vector[cell]; + seeds[p][count[p]] = it->seed(); + partitions[p][count[p]] = cell; + ++count[p]; + } + assert(count == sizes); + + // Create the domains. + for (int index = 0; index < num_domains; ++index) { + std::vector interior(partition_vector.size(), false); + for (int ix : partitions[index]) { + interior[ix] = true; + } + + Dune::SubGridPart view{grid, std::move(seeds[index])}; + + this->domains_.emplace_back(index, + std::move(partitions[index]), + std::move(interior), + std::move(view)); + } + + // Set up container for the local system matrices. + domain_matrices_.resize(num_domains); + + // Set up container for the local linear solvers. + for (int index = 0; index < num_domains; ++index) { + // TODO: The ISTLSolverEbos constructor will make + // parallel structures appropriate for the full grid + // only. This must be addressed before going parallel. + const auto& eclState = model_.ebosSimulator().vanguard().eclState(); + FlowLinearSolverParameters loc_param; + loc_param.template init(eclState.getSimulationConfig().useCPR()); + // Override solver type with umfpack if small domain. + // Otherwise hardcode to ILU0 + if (domains_[index].cells.size() < 200) { + loc_param.linsolver_ = "umfpack"; + } else { + loc_param.linsolver_ = "ilu0"; + loc_param.linear_solver_reduction_ = 1e-2; + } + loc_param.linear_solver_print_json_definition_ = false; + domain_linsolvers_.emplace_back(model_.ebosSimulator(), loc_param); + } + + assert(int(domains_.size()) == num_domains); + } + + //! \brief Called before starting a time step. + void prepareStep() + { + // Setup domain->well mapping. + model_.wellModel().setupDomains(domains_); + } + + //! \brief Do one non-linear NLDD iteration. + template + SimulatorReportSingle nonlinearIterationNldd(const int iteration, + const SimulatorTimerInterface& timer, + NonlinearSolverType& nonlinear_solver) + { + // ----------- Set up reports and timer ----------- + SimulatorReportSingle report; + Dune::Timer perfTimer; + + model_.initialLinearization(report, iteration, nonlinear_solver.minIter(), timer); + + if (report.converged) { + return report; + } + + // ----------- If not converged, do an NLDD iteration ----------- + + auto& solution = model_.ebosSimulator().model().solution(0); + auto initial_solution = solution; + auto locally_solved = initial_solution; + + // ----------- Decide on an ordering for the domains ----------- + const auto domain_order = this->getSubdomainOrder(); + + // ----------- Solve each domain separately ----------- + std::vector domain_reports(domains_.size()); + for (const int domain_index : domain_order) { + const auto& domain = domains_[domain_index]; + SimulatorReportSingle local_report; + switch (model_.param().local_solve_approach_) { + case DomainSolveApproach::Jacobi: + solveDomainJacobi(solution, locally_solved, local_report, + iteration, timer, domain); + break; + default: + case DomainSolveApproach::GaussSeidel: + solveDomainGaussSeidel(solution, locally_solved, local_report, + iteration, timer, domain); + break; + } + // This should have updated the global matrix to be + // dR_i/du_j evaluated at new local solutions for + // i == j, at old solution for i != j. + if (!local_report.converged) { + // TODO: more proper treatment, including in parallel. + OpmLog::debug("Convergence failure in domain " + std::to_string(domain.index)); + } + domain_reports[domain.index] = local_report; + } + + // Log summary of local solve convergence to DBG file. + { + int num_converged = 0; + SimulatorReportSingle rep; + for (const auto& dr : domain_reports) { + if (dr.converged) { + ++num_converged; + } + rep += dr; + } + std::ostringstream os; + os << fmt::format("Local solves finished. Converged for {}/{} domains.\n", + num_converged, domain_reports.size()); + rep.reportFullyImplicit(os, nullptr); + OpmLog::debug(os.str()); + local_reports_accumulated_ += rep; + } + + if (model_.param().local_solve_approach_ == DomainSolveApproach::Jacobi) { + solution = locally_solved; + model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + + // Finish with a Newton step. + // Note that the "iteration + 100" is a simple way to avoid entering + // "if (iteration == 0)" and similar blocks, and also makes it a little + // easier to spot the iteration residuals in the DBG file. A more sophisticated + // approach can be done later. + auto rep = model_.nonlinearIterationNewton(iteration + 100, timer, nonlinear_solver); + report += rep; + if (rep.converged) { + report.converged = true; + } + return report; + } + + /// return the statistics if the nonlinearIteration() method failed + const SimulatorReportSingle& localAccumulatedReports() const + { + return local_reports_accumulated_; + } + +private: + //! \brief Solve the equation system for a single domain. + std::pair + solveDomain(const Domain& domain, + const SimulatorTimerInterface& timer, + [[maybe_unused]] const int global_iteration, + const bool initial_assembly_required = false) + { + auto& ebosSimulator = model_.ebosSimulator(); + + SimulatorReportSingle report; + Dune::Timer solveTimer; + solveTimer.start(); + Dune::Timer detailTimer; + + ebosSimulator.model().newtonMethod().setIterationIndex(0); + + // When called, if assembly has already been performed + // with the initial values, we only need to check + // for local convergence. Otherwise, we must do a local + // assembly. + int iter = 0; + if (initial_assembly_required) { + detailTimer.start(); + ebosSimulator.model().newtonMethod().setIterationIndex(iter); + // TODO: we should have a beginIterationLocal function() + // only handling the well model for now + ebosSimulator.problem().wellModel().assembleDomain(ebosSimulator.model().newtonMethod().numIterations(), + ebosSimulator.timeStepSize(), + domain); + // Assemble reservoir locally. + report += this->assembleReservoirDomain(domain); + report.assemble_time += detailTimer.stop(); + } + detailTimer.reset(); + detailTimer.start(); + std::vector resnorms; + auto convreport = this->getDomainConvergence(domain, timer, 0, resnorms); + if (convreport.converged()) { + // TODO: set more info, timing etc. + report.converged = true; + return { report, convreport }; + } + + // We have already assembled for the first iteration, + // but not done the Schur complement for the wells yet. + detailTimer.reset(); + detailTimer.start(); + model_.wellModel().linearizeDomain(domain, + ebosSimulator.model().linearizer().jacobian(), + ebosSimulator.model().linearizer().residual()); + const double tt1 = detailTimer.stop(); + report.assemble_time += tt1; + report.assemble_time_well += tt1; + + // Local Newton loop. + const int max_iter = model_.param().max_local_solve_iterations_; + const auto& grid = ebosSimulator.vanguard().grid(); + do { + // Solve local linear system. + // Note that x has full size, we expect it to be nonzero only for in-domain cells. + const int nc = grid.size(0); + BVector x(nc); + detailTimer.reset(); + detailTimer.start(); + this->solveJacobianSystemDomain(domain, x); + model_.wellModel().postSolveDomain(x, domain); + report.linear_solve_time += detailTimer.stop(); + report.linear_solve_setup_time += model_.linearSolveSetupTime(); + report.total_linear_iterations = model_.linearIterationsLastSolve(); + + // Update local solution. // TODO: x is still full size, should we optimize it? + detailTimer.reset(); + detailTimer.start(); + this->updateDomainSolution(domain, x); + report.update_time += detailTimer.stop(); + + // Assemble well and reservoir. + detailTimer.reset(); + detailTimer.start(); + ++iter; + ebosSimulator.model().newtonMethod().setIterationIndex(iter); + // TODO: we should have a beginIterationLocal function() + // only handling the well model for now + // Assemble reservoir locally. + ebosSimulator.problem().wellModel().assembleDomain(ebosSimulator.model().newtonMethod().numIterations(), + ebosSimulator.timeStepSize(), + domain); + report += this->assembleReservoirDomain(domain); + report.assemble_time += detailTimer.stop(); + + // Check for local convergence. + detailTimer.reset(); + detailTimer.start(); + convreport = this->getDomainConvergence(domain, timer, iter, resnorms); + + // apply the Schur complement of the well model to the + // reservoir linearized equations + detailTimer.reset(); + detailTimer.start(); + model_.wellModel().linearizeDomain(domain, + ebosSimulator.model().linearizer().jacobian(), + ebosSimulator.model().linearizer().residual()); + const double tt2 = detailTimer.stop(); + report.assemble_time += tt2; + report.assemble_time_well += tt2; + } while (!convreport.converged() && iter <= max_iter); + + ebosSimulator.problem().endIteration(); + + report.converged = convreport.converged(); + report.total_newton_iterations = iter; + report.total_linearizations = iter; + report.total_time = solveTimer.stop(); + // TODO: set more info, timing etc. + return { report, convreport }; + } + + /// Assemble the residual and Jacobian of the nonlinear system. + SimulatorReportSingle assembleReservoirDomain(const Domain& domain) + { + // -------- Mass balance equations -------- + model_.ebosSimulator().model().linearizer().linearizeDomain(domain); + return model_.wellModel().lastReport(); + } + + //! \brief Solve the linearized system for a domain. + void solveJacobianSystemDomain(const Domain& domain, BVector& global_x) + { + const auto& ebosSimulator = model_.ebosSimulator(); + + Dune::Timer perfTimer; + perfTimer.start(); + + const Mat& main_matrix = ebosSimulator.model().linearizer().jacobian().istlMatrix(); + if (domain_matrices_[domain.index]) { + Details::copySubMatrix(main_matrix, domain.cells, *domain_matrices_[domain.index]); + } else { + domain_matrices_[domain.index] = std::make_unique(Details::extractMatrix(main_matrix, domain.cells)); + } + auto& jac = *domain_matrices_[domain.index]; + auto res = Details::extractVector(ebosSimulator.model().linearizer().residual(), + domain.cells); + auto x = res; + + // set initial guess + global_x = 0.0; + x = 0.0; + + auto& linsolver = domain_linsolvers_[domain.index]; + + linsolver.prepare(jac, res); + model_.linearSolveSetupTime() = perfTimer.stop(); + linsolver.setResidual(res); + linsolver.solve(x); + + Details::setGlobal(x, domain.cells, global_x); + } + + /// Apply an update to the primary variables. + void updateDomainSolution(const Domain& domain, const BVector& dx) + { + auto& ebosSimulator = model_.ebosSimulator(); + auto& ebosNewtonMethod = ebosSimulator.model().newtonMethod(); + SolutionVector& solution = ebosSimulator.model().solution(/*timeIdx=*/0); + + ebosNewtonMethod.update_(/*nextSolution=*/solution, + /*curSolution=*/solution, + /*update=*/dx, + /*resid=*/dx, + domain.cells); // the update routines of the black + // oil model do not care about the + // residual + + // if the solution is updated, the intensive quantities need to be recalculated + ebosSimulator.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } + + //! \brief Get reservoir quantities on this process needed for convergence calculations. + std::pair localDomainConvergenceData(const Domain& domain, + std::vector& R_sum, + std::vector& maxCoeff, + std::vector& B_avg, + std::vector& maxCoeffCell) + { + const auto& ebosSimulator = model_.ebosSimulator(); + const auto& grid = model_.ebosSimulator().vanguard().grid(); + + double pvSumLocal = 0.0; + double numAquiferPvSumLocal = 0.0; + const auto& ebosModel = ebosSimulator.model(); + const auto& ebosProblem = ebosSimulator.problem(); + + const auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + ElementContext elemCtx(ebosSimulator); + const auto& gridView = domain.view; + const auto& elemEndIt = gridView.template end(); + IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid()); + + OPM_BEGIN_PARALLEL_TRY_CATCH(); + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; + ++elemIt) + { + if (elemIt->partitionType() != Dune::InteriorEntity) { + continue; + } + const auto& elem = *elemIt; + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + + const auto pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * + ebosModel.dofTotalVolume(cell_idx); + pvSumLocal += pvValue; + + if (isNumericalAquiferCell(elem)) + { + numAquiferPvSumLocal += pvValue; + } + + model_.getMaxCoeff(cell_idx, intQuants, fs, ebosResid, pvValue, + B_avg, R_sum, maxCoeff, maxCoeffCell); + } + OPM_END_PARALLEL_TRY_CATCH("BlackoilModelEbos::localConvergenceData() failed: ", grid.comm()); + + // compute local average in terms of global number of elements + const int bSize = B_avg.size(); + for ( int i = 0; i& B_avg, + std::vector& residual_norms) + { + using Vector = std::vector; + + const int numComp = numEq; + Vector R_sum(numComp, 0.0 ); + Vector maxCoeff(numComp, std::numeric_limits::lowest() ); + std::vector maxCoeffCell(numComp, -1); + const auto [ pvSum, numAquiferPvSum] + = this->localDomainConvergenceData(domain, R_sum, maxCoeff, B_avg, maxCoeffCell); + + auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt); + cnvErrorPvFraction /= (pvSum - numAquiferPvSum); + + const double tol_mb = model_.param().local_tolerance_scaling_mb_ * + model_.param().tolerance_mb_; + // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0. + // For each iteration, we need to determine whether to use the relaxed CNV tolerance. + // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0. + const bool use_relaxed = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ && + iteration >= model_.param().min_strict_cnv_iter_; + // Tighter bound for local convergence should increase the + // likelyhood of: local convergence => global convergence + const double tol_cnv = model_.param().local_tolerance_scaling_cnv_ + * (use_relaxed ? model_.param().tolerance_cnv_relaxed_ + : model_.param().tolerance_cnv_); + + // Finish computation + std::vector CNV(numComp); + std::vector mass_balance_residual(numComp); + for (int compIdx = 0; compIdx < numComp; ++compIdx ) + { + CNV[compIdx] = B_avg[compIdx] * dt * maxCoeff[compIdx]; + mass_balance_residual[compIdx] = std::abs(B_avg[compIdx]*R_sum[compIdx]) * dt / pvSum; + residual_norms.push_back(CNV[compIdx]); + } + + // Create convergence report. + ConvergenceReport report{reportTime}; + using CR = ConvergenceReport; + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] }; + CR::ReservoirFailure::Type types[2] = { CR::ReservoirFailure::Type::MassBalance, + CR::ReservoirFailure::Type::Cnv }; + double tol[2] = { tol_mb, tol_cnv }; + for (int ii : {0, 1}) { + if (std::isnan(res[ii])) { + report.setReservoirFailed({types[ii], CR::Severity::NotANumber, compIdx}); + if (model_.terminalOutputEnabled()) { + OpmLog::debug("NaN residual for " + model_.compNames().name(compIdx) + " equation."); + } + } else if (res[ii] > model_.param().max_residual_allowed_) { + report.setReservoirFailed({types[ii], CR::Severity::TooLarge, compIdx}); + if (model_.terminalOutputEnabled()) { + OpmLog::debug("Too large residual for " + model_.compNames().name(compIdx) + " equation."); + } + } else if (res[ii] < 0.0) { + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + if (model_.terminalOutputEnabled()) { + OpmLog::debug("Negative residual for " + model_.compNames().name(compIdx) + " equation."); + } + } else if (res[ii] > tol[ii]) { + report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); + } + } + } + + // Output of residuals. + if (model_.terminalOutputEnabled()) + { + // Only rank 0 does print to std::cout + if (iteration == 0) { + std::string msg = fmt::format("Domain {}, size {}, containing cell {}\n| Iter", + domain.index, domain.cells.size(), domain.cells[0]); + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + msg += " MB("; + msg += model_.compNames().name(compIdx)[0]; + msg += ") "; + } + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + msg += " CNV("; + msg += model_.compNames().name(compIdx)[0]; + msg += ") "; + } + OpmLog::debug(msg); + } + std::ostringstream ss; + ss << "| "; + const std::streamsize oprec = ss.precision(3); + const std::ios::fmtflags oflags = ss.setf(std::ios::scientific); + ss << std::setw(4) << iteration; + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + ss << std::setw(11) << mass_balance_residual[compIdx]; + } + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + ss << std::setw(11) << CNV[compIdx]; + } + ss.precision(oprec); + ss.flags(oflags); + OpmLog::debug(ss.str()); + } + + return report; + } + + ConvergenceReport getDomainConvergence(const Domain& domain, + const SimulatorTimerInterface& timer, + const int iteration, + std::vector& residual_norms) + { + std::vector B_avg(numEq, 0.0); + auto report = this->getDomainReservoirConvergence(timer.simulationTimeElapsed(), + timer.currentStepLength(), + iteration, + domain, + B_avg, + residual_norms); + report += model_.wellModel().getDomainWellConvergence(domain, B_avg); + return report; + } + + //! \brief Returns subdomain ordered according to method and ordering measure. + std::vector getSubdomainOrder() + { + const auto& ebosSimulator = model_.ebosSimulator(); + const auto& solution = ebosSimulator.model().solution(0); + + std::vector domain_order(domains_.size()); + switch (model_.param().local_solve_approach_) { + case DomainSolveApproach::GaussSeidel: { + switch (model_.param().local_domain_ordering_) { + case DomainOrderingMeasure::AveragePressure: { + // Use average pressures to order domains. + std::vector> avgpress_per_domain(domains_.size()); + for (const auto& domain : domains_) { + double press_sum = 0.0; + for (const int c : domain.cells) { + press_sum += solution[c][Indices::pressureSwitchIdx]; + } + const double avgpress = press_sum / domain.cells.size(); + avgpress_per_domain[domain.index] = std::make_pair(avgpress, domain.index); + } + // Lexicographical sort by pressure, then index. + std::sort(avgpress_per_domain.begin(), avgpress_per_domain.end()); + // Reverse since we want high-pressure regions solved first. + std::reverse(avgpress_per_domain.begin(), avgpress_per_domain.end()); + for (std::size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = avgpress_per_domain[ii].second; + } + break; + } + case DomainOrderingMeasure::Residual: { + // Use maximum residual to order domains. + const auto& residual = ebosSimulator.model().linearizer().residual(); + const int num_vars = residual[0].size(); + std::vector> maxres_per_domain(domains_.size()); + for (const auto& domain : domains_) { + double maxres = 0.0; + for (const int c : domain.cells) { + for (int ii = 0; ii < num_vars; ++ii) { + maxres = std::max(maxres, std::fabs(residual[c][ii])); + } + } + maxres_per_domain[domain.index] = std::make_pair(maxres, domain.index); + } + // Lexicographical sort by pressure, then index. + std::sort(maxres_per_domain.begin(), maxres_per_domain.end()); + // Reverse since we want high-pressure regions solved first. + std::reverse(maxres_per_domain.begin(), maxres_per_domain.end()); + for (std::size_t ii = 0; ii < domains_.size(); ++ii) { + domain_order[ii] = maxres_per_domain[ii].second; + } + } + break; + } + break; + } + + case DomainSolveApproach::Jacobi: + default: + std::iota(domain_order.begin(), domain_order.end(), 0); + break; + } + + return domain_order; + } + + template + void solveDomainJacobi(GlobalEqVector& solution, + GlobalEqVector& locally_solved, + SimulatorReportSingle& local_report, + const int iteration, + const SimulatorTimerInterface& timer, + const Domain& domain) + { + auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveDomain(domain, timer, iteration); + local_report = res.first; + if (local_report.converged) { + auto local_solution = Details::extractVector(solution, domain.cells); + Details::setGlobal(local_solution, domain.cells, locally_solved); + Details::setGlobal(initial_local_solution, domain.cells, solution); + model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } else { + model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } + } + + template + void solveDomainGaussSeidel(GlobalEqVector& solution, + GlobalEqVector& locally_solved, + SimulatorReportSingle& local_report, + const int iteration, + const SimulatorTimerInterface& timer, + const Domain& domain) + { + auto initial_local_well_primary_vars = model_.wellModel().getPrimaryVarsDomain(domain); + auto initial_local_solution = Details::extractVector(solution, domain.cells); + auto res = solveDomain(domain, timer, iteration); + local_report = res.first; + if (!local_report.converged) { + // We look at the detailed convergence report to evaluate + // if we should accept the unconverged solution. + const auto& convrep = res.second; + // We do not accept a solution if the wells are unconverged. + if (!convrep.wellFailed()) { + // Calculare the sums of the mb and cnv failures. + double mb_sum = 0.0; + double cnv_sum = 0.0; + for (const auto& rc : convrep.reservoirConvergence()) { + if (rc.type() == ConvergenceReport::ReservoirFailure::Type::MassBalance) { + mb_sum += rc.value(); + } else if (rc.type() == ConvergenceReport::ReservoirFailure::Type::Cnv) { + cnv_sum += rc.value(); + } + } + // If not too high, we overrule the convergence failure. + const double acceptable_local_mb_sum = 1e-3; + const double acceptable_local_cnv_sum = 1.0; + if (mb_sum < acceptable_local_mb_sum && cnv_sum < acceptable_local_cnv_sum) { + local_report.converged = true; + OpmLog::debug("Accepting solution in unconverged domain " + std::to_string(domain.index)); + } + } + } + if (local_report.converged) { + auto local_solution = Details::extractVector(solution, domain.cells); + Details::setGlobal(local_solution, domain.cells, locally_solved); + } else { + model_.wellModel().setPrimaryVarsDomain(domain, initial_local_well_primary_vars); + Details::setGlobal(initial_local_solution, domain.cells, solution); + model_.ebosSimulator().model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0, domain.view); + } + } + + double computeCnvErrorPvLocal(const Domain& domain, + const std::vector& B_avg, double dt) const + { + double errorPV{}; + const auto& ebosSimulator = model_.ebosSimulator(); + const auto& ebosModel = ebosSimulator.model(); + const auto& ebosProblem = ebosSimulator.problem(); + const auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + for (const int cell_idx : domain.cells) { + const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * + ebosModel.dofTotalVolume(cell_idx); + const auto& cellResidual = ebosResid[cell_idx]; + bool cnvViolated = false; + + for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx) { + using std::fabs; + Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue; + cnvViolated = cnvViolated || (fabs(CNV) > model_.param().tolerance_cnv_); + } + + if (cnvViolated) { + errorPV += pvValue; + } + } + return errorPV; + } + + BlackoilModelEbos& model_; //!< Reference to model + 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 + SimulatorReportSingle local_reports_accumulated_; //!< Accumulated convergence report for subdomain solvers +}; + +} // namespace Opm + +#endif // OPM_BLACKOILMODELEBOS_NLDD_HEADER_INCLUDED